chracterize iLN all timepoints WT plus Ltbr floxed

Author

Mechthild Lütge

Published

July 1, 2023

load packages

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)
})

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/AllSamplesMerged_seurat.rds"))
seurat <- subset(seurat, location=="iLN")
table(seurat$EYFP, seurat$age)
     
         3w    8w E17to7wk   E18    P7
  neg   565  6041     2872 21202  1601
  pos  4398  8790      303   245   730
table(seurat$cond)

 LTbR    WT 
 3046 43701 
## rerunSeurat
seurat <- NormalizeData(object = seurat)
seurat <- FindVariableFeatures(object = seurat)
seurat <- ScaleData(object = seurat, verbose = FALSE)
seurat <- RunPCA(object = seurat, npcs = 30, verbose = FALSE)
seurat <- RunTSNE(object = seurat, reduction = "pca", dims = 1:20)
seurat <- RunUMAP(object = seurat, reduction = "pca", dims = 1:20)
seurat <- FindNeighbors(object = seurat, reduction = "pca", dims = 1:20)
res <- c(0.8,0.6,0.25,0.4)
for (i in 1:length(res)) {
  seurat <- FindClusters(object = seurat, resolution = res[i],
                         random.seed = 1234)
  }
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 46747
Number of edges: 1494535

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8911
Number of communities: 24
Elapsed time: 13 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 46747
Number of edges: 1494535

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9057
Number of communities: 19
Elapsed time: 12 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 46747
Number of edges: 1494535

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9411
Number of communities: 13
Elapsed time: 13 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 46747
Number of edges: 1494535

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9227
Number of communities: 14
Elapsed time: 12 seconds
colCond <- c("#446a7f", "#cb7457")
names(colCond) <- c("LTbR", "WT")

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",  "#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")

vis age

DimPlot(seurat, reduction = "umap", group.by = "age",
        cols = colAge)+
  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 = "age", pt.size=0.5,
        cols = colAge)+
  theme_void()

vis cond

DimPlot(seurat, reduction = "umap", group.by = "cond",
        cols = colCond, shuffle=T)+
  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 = "cond", pt.size=0.5,
        cols = colCond, shuffle=T)+
  theme_void()

vis age

DimPlot(seurat, reduction = "umap", group.by = "age", split.by = "EYFP",
        cols = colAge)+
  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 = "age", pt.size=0.5,
        cols = colAge, split.by = "EYFP",)+
  theme_void()

save seurat

saveRDS(seurat, file = paste0(basedir,
                              "/data/WT_allTime_iLNonly_WtplusLtbr",
                              "_seurat.rds"))

subset eyfp+ only

table(seurat$EYFP, seurat$age)
     
         3w    8w E17to7wk   E18    P7
  neg   565  6041     2872 21202  1601
  pos  4398  8790      303   245   730
seurat <- subset(seurat, EYFP=="pos")
table(seurat$cond, seurat$age)
      
         3w   8w E17to7wk  E18   P7
  LTbR    0  530      106    0    0
  WT   4398 8260      197  245  730
table(seurat$age)

      3w       8w E17to7wk      E18       P7 
    4398     8790      303      245      730 
## rerunSeurat
seurat <- NormalizeData(object = seurat)
seurat <- FindVariableFeatures(object = seurat)
seurat <- ScaleData(object = seurat, verbose = FALSE)
seurat <- RunPCA(object = seurat, npcs = 30, verbose = FALSE)
seurat <- RunTSNE(object = seurat, reduction = "pca", dims = 1:20)
seurat <- RunUMAP(object = seurat, reduction = "pca", dims = 1:20)
seurat <- FindNeighbors(object = seurat, reduction = "pca", dims = 1:20)
res <- c(0.8,0.6,0.25,0.4)
for (i in 1:length(res)) {
  seurat <- FindClusters(object = seurat, resolution = res[i],
                         random.seed = 1234)
  }
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 14466
Number of edges: 462776

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8248
Number of communities: 19
Elapsed time: 1 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 14466
Number of edges: 462776

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8431
Number of communities: 15
Elapsed time: 1 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 14466
Number of edges: 462776

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8959
Number of communities: 9
Elapsed time: 1 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 14466
Number of edges: 462776

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8699
Number of communities: 12
Elapsed time: 1 seconds

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")

vis age

DimPlot(seurat, reduction = "umap", group.by = "age",
        cols = colAge)+
  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 = "age", pt.size=0.5,
        cols = colAge)+
  theme_void()

vis cond

DimPlot(seurat, reduction = "umap", group.by = "cond",
        cols = colCond, shuffle=T)+
  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 = "cond", pt.size=0.5,
        cols = colCond, shuffle=T)+
  theme_void()

save seurat

saveRDS(seurat, file = paste0(basedir,
                              "/data/WT_allTime_iLNonly_WtplusLtbr_EYFPonly",
                              "_seurat.rds"))

transfer label adult

seuratLab <- readRDS(paste0(basedir,
                         "/data/WT_adultOnly_bothLabeled_integrated_",
                         "_seurat.rds"))
seuratLab <- subset(seuratLab, location=="iLN")
seuratLab <- subset(seuratLab, EYFP=="pos")
table(seuratLab$label)

  actMedRC    FDC/MRC      MedRC MedRC/IFRC    Pi16+RC        PRC       TBRC 
       156        113       1864       1079         10        505       1516 
       TRC       VSMC 
      2309         72 
