本发明涉及生物信息学领域，更特别地，涉及一种估算四倍体物种基因组的二倍化程度的方法。
背景技术：
四倍体基因组是正常二倍体基因组在自然情况或人工操作过程中，经自我加倍产生；或者两个亲缘关系相对较近物种通过自然或人工杂交并染色体加倍产生。前者往往称之为同源四倍体，后者称之为异源四倍体。四倍体物种的细胞核内包含组染色体，可以组成两套二倍体基因组。在进化过程中，四倍体基因组中的两套二倍体基因组会逐渐趋于形成更多的差异，向两个方向发展，使得同源四倍体异源化最终成为二倍体基因组。这个过程可以称之为四倍体基因组的二倍化。
当前基因组学研究中，并没有专门针对四倍体基因组二倍化程度进行直接分析的方法。在对四倍体基因组进行组装分析时，将四倍体基因组中两套基因组分别或同时组装出来，通过组装出来的contigs或染色体水平的基因组内互相比较估算两套基因组间的差异。如果两套基因组间差异相对较大，处于高度二倍化状态时，这种方法可以得到相对较好的结果。但由于四倍体基因组的特性，如果两套基因组间差异相对较小，那么往往不能得到理想的组装效果，不能分开两套基因组，因此这套方法对于同源四倍体效果较差。如果采用将两套基因组分别组装的方案，首先需要找到两套基因组对应的亲本物种，随后分别组装。由于异源四倍体是由两亲本物种杂交后又历经进化过程产生，亲本物种的基因组与异源四倍体中的基因组存在差异，这套方法会有一定程度的误差。
因此，需要一种直接估算四倍体基因组的二倍化程度并进行量化的方法。
技术实现要素：
为解决以上问题，本发明提供了一种估算四倍体物种基因组的二倍化程度的方法，其特征在于，包括以下步骤：
s1：对所述四倍体物种的基因组进行二代测序，得到四倍体基因组测序数据；
s2：将所述四倍体基因组测序数据与二倍体基因组测序数据进行比较，估算所述四倍体物种基因组的二倍化程度，所述二倍体基因组测序数据为所述四倍体物种的近缘二倍体物种基因组的二代测序数据。
本发明以分析四倍体基因组二倍化程度为目标的直接分析方法，首次提出了二倍化率的概念用于量化四倍体基因组的二倍化程度，不依赖于目标的基因组序列，具有成本低，速度快，成功率高等优势。
在一个实施方案中，s1和s2中所述二代测序为illumina测序。
在另一个实施方案中，所述二倍体基因组测序数据通过测序得到，或为已有的测序数据。
在另一个实施方案中，所述四倍体基因组测序数据的测序深度不小于100x。
在另一个实施方案中，所述二倍体基因组测序数据的测序深度不小于30x。
在另一个实施方案中，s2包括：
s21：获得所述二倍体基因组测序数据；
s22：对所述四倍体基因组测序数据和所述二倍体基因组测序数据进行分析处理，分别得到四倍体基因组k-mer集合和二倍体基因组k-mer集合；
s23：分别统计备所述四倍体基因组k-mer集合和所述二倍体基因组k-mer集合中的k-mer总数，并以k-mer的出现频数为横坐标，k-mer的种类数纵坐标分别制备所述四倍体基因组k-mer集合的k-mer种类数频数分布图和所述二倍体基因组k-mer集合的k-mer种类数频数分布图(例如分别统计在四倍体基因组和二倍体基因组中出现频率在1到1000次的k-mer的种类数)，并以所述k-mer种类数频数分布图中的第一波谷前k-mer为错误k-mer；
s24：根据所述四倍体基因组k-mer集合的k-mer种类数频数分布图和所述二倍体基因组k-mer集合的k-mer种类数频数分布图分别计算所述四倍体基因组的序列重复率和杂合度，以及所述二倍体基因组的序列重复率；
s25：根据所述四倍体基因组的序列重复率和杂合度，以及所述二倍体基因组的序列重复率计算所述四倍体基因组的二倍化率，计算公式如下：
公式ii：
d：四倍体基因组二倍化率
α：四倍体基因组序列重复率
β：二倍体基因组序列重复率
k：四倍体基因组杂合度。
基于基因组k-mer的频数以及种类数的分布估算基因组特征。该方法仅需要一定覆盖度的基本二代高通量数据就可以完成，且该过程对测序文库的种类及插入片段类型无要求且不需要进行基因组组装，所以受基因组复杂度影响极小。相比较，通过基因组序列估算方法受限于基因组序列，尤其是高度复杂的同源四倍体基因组序列的组装目前是一个世界性难题。往往需要构建多种测序文库，采用各种方法测序，并制定复杂的组装策略，而往往得不到质量足够好的基因组序列进行后续分析。
在另一个实施方案中，s24通过以下方法计算所述四倍体基因组的序列重复率和所述二倍体基因组的序列重复率：在所述四倍体基因组k-mer集合的k-mer种类数频数分布图中，以第一杂合峰2x处为主峰位置，以主峰1.8x处为界限，出现频率大于该界限的k-mer为所述四倍体基因组的重复k-mer；所述二倍体基因组k-mer集合的k-mer种类数频数分布图中，以主峰后1.8x处为界限，出现频率大于该界限的k-mers为重复k-mer，并根据公式i分别计算所述四倍体基因组的序列重复频率和所述二倍体基因组的序列重复频率
公式i：
r：基因组序列重复率
nkspecies：非重复k-mer种类数
nkfrequency：非重复k-mer频数
ekmer：错误k-mer数
akmer：总k-mer数；
并以所述四倍体基因组k-mer集合的k-mer种类数频数分布图中第一杂合峰计算其基因组杂合度。基因组杂合度的计算方法是现有技术，在本发明中本着突出重点的原则不做赘述。
在另一个实施方案中，s22中，通过以下方法处理所述四倍体基因组测序数据和所述二倍体基因组测序数据：
s221：过滤掉所述四倍体基因组测序数据和所述二倍体基因组测序数据中低质量碱基和/或短于一定长度的读序；
s222：将所述四倍体基因组测序数据和所述二倍体基因组测序数据分化成k-mer，分别得到四倍体基因组k-mer集合和二倍体基因组k-mer集合。
进一步地，s221中，所述低质量碱基读序为序列两端质量值小于20的读序，所述短于一定长度的读序为序列总长小于50的读序。
附图说明
图1为本发明方法的流程图。
具体实施方式
以下结合实例对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。
本发明的方法的流程示意图如图1所示。
以某种鱼类(因论文未发表，为保密起见，暂时不公开其种属)为例，该鱼类在自然种群中包括四倍体基因组类型和二倍体基因组类型。其四倍体基因组预估大小约2.4g，二倍体基因组预估大小约1.2g。我们基于本方法估算该四倍体基因组物种的两套二倍体基因组间是否发生了分化以及分化程度有多少。具体实施过程如下：
1)分别对两个物种进行建库测序，均建立插入片段300-350bp的illuminahiseq文库并进行pe150高通量测序。四倍体物种总共测得约280g数据，测序深度约117x，二倍体物种总共测得约52g数据，测序深度约43x。
2)使用htqc软件的ht-trim和ht-filter模块采用默认参数分别对两组数据进行碱基质量过滤和读序过滤，约过滤掉0.02％的数据，整体测序深度不变。
3)使用jellyfishcount首先以k-mer＝17计算两组数据的所有k-mer类型及频数；使用jellyfishstats统计并获得两数据的k-mer总数，四倍体数据有250,217,368,293个k-mer，二倍体数据有46,513,565,383个k-mer；使用jellyfishhisto绘制以出现频率为横坐标，k-mer的种类数为纵坐标统计k-mer种类数频数分布。
4)二倍体k-mer种类数频数分布中主峰位于k-mer频数为40处，因此频数大于40×1.8＝72的k-mer为重复k-mer，重复k-mer有2,3765,979,213个。其第一波谷为频数＝7处，因此错误k-mer有976,614,636个。通过公式(1)计算得其基因组序列重复约52％。
5)四倍体k-mer种类数频数分布中第一杂合峰位于k-mer频数46处，因此频数大于46×2×1.8≈166的k-mer为重复k-mer，重复k-mer有202,733,541,319个。其第一波谷为频数＝11处，因此错误k-mer有7,879,758,786个。通过公式计算的其基因组序列重复约84％。同过其第一杂合峰处的k-mer数计算可得杂合度约1.1％。
6)同过公式(2)代入以上计算得到的两基因组重复序列含量以及四倍体基因组杂合度可得该四倍体物种基因组二倍化率约32％。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。