关于k-mer与基因组（组装）的那些事

引言

随着越来越多物种的基因组被测序发表，极大丰富了我们对物种起源、进化等方面的认知，同时也为各种物种的研究打下了基础。但是在面对复杂基因组时或者参考基因组并未发表的情况下，我们通常会对这个物种的基因组进行De Novo组装，在得到一个相对完整的基因组图谱以后，接着进行接下来的分析如基因功能预测、物种间进化关系等。在这个过程中，我们经常会遇到k-mer这个名词，然而这个抽象的名词是什么意思呢？它又有什么用呢？接下来，就随着小编一起去探究这k-mer背后的含义吧！

k-mer是什么？

简单来说，k-mer是一段长度为k的序列，而后面的mer即为monomeric unit（单体单元），也就是每个碱基。因k-mer包含k个碱基，若一段核酸序列长度为L，以一个碱基为步长滑动，那么根据这个核酸序列就可以得到L-k+1个k-mer；由于每个位点的碱基可以为（A、T、C、G）中的任意一个，因此k-mer理论上说有个不同的序列。原本一条长片段，就变成了很多短的片段，因此计算机处理的碱基数量也会增加很多倍。而且，每次取k-mer是同一条reads正反取两次，这就是对这条reads的反向互补序列再取一次k-mer。下面的图就形象化了这一过程，长度为15的序列，选取k-mer为5，那么就会得到11（15-5+1=11）个5-mer。

那这个时候估计有人要问了，为什么非要用k-mer不可呢，直接用reads进行组装不行么？答案是可以的，当然是基于reads没有测序错误情况下。而在当前的技术条件下，虽然二代测序的准确率已经达到了99.9%，如果测序量达到10亿个碱基的时候，这些错误的碱基数量会达到1000万个，这些错误累加起来还是会有很非常大的影响。所以，无论是序列的组装还是SNP的calling，如果不对这些错误进行处理，而是直接用原始reads进行组装，这些错误就会累加到最终的序列中，导致结果的不准确性。通过将reads切割成以k为单位的k-mer，由于测序错误具有随机性，这些由于测序错误生成的k-mer绝大多数都是原测序物种中不存在的k-mer，因此都只出现了1次，要是将这些k-mer去掉，那么就会较大的可能除去测序错误，从而使得组装的结果更加可靠。

我们用k-mer做什么？

在了解了k-mer是什么以及通过去掉低频率的k-mer能够使得组装结果更加准确以后，k-mer就没有别的用途了吗？当然不是！要知道，全基因组De Novo组装还是没有那么容易的。由于二代测序测序长度较短以及组装软件算法的局限性，对于基因组中有较高杂合度以及重复序列的时候，结果往往显得不那么准确，此时就会借助其他方法如三代单分子测序、Hi-C测序等构建高质量的基因组序列。如果没有对该物种的基因组特征有着初步分析，直接上来就开始那些烧钱的分析，结果不能说不好，但是可能会因为测序深度较低或测序样本有污染等使得后期分析有着各种各样的麻烦问题。所以利用k-mer对物种基因组特征进行评估是很有必要的一步，何况这些进行评估的reads，还可以进行接下来的组装分析。何乐而不为呢？

1. 评估物种基因组大小

如果测序得到的reads能够覆盖整个基因组，那么从reads取出来的k-mer就能遍历整个基因组，因此我们就可以通过k-mer粗略计算该物种的基因组大小。根据Lander Waterman的算法，基因组大小（G）= 总reads数 ÷ 期望测序深度（在不考虑测序错误、序列重复性的条件下，k-mer的深度分布遵循泊松分布，可以将深度分布曲线的峰值作为期望测序深度）:

其中，L为reads的平均长度，K为k-mer的长度；为所有k-mer的总个数，为k-mer频数的期望值；为测序reads总个数，为期望覆盖度。

我们应该看到，k和的选择都会影响基因组大小的估计产生影响，因此我们在选择k和时都不应过高或过低。最好的办法是通过深度曲线，选出最佳的k-mer期望值，从而大致估算出该物种基因组大小。

仔细的小伙伴们，可能会注意到，在深度分布曲线的最前面，线条先是下降到最低点，然后上升。这就是我们刚刚说到的由于测序错误，会有很多k-mer会以孤立的状态出现，即frequency为较低频数如1、2等。因这种k-mer在所有的k-mers中所占比例不小，所以它们的比例还是会相对较高。

2. 杂合度、重复序列比例评估

由前面可以看到，如果没有测序错误或者重复序列的情况下，其k-mer深度曲线会符合泊松分布。如上面所展示的图片中，该物种的k-mer深度曲线就只存在一个明显的主峰，说明该物种为纯合体或者单倍体，且基因组构造简单几乎不存在重复序列。但是在实际样本中，会存在大量的杂合和重复序列的情况，其k-mer曲线往往不会呈现出良好的泊松分布，可能会有多个峰存在的情况。如下图，根据物种基因组杂合程度以及重复序列的比例，可以将基因组分为高重复、高杂合、高重复和高杂合三种类型。

杂合峰：假设纯合TAG在3-mer中出现了c次即深度为c，而杂合的TTG中的A变成了T，就会形成TAG/TTG两种k-mer，此时每种对应的k-mer平均出现的次数为c/2即深度为c/2，且杂合度越高，这个峰越高。

