当比较不同细胞的基因表达时，特定基因的read数被用作表达的代理指标。然而，这一指标会受到技术差异的影响，比如样本之间的测序深度差异，或者部分样本中未被检测到的表达基因。因此，要充分比较细胞间的基因表达水平，首先需要将它们归一化为相对的基因表达水平。最流行的归一化方法是Log归一化，然而，最近开发的SCTransform方法显示出能够更好的提高下游分析的表现。
单细胞RNA-seq数据的生物学异质性经常受到包括测序深度在内的技术因素的干扰。即使是在相同的细胞类型中，每个细胞中检测到的分子数量在不同的细胞之间可能有显著差异。scRNA-seq数据的解释需要有效的预处理和标准化，以消除这种技术的可变性。这里我们介绍了一个能从scRNA-seq实验中标准化和稳定统计分子数的模型框架，这个过程省去了包括伪计数加法或对数转换在内的启发式步骤的需要，并改进了常见的下游分析任务，如可变基因选择、降维和差异表达。
SCTransform 使用正则化负二项回归算法计算了一个技术噪音的模型，这个模型的输出是标准化后的值，可以是正值，也可以是负值。正值表示观察到了更多的UMI，比预期的该基因在群体和细胞测序深度中的平均表达要多，而负值则相反。
library(Seurat)
library(ggplot2)
library(sctransform)
导入数据，创建Seurat对象：
pbmc_data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc_data)
标准化：
以前的Seurat入门教程中都是NormalizeData()，ScaleData()，FindVariableFeatures()，但是现在要使用SCTransform()来代替。
使用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)
新版的SCTransform还支持glmGamPoi来加速处理过程。可以通过method="glmGamPoi"指定。
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("glmGamPoi")
pbmc <- SCTransform(pbmc, method = "glmGamPoi", vars.to.regress = "percent.mt", verbose = FALSE)
降维和可视化都使用常规Seurat流程：
# 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()
还可以使用一行代码：
pbmc <- CreateSeuratObject(pbmc_data) %>%
PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
SCTransform(vars.to.regress = "percent.mt") %>%
RunPCA() %>% FindNeighbors(dims = 1:30) %>%
RunUMAP(dims = 1:30) %>%
FindClusters()
标准化数据的获取：
- SCTransform()结果存储在 “SCT” assay中。就可以使用
pbmc[["SCT"]]@scale.data获取数据。 但是这个矩阵不是稀疏矩阵，因此如果要存所有基因的信息就会占很大内存。为了节省内存，可以设置return.only.var.genes = TRUE来返回variable基因的值。这也是默认设置。 - 为了协助可视化和解释。我们还将皮尔逊值转换回“corrected”的UMI counts。可以将这些解释为我们希望观察到的UMI counts，如果所有细胞被测序到相同的深度。
- The ‘corrected’ UMI counts are stored in pbmc[["SCT"]]@counts。它是一个log-normalized的值，可以很方便的可视化。
- 可以使用corrected log-normalized counts做差异表达和数据整合。但原则上最好是直接对scale.data slot进行计算。
关于注释结果的差异：
sctransform 标准化后的结果与标准Seurat 流程的结果有以下几个差异：
- 分群分的更开了，基于CD8A，GZMK，CCL5这些marker显示分出了至少3个CD8T 细胞的群（naive，memory, e）
- 基于S100A4、CCR7、IL32和ISG15，明确分离三种CD4 T细胞群体(naive, memory, ifn activated)
- 基于TCL1A, FCER2的B细胞簇的额外发育亚结构
- 基于XCL1和FCGR3A, NK细胞进一步分离为CD56dim vs bright cluster.
# 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)
参考：
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1874-1
https://satijalab.org/seurat/articles/sctransform_vignette.html
