【单细胞分析】换个姿势标准化单细胞数据

简介

由于单细胞数据的高维度，基因长度差异，覆盖度差异及实验过程中的偏好性等因素，对前期数据进行有质量的标准化，对后续分析结果解读至关重要。标准化的分析比较多，对于常规的群体RNA分析而言，各种软件包也有相应的标准化方式，基因的定量也有RPKM ，TPM等标准化指标，单细胞数据对标准化的要求更高，本篇文章主要介绍新的标准化方法sctransform，他可以应用到Seurat中替换其原有的标准化函数。seuratV3简介及实操有对其相应的介绍。

安装R包及数据导入创建对象

#打开R环境直接安装
devtools::install_github(repo = 'ChristophH/sctransform', ref = 'develop')
devtools::install_github(repo = 'satijalab/seurat', ref = 'release/3.0')
#也可以直接安装SeuratV3包
install.packages("seurat")
library(Seurat)
library(ggplot2)
library(sctransform)
#测试数据在10X官网下载
pbmc_data <- Read10X(data.dir = "pbmc3k/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc_data)

细节介绍

为什么在使用sctransform时可以选择更多的PC？
答：

在标准的Seurat工作流程中，我们专注于此数据集的10PC，但我们强调结果与此参数的更高设置相似。有趣的是，我们发现在使用sctransform时，我们通常会将此参数推得更高。我们相信这是因为sctransform工作流程执行更有效的规范化，强烈消除了数据的技术效果。
即使在标准对数归一化之后，测序深度的变化仍然是一个混淆因素（见图1），这种影响可以巧妙地影响更高的PC。在sctransform中，这种影响大大减轻了（见图3）。这意味着更高的PC更有可能代表微妙的，但与生物学相关的异质性来源 - 因此包括它们可能会改善下游分析。
此外，sctransform默认返回3,000个变量功能，而不是2,000。基本原理类似，额外的变量特征不太可能由细胞间的技术差异驱动，而是可能代表更微妙的生物波动。一般来说，我们发现用sctransform产生的结果较少依赖于这些参数（实际上，当在转录组中使用所有基因时，我们获得了几乎相同的结果，尽管这确实降低了计算效率）。这可以帮助用户生成更强大的结果，此外，还可以应用具有相同参数设置的标准anlaysis管道，这些参数设置可以快速应用于新数据集：

例如，以下代码在单个命令中复制完整的端到端工作流：

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存储的归一化值在哪里？
答：
如的预印本中所述，sctransform使用“正则化负二项式回归”计算scRNA-seq数据中的技术噪声模型。该模型的残差是标准化值，可以是正数或负数。给定基因中给定基因的阳性残基表明，鉴于基因在群体和细胞测序深度中的平均表达，观察到比预期更多的UMI，而负残差表明相反。
sctransfrom的结果存储在“SCT”测定中。您可以在插图，命令备忘单或开发者指南中了解有关Seurat中多种化验数据和命令的更多信息。

pbmc[["SCT"]]@scale.data包含残差（标准化值），并直接用作PCA的输入。请注意，这个矩阵是非稀疏的，因此如果为所有基因存储，则会占用大量内存。为了节省内存，我们通过在SCTransform函数调用中默认设置return.only.var.genes = TRUE来仅为变量基因存储这些值。
协助可视化和解释。我们还将Pearson残差转换回'更正'的UMI计数。您可以将这些解释为我们希望观察到的所有细胞都被测序到相同深度的UMI计数。如果您想确切了解我们如何做到这一点，请在此处查看正确的功能。
“更正的”UMI计数存储在pbmc[["SCT"]]@counts。我们存储这些校正计数的日志标准化版本，这pbmc[["SCT"]]@data对可视化非常有用。
您可以使用更正的对数标准化计数进行差异表达和积分。然而，原则上，直接对残差（存储在scale.data槽中）本身执行这些计算是最佳的。Seurat v3目前不支持此功能，但很快就会支持。

两种标准化方式比较

作为参考，首先应用标准的Seurat工作流程，并使用对数标准化

pbmc_logtransform <- pbmc
pbmc_logtransform <- NormalizeData(pbmc_logtransform, verbose = FALSE)
pbmc_logtransform <- FindVariableFeatures(pbmc_logtransform, verbose = FALSE) 
pbmc_logtransform <- ScaleData(pbmc_logtransform, verbose = FALSE) 
pbmc_logtransform <- RunPCA(pbmc_logtransform, verbose = FALSE) 
pbmc_logtransform <- RunUMAP(pbmc_logtransform,dims = 1:20, verbose = FALSE)

为了比较，现在应用sctransform规范化

# Note that this single command replaces NormalizeData, ScaleData, and FindVariableFeatures.
# Transformed data will be available in the SCT assay, which is set as the default after running sctransform
pbmc <- SCTransform(object = pbmc, verbose = FALSE)

通过PCA和UMAP嵌入执行降维

# These are now standard steps in the Seurat workflow for visualization and clustering
pbmc <- RunPCA(object = pbmc, verbose = FALSE)
pbmc <- RunUMAP(object = pbmc, dims = 1:20, verbose = FALSE)
pbmc <- FindNeighbors(object = pbmc, dims = 1:20, verbose = FALSE)
pbmc <- FindClusters(object = pbmc, verbose = FALSE)

在sctransform和log-normalized嵌入上可视化聚类结果。

pbmc_logtransform$clusterID <- Idents(pbmc)
Idents(pbmc_logtransform) <- 'clusterID'
plot1 <- DimPlot(object = pbmc, label = TRUE) + NoLegend() + ggtitle('sctransform') 
plot2 <- DimPlot(object = pbmc_logtransform, label = TRUE)
plot2 <- plot2 + NoLegend() + ggtitle('Log-normalization') 
CombinePlots(list(plot1,plot2))

用户可以基于规范标记单独注释群集。然而，与对数归一化分析相比，sctransform归一化显示出更明显的生物学区别。例如，注意在对数归一化分析中如何将簇0,1,9和11（所有不同的T细胞簇）混合在一起。sctransform分析显示：

基于CD8A，GZMK，CCL5，GZMK表达，清楚分离三个CD8 T细胞簇（幼稚，记忆，效应器）
基于S100A4，CCR7，IL32和ISG15清楚分离三个CD4 T细胞簇（幼稚，记忆，IFN激活）
B细胞群中的附加发育子结构，基于TCL1A，FCER2
基于XCL1和FCGR3A的NK细胞簇内的额外亚结构（CD56dim对比亮）
将规范标记基因可视化为小提琴图。

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

在sctransform嵌入上可视化规范标记基因。

FeaturePlot(object = pbmc, features = c("CD8A","GZMK","CCL5","S100A4","S100A4","CCR7"), pt.size = 0.2,ncol = 3)

FeaturePlot(object = pbmc, features = c("CD3D","ISG15","TCL1A","FCER2","XCL1","FCGR3A"), pt.size = 0.2,ncol = 3)

参考材料:

https://rawgit.com/ChristophH/sctransform/master/inst/doc/seurat.html
https://satijalab.org/seurat/v3.0/sctransform_vignette.html