suppressPackageStartupMessages({
library(tidyverse)
library(Seurat)
library(magrittr)
library(dplyr)
library(purrr)
library(ggplot2)
library(here)
library(runSeurat3)
library(SingleCellExperiment)
library(RColorBrewer)
library(pheatmap)
library(scater)
library(scran)
library(ggsci)
})label and characterize adult iLN
load packages
heatmap function
avgHeatmap <- function(seurat, selGenes, colVecIdent, colVecCond=NULL,
ordVec=NULL, gapVecR=NULL, gapVecC=NULL,cc=FALSE,
cr=FALSE, condCol=FALSE){
selGenes <- selGenes$gene
## assay data
clusterAssigned <- as.data.frame(Idents(seurat)) %>%
dplyr::mutate(cell=rownames(.))
colnames(clusterAssigned)[1] <- "ident"
seuratDat <- GetAssayData(seurat)
## genes of interest
genes <- data.frame(gene=rownames(seurat)) %>%
mutate(geneID=gsub("^.*\\.", "", gene)) %>% filter(geneID %in% selGenes)
## matrix with averaged cnts per ident
logNormExpres <- as.data.frame(t(as.matrix(
seuratDat[which(rownames(seuratDat) %in% genes$gene),])))
logNormExpres <- logNormExpres %>% dplyr::mutate(cell=rownames(.)) %>%
dplyr::left_join(.,clusterAssigned, by=c("cell")) %>%
dplyr::select(-cell) %>% dplyr::group_by(ident) %>%
dplyr::summarise_all(mean)
logNormExpresMa <- logNormExpres %>% dplyr::select(-ident) %>% as.matrix()
rownames(logNormExpresMa) <- logNormExpres$ident
logNormExpresMa <- t(logNormExpresMa)
rownames(logNormExpresMa) <- gsub("^.*?\\.","",rownames(logNormExpresMa))
## remove genes if they are all the same in all groups
ind <- apply(logNormExpresMa, 1, sd) == 0
logNormExpresMa <- logNormExpresMa[!ind,]
genes <- genes[!ind,]
## color columns according to cluster
annotation_col <- as.data.frame(gsub("(^.*?_)","",
colnames(logNormExpresMa)))%>%
dplyr::mutate(celltype=gsub("(_.*$)","",colnames(logNormExpresMa)))
colnames(annotation_col)[1] <- "col1"
annotation_col <- annotation_col %>%
dplyr::mutate(cond = gsub(".*_","",col1)) %>%
dplyr::select(cond, celltype)
rownames(annotation_col) <- colnames(logNormExpresMa)
ann_colors = list(
cond = colVecCond,
celltype=colVecIdent)
if(is.null(ann_colors$cond)){
annotation_col$cond <- NULL
}
## adjust order
logNormExpresMa <- logNormExpresMa[selGenes,]
if(is.null(ordVec)){
ordVec <- levels(seurat)
}
logNormExpresMa <- logNormExpresMa[,ordVec]
## scaled row-wise
pheatmap(logNormExpresMa, scale="row" ,treeheight_row = 0, cluster_rows = cr,
cluster_cols = cc,
color = colorRampPalette(c("#2166AC", "#F7F7F7", "#B2182B"))(50),
annotation_col = annotation_col, cellwidth=15, cellheight=10,
annotation_colors = ann_colors, gaps_row = gapVecR, gaps_col = gapVecC)
}set dir and load sample
basedir <- here()
seurat <- readRDS(paste0(basedir,
"/data/WT_adultOnly_iLNonly_filtered_seurat.rds"))
colLoc <- c("#61baba", "#ba6161")
names(colLoc) <- unique(seurat$location)
colAge <- c("#440154FF", "#3B528BFF", "#21908CFF", "#5DC863FF", "#FDE725FF")
names(colAge) <- c("E18" , "P7", "3w", "8w","E17to7wk")
colPal <- c("#DAF7A6", "#FFC300", "#FF5733", "#C70039", "#900C3F", "#b66e8d",
"#61a4ba", "#6178ba", "#54a87f", "#5468a8", "#25328a",
"#b6856e", "#0073C2FF", "#EFC000FF", "#868686FF", "#CD534CFF",
"#7AA6DCFF", "#003C67FF", "#8F7700FF", "#3B3B3BFF", "#A73030FF",
"#4A6990FF")[1:length(unique(seurat$RNA_snn_res.0.4))]
names(colPal) <- unique(seurat$RNA_snn_res.0.4)
colDat <- colDat <- c(pal_npg()(10),pal_futurama()(12), pal_aaas()(10),
pal_jama()(8))[1:length(unique(seurat$dataset))]
names(colDat) <- unique(seurat$dataset)DimPlot all
clustering
DimPlot(seurat, reduction = "umap", group.by = "RNA_snn_res.0.4",
cols = colPal)+
theme_bw() +
theme(axis.text = element_blank(), axis.ticks = element_blank(),
panel.grid.minor = element_blank()) +
xlab("UMAP1") +
ylab("UMAP2")
add label
seurat$label <- "TRC"
seurat$label[which(seurat$RNA_snn_res.0.4 == "1")] <- "medRC1"
seurat$label[which(seurat$RNA_snn_res.0.4 == "0")] <- "IFRC"
seurat$label[which(seurat$RNA_snn_res.0.4 == "5")] <- "medRC2"
seurat$label[which(seurat$RNA_snn_res.0.4 == "3")] <- "TBRC"
seurat$label[which(seurat$RNA_snn_res.0.4 == "4")] <- "PRC1"
seurat$label[which(seurat$RNA_snn_res.0.4 == "6")] <- "Pi16+RC"
seurat$label[which(seurat$RNA_snn_res.0.4 == "9")] <- "PRC3"
seurat$label[which(seurat$RNA_snn_res.0.4 == "7")] <- "VSMC"
seurat$label[which(seurat$RNA_snn_res.0.4 == "8")] <- "FDC"
colLab <- c("#42a071", "#900C3F","#b66e8d", "#8F7700FF", "#61a4ba","#003C67FF",
"#e3953d","#714542", "#b6856e", "#FFC300")
names(colLab) <- c("FDC", "TRC", "TBRC", "IFRC", "medRC1" , "medRC2",
"PRC1", "Pi16+RC", "PRC3", "VSMC")vis label
DimPlot(seurat, reduction = "umap", group.by = "label",
cols = colLab)+
theme_bw() +
theme(axis.text = element_blank(), axis.ticks = element_blank(),
panel.grid.minor = element_blank()) +
xlab("UMAP1") +
ylab("UMAP2")
DimPlot(seurat, reduction = "umap", group.by = "label", pt.size=0.5,
cols = colLab)+
theme_void()
cluster characterization
Idents(seurat) <- seurat$label
seurat_markers <- data.frame(gene=c("Fcgr2b","Fcer2a","Cr2","Cxcl13",
"Sulf1", "Thbs4", "Glycam1", "Ccl19",
"Ccl21a", "Fmod", "Grem1", "Bmp4",
"Hamp2", "Tnfsf11", "Stc1", "Fbn2", "Ptn",
"Lepr", "Ccl2", "Ccl7", "Cxcl10",
"Nr4a1","F3", "Fbln1","Ly6a", "Gdf10",
"Cd34","Ly6c1","Igfbp6", "Adamtsl1",
"Flrt2",
"Pi16", "Thy1", "Lsp1", "Npr1", "Mfap5",
"Fndc1", "Cnn1", "Myh1", "Myh2", "Flt1",
"Acta2", "Tinagl1", "Mcam", "Itga7", "Esam",
"Fabp4", "Nrarp", "Cox4i2", "Rgs4"
))
genes <- data.frame(geneID=rownames(seurat)) %>%
mutate(gene=gsub(".*\\.", "", geneID))
markerAll <- seurat_markers %>% left_join(., genes, by="gene")
ordVec <- names(colLab)
pOut <- avgHeatmap(seurat = seurat, selGenes = markerAll,
colVecIdent = colLab,
ordVec=ordVec,
gapVecR=NULL, gapVecC=NULL,cc=F,
cr=F, condCol=F)
vis selected PRC marker
genes <- data.frame(gene=rownames(seurat)) %>%
mutate(geneID=gsub("^.*\\.", "", gene))
selGenesAll <- data.frame(geneID=c( "Acta2", "Pi16", "Cd34", "Fbln1", "Gdf10",
"Adamtsl1", "Fndc1", "Dpp4", "F3", "Mfap5",
"Gsc", "Fgf10", "Col15a1", "Bmper",
"Sema3c", "Igfbp6", "Thy1"
)) %>%
left_join(., genes, by = "geneID")
pList <- sapply(selGenesAll$gene, function(x){
p <- FeaturePlot(seurat, reduction = "umap",
features = x,
cols=c("lightgrey", "darkred"),
order = F)+
theme(legend.position="right")
plot(p)
})
fraction EYFP+ cells
Barplot
eyfpCnt <- data.frame(table(seurat$label, seurat$EYFP)) %>%
spread(.,Var2 , Freq) %>% mutate(tot=pos+neg) %>%
mutate(freqPos=pos*100/tot) %>% mutate(freqNeg=neg*100/tot)
eyfpCntDat <- eyfpCnt %>% select(Var1, freqPos, freqNeg) %>%
gather(., eyfp, freq, freqPos:freqNeg)
p <- ggpubr::ggbarplot(eyfpCntDat, x="Var1", y="freq", fill="eyfp",
palette = c("#9d9f9e","#09983f"),
order = rev(names(colLab)),
xlab = "", ylab = "Frequency",
orientation = "horizontal") +
theme(legend.position = "right")
p
DotPlot
seurat$label <- factor(seurat$label, levels = names(colLab))
Idents(seurat) <- seurat$label
DotPlot(seurat, assay="RNA", features = "Rosa26eyfp.Rosa26eyfp", scale =F,
dot.min = 0, dot.scale = 4, scale.by = "size") +
scale_color_gradient(low="#bbd2c6", high="#017e40") +
coord_flip() +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
scale_x_discrete(breaks="Rosa26eyfp.Rosa26eyfp", labels="Rosa26eyfp") +
xlab("") + ylab("")
save seurat
saveRDS(seurat, file = paste0(basedir,
"/data/WT_adultOnly_iLNonly_filtered",
"_labeled_seurat.rds"))session info
sessionInfo()R version 4.3.0 (2023-04-21)
Platform: x86_64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] ggsci_3.0.1 scran_1.28.2
[3] scater_1.28.0 scuttle_1.10.3
[5] pheatmap_1.0.12 RColorBrewer_1.1-3
[7] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.2
[9] Biobase_2.60.0 GenomicRanges_1.52.1
[11] GenomeInfoDb_1.36.4 IRanges_2.36.0
[13] S4Vectors_0.40.1 BiocGenerics_0.48.0
[15] MatrixGenerics_1.12.3 matrixStats_1.2.0
[17] runSeurat3_0.1.0 here_1.0.1
[19] magrittr_2.0.3 Seurat_5.0.2
[21] SeuratObject_5.0.1 sp_2.1-3
[23] lubridate_1.9.3 forcats_1.0.0
[25] stringr_1.5.1 dplyr_1.1.4
[27] purrr_1.0.2 readr_2.1.5
[29] tidyr_1.3.1 tibble_3.2.1
[31] ggplot2_3.5.0 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] RcppAnnoy_0.0.22 splines_4.3.0
[3] later_1.3.2 bitops_1.0-7
[5] polyclip_1.10-6 fastDummies_1.7.3
[7] lifecycle_1.0.4 rstatix_0.7.2
[9] edgeR_3.42.4 rprojroot_2.0.4
[11] globals_0.16.2 lattice_0.22-5
[13] MASS_7.3-60.0.1 backports_1.4.1
[15] limma_3.56.2 plotly_4.10.4
[17] rmarkdown_2.26 yaml_2.3.8
[19] metapod_1.8.0 httpuv_1.6.14
[21] sctransform_0.4.1 spam_2.10-0
[23] spatstat.sparse_3.0-3 reticulate_1.35.0
[25] cowplot_1.1.3 pbapply_1.7-2
[27] abind_1.4-5 zlibbioc_1.46.0
[29] Rtsne_0.17 RCurl_1.98-1.14
[31] GenomeInfoDbData_1.2.10 ggrepel_0.9.5
[33] irlba_2.3.5.1 listenv_0.9.1
[35] spatstat.utils_3.0-4 goftest_1.2-3
[37] RSpectra_0.16-1 dqrng_0.3.2
[39] spatstat.random_3.2-3 fitdistrplus_1.1-11
[41] parallelly_1.37.1 DelayedMatrixStats_1.22.6
[43] leiden_0.4.3.1 codetools_0.2-19
[45] DelayedArray_0.26.7 tidyselect_1.2.0
[47] farver_2.1.1 viridis_0.6.5
[49] ScaledMatrix_1.8.1 spatstat.explore_3.2-6
[51] jsonlite_1.8.8 BiocNeighbors_1.18.0
[53] ellipsis_0.3.2 progressr_0.14.0
[55] ggridges_0.5.6 survival_3.5-8
[57] tools_4.3.0 ica_1.0-3
[59] Rcpp_1.0.12 glue_1.7.0
[61] gridExtra_2.3 xfun_0.42
[63] withr_3.0.0 fastmap_1.1.1
[65] bluster_1.10.0 fansi_1.0.6
[67] digest_0.6.34 rsvd_1.0.5
[69] timechange_0.3.0 R6_2.5.1
[71] mime_0.12 colorspace_2.1-0
[73] scattermore_1.2 tensor_1.5
[75] spatstat.data_3.0-4 utf8_1.2.4
[77] generics_0.1.3 data.table_1.15.2
[79] httr_1.4.7 htmlwidgets_1.6.4
[81] S4Arrays_1.0.6 uwot_0.1.16
[83] pkgconfig_2.0.3 gtable_0.3.4
[85] lmtest_0.9-40 XVector_0.40.0
[87] htmltools_0.5.7 carData_3.0-5
[89] dotCall64_1.1-1 scales_1.3.0
[91] png_0.1-8 knitr_1.45
[93] rstudioapi_0.15.0 tzdb_0.4.0
[95] reshape2_1.4.4 nlme_3.1-164
[97] zoo_1.8-12 KernSmooth_2.23-22
[99] vipor_0.4.7 parallel_4.3.0
[101] miniUI_0.1.1.1 pillar_1.9.0
[103] grid_4.3.0 vctrs_0.6.5
[105] RANN_2.6.1 ggpubr_0.6.0
[107] promises_1.2.1 car_3.1-2
[109] BiocSingular_1.16.0 beachmat_2.16.0
[111] xtable_1.8-4 cluster_2.1.6
[113] beeswarm_0.4.0 evaluate_0.23
[115] locfit_1.5-9.9 cli_3.6.2
[117] compiler_4.3.0 rlang_1.1.3
[119] crayon_1.5.2 ggsignif_0.6.4
[121] future.apply_1.11.1 labeling_0.4.3
[123] plyr_1.8.9 ggbeeswarm_0.7.2
[125] stringi_1.8.3 viridisLite_0.4.2
[127] deldir_2.0-4 BiocParallel_1.34.2
[129] munsell_0.5.0 lazyeval_0.2.2
[131] spatstat.geom_3.2-9 Matrix_1.6-5
[133] RcppHNSW_0.6.0 hms_1.1.3
[135] patchwork_1.2.0 sparseMatrixStats_1.12.2
[137] future_1.33.1 statmod_1.5.0
[139] shiny_1.8.0 ROCR_1.0-11
[141] broom_1.0.5 igraph_2.0.2 date()[1] "Wed Apr 3 13:57:54 2024"