RNA-seq中的那些统计学问题（二）FPKM/RPKM之外的那些标准化方法

Posted 2023-02-28

参考技术A

目录

由于不同文库测序深度不同，比较前当然要进行均一化！用总reads进行均一化可能最简单，其基于以下两个基本假设：

但在转录组中，通常一小部分极高丰度基因往往会贡献很多reads，如果这些“位高权重”的基因还是差异表达的，则会影响所有其它基因分配到的reads数，而且，两个样本总mRNA量完全相同的前提假设也过于理想了。那如何比较呢，各个方家使出浑身解数，有用中位数的，有用75分位数的，有用几何平均数的，有用TMM(trimmed mean of Mvalues)的等等，总之要 找一个更稳定的参考 值。

矫正的思路很简单，就是在变化的样本中寻找不变的量

那么在不同RNA-seq样本中，那些是不变的量呢？一个很容易想到的就是 管家基因 (House-keeping gene(s))

那么 Human 常用的 House-keeping gene 怎么确定？

目前大家用的比较多的一个human housekeeping gene list 来源于下面这篇文章，是2013年发表在 Cell系列的 Trends in Genetics 部分的一篇文章

使用Housekeeping gene的办法来进行相对定量，这种办法在一定程度上能够解决我们遇到的问题。但其实这种办法有一个 非常强的先验假设 ：housekeeping gene的表达量不怎么发生变化。其实housekeeping gene list有几千个，这几千个基因有一定程度上的变化是有可能的

spike-in方法 ：在RNA-Seq建库的过程中掺入一些预先知道序列信息以及序列绝对数量的内参。这样在进行RNA-Seq测序的时候就可以通过不同样本之间内参（spike-in）的量来做一条标准曲线，就可以非常准确地对不同样本之间的表达量进行矫正

比较常用的spike-in类型：ERCC Control RNA

CPM(count-per-million)

简单来说，就是找出多个样本中library size为中位数的样本，作为参考样本，将所有的样本的library size按比例缩放到参考样本的水平

选择一个library size为中位数的sample，以它为baseline，计算出其它sample对于baseline的normalization factor，即一个缩放因子：

然后基于该缩放因子对特定的sample中的每个基因的read count进行标准化（缩放）：

简单来说，就是排序后求平均，然后再回序

在R里面，推荐用preprocessCore 包来做quantile normalization，不需要自己造轮子啦！
但是需要明白什么时候该用quantile normalization，什么时候不应该用，就复杂很多了

该方法基于的假设是，即使处在不同条件下的不同个样本，大多数基因的表达是不存在差异的，实际存在差异的基因只占很小的部分那么我们只需要将这些稳定的部分找出来，作为标准化的内参，依据内参算出各个样本的标准化因子

（1）对每个基因计算几何平均数，得到一个假设的参考样本(pseudo-reference sample)

（2）对每个样本的每个基因对于参考样本计算Fold Change

（3）获取每个样本中Fold Change的中位数，我们就得到了非DE基因代表的Fold Change，该基因就是我们选择的该样本的内参基因，它的Fold Change就是该样本的标准化因子

该方法的思想与DESeq2的Median of Ratio相同，假设前提都是：大多数基因的表达是不存在差异的

它与DESeq2的不同之处在于对内参的选择上：

（1）移除所有未表达基因

（2）从众多样本中找出一个数据趋势较为平均的样本作为参考样本

（3）找出每个样本中的代表基因集，参考这些代表基因集的fold change，计算出该样本的标准化因子

寻找样本的代表基因集：依据基因的 偏倚程度 和 Reads数 大小选出——偏倚程度小、reads数居中的基因

由参考代表基因集计算样本的标准化因子：

FPKM和RPKM分别是什么

这两个量的计算方式的目的是为了解决计算RNA-seq转录本丰度时的两个bias：

首先，到底什么是RNA转录本的表达丰度这个问题

对于样本X，其有一个基因g被转录了mRNA_g次，同时样本X中所有基因的转录总次数假定是mRNA_total, 那么正确描述基因g转录丰度的值应该是：

则一个样本中基因表达丰度的均值为

而

所以

这个期望值竟然和测序状态无关！仅仅由样本中基因的总数所决定的

也就是说，对于同一个物种，不管它的样本是哪种组织（正常的或病变的），也不管有多少个不同的样本，只要它们都拥有相同数量的基因，那么它们的r_mean都将是一致的

由于上面的结果是在理论情况下推导出来的，实际上我们无法直接计算这个r，那么我们可以尝试通过其他方法来近似估计r， 只要这些近似统计量可以隐式地包含这一恒等关系即可

实际数据来证明

究竟为什么会有如此之大的差异？

可以从其公式上找到答案

这个统计量在2012年所发表的一篇讨论RPKM的文章（RPKM measure is inconsistent among samples. Wagner GP, Kin K, Lynch VJ. Theory Biosci. 2012.）中就被提出来了，称之为TPM —— Transcripts Per Million，它的计算是：

简单计算之后我们就可以发现TPM的均值是一个独立于样本之外的恒定值，它等于：

这个值刚好是r_mean的一百万倍，满足等价描述的关系。

参考资料：

(1) 孟浩巍《生物信息学100个基础问题 —— 第38题 当转录组普遍变化时RNA-Seq怎么进行分析(1)？》

(2) 孟浩巍《生物信息学100个基础问题 —— 第38题 当转录组普遍变化时RNA-Seq怎么进行分析(2)？》

(3) 【生信菜鸟团】quantile normalization到底对数据做了什么？

(4) Introduction to DGE

(5) 生信菜鸟团：StatQuest生物统计学专题 - library normalization进阶之edgeR的标准化方法

(6) 【】为什么说FPKM和RPKM都错了？