从metaWRAP quant_bins计算模块理解宏基因组分箱bin的丰度计算

背景

在进行扩增子分析时，我们拿到的最关键的一个中间数据就是OTU/ASV表，在这个矩阵中，我们能获得我们的分析对象（OTU/ASV）在样本间的分布规律，并通过微生物群落的结构在样本之间的差异来解决一系列的科学问题。其中，我们常常可以通过OTU/ASV在不同样本间的共现关系得到它们之间的相关性关系。

在WGS支持下的宏基因组测序分析，通过组装、分箱等步骤可以获得metagenome assembled genomes（MAGs），每一个MAG常被认作一个单菌，我们选取高质量（完整度高、污染率低）的MAG可以进行很多分析，其中也包括不同的MAG的共现或相关性规律的探寻，为此，需要研究一种MAG（bin）的丰度计算方法。

用于宏基因组数据分箱的metaWRAP pipeline在分箱后提供了计算相对丰度的模块quant_bins，我也试过一次，不过一直还对这个里面给出来的结果有一些疑问。这个模块使用的软件是Salmon，一般来说在转录组定量分析中经常使用。为了记录这一部分的计算方法，写一篇学习笔记以供和大家交流，可能有不对的地方可以留言指出，一起讨论。

从RNA-seq常用定量指标聊起

RNA-seq：Read Count与定量

RNA-seq（转录组测序技术）常用于通过对各种RNA（如mRNA、小型RNA等）定量进而对基因表达水平或差异表达等进行具体的分析，这一部分我没有过多的研究，不细说。主要说一下关于定量方法的问题。

我们回到OTU表的计算，实际上常见的基于NGS的扩增子分析得到的OTU的表中的数据就是Read Count，即测到的序列（Read）我将它分到哪个OTU，该OTU在该样本下的值就+1。可以这么做的前提是，每一条Read都能被assign到单独的一个分析单元（即Operational Taxonomic Units）。为了进行样本间的比较与统计检验，我们通常需要对OTU表进行标准化，这种标准化是按照样本测序深度进行的，也就是说我需要保证样本与样本之间的总Read Count是一致的，这样其比较才是一致的，这非常重要。这一步通常是通过稀释化方法（Rarefaction）与按标准重抽（Resample by Standard）进行的。

那么，在RNA-seq或者在MAG的定量中，可以直接使用Read Count进行操作吗？一般来说，需要思考下面两个问题：

1. 假定我们认为RNA-seq中，测到能对比上gene A的碱基越多，则该片段对应的基因的表达量越高，那么我们可以认为mapping到gene A上的Reads越多，该基因的表达量越高吗？
2. 假设有两组样品（比如是wild type和treatment）里都需要检测同一基因gene B，那么gene B在两个样品中的mapped read count可以用来比较这个基因在两个样本间的表达量差异吗？

以上两个问题的答案当然都是否定的。对于第1个问题，与扩增子测序不同，shotgun技术将片段打断后，对于指定片段的测序概率和片段在总体中的分布频率是有关的，简单来说就是长度较长的gene，测序时测到属于它的片段的概率就更高，因此对于这一类gene，它的mapped read count多并不等同于它的表达量高；对于第2个问题，实际上就是样本测序深度的问题，即两个样本的深度需要保证一致才能将两组数据的同一基因的丰度进行比较。

因此，我们常需要对得到的Read Count进行一些计算，来方便更科学地比较基因表达差异。在我们比较MAGs的时候，类似地，也需要考虑这些因素——MAG size与sample size。

标准化参数 RPKM/FPKM

RPKM：Reads Per Kilobase per Million mapped reads
FPKM：Fragments Per Kilobase per Million mapped reads

在双端测序中，两端Reads（不管测没测通，联合insert size一起）组成了一个fragment，有时也用fragment作为计数指标计算，得到FPKM。

对于单端（Single-End, SE）测序来说，常用RPKM；对于双端（Paired-End, PE）测序来说，常用FPKM。

有的地方会说RPKM是FPKM值的两倍。这个说法我认为并不完全正确。对于PE测序来说，如果一对reads都比对上一个gene，那么该对reads算一个fragment，但如果仅有一端的read比对上了gene，那么应该仅将这一个read算作fragment，另外还应考虑双端read总长与文库长度的关系，因此我认为RPKM与FPKM的数量关系并不是绝对的2倍。
另外，很多软件对这两种参数采取的计算方法也不一样。

TPM

RPKM解决了之前提出的问题中的第一项，即对基因长度进行标准化。如果我们仅想做单样本内genes之间的比较，那么RPKM就已经达到目的了。但是很多情况下，我们还需要对同一基因的表达量在不同样本中的情况进行分析，那么久也需要对sample size进行标准化，于是引入了TPM的概念：

