单细胞分析之标准化处理

书籍原始链接:https://www.bioconductor.org/books/release/OSCA/normalization.html

本篇笔记根据书中介绍的内容,按照自己的实际操作情况,进行复现和重新整理。

为什么要做标准化?要解决什么样的问题?进一步消除细胞之间的系统差异(我们前面做质控分析,其实也涉及到这一点),使得后面在比较细胞之间的表达量的时候,这种差异主要来源于生物条件的不同,而非是一些实验过程中的系统误差。

补充小知识点1:标准化与批次校正的区别

标准化:不管批次结构,仅考虑技术偏差(文库之间cDNA捕获与PCR扩增的技术差异,导致覆盖率的不同)

批次校正:仅在批次之间发生,而且必须同时考虑技术偏差和生物学差异,技术偏差影响基因以相似的方式,或者至少是一种与他们的生物物理特性相关的方式,如长度,GC含量,而生物学差异则高度不可预测。

因此,标准化与批次校正是不同的东西,需要使用不同的计算方法。切忌将两者混为一谈。

在众多的标准化策略中,缩放标准化(scaling normalization)是最常用也是最简单的标准化策略。

补充知识点2:缩放标准化

计算方法:将每一个细胞的所有计数值除以细胞特异性的scaling factor。

基本假设:任何细胞特异性偏倚(例如捕获或扩增效率)都会通过缩放该细胞的预期平均数来同等地影响所有基因。

标准化的结果可以用于下游的聚类与降维的分析。

1、加载测试数据集

library(scRNAseq)
example_sce <- ZeiselBrainData()

example_sce数据集的详细信息:

example_sce

 

class: SingleCellExperiment 
dim: 20006 3005 
metadata(0):
assays(1): counts
rownames(20006): Tspan12 Tshz1 ... mt-Rnr1 mt-Nd4l
rowData names(1): featureType
colnames(3005): 1772071015_C02 1772071017_G12 ... 1772066098_A12 1772058148_F03
colData names(22): tissue group # ... altexps_repeat_percent total
reducedDimNames(0):
altExpNames(2): ERCC repeat

2、方案一:文库大小的标准化

文库大小:每一个细胞所有基因表达量的总和。

每个细胞的“library size factor”直接与其库大小成正比,其中定义了比例常数,以便所有单元的平均大小因子等于1。此种标准化可以确保标准化后的表达值与原始数据的范围相同。

library(scater)
#逐个计算3005个细胞的factor,每个细胞有一个对应的factor
lib.sf.zeisel <- librarySizeFactors(example_sce)  
summary(lib.sf.zeisel)

 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1721  0.5437  0.8635  1.0000  1.2895  4.2466 

在这个数据集中,最大的factor值(4.2466)与最小的factor值(0.1721)之间差40倍之多。这是scRNA-seq数据覆盖率变异性的典型特征。

绘制直方图,统计3005个细胞的library size factor的基本分布情况。

hist(log10(lib.sf.zeisel), xlab="Log10[Size factor]", col='lightblue')
#注意加log10()处理之后,使图像整体更加美观
每一个细胞的library size factor的具体分布
(我感觉通过上述步骤只是计算了一个factor,并没有把数据矩阵中的每一个数据都最终标准化。是做的还不够的)

这种归一化方法的限制条件:这种方法在使用的过程中,存在一个基本的假设,即两个细胞内的上调的基因集与其下调的基因集幅度相同,可以互相抵消(称之为平衡)。但是实际的情况是,这种情况在scRNA-seq的数据中几乎不存在,因此,这种归一化方法不能为下游的分析提供可靠的归一化数据。

由于,归一化的可靠性,并不是探索scRNA-seq数据的主要考虑的因素(所以,在某些情况下,可以睁一只眼闭一只眼?)。在聚类的过程中,这种方法所引起的归一化偏差不会影响簇的分离,而只会影响细胞类型或者簇之间的log10()值的幅度或者方向(方向?)。

因此,这种归一化方法,适用于下游的聚类分析

3、方案二:去卷积标准化

我们对上文提到的成分偏见(composition biases)这一现象进行重新的分析。假设有两个细胞,细胞A中的geneX相较于细胞B而言是上调的。

(1)这种上调意味着,在A细胞中更多的资源被分配给了geneX,而其他非差异表达基因的覆盖率则相应的降低(当每个细胞中测序的总counts数是固定的)。

(2)当geneX被分配有更多的reads或者UMIs,A细胞的测序总counts是上升的。增加了上述方法中的library size factor的大小,而使得所有正常表达的非差异表达基因的值减少。

在这以上的两种情况下,A细胞中的非差异表达的基因,相较于B细胞而言,会被错误的认为下调了(也就是存在一个超异常值,从而拉高了整体的平均值,在归一化的过程中,一般值相对平均水平低了,如果除去异常值,则是不一样的效果)。

这时,一种标准化的策略是,汇总了许多细胞中的counts以增加counts的大小(单细胞数据中普遍的0或者很低的值,所以总counts很小?),以进行准确的size factor估算。然后,将基于库的size factor“分解”为cell-based factors,以标准化每个细胞的表达谱。(这里我是没怎么明白的?许多撮面疙瘩,揉揉,变成一个面团,然后再重新搓成许多的面疙瘩?)

