Seurat从3.0版本引进了SCTransform这个函数用来对数据做标准化，并且这一个函数可以代替三个函数（NormalizeData, ScaleData, FindVariableFeatures）的运行。 且其对测序深度的校正效果要好于log标准化。（10万以内的细胞都建议使用SCT标准化）
️SCTransform对测序深度的校正效果很好，也可用于矫正线粒体等因素的影响，但不能用于批次矫正。
官网：
https://cran.r-project.org/web/packages/sctransform/index.html
https://satijalab.org/seurat/archive/v3.0/sctransform_vignette.html
官网对这个函数的描述：A normalization method for single-cell UMI count data using a variance stabilizing transformation. The transformation is based on a negative binomial regression model with regularized parameters. As part of the same regression framework, this package also provides functions for batch correction, and data correction. See Hafemeister and Satija 2019 <doi:10.1186/s13059-019-1874-1> for more details. （这是一个用方差稳定变换对单细胞UMI count 数据标准化的方法，方差稳定变换是基于负二项回归。这个函数在对数据进行均一化的同时还可以去除线粒体红细胞等混杂因素的影响。）
1. NormalizeData, ScaleData, FindVariableFeatures和SCTransform结果的存放位置
示例数据下载：pbmc3k
#读入数据，创建seurat对象
pbmc <- Read10X('./filtered_gene_bc_matrices/hg19/')
pbmc <- CreateSeuratObject(pbmc,project = 'pbmc3k',min.cells = 3,min.features = 200)
## 质控
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
常规NormalizeData, ScaleData, FindVariableFeatures操作
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
View(pbmc)
这一步的结果都储存在pbmc@assays[["RNA"]]中。
再使用SCTransform进行标准化
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)
View(pbmc)
再次查看pbmc的数据结构，新出现了一个pbmc@assays[["SCT"]]的槽，并将SCTransform的结果储存在中，结果中也有对应的原始矩阵，Normalize矩阵，scale后的矩阵和找到的var.features（SCTransform默认是找3000个var.features，上面的方法 默认是2000个）。
在常规分析中，使用少量的PC既能关注到关键的生物学差异，又能够不引入更多的技术差异，相当于一种保守性的做法。它会失去一些生物差异信息但是同时又在常规手段中比较安全。
但sctransform的归一化、标准化都做得不错，多输入一些PCs能提取更多的生物差异，并且兼顾不引入技术误差。sctransform认为：新增加的这1000个基因就包含了之前没有检测到的微弱的生物学差异。而且，即使使用全部的全部的基因去做下游分析，得到的结果也是和sctransform这3000个基因的结果相似
因此NormalizeData, ScaleData, FindVariableFeatures和SCTransform的结果存放在不同的slot中，不必担心进行了Normalize和scale再进行SCTransform会对数据过度处理，在后续处理数据的时候注意选择需要的矩阵即可。
SCTransform得到的结果都保存在pbmc[["SCT"]]中：
值得注意的是，SCTransform分析结果是@scale.data，@counts和@data都是用@scale.data反向计算出来的。
pbmc[["SCT"]]@scale.data：包含了残差数据，用作PCA的输入。这个数据不是稀疏矩阵，因此会占用大量内存。不过SCTransform函数计算的时候，为了节省内存，默认使用了return.only.var.genes = TRUE，只对差异基因进行scale的结果。设置为FALSE相当于三步法的ScaleData(features=all.genes)。- 关于为什么SCT运行完基因数会减少
Christoph作了回复：sctransform的min_cells参数会忽略一部分细胞。如果想要所有细胞，调整参数即可。
- RNA和SCT的区别
pbmc[["SCT"]]@counts：相当于细胞数据量相同情况下的值（也就是均一化了的UMI count值）。
pbmc[["SCT"]]@data：包含了上面count值的log-normalized结果，有利于后面可视化
目前可以使用pbmc[["SCT"]]@data结果进行差异分析和数据整合，实际上，官方更推荐直接使用残差值pbmc[["SCT"]]@scale.data。 但是seurat目前不支持，因此还得用@data的数据。
⚠️在运行完三步法标准化再运行SCT标准化，后续分析默认使用SCT矩阵。可以使用DefalutAssay函数可以查看目前默认用的是哪个矩阵。
2. 两种方法比较
library(VennDiagram)
library(gplots)
v1 <- pbmc@assays[['RNA']]@var.features
v2 <- pbmc@assays[['SCT']]@var.features
p <- venn.diagram(list(v1,v2), col=c("red","green"), alpha = c(0.5, 0.5),
category.names=c('pbmc','pbmc.sc'), filename=NULL, margin=0.15)
grid.draw(p)
可以看到两种方法找到的高变基因还是有一些差别
3. SCTransform做了些什么
SCTransform利用了正则化负二项分布（regularized negative binomial regression）计算了技术噪音模型，得到的残差是归一化值，有正有负。正值表示：考虑到细胞群体中基因的平均表达量和细胞测序深度，某个细胞的某个基因所包含的UMIs比预测值要高。
首先提取出表达矩阵并对每个基因做简单的统计包括表达量均值，方差，以及可检测率（detection rate）。从一张均值和方差的散点图可以看到在均值较小的一定范围内，均值和方差有线性关系。均值再变大时，方差呈现过度分散（overdispersion）。
pbmc_data <- as.matrix(pbmc@assays$RNA@counts)
gene_attr <- data.frame(mean = rowMeans(pbmc_data),
detection_rate = rowMeans(pbmc_data > 0),
var = apply(pbmc_data, 1, var))
gene_attr$log_mean <- log10(gene_attr$mean)
gene_attr$log_var <- log10(gene_attr$var)
rownames(gene_attr) <- rownames(pbmc_data)
ggplot(gene_attr, aes(log_mean, log_var)) + geom_point(alpha = 0.3, shape = 16) +
geom_density_2d(size = 0.3) + geom_abline(intercept = 0, slope = 1, color = "red")
接下来我们画一张基因表达均值和基因可检测率的关系图，图中红线是假设基因表达均值符合柏松分布的期望可检测率。可以看到高表达和低表达的基因符合预期值，基因表达量局中的基因的可检测率低于预期。
x = seq(from = -3, to = 2, length.out = 1000)
poisson_model <- data.frame(log_mean = x, detection_rate = 1 - dpois(0, lambda = 10^x))
ggplot(gene_attr, aes(log_mean, detection_rate)) + geom_point(alpha = 0.3, shape = 16) +
geom_line(data = poisson_model, color = "red") + theme_gray(base_size = 8)
第三步来看对每个细胞，其总UMI计数和能检测到的基因数量是正相关的。
cell_attr <- data.frame(n_umi = colSums(pbmc_data),
n_gene = colSums(pbmc_data > 0))
ggplot(cell_attr, aes(n_umi, n_gene)) + geom_point(alpha = 0.3, shape = 16) + geom_density_2d(size = 0.3)
综合以上观察，作者认为每个基因的表达量符合负二项分布且其在每个细胞内的表达量和该细胞的总基因表达量（total UMI）呈正相关。公式更为细节的地方包括复杂度调整（regularization）和对数联系函数（link function）。最终对每个基因在每个细胞的表达量是真实表达量和拟合后的表达预期值的差值。这就是为什么用SCTransform省略了ScaleData那一步，因为是先得到scaled data（皮尔逊残差），再去反推得到normalized data 和count data的。
我们跑一下SCTransform这个函数并看一下基因表达量的改变。我们挑选了三个基因MALAT1，RPL10和FTL。下图上半部显示矫正前，下半部显示矫正后基因表达量和测序深度的关系。可以看到经过SCTransform我们达到了目的，基因表达量不再和每个细胞的测序深度有相关性。
# run sctransform
set.seed(44)
vst_out <- sctransform::vst(pbmc_data, latent_var = c("log_umi"), return_gene_attr = TRUE, return_cell_attr = TRUE)
sctransform::plot_model(vst_out, pbmc_data, c("MALAT1", "RPL10", "FTL"), plot_residual = TRUE)
最后我们跑一下Seurat标准scale data流程，为保证可比性，我们regress out总基因表达量（total UMI）。同样我们画一下标准流程出来的残差和细胞总表达量的关系，可以看到同样的基因表达量不再和每个细胞的测序深度有相关性。SCTransform对比标准流程的优势可以看到对于像FTL这样的基因，用SCTransform后能明显看到其表达分为高和低两组细胞，而用标准流程此现象基本看不到。（pbmc3k的数据使用标准scale流程也可以分开）
pbmc <- NormalizeData(pbmc)
pbmc <- ScaleData(pbmc,vars.to.regress = c("nCount_RNA","percent.mt"), features = rownames(pbmc@assays$RNA@counts))
regular_df <- data.frame()
for(i in c("MALAT1", "RPL10", "FTL")){
tmp <- data.frame(cell_log_umi = log10(pbmc@meta.data$nCount_RNA + 1),
scale.residual = pbmc@assays$RNA@scale.data[rownames(pbmc@assays$RNA@scale.data) == i , drop=FALSE],
gene = rep(i, length(log10(pbmc@meta.data$nCount_RNA + 1)))
)
regular_df <- rbind(regular_df, tmp)
}
ggplot(regular_df, aes(cell_log_umi, scale.residual)) + geom_point(alpha = 0.3, shape = 16) +
geom_density_2d(size = 0.3, colour="cadetblue1")+
facet_wrap(.~factor(gene,levels=c("MALAT1", "RPL10", "FTL")))
