count normalization的方法
这应该会是一个小系列，我会简单介绍下现在的一些count normaliztion的方法，以及对应的代码解释、思考，从而最后得出我所认为比较适合的count normaliztion的方法。
这个部分我会介绍下Count nromaliztion需要考虑的因素、以及一些方法
Count normaliztion需要考虑的因素
首先我们需要明确的一点是如果你的Seq count如果没有经过矫正的话，那么不管是组内还是组间的基因比较，都是不可取的。所以，为了让Gene count能够在组间或者组内进行比较，我们需要考虑一些因素，然后根据这些因素来对我们的count进行normalize，从而来让组间或者组内的基因比较变得更加地合理。
Sequencing depth（测序深度）
image
我们在图中可以看到，Sample A中的Gene reads counts都基本上是Sample B中的reads counts的两倍，但这并不意味着Sample A中的Gene表达量就要比Sample B中的高，因为Sample A自身的测序深度就要比Sample B高（通俗来说，就是Sample A可能测了10G，而Sample B测了5G）
Gene length（基因长度）
image
Gene length在进行组内比较的时候，也是一个重要的因素。可能实际上Gene X跟Gene Y的分子数是差不多的，但因为 Gene X的长度更长，所以其对应的counts就越多。
RNA composition
image
RNA composition也是一个重要的因素。以上图为例，我们可以看到Sample A的total reads counts是远大于Sample B的，所以我们就可以考虑Sample A和B都除以其total reads counts，这样得到的就是每个基因所占全部基因的比例了（其实就是CPM，还是比较合理的 ：））。
就像你的资产是1w人民币，而你美国朋友的资产是2k美元， 你想比较下谁更有钱，但直接比是不太恰当的。这时候就可以考虑你的资产除以中国全部的资产，你美国朋友的资产除以美国全部的资产。比如你是1%，而你美国朋友是5%，那我们就有理由相信他更有钱一点。
但这样其实有个很大的问题就是，Sample A和Sample B的RNA composition是非常不一样的，Sample A中的Gene DE占据了total counts中的大部分，这样就会压缩Gene X、Y、Z的比例，从而让X、Y、Z无辜地 “被” 降低了表达量。因为我们算的是"the proportion of the total count for each gene"，proportion总和加起来肯定是1，DE如果占据了大部分，那么对应的剩下就肯定会“少了”。
这个问题其实也引出了TPM、FPKM、CPM等根据proportion矫正的方法的一些局限性，这一点我会在后面思考的部分提到。
当然，我还见过有考虑测序文库的reads length因素的……但我不知道为啥要考虑这个
三张图的来源是：Introduction to DGE
其实在考虑normaliztion factor的时候，还有一个很有意思的点，即我们其实真正关心的是molecule counts，但我们几乎是无法知道Gene真正的molecule counts的。对应的，我们只能对Gene molecular进行破碎建库，然后进行抽样测序，用reads count来表征对应的Gene molecular counts（其实这也就是用样本分布来代表总体分布）。在这其中，除了我们上述提到的factor之外，还会有诸多的问题，包括Gene与isoform的关系，Reads multi-mapped，library sequence complexity等等问题。
Count nromaliztion所用的方法
这里我用一张表来先来总结下
| Normalization method | Description | Accounted factors | Recommendations for use |
|---|---|---|---|
| CPM (counts per million) | counts scaled by total number of reads | sequencing depth | gene count comparisons between replicates of the same sample group; NOT for within sample comparisons or DE analysis |
| TPM (transcripts per kilobase million) | counts per length of transcript (kb) per million reads mapped | sequencing depth and gene length | gene count comparisons within a sample or between samples of the same sample group; NOT for DE analysis |
| RPKM/FPKM (reads/fragments per kilobase of exon per million reads/fragments mapped) | similar to TPM | sequencing depth and gene length | gene count comparisons between genes within a sample; NOT for between sample comparisons or DE analysis |
| DESeq2’s median of ratios [1] | counts divided by sample-specific size factors determined by median ratio of gene counts relative to geometric mean per gene | sequencing depth and RNA composition | gene count comparisons between samples and for DE analysis; NOT for within sample comparisons |
| EdgeR’s trimmed mean of M values (TMM) [2] | uses a weighted trimmed mean of the log expression ratios between samples | sequencing depth, RNA composition, and gene length | gene count comparisons between and within samples and for DE analysis |
关于TMM能根据gene length矫正，从而进行组内基因的比较，我好像还没听说过……不过可能是我没怎么细看吧
CPM
公式：
是对应基因的count值，而N则是你整个文库count。之所以还有乘上
只是单纯的有些比例太小了，乘上
CPM (counts per million，有时候也叫RPM, Reads per million mapped reads) 其实就是单纯根据测序深度来进行矫正，但我个人觉得“根据测序深度矫正“解释不好。我更偏爱的解释是CPM展示的是这个基因在这个样本测序文库中的比例，也只有这个解释，才可以帮助我们更好地理解后面的TPM、FPKM。
CPM 这里还有一个问题就是，你如何定义N这个值。在普通的RNA-Seq中，我们可以认为 N 是你基因count数目的总和，亦或是mapped reads，（虽然这也是值得商榷的，因为如果你的基因组组装的稀烂，或者注释稀烂的话，那么不管是基因count数目的总和，或者mapped reads的总和，都不能充分定量你文库中的分子数目）。
在 ChIP-Seq 或者 ATAC-Seq 中，这个N更是一个很难定义的问题，你是数你Peak所含的count呢，还是数整个bam文件里面的count呢。
TPM
公式：
则是基因的长度。
这个则是 每个基因count除对应基因长度 的值的总和。
我上面在CPM中说过，我们可以通过 基因在这个样本测序文库中的比例 来对不同文库的测序深度进行矫正，但有一个问题是 如果基因A是30k，count是100，而基因B是300bp，count是50，我们能说基因A比基因B表达量高么，那肯定是不能的。因为基因A的长度是基因B的100倍，理论上来说count应该也是100倍。但现在 基因A的count 只是 基因B 的2倍而已。所以我们在考虑 基因在这个样本测序文库中的比例 也应该考虑基因长度这个问题。
用我之前那个中国和美国人资产比较例子的话，就是 我们家有100个人，年薪 20w RMB，而美国小伙伴家里是2个人，年薪2k美元。我如果想要比较 我们家 和 美国小伙伴家 的赚钱能力的话，应该是 用我们家的人均年薪 除以 中国每个家庭的人均年薪的总和，比上 对应美国小伙伴的。
但这里有一个问题就是我们怎么界定基因长度，是整个基因的长度，还是外显子的长度，还是各种isoform长度的累加。其实这也是不同 数count软件 在算法上的区别。
hhhh，对应到家庭那个例子的话，就是 我们家有些人是不赚钱的（有些isoform是不表达的）……
FPKM
公式：
这里之所以要乘
是因为
这个值太小了而已……之所以是
可以看到，FPKM和TPM其实是差不多的，唯一的区别就是一个除的是N，另一个除的是 。所以转换也比较方便
从 基因在这个样本测序文库中的比例 这个角度来看，自然是TPM更加的合理一点，因为TPM才真正做到从百分比的角度来进行矫正。
但其实还有一个问题，即如果两个样本的RNA composition很不一样的话，TPM 真的是一个合理的样本间比较方法么？
参考：
DESeq2 normalized count
这个是DESeq2自己的count矫正方法，主要是为了矫正不同文库的深度以及RNA组成，从而使得大部分基因在样本之间保持不变，本质上就是为每个样本计算一个size Factor，从而得到normalize count，进行后续的差异分析。就像下图一样，如果我们只是根据样本文库深度进行矫正的话（蓝线），那么所有基因都会有差异，如果我们根据 “大部分基因保持不变” 的原则来进行矫正（红线），那么就只会是 点C 有差异了。
image
尽管常见的normalized方法是每个样本得到一个size factor，但其实你也可以对每个样本里的每个基因都指定一个size factor
大部分基因在样本之间保持不变这个假设前提在大部分实验中都是成立的，但有时候也会有例外。比如说在正常组织和肿瘤样本中，肿瘤样本的基因就会有极大的改变，从而使得不符合这个假设前提。其实哪怕你是不同组织之间比较，也是不符合这个假设的，因为每个组织是各司其职，都有其主要表达的基因簇，你拿心脏和脑组织的RNA做差异分析，结果肯定是有很多很多的差异基因的。还有就是你拿DESeq2去做ChIP-Seq的差异分析的时候，可能你的对照组是野生型，而你的处理组是组蛋白酶缺失的突变体，那么其实组蛋白水平是整体降低的，这也不符合大部分保持不变的原则。
但话说回来，怎么样才算 大部分呢？超过50%就算大部分么？回答就是我也不知道╮（╯＿╰）╭
其实DESeq2对于这种情况应该还是比较robust的……不然也不会广泛使用了，尤其是在各种癌症组织中。
图来自：
edgeR的TMM
TMM应该是也是基于大部分基因保持不变的原则，但TMM应该是更加的复杂，有更多的参数，我也更加地不懂……
网友评论