BiocManager::install("scran")
library(scran)
set.seed(100)

#我们使用了一个预集群步骤quickCluster(),具体的原理下面会介绍。
clust.zeisel <- quickCluster(example_sce) 

#查看cluster的分群
table(clust.zeisel)

#按照分群,计算factor
deconv.sf.zeisel <- calculateSumFactors(example_sce, cluster=clust.zeisel)
summary(deconv.sf.zeisel)

clust.zeisel的具体的数值:

clust.zeisel
  1   2   3   4   5   6   7   8   9  10  11  12  13  14 
139 164 251 234 372 272 151 428 117 148 289 225 111 104 

deconv.sf.zeisel的具体数值(可以与上面的文库大小的factor进行比较,同一套数据):

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1209  0.4621  0.8082  1.0000  1.3177  4.8508 
  • 使用quickCluster()进行聚类,每一个类进行分组的标准化。重新调整分组因子,以便群集之间具有可比性。

比较文库大小标准化得到的每一个细胞的factor与去卷积化得到的factor之间的关系。

plot(lib.sf.zeisel, deconv.sf.zeisel, xlab="Library size factor",
     ylab="Deconvolution size factor", pch=16,
     col=as.integer(factor(example_sce$tissue)))
abline(a=0, b=1, col="green")
绿色的那条线为y=x的渐近线。这张图展示的是两组标准化得到的factor的比较。

 

标准化的方法,对于差异表达基因的分析影响还是挺大的,所以在操作的时候要尽量回避能够产生composition biases的策略(如文库大小标准化策略)。

4、方案三、通过spike-in标准化(在本节将重点介绍spike-in的用途)

基于的基本假设:每一个细胞中添加有等量的spike-inRNA。spike-in转录本的覆盖率的系统差异(systematic difference)只能归因于细胞特异性偏差,如捕获的效率,测序的深度。

为了移除这些偏差,我们通过“spike-in size factor”来均衡细胞之间的spike-in覆盖率。相较于之前的标准化的方法,spike-in标准化方法对系统的生物学条件无要求(比如,缺少许多差异表达基因)。

在实际的处理条件中,spike-in标准化通常被用在如果每一个细胞的总RNA含量的差异是我们感兴趣的内容,并且这必需保留在下游的分析中。对于给定的细胞,其内源RNA总量的上升不会引起其spike-in size factor的变化。这确保着,总RNA含量对整个细胞群表达的影响,不会在scaling的过程中除去。相比之下,上述其他标准化方法只是将总RNA含量的任何变化解释为偏差的一部分,并将其消除。

#加载scran包,librarySizeFactors()函数位于这个包中,如果之前已经加载过不必重复加载了
library(scran)
#计算spike-in factor
spike.sf.zeisel<-librarySizeFactors(altExp(example_sce,"ERCC"))
summary(spike.sf.zeisel)

关于函数computeSpikeFactors(),librarySizeFactor(altExp())等函数的具体用法,可参考链接

对于同样一套数据集,spike-in factor的值分布情况。

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1565  0.7697  0.9844  1.0000  1.1808  3.4171 

5、应用size factors。

(1)缩放以及对数转换

补充知识点3:为什么要进行log转换?

对数值的变化,代表着表达值的对数倍的变化。这对于下游的基于欧氏距离的处理步骤如聚类和降维,非常必要。通过对对数转换后的数据进行操作,我们确保这些过程是基于细胞表达的对数倍来计算距离。

在进行对数转换时,我们通常会添加一个伪计数以避免未定义的零值。对于低丰度基因,较大的伪计数将有效地将细胞之间的对数倍变化缩小至零,这意味着下游的高维分析将更多地由高丰度基因的表达差异来驱动。相反,较小的伪计数将增加低丰度基因的相对贡献。常见的做法是使用1作为伪计数,原因很简单,即它保留了原始矩阵中的稀疏性(即输入中的零在变换后仍为零)。除大多数病理情况外,此方法在所有情况下均适用。

set.seed(100)
clust.zeisel <- quickCluster(example_sce) 
#这里和前面计算factor的操作不太一样的是,将其赋值给了原始数据集
example_sce<- computeSumFactors(example_sce, cluster=clust.zeisel)
example_sce<- logNormCounts(example_sce)
assayNames(example_sce)

 查看结果。

 colnames(colData(example_sce))

下面是colnames()的结果,可以看到sizeFactor被添加到了数据集中。

[1] "tissue"         "group #"        "total mRNA mol" "well"           "sex"           
 [6] "age"            "diameter"       "cell_id"        "level1class"    "level2class"   
[11] "sizeFactor"   

6、还有一种情况,没怎么看懂(这里不写了)。

 

后记:

其实按照操作流程一步一步将代码跑出来,倒不是难的。

我觉得这里非常难的地方,在于理解为什么要这样去处理?应该如何选择适合的方法?这个方法是否适合于我的数据?处理之后会对我的数据下游的分析产生什么样的影响?这个影响是我们想要的吗?

我们要从中做出选择,并尽可能做出优化。很多函数是直接内嵌在一体化操作的R包里,运行,可以直接得到结果。但是,你需要明白,这个函数内部进行着什么操作?是我们想要的吗?

 

 

 

 

 

 

 

 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值