lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

基因组分析简介之K-mer分析

已有 17761 次阅读 2019-2-15 21:15 |个人分类:概述|系统分类:科研笔记| _K-mer分析

基因组分析简介之K-mer分析

 


为什么要进行K-mer分析


在对某物种进行全基因组测序时,若该物种没有已知的参考基因组信息(或者说该物种第一次被测序),那么我们若想详细得知该物种的基因组特征,就必须首先对测序序列进行拼接组装(基因组De novo),最终得到该物种的全基因组图谱后才可进行更精细的分析,如基因结构预测及功能注释等。

虽然说当今的测序成本也在逐年下降,全基因组测序的价格也越来越容易被广大实验室所接受,伴随着越来越多物种的全基因组图谱被公布。但是对于高等真核生物(特别是高等植物)来讲,它们的基因组一般比较复杂,通常伴随着较高的杂合度或者重复序列等,进行基因组De novo也并非一件容易的事情。特别是进行高精细基因组图谱的构建时,只基于常规DNA文库的二代测序根本无法满足需要,还需结合三代测序、BioNanoHi-C数据等才可获得高质量的基因组。相对而言,这时的测序成本也会是相当大的一笔开销,需要谨慎对待。因此在对高等真核生物进行全基因组De novo之前,我们需要设法提前获知该物种基因组特征的一些信息,为后续的测序方案、基因组组装方案、基因组结构注释等提供参考依据。

这种情况下,我们一般会在基因组大规模测序或者正式组装之前,首先构建DNA小片段文库进行中低深度的二代测序,使用PE文库测序所得的reads信息进行基因组Survey分析以初步评估基因组特征,这是很重要的一步(除非你们实验室很土豪不在乎,直接上三代测序,就可以当我没说了......)。基因组Survey分析的核心内容即为k-mer分析,使用k-mer分析可直接在相对较小数据量的二代测序reads水平上(直接使用原始测序reads数据,无需组装拼接),对该物种基因组大小、GC含量、杂合度高低以及重复序列比例等进行评估。

 


二代测序中的K-mer分析及基于K-mer的基因组特征评估简介


基因组二代测序中,首先对原始下机数据Raw reads进行过滤质控,得到Clean reads;之后再进行一系列的readsDuplication和纠错环节,以去除PCR扩增中带来的重复reads以及测序错误所引入的错误碱基后,这时的测序reads数据即可进行k-mer分析了。

那么,什么是k-mer分析?通过基因组二代测序的k-mer分析,我们可以获得哪些有用的信息呢?

 

何为K-mer


这里首先需要知道一个专有名词的概念,mer,其在分子生物学领域中意义为单体单元 monomeric unitmer)。通常用于核酸序列中的单位,代表nt或者bp,例如,100 mer DNA代表这段DNA序列单链长度100nt,或者双链长度100bp

k-mer则是指将核酸序列分成包含k个碱基的字符串,即从一段连续的核酸序列中迭代地选取长度为K个碱基的序列,若核酸序列长度为Lk-mer长度为K,那么可以得到L-K+1k-mers。如下图所示,假设这里存在某序列长度为21,设定选取的k-mer长度为7,则得到(21-7+1=15)个7-mers

1.png

 

二代测序reads中的K-mer统计计数及K-mer频数分布曲线


在测序序列中迭代截取k-mer,统计所出现的k-mer类型及各类型k-mer的出现频数,是进行k-mer分析的第一步。

对于二代高通量测序数据来讲,会得到许许多多的测序reads。每一条测序reads视为一条序列,并且测序reads也有一定的长度(如使用Illumina PE150测序平台测序,原始测序raw reads长度为150bp,过滤后的clean reads平均长度一般也有140bp左右),那么当我们在该二代测序数据所有reads的基础上依据特定长度截取k-mer后,会获得数量众多的k-mers

此时统计各k-mer及其出现的频数,即有多少k-mer片段只出现1次,有多少出现了2次,有多少出现了3次、4次、5......即可得到一个统计表,然后以出现频数为横坐标,以出现该频数的k-mer片段总数(或该数占总k-mer片段数的百分比)为纵坐标作图,即可得到某特定k-mer长度设定值下的所有k-mer的频数分布图。

 

