SCTransform到底做了什么

大名鼎鼎的单细胞分析工具Seurat从3.0版本开始引进了SCTransform这个函数用来对数据做标准化，并且这一个函数可以代替三个函数（NormalizeData, ScaleData, FindVariableFeatures）的运行。今天想来看一下SCTransform这个函数到底做了什么，以及和标准的流程出来的结果有什么不一样的地方。

首先读入10X的PBMC 4K数据（数据下载地址：https://support.10xgenomics.com/single-cell-gene-expression/datasets/2.1.0/pbmc4k）。

tmp.data <- Read10X(data.dir = "./PBMC_4k/filtered_gene_bc_matrices/GRCh38/")

pbmc <- CreateSeuratObject(counts = tmp.data, project = "PBMC_4K", min.cells = 3, min.features = 200)

首先提取出表达矩阵并对每个基因做简单的统计包括表达量均值，方差，以及可检测率（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")

image.png

接下来我们画一张基因表达均值和基因可检测率的关系图，图中红线是假设基因表达均值符合柏松分布的期望可检测率。可以看到高表达和低表达的基因符合预期值，基因表达量局中的基因的可检测率低于预期。

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)

image.png

第三步来看对每个细胞，其总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)

image.png

综合以上观察，作者认为每个基因的表达量符合负二项分布且其在每个细胞内的表达量和该细胞的总基因表达量（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)

image.png

最后我们跑一下Seurat标准scale data流程，为保证可比性，我们regress out总基因表达量（total UMI）。同样我们画一下标准流程出来的残差和细胞总表达量的关系，可以看到同样的基因表达量不再和每个细胞的测序深度有相关性。SCTransform对比标准流程的优势可以看到对于像FTL这样的基因，用SCTransform后能明显看到其表达分为高和低两组细胞，而用标准流程此现象基本看不到。

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")))

image.png

Reference：
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1874-1