去除细胞效应和基因效应

Posted on by

前面我们一起学习了单细胞转录组数据的上游分析，而且了解了自己的项目的样本数量和测序量，还过滤了不合格的细胞和基因， 教程目录如下：

其实质量控制三部曲，还有一个很关键的点没有讲解，就是多个样本整合，并且区分批次效应和生物学差异。但是这个点很大程度是依赖于经验，就是说，要想搞清楚，需要写很多自定义的代码去一点一滴的探索，而不仅仅是流程。所以我们就只能介绍到这里，假设大家都拿到了干净的表达矩阵，而且可以很肯定的说这个表达矩阵做下游分析是ok的。

那么我们就来看看，有了干净的表达矩阵后下游分析的第一个分析要点就是：归一化和标准化。

归一化和标准化的区别

实际上口语里面通常是没办法很便捷的区分这两个概念，我查阅了大家的资料，发现基本上都是混淆在z-score转换和0-1转换这两个上面。所以我这里把归一化和标准化替换成为去除样本/细胞效应或者去除基因效应：

  • 首先去除样本/细胞效应：因为不同样本或者细胞的测序数据量不一样，所以同样的一个基因在不同细胞，哪怕你看到的表达量是一样的，但是背后细胞整体测序数据量的差异其实反而说明了这个基因在不同细胞表达量其实是有差异的。在seurat3里面的代码是：
 sce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000)

默认文库大小是1万，默认会取log值。

  • 然后去除基因效应，这个主要是在绘制热图的时候会需要使用，因为个别基因表达量超级高，在热图里面一枝独秀，实际上我们并不会关心不同基因的表达量高低，我们仅仅是想看指定基因在不同细胞的高低而已，这样的话，就把该基因的表达量在不同细胞的数值，进行z-score这样的转换。在seurat3里面的代码是：
 sce <- ScaleData(sce)

这样处理后的表达矩阵，就可以进行后续的降维聚类分群啦，我们下期再讲。

取log让表达量离散程度降低

这个时候眼尖的朋友其实看到了，在使用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下载到的表达矩阵的形式。

使用GetAssayData函数查看不同形式的表达矩阵

在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函数。

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

​

Comments are closed.