如下所示,使用某物种的二代测序数据计算k-mer(选取k-mer长度17),最后可得到一个k-mer频数分布表(下图左图),第一列为k-mer深度,即各k-mer的出现频数;第二列为出现该频数的k-mer片段总数。下图右图为k-mer频数分布图,使用左图的统计表数据所绘制,图中横坐标为各k-mer的出现频数(Frequency),纵坐标为出现该频数的k-mer片段总数(Number)。

2.png                3.png

可以发现原始图中,最左侧(Frequency = 12等起始位置处)出现了很高的值,表明测序结果中存在大量的k-mer仅出现了1-2次,这个在k-mer频数统计表中也可轻易发现。这是因为在实际的二代测序数据中,由于测序错误(如Illumina测序平台的平均错误率约1%)的存在会引入许多带有错误碱基的reads,将这些reads打断成长度Kk-mer后,会产生许多错误的k-mer。由于测序错误带来的碱基类型是随机的,因此可知这些错误k-mer的出现频数很低,但总数目却非常的多。因此在上图中,低频数的k-mer数目占很大的比例,即在Frequency = 12等起始位置处出现很高的k-mer数目,使得图中曲线峰值很难分辨;为了增强曲线的可读性,可选择在作图时屏蔽掉曲线最左侧区域。当然也不排除一些真实的核酸序列,由于其碱基组成具有特异性且其只被测序测到了一次,将该序列截断为一定长度的k-mer之后这些k-mer只出现了唯一一次。但是相较于测序错误所产生的k-mer数量,后面这种情况所产生的k-mer数量基本上可忽略了,除非在很低深度的测序模式下。

此外,我们也可轻易看到,出现次数为几百上千次的k-mer数量其实很少。尽管在统计时不可丢弃这些出现频数很高但总体数量很少k-mer,但只是作图展示k-mer频数分布的话,是无需展示这些高频数深度的k-mer的,以便增强曲线的可读性(一些k-mer分析软件会统计至很高的k-mer频数深度,如10000,事实上在绘制k-mer曲线图时用不到这么多,视情况加以取舍)。

通常情况下,会考虑将低频数和高频数的数据屏蔽掉,屏蔽频数区间根据实际情况而定。屏蔽Frequency = 12等起始位置处以及Frequency > 5001000等高频深度的数据后,峰值即可呈现出,结果示例如下图所示(使用数据同上,只展示5 Frequency 500的区域)。此时,在不考虑测序错误率、基因组的杂合度和重复度的情况下,逐碱基取k-mer,则k-mer曲线在理想状态下服从泊松分布。

4.png 

 

基因组大小估算


上述我们获得了k-mer频数统计结果,接下来可以根据这个统计结果初步估算测序物种基因组特征。其中,k-mer分析估算基因组大小的原理如下。

reads中逐碱基取出的所有k-mer能够遍历整个基因组。根据Lander waterman算法,基因组大小(G)满足如下公式:

5.png

其中,Lreads平均长度,Kk-mer长度;knum为所有的k-mer总个数,kdepthk-mer频数的期望深度(即k-mer曲线中主峰对应的横坐标位置);bnum为测序reads覆盖碱基的总个数,bdepth为覆盖碱基的期望深度。

在这里,我们即可根据测序数据中的k-mer频数分布统计结果,大致估算出物种基因组大小了。

 

前面提到,因二代测序错误的存在会带来许多低出现频数的k-mer。这些由测序错误所引入的错误k-mer序列绝大多数是原测序物种基因组中所不存在的碱基组合,因此只出现了一次,也就是说它们基本都分布在k-mer频数Frequency = 1的位置。有的错误k-mer序列也可能会与原测序物种中已存在的一些k-mer序列恰好一致,或者导致出现了两个或多个相同的错误k-mer,使得k-mer计数结果与实际值有所偏差,如Frequency = 2等的位置处也出现了很高的数量等。

此时若我们忽略后面这种情况,只将k-mer频数Frequency = 1的情况认为是错误k-mer,并且忽略测序数据中确实只出现了唯一一次的真实k-mer(相较于测序错误所带来的只出现一次的k-mer,这些真实k-mer的数量可忽略,除了在很低深度的测序模式下),可据此大致计算测序错误率(Rerr),并用于修正基因组大小。

7.png

其中,Gcor为修正后的物种基因组大小,G为初始使用k-mer频数统计结果所估算出的物种基因组大小,knum为所有的k-mer总个数,knum_depth1为出现次数为1k-mer片段总数。

 

