Tutorial: scRNA-seq data
Resolution of curse of dimensionality (RECODE) is a noise reduction method for single-cell sequencing data. This vignette briefly demonstrates how to run RECODE on single-cell RNA-seq data and integrate the result into Seurat analysis pipeline.
1 Set up
This tutorial utilizes the screcode Python package via the reticulate R package, and requires Python, the screcode Python package, and the reticulate R package installed in advance. You can install the screcode into your preferred Python environment using pip:
# pip install screcodeYou can also install the R packages as follows.
# install.packages("reticulate") # For interfacing with Python
# install.packages("Seurat") # For the downstream analysis
# install.packages("Matrix") # For processing sparse matrices
# install.packages("ggplot2") # For data visualization1.1 Import libraries
Though we use Seurat v4 in this tutorial, we have checked Seurat v3 also works.
library(Seurat)
library(Matrix)
library(ggplot2)
library(reticulate)
set.seed(0)You can select the Python environment to work in using the use_python function in the reticulate R package.
reticulate::use_python("C:/Users/ashbi/anaconda3/envs/r", # Path to python.exe
required=TRUE)
py_config() # Confirm which Python to be usedtheme_set(theme(text = element_text(size = 18),
panel.background = element_rect(fill = "white", color = "gray20"),
panel.grid = element_line(color = "gray92"),
legend.key = element_rect(colour = NA, fill = NA)))1.2 Import data
Sample data： For this tutorial, we use sample 10k Human PBMCs, 3’ v3.1, Chromium Controller (11,485 cells and 36,601 genes) in 10X Genomics Datasets. The test data is deposited as Feature / cell matrix HDF5 (filterd) here (registration required).
seurat <- Read10X_h5("10k_PBMC_3p_nextgem_Chromium_Controller_filtered_feature_bc_matrix.h5")
seurat <- CreateSeuratObject(seurat)
seurat[["percent.mt"]] <- PercentageFeatureSet(seurat, pattern = "^MT-")2 Pre-processing
2.1 Run RECODE
We apply RECODE on a count matrix where each row is a gene and each column is a cell (gene x cell). The output is a denoised count matrix (gene x cell).
plt <- reticulate::import(module="matplotlib.pyplot")
screcode <- reticulate::import(module="screcode.screcode")
recode<-screcode$RECODE()
data<-t(as.matrix(seurat[["RNA"]]@counts))
data_RECODE <- recode$fit_transform(data)
rownames(data_RECODE)<-colnames(seurat)
colnames(data_RECODE)<-rownames(seurat)recode$report()
plt$show()
We store the denoised count matrix in seurat[["RECODE"]]@counts.
seurat[["RECODE"]] <- CreateAssayObject(counts = Matrix(t(data_RECODE), sparse = TRUE))
DefaultAssay(seurat) <- "RECODE"2.2 Quality control
We remove low-quality cells using the following thresholds.
700 \(\leq\) No. of detected genes \(\leq\) 6000
% Mitochondrial transcript \(\leq\) 15
1000 \(\leq\) Total UMI \(\leq\) 25000
ngene_low <- 700
ngene_high <- 6000
mito_high <- 15
numi_low <- 1000
numi_high <- 25000meta.data <- seurat@meta.data
meta.data$quality <- rep("High", nrow(meta.data))
meta.data$quality[meta.data$nFeature_RNA<ngene_low|ngene_high<meta.data$nFeature_RNA] <- "Low"
meta.data$quality[meta.data$nCount_RNA<numi_low|numi_high<meta.data$nCount_RNA] <- "Low"
meta.data$quality[meta.data$percent.mt>mito_high] <- "Low"
meta.data$quality <- factor(meta.data$quality)p1 <- ggplot(meta.data, aes(x = nFeature_RNA, y = percent.mt, color = quality))+
geom_point(size = 0.1, show.legend = F)+
scale_color_manual(values = c("High"="black", "Low"="red"))+
xlab("nGene")
p2 <- ggplot(meta.data, aes(x = nCount_RNA, y = nFeature_RNA, color = quality))+
geom_point(size = 0.1)+
scale_color_manual(values = c("High"="black", "Low"="red"))+
xlab("nUMI")+
ylab("nGene")+
guides(color = guide_legend(override.aes = list(size = 5)))
p1+p2
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
ncol = 3,
pt.size = 0)
seurat <- subset(seurat, subset = nFeature_RNA > ngene_low & nFeature_RNA < ngene_high & percent.mt < mito_high & nCount_RNA > numi_low & nCount_RNA < numi_high)
dim(seurat)## [1] 36601 110052.3 Normalization, centering, and PCA
Then, we perform library size normalization and run PCA on 2000 highly variable genes. To retain as much biological information as possible, We do not scale gene expression variance to 1 in this tutorial.
seurat <- NormalizeData(seurat,scale.factor=10^5)
seurat <- FindVariableFeatures(seurat, selection.method = "vst", nfeatures = 2000)
seurat <- ScaleData(seurat, features = rownames(seurat), do.scale = FALSE, do.center = TRUE)## Centering data matrixseurat <- RunPCA(seurat, features = VariableFeatures(seurat))## PC_ 1
## Positive: LYZ, S100A9, S100A8, CST3, IFI30, SERPINA1, FOS, VCAN, AIF1, LST1
## HLA-DRA, AC020656.1, MARCKS, FCER1G, S100A12, WARS, CFD, CPVL, MAFB, HLA-DRB1
## IFITM3, SAT1, TNFSF10, PLBD1, LILRB2, LGALS2, HLA-DPA1, LRRC25, AC020916.1, CD74
## Negative: CD247, GZMM, CTSW, KLRK1, CD8B, CCL5, CD8A, GZMA, SAMD3, GRAP2
## NKG7, STMN1, C12orf75, HOPX, KLRB1, CST7, VEGFB, LINC02446, PYHIN1, NCR3
## GNLY, MATK, PRF1, KLRG1, LINC01871, LAIR2, KLRD1, PLEKHF1, TRGC2, MYBL1
## PC_ 2
## Positive: GZMM, CD247, CTSW, NKG7, GZMA, CCL5, KLRK1, S100A8, CST7, S100A9
## GNLY, PRF1, SAMD3, KLRD1, CD8A, FCER1G, AIF1, CD8B, KLRG1, GBP1
## S100A12, FOS, LYZ, LAIR2, GZMH, TRDC, VCAN, KLRB1, LINC01871, MATK
## Negative: IGHM, HLA-DRA, MS4A1, IGKC, CD79A, HLA-DQA1, CD74, HLA-DRB1, BANK1, CD79B
## HLA-DPB1, HLA-DPA1, IGHD, RALGPS2, LINC00926, TNFRSF13C, BCL11A, HLA-DRB5, VPREB3, CD22
## AFF3, LINC02397, SPIB, NIBAN3, HLA-DQB1, FCER2, TCL1A, PAX5, BLK, TCF4
## PC_ 3
## Positive: NKG7, GNLY, GZMA, CCL5, CST7, PRF1, KLRD1, FCGR3A, HOPX, GZMB
## GZMH, FGFBP2, KLRG1, KLRF1, TRDC, RHOC, SPON2, CLIC3, ABI3, CCL4
## CD160, PTGDR, HLA-DPA1, ADGRG1, HLA-DPB1, FCER1G, MATK, SLAMF7, LAIR2, XCL2
## Negative: BEX3, SLC40A1, LINC02446, PADI4, CD8B, CKAP4, F5, APP, S100A12, S100A8
## VCAN, HEMGN, GRAP2, FLNB, S100A9, RETN, CYP1B1, AIF1, RTKN2, QPCT
## FHL1, MCEMP1, S100B, CRISPLD2, AEBP1, LAPTM4B, MARC1, NET1, STRBP, MGST1
## PC_ 4
## Positive: FCGR3A, IFITM3, CDKN1C, HES4, HLA-DQA1, CTSL, TCF7L2, BATF3, MS4A7, RHOC
## CSF1R, RRAS, ABI3, MTSS1, APOBEC3A, CASP5, HLA-DPA1, AIF1, SIGLEC10, HLA-DPB1
## PTP4A3, FAM110A, CALHM6, WARS, ISG15, TUBA1B, NEURL1, CKS1B, GBP1, ADA
## Negative: S100A8, S100A9, S100A12, VCAN, FOS, CYP1B1, ALOX5AP, LYZ, AC020916.1, THBS1
## AC020656.1, CES1, PADI4, MEGF9, RNASE6, MGST1, MCEMP1, QPCT, CKAP4, CLEC4E
## RNASE2, STEAP4, PLBD1, F5, CRISPLD2, RBP7, MARC1, CD163, ARHGAP24, S1PR3
## PC_ 5
## Positive: CD79B, CD79A, MS4A1, TNFSF10, LINC00926, AIM2, MARCKS, MTSS1, IFIT3, VPREB3
## TNFRSF13C, IFIT2, FCRL3, IGHM, BANK1, MS4A7, IGHD, FCGR3A, IFIT1, RALGPS2
## SWAP70, FCRL1, CD40, C3AR1, PAX5, CCL5, ISG15, FCER2, GBP1, IFITM3
## Negative: GZMB, PLD4, LILRA4, ITM2C, SERPINF1, PPP1R14B, CLEC4C, TPM2, CCDC50, CST3
## C12orf75, SCT, UGCG, IL3RA, CLIC3, FCER1A, RUNX2, SMPD3, DAB2, FABP5
## LRRC26, MZB1, GAS6, DNASE1L3, APP, PTCRA, LAMP5, FCER1G, GSN, LILRB43 Analysis
3.1 Clustering
seurat <- FindNeighbors(seurat, dims = 1:50)## Computing nearest neighbor graph## Computing SNNseurat <- FindClusters(seurat, resolution = 1.2)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11005
## Number of edges: 383227
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8413
## Number of communities: 19
## Elapsed time: 1 seconds3.2 Run UMAP
seurat <- RunUMAP(seurat, dims = 1:50, seed.use = 42)## 16:53:34 UMAP embedding parameters a = 0.9922 b = 1.112## 16:53:34 Read 11005 rows and found 50 numeric columns## 16:53:34 Using Annoy for neighbor search, n_neighbors = 30## 16:53:34 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 16:53:35 Writing NN index file to temp file C:\Users\ashbi\AppData\Local\Temp\RtmpgBbf3N\filea5fc34a6136e
## 16:53:35 Searching Annoy index using 1 thread, search_k = 3000
## 16:53:37 Annoy recall = 100%
## 16:53:37 Commencing smooth kNN distance calibration using 1 thread
## 16:53:38 Initializing from normalized Laplacian + noise
## 16:53:38 Commencing optimization for 200 epochs, with 489256 positive edges
## 16:53:48 Optimization finishedDimPlot(seurat, reduction = "umap", label = TRUE, label.size = 5)
3.3 Marker gene expression
DefaultAssay(seurat) <- "RECODE"
FeaturePlot(seurat, features = c("CD3D", "CD4", "CD8B", "CCR7", "GZMK", "PRF1", "KLRB1", "TRDC", "NKG7", "CD14", "FCGR3A", "CD19", "TCL1A", "FCER1A", "LILRA4"), ncol = 4)
3.4 Annotation
use.switch <- function(x)
{
sapply(x, function(y){
switch(y,
"2" = "Naive CD8+ T cells",
"1" = "Naive CD4+ T cells",
"3" = "Memory CD4+ T cells",
"16" = "Memory CD4+ T cells",
"11" = "Central memory CD8+ T cells",
"7" = "Effector memory CD8+ T cells",
"14" = "gd T cells",
"10" = "Natural Killer cells",
"6" = "Memory B cells",
"9" = "Naive B cells",
"0" = "CD14+ monocytes",
"4" = "CD14+ monocytes",
"5" = "CD14+ monocytes",
"8" = "FCGR3A+ monocytes",
"13" = "Myeloid dendritic cells",
"15" = "Plasmacytoid dendritic cells",
"18" = "Plasmacytoid dendritic cells",
"12" = "Doublets",
"17" = "Doublets",
y
)
})
}
seurat$cluster <- as.character(seurat$seurat_clusters)
seurat$cluster <- use.switch(seurat$cluster)
seurat$cluster <- factor(seurat$cluster, levels = c("Naive CD8+ T cells", "Naive CD4+ T cells", "Memory CD4+ T cells", "Central memory CD8+ T cells", "Effector memory CD8+ T cells", "gd T cells", "Natural Killer cells", "Memory B cells", "Naive B cells", "CD14+ monocytes", "FCGR3A+ monocytes", "Myeloid dendritic cells", "Plasmacytoid dendritic cells", "Doublets"))DimPlot(seurat, group.by = "cluster")
3.5 Check noise reduction effect of RECODE
seurat <- NormalizeData(seurat, assay = "RNA",scale.factor=10^5)genes <- c("CD14", "FCGR3A")
# Check gene expression value before RECODE
dat <- data.frame(A = seurat[["RNA"]]@data[genes[1], ],
B = seurat[["RNA"]]@data[genes[2], ],
cluster = seurat$cluster)
dat <- dat[(grep("monocytes", dat$cluster)), ]
p1 <- ggplot(dat, aes(x = A, y = B, color = cluster))+
geom_point(size = 0.7)+
xlab(genes[1])+
ylab(genes[2])+
guides(color = guide_legend(override.aes = list(size = 5)))+
ggtitle("Original")
# Check gene expression value after RECODE
dat <- data.frame(A = seurat[["RECODE"]]@data[genes[1], ],
B = seurat[["RECODE"]]@data[genes[2], ],
cluster = seurat$cluster)
dat <- dat[(grep("monocytes", dat$cluster)), ]
p2 <- ggplot(dat, aes(x = A, y = B, color = cluster))+
geom_point(size = 0.7)+
xlab(genes[1])+
ylab(genes[2])+
guides(color = guide_legend(override.aes = list(size = 5)))+
ggtitle("RECODE")
p1/p2
The distribution of two monocyte populations (CD14+ cells and FCGR3A+ cells) is clearly captured with RECODE.
genes <- c("CCR7", "S100A4")
# Check gene expression value before RECODE
dat <- data.frame(A = seurat[["RNA"]]@data[genes[1], ],
B = seurat[["RNA"]]@data[genes[2], ],
cluster = seurat$cluster)
dat <- dat[(grep("CD8", dat$cluster)), ]
p1 <- ggplot(dat, aes(x = A, y = B, color = cluster))+
geom_point(size = 0.7)+
xlab(genes[1])+
ylab(genes[2])+
guides(color = guide_legend(override.aes = list(size = 5)))+
ggtitle("Original")
# Check gene expression value after RECODE
dat <- data.frame(A = seurat[["RECODE"]]@data[genes[1], ],
B = seurat[["RECODE"]]@data[genes[2], ],
cluster = seurat$cluster)
dat <- dat[(grep("CD8", dat$cluster)), ]
p2 <- ggplot(dat, aes(x = A, y = B, color = cluster))+
geom_point(size = 0.7)+
xlab(genes[1])+
ylab(genes[2])+
guides(color = guide_legend(override.aes = list(size = 5)))+
ggtitle("RECODE")
p1/p2
This data includes three CD8+ T cell subsets: naive cells (CCR7+, S100A4-), central memory cells (CCR7+, S100A4+), and effector memory cells (CCR7-, S100A4+). This expression pattern is visible with RECODE.
While the original expression values are sparse and discontinuous due to random sampling effects, RECODE can reconstruct the continuous denoised data, which more accurately reflects the biological information. These two simple examples demonstrate the successful noise reduction by RECODE.
4 Session information
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Japanese_Japan.932 LC_CTYPE=Japanese_Japan.932 LC_MONETARY=Japanese_Japan.932 LC_NUMERIC=C
## [5] LC_TIME=Japanese_Japan.932
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rmarkdown_2.8 reticulate_1.16 ggplot2_3.3.5 Matrix_1.3-4 SeuratObject_4.0.2 Seurat_4.0.4 RevoUtils_11.0.2
## [8] RevoUtilsMath_11.0.0
##
## loaded via a namespace (and not attached):
## [1] plyr_1.8.6 igraph_1.2.5 lazyeval_0.2.2 splines_4.0.2 listenv_0.8.0 scattermore_0.7
## [7] usethis_2.0.1 digest_0.6.25 htmltools_0.5.2 fansi_0.4.1 magrittr_2.0.1 memoise_2.0.0
## [13] tensor_1.5 cluster_2.1.0 ROCR_1.0-11 remotes_2.4.0 globals_0.14.0 matrixStats_0.61.0
## [19] spatstat.sparse_2.0-0 prettyunits_1.1.1 colorspace_1.4-1 ggrepel_0.8.2 xfun_0.23 dplyr_1.0.7
## [25] callr_3.7.0 crayon_1.4.1 jsonlite_1.7.0 spatstat.data_2.1-0 survival_3.1-12 zoo_1.8-8
## [31] glue_1.4.1 polyclip_1.10-0 gtable_0.3.0 leiden_0.3.9 pkgbuild_1.2.0 future.apply_1.8.1
## [37] abind_1.4-5 scales_1.1.1 DBI_1.1.1 miniUI_0.1.1.1 Rcpp_1.0.7 viridisLite_0.3.0
## [43] xtable_1.8-4 spatstat.core_2.3-0 bit_4.0.4 htmlwidgets_1.5.4 httr_1.4.2 RColorBrewer_1.1-2
## [49] ellipsis_0.3.2 ica_1.0-2 farver_2.0.3 pkgconfig_2.0.3 sass_0.4.0 uwot_0.1.10
## [55] deldir_0.2-10 utf8_1.2.2 labeling_0.4.2 tidyselect_1.1.1 rlang_0.4.11 reshape2_1.4.4
## [61] later_1.3.0 munsell_0.5.0 tools_4.0.2 cachem_1.0.6 cli_3.0.1 generics_0.1.0
## [67] devtools_2.4.2 ggridges_0.5.3 evaluate_0.14 stringr_1.4.0 fastmap_1.1.0 yaml_2.2.1
## [73] goftest_1.2-2 processx_3.5.2 knitr_1.34 bit64_4.0.5 fs_1.5.0 fitdistrplus_1.1-5
## [79] purrr_0.3.4 RANN_2.6.1 pbapply_1.5-0 future_1.22.1 nlme_3.1-148 mime_0.11
## [85] hdf5r_1.3.4 compiler_4.0.2 plotly_4.9.4.1 png_0.1-7 testthat_3.0.4 spatstat.utils_2.2-0
## [91] tibble_3.0.3 bslib_0.3.0 stringi_1.4.6 highr_0.9 ps_1.6.0 RSpectra_0.16-0
## [97] desc_1.3.0 lattice_0.20-41 vctrs_0.3.8 pillar_1.6.2 lifecycle_1.0.0 spatstat.geom_2.2-2
## [103] lmtest_0.9-37 jquerylib_0.1.4 RcppAnnoy_0.0.18 data.table_1.12.8 cowplot_1.1.1 irlba_2.3.3
## [109] httpuv_1.6.3 patchwork_1.0.1 R6_2.4.1 promises_1.2.0.1 KernSmooth_2.23-17 gridExtra_2.3
## [115] parallelly_1.28.1 sessioninfo_1.1.1 codetools_0.2-16 MASS_7.3-51.6 assertthat_0.2.1 pkgload_1.2.2
## [121] rprojroot_2.0.2 withr_2.4.2 sctransform_0.3.2 mgcv_1.8-31 parallel_4.0.2 grid_4.0.2
## [127] rpart_4.1-15 tidyr_1.1.0 Rtsne_0.15 shiny_1.5.0