Test Seurat on SciServer

Summary

• Based on https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

Load packages

Location

.libPaths( "/home/idies/workspace/c_moor_data/R/4.0.3/scRNA-seq/" )
.libPaths()

[1] "/home/idies/workspace/c_moor_data/R/4.0.3/scRNA-seq"
[2] "/home/idies/R/lib64/R/library"                      

Packages

# reticulate::use_virtualenv( "/home/idies/workspace/c_moor_data/R/4.0.3/Seurat/4.0.1/env" )
# reticulate::py_config()

library( "dplyr" )

Prepare data

Load counts

pbmc.data <- Read10X( "/home/idies/workspace/c_moor_data/pbmc3k/filtered_gene_bc_matrices/hg19/" )
pbmc <- CreateSeuratObject( pbmc.data, min.cells = 3, min.features = 200 )

Feature names cannot have underscores ('_'), replacing with dashes ('-')

pbmc

An object of class Seurat 
13714 features across 2700 samples within 1 assay 
Active assay: RNA (13714 features, 0 variable features)

Normalize

pbmc <- NormalizeData( pbmc )

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Find variable genes

pbmc <- FindVariableFeatures( pbmc )

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

top10 <- head( VariableFeatures(pbmc), 10 )
plot1 <- VariableFeaturePlot( pbmc )
LabelPoints( plot1, points = top10, repel=TRUE )

When using repel, set xnudge and ynudge to 0 for optimal results

Scale

all.genes <- rownames( pbmc )
pbmc <- ScaleData( pbmc, features=all.genes, verbose=FALSE )

PCA

pbmc <- RunPCA( pbmc, verbose=FALSE )

VizDimLoadings( pbmc, dims=1:2 )

Cluster

pbmc <- FindNeighbors( pbmc, dims=1:10 )

Computing nearest neighbor graph

Computing SNN

pbmc <- FindClusters( pbmc, resolution=0.5 )

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 2700
Number of edges: 97892

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8719
Number of communities: 9
Elapsed time: 0 seconds

UMAP

pbmc <- RunUMAP( pbmc, dims=1:10 )

The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session21:31:41 UMAP embedding parameters a = 0.9922 b = 1.112
21:31:41 Read 2700 rows and found 10 numeric columns
21:31:41 Using Annoy for neighbor search, n_neighbors = 30
21:31:41 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
21:31:42 Writing NN index file to temp file /tmp/RtmprNLYRG/file14275bf786a
21:31:42 Searching Annoy index using 1 thread, search_k = 3000
21:31:44 Annoy recall = 100%
21:31:45 Commencing smooth kNN distance calibration using 1 thread
21:31:47 Initializing from normalized Laplacian + noise
21:31:47 Commencing optimization for 500 epochs, with 107868 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
21:31:54 Optimization finished

Visualize cells

DimPlot( pbmc, reduction="umap" )

Find markers

Detect genes

pbmc.markers <- FindAllMarkers( 
  pbmc, only.pos=TRUE, 
  min.pct=0.25, logfc.threshold=0.25,
  verbose=FALSE
)

Summarize

pbmc.markers %>% 
  group_by( cluster ) %>% 
  top_n( n=2, wt=avg_log2FC
)

Violin plot

VlnPlot( pbmc, features=c( "MS4A1", "CD79A" ) )

Feature plot

FeaturePlot( pbmc, features=c( "MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A" ) )

Label clusters

new.cluster.ids <- c( "Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet" )
names( new.cluster.ids ) <- levels( pbmc )
pbmc <- RenameIdents( pbmc, new.cluster.ids )
DimPlot( pbmc, reduction="umap", label=TRUE, pt.size=0.5 ) + NoLegend()

Document software

sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 8

Matrix products: default
BLAS:   /home/idies/R/lib64/R/lib/libRblas.so
LAPACK: /home/idies/R/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] dplyr_1.0.5        SeuratObject_4.0.0 Seurat_4.0.1      

loaded via a namespace (and not attached):
  [1] Rtsne_0.15            colorspace_2.0-0      deldir_0.2-10        
  [4] ellipsis_0.3.1        ggridges_0.5.3        rstudioapi_0.13      
  [7] spatstat.data_2.1-0   leiden_0.3.7          listenv_0.8.0        
 [10] farver_2.1.0          ggrepel_0.9.1         RSpectra_0.16-0      
 [13] fansi_0.4.2           codetools_0.2-16      splines_4.0.3        
 [16] knitr_1.31            polyclip_1.10-0       jsonlite_1.7.2       
 [19] ica_1.0-2             cluster_2.1.0         png_0.1-7            
 [22] uwot_0.1.10           shiny_1.6.0           sctransform_0.3.2    
 [25] spatstat.sparse_2.0-0 compiler_4.0.3        httr_1.4.2           
 [28] assertthat_0.2.1      Matrix_1.2-18         fastmap_1.1.0        
 [31] lazyeval_0.2.2        limma_3.46.0          later_1.1.0.1        
 [34] htmltools_0.5.1.1     tools_4.0.3           igraph_1.2.6         
 [37] gtable_0.3.0          glue_1.4.2            RANN_2.6.1           
 [40] reshape2_1.4.4        Rcpp_1.0.6            scattermore_0.7      
 [43] jquerylib_0.1.3       vctrs_0.3.6           nlme_3.1-149         
 [46] lmtest_0.9-38         xfun_0.22             stringr_1.4.0        
 [49] globals_0.14.0        mime_0.10             miniUI_0.1.1.1       
 [52] lifecycle_1.0.0       irlba_2.3.3           goftest_1.2-2        
 [55] future_1.21.0         MASS_7.3-53           zoo_1.8-9            
 [58] scales_1.1.1          spatstat.core_2.1-2   promises_1.2.0.1     
 [61] spatstat.utils_2.1-0  parallel_4.0.3        RColorBrewer_1.1-2   
 [64] yaml_2.2.1            reticulate_1.19       pbapply_1.4-3        
 [67] gridExtra_2.3         ggplot2_3.3.3         sass_0.3.1           
 [70] rpart_4.1-15          stringi_1.5.3         rlang_0.4.10         
 [73] pkgconfig_2.0.3       matrixStats_0.58.0    evaluate_0.14        
 [76] lattice_0.20-41       ROCR_1.0-11           purrr_0.3.4          
 [79] tensor_1.5            patchwork_1.1.1       htmlwidgets_1.5.3    
 [82] labeling_0.4.2        cowplot_1.1.1         tidyselect_1.1.0     
 [85] parallelly_1.24.0     RcppAnnoy_0.0.18      plyr_1.8.6           
 [88] magrittr_2.0.1        R6_2.5.0              generics_0.1.0       
 [91] DBI_1.1.1             withr_2.4.1           pillar_1.5.1         
 [94] mgcv_1.8-33           fitdistrplus_1.1-3    survival_3.2-7       
 [97] abind_1.4-5           tibble_3.1.0          future.apply_1.7.0   
[100] crayon_1.4.1          KernSmooth_2.23-17    utf8_1.1.4           
[103] spatstat.geom_2.1-0   plotly_4.9.3          rmarkdown_2.7        
[106] grid_4.0.3            data.table_1.14.0     digest_0.6.27        
[109] xtable_1.8-4          tidyr_1.1.3           httpuv_1.5.5         
[112] munsell_0.5.0         viridisLite_0.4.0     bslib_0.2.4