测序深度估算


估算出基因组大小后,即可使用测序总数据量(测序总碱基数)与估算出的基因组大小(基因组总碱基数)的比值来代表测序深度。

 

评估基因组杂合度及重复序列


首先可通过k-mer频数分布曲线图,直观地查看测序物种基因组的杂合度及重复序列含量情况。

上文提到,若忽略测序错误,且测序物种基因组中不存在杂合区或者重复序列的情况下,则其k-mer曲线在理想状态下服从泊松分布。例如上文中所展示的某物种,该物种的二代测序reads数据的k-mer频数分布曲线图中只存在一个明显的主峰,说明该物种个体为纯合体或单倍体,且其基因组构造简单几乎不存在重复序列。

而对于某些物种来讲,它们的基因组可能高度杂合,或者具有较高的重复序列比例。因此,其k-mer曲线往往不会呈现出良好的泊松分布。由于一定杂合度,会导致在主峰对应的横坐标的二分之一处出现杂合峰(如下图A所示);而一定的重复度,则会在主峰对应的横坐标的整数倍处出现重复峰(如下图B所示)。

9.png     10.png

11.png

A:在在x = a处出现主峰,x = 1/2 a处有一个次峰,说明部分片段出现的期望值是大部分的1/2,当序列有杂合时,包含杂合位点的k-mer因为分成了两部分,所以出现频率变为一半,次峰为杂合峰。

B:在x = a处出现主峰,x = 2a处有一个次峰,说明一部分片段出现的期望值是大部分的2倍,这些片段为重复片段,次峰为重复峰。

C:分别在x = ax = 2a处各出现一个峰,因两个峰高相差不大,两峰横坐标又是2倍关系,说明该个体高杂合或高重复。

 

那么这个时候,对于这些基因组构造较为复杂的物种来讲,我们还需额外对其基因组杂合度及重复序列含量进行评估。在不考虑测序错误、PCR扩增偏好性、测序深度等的情况下,一个简单的方法可基于k-mer曲线,计算杂合峰值与纯合峰值的比值,从而得到杂合率;计算标准泊松分布和实际数据曲线峰值后的面积差值,可得到重复序列百分比。使用软件计算时,由于不同的分析软件所涉及的算法可能更为复杂(如某些算法还需统计k-mer种类数,这个信息在k-mer频数分布表中是没有体现的,此时还需输入包含各k-mer类型及其数量的更原始信息),会导致不同软件的计算结果会略有差异,不再多说。

 

评估测序质量(错误率)


前文在进行基因组大小估算时提到,可将出现次数为1k-mer大致认为是错误k-mer,并且可据此大致估算测序错误率(详见上文),以评估测序质量。

 

此外,根据k-mer曲线也可以初步评估测序质量。

由于测序错误的存在产生了许多只出现了1次、2次等的低频数k-mer片段,因此可知正常情况下的k-mer曲线,其左侧区间应迅速下降至低谷后,再逐渐上升达到峰值(可参见上述任意一张展示图)。若k-mer曲线左侧未完全下降至低谷即出现上升幅度,此时需要考虑两个原因。

1)选取的k-mer长度过长。k-mer长度越长,则截取后所得k-mer片段中出现错误碱基的概率也越高,即包含错误碱基的低频数的k-mer数量越多,特别是在较低深度的测序模式下。为了避免这个情况,我们一般不会选择较长的k-mer进行计算,通常选择长度大小为1517192123等的短k-mer即可满足需求(简单的微生物基因组可选1517等,重复序列较多的动植物基因组可选2123等)。当然,不排除测序数据本身质量较差的情况存在。

2)若在k-mer = 17等短k-mer长度的分析中,仍然出现这种情形,则说明本次测序质量确实较差,测序碱基平均错误率偏高,此时的测序数据不建议使用。

 

如下所示,此处存在一个测序物种(已确定该物种基因组无杂合区,极低重复区,即为简单基因组)的DNA样本,已知该样本的DNA质量不合格,因此导致了其测序数据的质量很差。使用该低质量的测序数据进行k-mer分析后,我们可以较直观地在k-mer曲线(下图左)中看到,本次测序数据的错误率较高。主要依据大致有三点:

1)其k-mer曲线左侧未下降至平缓便出现上升趋势,原因已在上述说明。

