Seurat4.0系列教程11：使用sctransform
单细胞RNA-seq数据中的生物异质性经常受技术因素（包括测序深度）影响。每个细胞中检测到的分子数量在细胞之间可能显著变化，即使在相同的细胞类型内也是如此。对 scRNA-seq 数据的解释需要有效的预处理和标准化才能消除这种技术差异。
在此教程中，我们演示了使用基于标准化的sctransform ，能够得到log后较清晰的生物学差异。
library(Seurat)
library(ggplot2)
library(sctransform)
加载数据并创建seurat对象
pbmc_data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc_data)
应用sctransform标准化
请注意，此单个命令替换了 NormalizeData()，ScaleData()，FindVariableFeatures()sctransform后的数据默认存贮在 SCT assay 在标准化过程中，我们还可以删除干扰变异的因素，例如线粒体
# store mitochondrial percentage in object meta data
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")
# run sctransform
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)
最新版本还支持使用glmGamPoi包，这大大提高了运行的速度。它可以通过指定sctransformmethod=glmGamPoi来调用。
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("glmGamPoi")
pbmc <- SCTransform(pbmc, method = "glmGamPoi", vars.to.regress = "percent.mt", verbose = FALSE)
PCA 和 UMAP 降维
# These are now standard steps in the Seurat workflow for visualization and clustering
pbmc <- RunPCA(pbmc, verbose = FALSE)
pbmc <- RunUMAP(pbmc, dims = 1:30, verbose = FALSE)
pbmc <- FindNeighbors(pbmc, dims = 1:30, verbose = FALSE)
pbmc <- FindClusters(pbmc, verbose = FALSE)
DimPlot(pbmc, label = TRUE) + NoLegend()
用户可以根据传统marker单独注释亚群。与标准的 Seurat 工作流程相比，sctransform在几个方面揭示了更清晰的生物学区别，如：
基于CD8A、GZMK、CCL5、GZMK表达的至少3个CD8+T细胞群（naive, memory, effector）的分离 基于S100A4、CCR7、IL32和ISG15的三个CD4+T细胞群（naive, memory, IFN-activated）的清晰分离 基于 TCL1A、FCER2 的 B 细胞群中的其他亚型 基于 XCL1 和 FCGR3A 的 NK 细胞CD56dim vs. bright亚群的进一步分离
# These are now standard steps in the Seurat workflow for visualization and clustering Visualize
# canonical marker genes as violin plots.
VlnPlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7", "ISG15", "CD3D"),
pt.size = 0.2, ncol = 4)
# Visualize canonical marker genes on the sctransform embedding.
FeaturePlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7"), pt.size = 0.2,
ncol = 3)
FeaturePlot(pbmc, features = c("CD3D", "ISG15", "TCL1A", "FCER2", "XCL1", "FCGR3A"), pt.size = 0.2,
ncol = 3)