TPM：Transcripts Per Million

其计算公式可以写为：

从公式可以看出，TPM是下对基因长度按照kbp进行标准化，然后对所有样本的总Reads数求指定样本所占的比例。如果将RPKM公式中的Total No. of mapped reads按求和符号提公因式的话，就可以把RPKM的计算式转化为TPM的公式：

可以知道，TPM实际上是在1Mbp的标准化前提下，RPKM的比例。TPM表中，所有样本中每一样本的所有gene的TPM值加和应该都是10^6，即完成了两项标准化。

将RNA-seq问题转化宏基因组bin quantify问题

metaWRAP pipeline提供了quant_bins模块计算bin的相对丰度（relative abundance）。作者将其计算出的数值称为CPM（Copies Per Million），与RNA-seq中的TPM类似。简单来说，就是每个bin的Read平均覆盖度按照10^6bp（1Mbp）的标准化。

不过，这种计算方法与TPM也并不完全相同。在Github上有几个和这个模块有关的issues，作者也参与了一些讨论，我节选了一些作者对这个模块的解释说明，然后自己理解一下。

The bin.9 10.50602416 value means that on average, the contigs in bin.9 have a read coverage of ~10. So in your example, the total estimated number of reads mapping to your bin is 1073567*10.50602416/100 = 113,755 reads (I used 100bp for read length because I assumed it was a paired-end library). So of your total library size of 23,887,914 reads, at least 0.5% came from this organism. However, if you goal is to know the bin’s “abundance”, using the original read depth as a representation of the abundance (10.50) is probably more accurate. Because you may have only a part of the whole genome, read depth is a better metric. Just dont forget to standardize them if you start comparing samples. (https://github.com/bxlab/metaWRAP/issues/40)

这里是一个例子，即某一个bin经过quant_bins计算后得到其abundance是10点几，即说明这个bin中的contig的平均read覆盖度约为10。
因此，可以知道能够map到该bin上所有的reads的估计应该是bin size * abundance / read length，如这里例子的bin的长度就是1073567bp，然后用的双端50bp的测序。因此，如果总的Read数是23,887,914的话，那么就至少有113,755/23,887,914约等于0.47%的Read是来自这个bin的。从这个角度上来说，这个百分比数据也是可以用的（roughly）。
但是后续作者也指出这个数据的问题，即例子中这个bin的完整度仅有78%，也就是说，理论上还有更多的read也是来自这个bin（那22%未知的部分），只是我们没法知道了。所以我认为，仅有refine bins的过程中，仅选取高质量，尤其是完整度高的bins作为后续分析的基础时，使用这个百分比才具有一定的参考意义。

The values reported by the Quant_bins module are essentially estimated average read coverage values for each bin in each sample, standardized to the number of reads in each sample. So they do not have to add up to any particular number. (https://github.com/bxlab/metaWRAP/issues/67)

实际使用quant_bins模块的时候，很容易发现同一批bin在不同样品间计算丰度时，样本下总数并不是加起来等于某一固定数值，这和TPM在样本间的总和均加起来等于1M不同。因此有人对这件事提出了质疑，作者给出了如下的解释：

To put it simply, the total abundance of the bins absolutely do NOT have to be the same in each sample. This is very different from something like RNAseq gene expression values, because we cannot reliably reconstruct all the bins from all the samples. Because assembly and binning biases vary between samples, the total of bin (and contig) abundances can be different.

这件事是非常合理的，因为得到bin的过程是非常不固定的。在扩增子分析中，我们必须保证测序是完全的，体现为稀释曲线需要达到平台，也就是测到的所有序列能够代表整体。但是bin的分析里，Reads是无法代表所有bin的，或者说我们得到的bin set加起来并不是一个整体，这么来说，它们的abundance加和得到定值这件事是没有意义的。
我们称这里的值为Relative Abundance，因为我们并不能知道真正的abundance，因为我们也不知道样本之间的biomass difference。

参考

https://github.com/bxlab/metaWRAP/issues

生物信息学100个基础问题 —— 第35题 RNA-Seq 数据的定量之RPKM和FPKM

Garber M, Grabherr MG, Guttman M, Trapnell C. Computational methods for transcriptome annotation and quantification using RNA-seq. Nat Methods. 2011 Jun;8(6):469-77. doi: 10.1038/nmeth.1613. Epub 2011 May 27. Erratum in: Nat Methods. 2011 Jun;8(6):following 477. PMID: 21623353.