其实质量控制三部曲,还有一个很关键的点没有讲解,就是多个样本整合,并且区分批次效应和生物学差异。但是这个点很大程度是依赖于经验,就是说,要想搞清楚,需要写很多自定义的代码去一点一滴的探索,而不仅仅是流程。所以我们就只能介绍到这里,假设大家都拿到了干净的表达矩阵,而且可以很肯定的说这个表达矩阵做下游分析是ok的。
那么我们就来看看,有了干净的表达矩阵后下游分析的第一个分析要点就是:归一化和标准化。
实际上口语里面通常是没办法很便捷的区分这两个概念,我查阅了大家的资料,发现基本上都是混淆在z-score转换和0-1转换这两个上面。所以我这里把归一化和标准化替换成为去除样本/细胞效应或者去除基因效应:
sce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000)
sce <- ScaleData(sce)
这样处理后的表达矩阵,就可以进行后续的降维聚类分群啦,我们下期再讲。
这个时候眼尖的朋友其实看到了,在使用NormalizeData函数的时候,有一个 normalization.method = "LogNormalize" 参数被设置了,这个是为什么呢?
其实很简单,原始的raw counts矩阵是一个离散型的变量,离散程度很高。有的基因表达丰度比较高,counts数为10000,有些低表达的基因counts可能10,甚至在有些样本中为0。因为是单细胞转录组,drop-out现象很严重,其实大量的基因在很多细胞的表达量都是0,如果是10x数据,甚至会出现一个表达矩阵里面97%的数值都是0 的现象。
如果对表达量取一下log10,发现10000变成了4,10变成了1,这样之前离散程度很大的数据就被集中了。
有时当表达量为0时,取log会出现错误,可以log(counts+1)来取log值。当x=1时,所有的log系列函数值都为0。这样原本表达量为0的值,取log后仍为0。
这也就是UCSC的XENA下载到的表达矩阵的形式。
在seurat3里面,很方便进行各种形式的归一化或者标准化,同样的,也很容易查看处理前后的表达矩阵,使用GetAssayData函数即可。
This function can be used to pull information from any of the slots in the Assay class. For example, pull one of the data matrices("counts", "data", or "scale.data").
如下代码:
# 最原始数据
GetAssayData(sce,'counts')[1:6,1:6]
# 去除了细胞测序数据量后
1/(colSums(as.data.frame(GetAssayData(sce,'counts')))[2]/10000)
GetAssayData(sce,'data')[1:6,1:6]
# z-score后
GetAssayData(sce,'scale.data')[1:6,1:6]
如下结果:
> # 最原始数据
> GetAssayData(sce,'counts')[1:6,1:6]
6 x 6 sparse Matrix of class "dgCMatrix"
6A-11 6A-13 6A-14 6A-15 6A-16 6A-17
DDX11L1 . . . . . .
WASH7P . 1 . . . .
MIR6859-1 . . . . . .
RP11-34P13.3 . . . . . .
OR4G11P . . . . . .
OR4F5 . . . . . .
> # 去除了细胞测序数据量后
> 1/(colSums(as.data.frame(GetAssayData(sce,'counts')))[2]/10000)
6A-13
0.009725971
> GetAssayData(sce,'data')[1:6,1:6]
6 x 6 sparse Matrix of class "dgCMatrix"
6A-11 6A-13 6A-14 6A-15 6A-16 6A-17
DDX11L1 . . . . . .
WASH7P . 0.009678978 . . . .
MIR6859-1 . . . . . .
RP11-34P13.3 . . . . . .
OR4G11P . . . . . .
OR4F5 . . . . . .
> # z-score后
> GetAssayData(sce,'scale.data')[1:6,1:6]
6A-11 6A-13 6A-14 6A-15 6A-16 6A-17
RP5-857K21.9 -0.8302953 -0.8302953 -0.7327066 -0.7224892 -0.7683106 -0.7692732
RP5-857K21.8 -0.5066862 -0.5688591 -0.5349974 -0.4881690 -0.3650311 -0.4524682
PRKCZ -0.8350477 -0.5996852 -0.8350477 -0.8350477 -0.5988844 -0.8350477
TNFRSF14 -0.2550752 -0.2550752 -0.2550752 -0.2550752 -0.2550752 -0.2550752
TP73 -0.2761133 -0.2761133 -0.2761133 -0.2761133 -0.2761133 -0.2337521
HES2 -0.4569104 -0.4569104 -0.4569104 -0.4569104 -0.4569104
每个细胞测序数据量的不一致是很容易理解的,但其实细胞之间还有很多其它效应,比如线粒体基因含量,ERCC含量等等,那些处理起来,其实就是深入了解我们讲解seurat里面的NormalizeData和ScaleData函数。