Merge human LN imm cells

Author

Mechthild Lütge

Published

April 5, 2023

load packages

suppressPackageStartupMessages({
  library(tidyverse)
  library(Seurat)
  library(magrittr)
  library(dplyr)
  library(purrr)
  library(ggplot2)
  library(here)
  library(runSeurat3)
  library(ggsci)
})

set dir

basedir <- here()
metaDat <- read_tsv(paste0(basedir, "/data/metadataImm.txt"), col_names = T)

load and assign samples

assignSamples <- function(smpNam, basedirSmp, smpPat, smpGrp, smpOri,
                          smpSort){
  smpNamFull <- list.files(path = paste0(basedirSmp, "/data/immune/"),
                 pattern = paste0(smpNam, ".*_seurat.rds"))
  seuratSmp <- readRDS(paste0(basedirSmp, "/data/immune/", smpNamFull))
  seuratSmp$patient <- smpPat
  seuratSmp$cond <- smpGrp
  seuratSmp$origin <- smpOri
  seuratSmp$sort <- smpSort
  return(seuratSmp)
}


####################################################################

for(i in 1:length(metaDat$dataset)){
  seuratX <- assignSamples(smpNam = metaDat$dataset[i],
                           basedirSmp = basedir,
                           smpPat =  metaDat$patient[i],
                           smpGrp = metaDat$cond[i],
                           smpOri = metaDat$origin[i],
                           smpSort = metaDat$sort[i])
  if(exists("seurat")){
    seurat <- merge(x = seurat, y = seuratX, project = "LymphNode_IMM")
  }else{
    seurat <- seuratX
  }
}

remove(seuratX)

filter cells

dim(seurat)

[1]  36180 170115

seurat$cond2 <- seurat$cond
seurat$cond2[which(seurat$cond %in% c("chronic", "acute"))] <- "activated"

## reprocess
res <- c(0.8,0.6,0.4,0.25)

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)
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: 170115
Number of edges: 4993624

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

Number of nodes: 170115
Number of edges: 4993624

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

Number of nodes: 170115
Number of edges: 4993624

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

Number of nodes: 170115
Number of edges: 4993624

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9481
Number of communities: 14
Elapsed time: 75 seconds

## remove cluster with high expression of COL1A1, PDPN, LTF, KRTs, PROX1, ..
## PECAM1 + no PTPRC (stromal cells)

seurat <- subset(seurat, RNA_snn_res.0.25 %in% c("10", "11"), invert = T)

dim(seurat)

[1]  36180 168298

## reprocess
res <- c(0.8,0.6,0.4,0.25)

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)
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: 168298
Number of edges: 4986183

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

Number of nodes: 168298
Number of edges: 4986183

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

Number of nodes: 168298
Number of edges: 4986183

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

Number of nodes: 168298
Number of edges: 4986183

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9495
Number of communities: 18
Elapsed time: 72 seconds

visualize data

clustering

colPal <- c(pal_nejm()(8),pal_futurama()(12))[1:length(levels(seurat))]
names(colPal) <- levels(seurat)

## visualize input data
DimPlot(seurat, reduction = "umap", cols=colPal)+
  theme_bw() +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), 
        panel.grid.minor = element_blank()) +
  xlab("UMAP1") +
  ylab("UMAP2")

patient

colPat <- c(pal_jco()(10),pal_futurama()(12))[1:length(unique(seurat$patient))]
names(colPat) <- unique(seurat$patient)

## visualize input data
DimPlot(seurat, reduction = "umap", cols=colPat, group.by = "patient")+
  theme_bw() +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), 
        panel.grid.minor = element_blank()) +
  xlab("UMAP1") +
  ylab("UMAP2")

condition

colCond <- pal_igv()(length(unique(seurat$cond)))
names(colCond) <- unique(seurat$cond)

## visualize input data
DimPlot(seurat, reduction = "umap", cols=colCond, group.by = "cond")+
  theme_bw() +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), 
        panel.grid.minor = element_blank()) +
  xlab("UMAP1") +
  ylab("UMAP2")

origin

colOrig <- pal_aaas()(length(unique(seurat$origin)))
names(colCond) <- unique(seurat$cond)

## visualize input data
DimPlot(seurat, reduction = "umap", cols=colOrig, group.by = "origin")+
  theme_bw() +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), 
        panel.grid.minor = element_blank()) +
  xlab("UMAP1") +
  ylab("UMAP2")

save seurat object