2)事实上,该k-mer曲线为屏蔽了k-mer频数Frequency < 20的区域后所绘制。正常情况下是无需屏蔽这么多的区域即可观测到明显峰的存在的,但这里情况特殊,比方说,在只屏蔽了k-mer频数Frequency < 10的区域时根本无法在图中看到曲线“弧形突起”(下图右,使用与左图相同的数据),因此屏蔽了更多的低频数k-mer区域进行作图展示。也就是说,该k-mer曲线左侧下降较为缓慢,即统计结果中存在极其多数量的低频数k-mer。这种情况也有可能是测序数据质量偏低所致(注,只是有可能,也可能是该物种基因组本身特征所致,因此仅作为参考)。

3k-mer曲线显得“怪异”,如“弧形突起”区域的整体跨度较大且略为平坦,曲线峰值不明显;且在细节上曲线存在上下波动,不平滑。这种情况也有可能是测序数据质量偏低所致(注,只是有可能,也可能是该物种基因组本身特征所致,因此仅作为参考)。

12.png     13.png

 

检查样本纯度(是否存在较多其他物种DNA的污染)


使用k-mer频数分布曲线图还可初步判定样本DNA是否受到了其它物种DNA的污染,若k-mer曲线显得不正常(如多峰情况),这时候就需要考虑这种情况了。特别是在微生物物种的测序样本中,因为它们在培养时很容易被污染,特别是有时因菌落外形几乎一致还不易区分,这时需多加注意。

如下所示,使用某细菌物种的二代测序数据所计算(选取k-mer = 151719分别用于分析,已屏蔽Frequency < 5的区域)。已知该测序样本不纯,归因于单菌培养时不慎混入了另一杂菌,使得提取的基因组DNA并非来自于单一细菌物种,而是来自两个细菌物种基因组DNA的混合样本。可以看到该测序样本的k-mer曲线呈现双峰,且两个峰值的横坐标的数值关系也非“2的关系(上述提到,杂合峰横坐标数值约为主峰的1/2重复峰横坐标数值约为主峰的2倍,此处均不对应)。况且对于细菌物种来讲,由于其基因组结构简单,不可能存在杂合或高重复的情况,其k-mer曲线中理论上也只应该存在一个明显的主峰。因此可初步判断该测序样本存在污染,测序reads包含了两个物种的DNA序列。

14.png

 

区分物种或判断物种间亲缘关系远近


我们知道,每个物种的基因组各具特色。在k-mer分析中,我们将测序物种的基因组截成一段段k-mer,因此可知不同物种的k-mer是很不同的,而同一物种的不同个体之间的k-mer是差别不大的,基于这个原理可以基于k-mer区分物种。当然,上述根据k-mer曲线判断样本是否受到污染也差不多可以认为是基于这个原理。

在二代测序中,若不考虑测序错误率、PCR扩增偏好性等情况,即认为基因组中每个碱基被测到的概率是等同的,那么可以得知,对于每个物种的基因组测序数据,在某特定长度的k-mer下,其k-mer曲线肯定都是在一段特定区间内出现峰值。

同理,不同物种间进化关系越近,则它们的基因组组成就越为相似,反之越远。若两个物种具有相近的亲缘关系,则两个物种的k-mer曲线肯定也是较为一致的;反正,若两个物种亲缘关系较远,则两个物种的k-mer曲线肯定相差很大。

 

越长的k-mer片段越具很强的物种特异性,这也是毫无疑问的。例如,“ATCGGTCTCAGCGCGCAAACT”肯定比“GCGC”更具物种基因组特异性。虽然我们在k-mer分析中使用了短k-mer进行了分析,但接下来在基于k-mer原理的基因组二代组装中,为了得到更完整的基因组,会尽可能使用较长的k-mer用于组装。尽管由于二代测序错误率的存在,选择较长的k-mer会带来较高的错误率,但可以使用加大测序深度来弥补。此时,在高深度的二代测序模式下使用相对较长的k-mer进行组装,可以得到更准确更完整的基因组草图。

 


附:K-mer分析常用工具


JELLYFISH

KmerGenie

GCE




http://blog.sciencenet.cn/blog-3406804-1162384.html

上一篇:K-mer分析及基因组特征评估软件GCE的安装及使用说明
下一篇:R包Tax4Fun基于16S物种丰度预测细菌群落代谢功能及统计作图的一个简单实例

2 熊星辰 张成岗

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-12 18:34

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部