重复峰：如果在x=c处出现主峰，x=2c处有一个次峰，说明有一部分片段出现的期望值是大部分的2倍，这些片段为重复片段，次峰为重复峰

3. 物种样品污染评估

k-mer还能被用于评估样品是否被污染。由于提取DNA的时候会受到很多方面的影响，如提取时的实验环境是否干净、提取物种的复杂性（微生物容易被外源微生物污染）等，导致提取到的DNA是不纯的。如不对该DNA进行纯性分析，后期也会遇到各种不知名的麻烦。

通常评估物种DNA样品的污染程度有两种方法，一种是通过k-mer深度分布曲线，主要观看是否有着峰异常的情况；另一种是通过GC含量分布图，查看图中是否存在多个密度集中的类群。

下图是在k-mer=15、17、19时分别作的k-mer深度分布曲线。我们可以看到，该测序样本在每一种k-mer选择下都呈现相同的趋势，即呈现双峰，且两个峰的横坐标之间也并非2倍的关系，由此可以初步判定该测序样本存在污染，至于是因为什么物种造成污染，则可以通过blast比对nr库进行简单的判断。

下图是对另外某物种进行GC含量分析的分布图。横坐标表示GC含量，纵坐标表示测序深度，右方是contig覆盖深度分布，上方是GC含量分布。可以看到在红色区域部分，该物种的GC分布呈现出两部分的区域分布，两部分的重心分别在深度20左右和40左右，即对应右边的contig覆盖分布，在20位置的小峰为杂合峰，在40位置的为纯合峰。对应红色的散点图，深度在20左右的红色区域是杂合contig的GC分布区域；再看上方的GC含量分布，主峰在35的位置，和我们计算得到的基因组的GC含量基本一致，且红色的散点也分布在GC含量在35%附近，说明该基因组没有受到其他外源物种的污染[1]。

4. 组装基因组

由于NGS测序技术的限制，其短读长不能使得NGS一次性测通一条染色体，折中的办法就是通过将长序列打断成100bp-300bp的reads，为了保证准确性，测序的时候基因组序列会被反复测很多遍，从而会生成以亿为单位reads数量。接着通过算法，利用reads之间的重复部分对基因组进行组装。然而就跟上文所说的一样，由于reads的测序错误会导致组装的结果不准确，所以许多科学家就想了很多算法来优化这一过程，De Bruijn graph就是最常用的算法之一（其他还有Overlap-Layout-Consensus，String Graph等）。De Bruijn graph就是通过将reads打断成k-mer后，利用k-mer之间的重复部分构建图，得到最优化路径从而拼接contig。

小结：通过以上对k-mer相关信息的介绍，可以看出在进行大规模测序前，简要的k-mer分析对研究物种基因组信息以及测序质量评估有着重要的作用。

k-mer为什么只能是奇数？

说了那么多使用k-mer分析的优点，好像忘了一个重要的点：k-mer怎么好像只有奇数呢？

是的，k-mer只能是奇数，就是为了防止通过k-mer组装时，正反链混淆。因为，当我们构建De Bruijn graph的时候，我们得准确把属于同一条read上的k-mer连接起来，我们不能一会儿把A kmer正确地连到它自己所在的read，一会儿又连到它互补链的read上去了！例如，4mer的GTAC反向互补后仍是GTAC，而3mer的GTA反向互补是什么？TAC。奇数的k-mer已经被证明是不能够匹配其反义互补链的，有兴趣的同学可以试着去证明哦！因此，我们得用奇数的k-mer。

一般来说，我们一般选用17-mer来估算基因组大小，因为ATCG四种不同的碱基组成长度为17的核苷酸有，足以覆盖一般大小的基因组；如果选择15的话，则只有1G，这就对正常的基因组来说可能覆盖度不够，导致估计不准。越长的k-mer片段越具很强的物种特异性，这也是毫无疑问的。例如，"ATGGACTTCGAGCTCCAA"肯定比"CGAGC"更具物种基因组特异性。当基因组中有较多的重复序列时，这时就可以使用较大的k-mer来跨过高重复的区域，从而获得更加准确及完整的基因组草图；由于reads上的碱基错误率的存在，选择较长的k-mer会带来较高的错误率，但这也可以加大测序深度来弥补。所以为了得到更加完整的基因组，要尽可能使用较长的k-mer用于组装。

k-mer常用的工具

• JELLYFISH[2]：该软件运用Hash表来存储数据，同时能多线程运行，具有速度快，内存消耗小的优点。
• KAT(The K-mer Analysis Toolkit)[3]: 包含多个工具来帮助用户通过使用k-mer对测序数据进行简单分析，如组装完整性、测序错误、是否有污染等。
• KmerGenie[4]：KmerGenie的最大优点在于可以实现在多个预设k-mer下的自动分析，除了进行常规的k-mer频数统计之外，还能够基于不同k-mer自动计算基因组大小，并为基因组组装评估一个最佳组装k-mer数值作为备选。

参考资料

[1]

诺禾致源survey分析: http://www.novogene.com/tech/service/Survey/result/

[2]

JELLYFISH: https://www.cbcb.umd.edu/software/jellyfish/

[3]

KAT: https://github.com/TGAC/KAT

[4]

KmerGenie: http://kmergenie.bx.psu.edu/