Seurat從3.0版本引進(jìn)了SCTransform這個(gè)函數(shù)用來對數(shù)據(jù)做標(biāo)準(zhǔn)化，并且這一個(gè)函數(shù)可以代替三個(gè)函數(shù)（NormalizeData, ScaleData, FindVariableFeatures）的運(yùn)行。 且其對測序深度的校正效果要好于log標(biāo)準(zhǔn)化。（10萬以內(nèi)的細(xì)胞都建議使用SCT標(biāo)準(zhǔn)化）
?SCTransform對測序深度的校正效果很好，也可用于矯正線粒體等因素的影響，但不能用于批次矯正。
官網(wǎng)：
https://cran.r-project.org/web/packages/sctransform/index.html
https://satijalab.org/seurat/archive/v3.0/sctransform_vignette.html
官網(wǎng)對這個(gè)函數(shù)的描述：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. （這是一個(gè)用方差穩(wěn)定變換對單細(xì)胞UMI count 數(shù)據(jù)標(biāo)準(zhǔn)化的方法，方差穩(wěn)定變換是基于負(fù)二項(xiàng)回歸。這個(gè)函數(shù)在對數(shù)據(jù)進(jìn)行均一化的同時(shí)還可以去除線粒體紅細(xì)胞等混雜因素的影響。）
1. NormalizeData, ScaleData, FindVariableFeatures和SCTransform結(jié)果的存放位置
示例數(shù)據(jù)下載：pbmc3k
#讀入數(shù)據(jù)，創(chuàng)建seurat對象
pbmc <- Read10X('./filtered_gene_bc_matrices/hg19/')
pbmc <- CreateSeuratObject(pbmc,project = 'pbmc3k',min.cells = 3,min.features = 200)
## 質(zhì)控
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
常規(guī)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)
這一步的結(jié)果都儲存在pbmc@assays[["RNA"]]中。
再使用SCTransform進(jìn)行標(biāo)準(zhǔn)化
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)
View(pbmc)
再次查看pbmc的數(shù)據(jù)結(jié)構(gòu)，新出現(xiàn)了一個(gè)pbmc@assays[["SCT"]]的槽，并將SCTransform的結(jié)果儲存在中，結(jié)果中也有對應(yīng)的原始矩陣，Normalize矩陣，scale后的矩陣和找到的var.features（SCTransform默認(rèn)是找3000個(gè)var.features，上面的方法 默認(rèn)是2000個(gè)）。
在常規(guī)分析中，使用少量的PC既能關(guān)注到關(guān)鍵的生物學(xué)差異，又能夠不引入更多的技術(shù)差異，相當(dāng)于一種保守性的做法。它會失去一些生物差異信息但是同時(shí)又在常規(guī)手段中比較安全。
但sctransform的歸一化、標(biāo)準(zhǔn)化都做得不錯，多輸入一些PCs能提取更多的生物差異，并且兼顧不引入技術(shù)誤差。sctransform認(rèn)為：新增加的這1000個(gè)基因就包含了之前沒有檢測到的微弱的生物學(xué)差異。而且，即使使用全部的全部的基因去做下游分析，得到的結(jié)果也是和sctransform這3000個(gè)基因的結(jié)果相似
因此NormalizeData, ScaleData, FindVariableFeatures和SCTransform的結(jié)果存放在不同的slot中，不必?fù)?dān)心進(jìn)行了Normalize和scale再進(jìn)行SCTransform會對數(shù)據(jù)過度處理，在后續(xù)處理數(shù)據(jù)的時(shí)候注意選擇需要的矩陣即可。
SCTransform得到的結(jié)果都保存在pbmc[["SCT"]]中：
值得注意的是，SCTransform分析結(jié)果是@scale.data，@counts和@data都是用@scale.data反向計(jì)算出來的。
pbmc[["SCT"]]@scale.data：包含了殘差數(shù)據(jù)，用作PCA的輸入。這個(gè)數(shù)據(jù)不是稀疏矩陣，因此會占用大量內(nèi)存。不過SCTransform函數(shù)計(jì)算的時(shí)候，為了節(jié)省內(nèi)存，默認(rèn)使用了return.only.var.genes = TRUE，只對差異基因進(jìn)行scale的結(jié)果。設(shè)置為FALSE相當(dāng)于三步法的ScaleData(features=all.genes)。- 關(guān)于為什么SCT運(yùn)行完基因數(shù)會減少
Christoph作了回復(fù)：sctransform的min_cells參數(shù)會忽略一部分細(xì)胞。如果想要所有細(xì)胞，調(diào)整參數(shù)即可。
- RNA和SCT的區(qū)別
pbmc[["SCT"]]@counts：相當(dāng)于細(xì)胞數(shù)據(jù)量相同情況下的值（也就是均一化了的UMI count值）。
pbmc[["SCT"]]@data：包含了上面count值的log-normalized結(jié)果，有利于后面可視化
目前可以使用pbmc[["SCT"]]@data結(jié)果進(jìn)行差異分析和數(shù)據(jù)整合，實(shí)際上，官方更推薦直接使用殘差值pbmc[["SCT"]]@scale.data。 但是seurat目前不支持，因此還得用@data的數(shù)據(jù)。
??在運(yùn)行完三步法標(biāo)準(zhǔn)化再運(yùn)行SCT標(biāo)準(zhǔn)化，后續(xù)分析默認(rèn)使用SCT矩陣?？梢允褂肈efalutAssay函數(shù)可以查看目前默認(rèn)用的是哪個(gè)矩陣。
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利用了正則化負(fù)二項(xiàng)分布（regularized negative binomial regression）計(jì)算了技術(shù)噪音模型，得到的殘差是歸一化值，有正有負(fù)。正值表示：考慮到細(xì)胞群體中基因的平均表達(dá)量和細(xì)胞測序深度，某個(gè)細(xì)胞的某個(gè)基因所包含的UMIs比預(yù)測值要高。
首先提取出表達(dá)矩陣并對每個(gè)基因做簡單的統(tǒng)計(jì)包括表達(dá)量均值，方差，以及可檢測率（detection rate）。從一張均值和方差的散點(diǎn)圖可以看到在均值較小的一定范圍內(nèi)，均值和方差有線性關(guān)系。均值再變大時(shí)，方差呈現(xiàn)過度分散（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")
接下來我們畫一張基因表達(dá)均值和基因可檢測率的關(guān)系圖，圖中紅線是假設(shè)基因表達(dá)均值符合柏松分布的期望可檢測率?？梢钥吹礁弑磉_(dá)和低表達(dá)的基因符合預(yù)期值，基因表達(dá)量局中的基因的可檢測率低于預(yù)期。
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)
第三步來看對每個(gè)細(xì)胞，其總UMI計(jì)數(shù)和能檢測到的基因數(shù)量是正相關(guān)的。
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)
綜合以上觀察，作者認(rèn)為每個(gè)基因的表達(dá)量符合負(fù)二項(xiàng)分布且其在每個(gè)細(xì)胞內(nèi)的表達(dá)量和該細(xì)胞的總基因表達(dá)量（total UMI）呈正相關(guān)。公式更為細(xì)節(jié)的地方包括復(fù)雜度調(diào)整（regularization）和對數(shù)聯(lián)系函數(shù)（link function）。最終對每個(gè)基因在每個(gè)細(xì)胞的表達(dá)量是真實(shí)表達(dá)量和擬合后的表達(dá)預(yù)期值的差值。這就是為什么用SCTransform省略了ScaleData那一步，因?yàn)?strong>是先得到scaled data（皮爾遜殘差），再去反推得到normalized data 和count data的。
我們跑一下SCTransform這個(gè)函數(shù)并看一下基因表達(dá)量的改變。我們挑選了三個(gè)基因MALAT1，RPL10和FTL。下圖上半部顯示矯正前，下半部顯示矯正后基因表達(dá)量和測序深度的關(guān)系?？梢钥吹浇?jīng)過SCTransform我們達(dá)到了目的，基因表達(dá)量不再和每個(gè)細(xì)胞的測序深度有相關(guān)性。
# 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標(biāo)準(zhǔn)scale data流程，為保證可比性，我們r(jià)egress out總基因表達(dá)量（total UMI）。同樣我們畫一下標(biāo)準(zhǔn)流程出來的殘差和細(xì)胞總表達(dá)量的關(guān)系，可以看到同樣的基因表達(dá)量不再和每個(gè)細(xì)胞的測序深度有相關(guān)性。SCTransform對比標(biāo)準(zhǔn)流程的優(yōu)勢可以看到對于像FTL這樣的基因，用SCTransform后能明顯看到其表達(dá)分為高和低兩組細(xì)胞，而用標(biāo)準(zhǔn)流程此現(xiàn)象基本看不到。（pbmc3k的數(shù)據(jù)使用標(biāo)準(zhǔn)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")))