labCells <- data.frame(label=seuratLab$label) %>% rownames_to_column(., "cell")
allCell <- data.frame(cell=colnames(seurat)) %>% 
  left_join(., labCells, by= "cell")
allCell$label[which(is.na(allCell$label))] <- "unassigned"
seurat$label <- allCell$label

colLab <- c("#42a071", "#900C3F","#b66e8d", "#61a4ba", "#424671", "#003C67FF",
            "#e3953d", "#714542", "#b6856e", "#a4a4a4")

names(colLab) <- c("FDC/MRC", "TRC", "TBRC", "MedRC/IFRC", "MedRC" , "actMedRC",
                   "PRC", "Pi16+RC", "VSMC", "unassigned")

DimPlot(seurat, reduction = "umap", group.by = "label",
        cols = colLab, shuffle=T)+
  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, shuffle=F)+
  theme_void()

save seurat

saveRDS(seurat, file = paste0(basedir,
                              "/data/WT_allTime_iLNonly_WtplusLtbr_EYFPonly_labelTrans",
                              "_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           edgeR_3.42.4             
  [9] rprojroot_2.0.4           globals_0.16.2           
 [11] lattice_0.22-5            MASS_7.3-60.0.1          
 [13] limma_3.56.2              plotly_4.10.4            
 [15] rmarkdown_2.26            yaml_2.3.8               
 [17] metapod_1.8.0             httpuv_1.6.14            
 [19] sctransform_0.4.1         spam_2.10-0              
 [21] spatstat.sparse_3.0-3     reticulate_1.35.0        
 [23] cowplot_1.1.3             pbapply_1.7-2            
 [25] abind_1.4-5               zlibbioc_1.46.0          
 [27] Rtsne_0.17                RCurl_1.98-1.14          
 [29] GenomeInfoDbData_1.2.10   ggrepel_0.9.5            
 [31] irlba_2.3.5.1             listenv_0.9.1            
 [33] spatstat.utils_3.0-4      goftest_1.2-3            
 [35] RSpectra_0.16-1           dqrng_0.3.2              
 [37] spatstat.random_3.2-3     fitdistrplus_1.1-11      
 [39] parallelly_1.37.1         DelayedMatrixStats_1.22.6
 [41] leiden_0.4.3.1            codetools_0.2-19         
 [43] DelayedArray_0.26.7       tidyselect_1.2.0         
 [45] farver_2.1.1              viridis_0.6.5            
 [47] ScaledMatrix_1.8.1        spatstat.explore_3.2-6   
 [49] jsonlite_1.8.8            BiocNeighbors_1.18.0     
 [51] ellipsis_0.3.2            progressr_0.14.0         
 [53] ggridges_0.5.6            survival_3.5-8           
 [55] tools_4.3.0               ica_1.0-3                
 [57] Rcpp_1.0.12               glue_1.7.0               
 [59] gridExtra_2.3             xfun_0.42                
 [61] withr_3.0.0               fastmap_1.1.1            
 [63] bluster_1.10.0            fansi_1.0.6              
 [65] digest_0.6.34             rsvd_1.0.5               
 [67] timechange_0.3.0          R6_2.5.1                 
 [69] mime_0.12                 colorspace_2.1-0         
 [71] scattermore_1.2           tensor_1.5               
 [73] spatstat.data_3.0-4       utf8_1.2.4               
 [75] generics_0.1.3            data.table_1.15.2        
 [77] httr_1.4.7                htmlwidgets_1.6.4        
 [79] S4Arrays_1.0.6            uwot_0.1.16              
 [81] pkgconfig_2.0.3           gtable_0.3.4             
 [83] lmtest_0.9-40             XVector_0.40.0           
 [85] htmltools_0.5.7           dotCall64_1.1-1          
 [87] scales_1.3.0              png_0.1-8                
 [89] knitr_1.45                rstudioapi_0.15.0        
 [91] tzdb_0.4.0                reshape2_1.4.4           
 [93] nlme_3.1-164              zoo_1.8-12               
 [95] KernSmooth_2.23-22        vipor_0.4.7              
 [97] parallel_4.3.0            miniUI_0.1.1.1           
 [99] pillar_1.9.0              grid_4.3.0               
[101] vctrs_0.6.5               RANN_2.6.1               
[103] promises_1.2.1            BiocSingular_1.16.0      
[105] beachmat_2.16.0           xtable_1.8-4             
[107] cluster_2.1.6             beeswarm_0.4.0           
[109] evaluate_0.23             locfit_1.5-9.9           
[111] cli_3.6.2                 compiler_4.3.0           
[113] rlang_1.1.3               crayon_1.5.2             
[115] future.apply_1.11.1       labeling_0.4.3           
[117] plyr_1.8.9                ggbeeswarm_0.7.2         
[119] stringi_1.8.3             viridisLite_0.4.2        
[121] deldir_2.0-4              BiocParallel_1.34.2      
[123] munsell_0.5.0             lazyeval_0.2.2           
[125] spatstat.geom_3.2-9       Matrix_1.6-5             
[127] RcppHNSW_0.6.0            hms_1.1.3                
[129] patchwork_1.2.0           sparseMatrixStats_1.12.2 
[131] future_1.33.1             statmod_1.5.0            
[133] shiny_1.8.0               ROCR_1.0-11              
[135] igraph_2.0.2             
date()
[1] "Wed Apr 24 10:03:47 2024"