saveRDS(seurat, file=paste0(basedir,
                            "/data/AllPatWithoutCM_IMMMerged_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] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggsci_3.0.1        runSeurat3_0.1.0   here_1.0.1         magrittr_2.0.3     Seurat_5.0.2      
 [6] SeuratObject_5.0.1 sp_2.1-3           lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1     
[11] dplyr_1.1.4        purrr_1.0.2        readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
[16] ggplot2_3.5.0      tidyverse_2.0.0   

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3          rstudioapi_0.15.0           jsonlite_1.8.8             
  [4] spatstat.utils_3.0-4        farver_2.1.1                rmarkdown_2.26             
  [7] zlibbioc_1.46.0             vctrs_0.6.5                 ROCR_1.0-11                
 [10] spatstat.explore_3.2-6      RCurl_1.98-1.14             S4Arrays_1.0.6             
 [13] htmltools_0.5.7             sctransform_0.4.1           parallelly_1.37.1          
 [16] KernSmooth_2.23-22          htmlwidgets_1.6.4           ica_1.0-3                  
 [19] plyr_1.8.9                  plotly_4.10.4               zoo_1.8-12                 
 [22] igraph_2.0.2                mime_0.12                   lifecycle_1.0.4            
 [25] pkgconfig_2.0.3             Matrix_1.6-5                R6_2.5.1                   
 [28] fastmap_1.1.1               GenomeInfoDbData_1.2.10     MatrixGenerics_1.12.3      
 [31] fitdistrplus_1.1-11         future_1.33.1               shiny_1.8.0                
 [34] digest_0.6.34               colorspace_2.1-0            S4Vectors_0.40.1           
 [37] patchwork_1.2.0             rprojroot_2.0.4             tensor_1.5                 
 [40] RSpectra_0.16-1             irlba_2.3.5.1               GenomicRanges_1.52.1       
 [43] labeling_0.4.3              progressr_0.14.0            fansi_1.0.6                
 [46] spatstat.sparse_3.0-3       timechange_0.3.0            httr_1.4.7                 
 [49] polyclip_1.10-6             abind_1.4-5                 compiler_4.3.0             
 [52] bit64_4.0.5                 withr_3.0.0                 fastDummies_1.7.3          
 [55] MASS_7.3-60.0.1             DelayedArray_0.26.7         tools_4.3.0                
 [58] lmtest_0.9-40               httpuv_1.6.14               future.apply_1.11.1        
 [61] goftest_1.2-3               glue_1.7.0                  nlme_3.1-164               
 [64] promises_1.2.1              grid_4.3.0                  Rtsne_0.17                 
 [67] cluster_2.1.6               reshape2_1.4.4              generics_0.1.3             
 [70] gtable_0.3.4                spatstat.data_3.0-4         tzdb_0.4.0                 
 [73] data.table_1.15.2           hms_1.1.3                   XVector_0.40.0             
 [76] utf8_1.2.4                  BiocGenerics_0.48.0         spatstat.geom_3.2-9        
 [79] RcppAnnoy_0.0.22            ggrepel_0.9.5               RANN_2.6.1                 
 [82] pillar_1.9.0                vroom_1.6.5                 spam_2.10-0                
 [85] RcppHNSW_0.6.0              later_1.3.2                 splines_4.3.0              
 [88] lattice_0.22-5              bit_4.0.5                   survival_3.5-8             
 [91] deldir_2.0-4                tidyselect_1.2.0            SingleCellExperiment_1.22.0
 [94] miniUI_0.1.1.1              pbapply_1.7-2               knitr_1.45                 
 [97] gridExtra_2.3               IRanges_2.36.0              SummarizedExperiment_1.30.2
[100] scattermore_1.2             stats4_4.3.0                xfun_0.42                  
[103] Biobase_2.60.0              matrixStats_1.2.0           pheatmap_1.0.12            
[106] stringi_1.8.3               lazyeval_0.2.2              yaml_2.3.8                 
[109] evaluate_0.23               codetools_0.2-19            cli_3.6.2                  
[112] uwot_0.1.16                 xtable_1.8-4                reticulate_1.35.0          
[115] munsell_0.5.0               GenomeInfoDb_1.36.4         Rcpp_1.0.12                
[118] globals_0.16.2              spatstat.random_3.2-3       png_0.1-8                  
[121] parallel_4.3.0              ellipsis_0.3.2              dotCall64_1.1-1            
[124] bitops_1.0-7                listenv_0.9.1               viridisLite_0.4.2          
[127] scales_1.3.0                ggridges_0.5.6              crayon_1.5.2               
[130] leiden_0.4.3.1              rlang_1.1.3                 cowplot_1.1.3

date()

[1] "Wed Mar 13 21:29:08 2024"