本發明涉及生物信息學領域，更特別地，涉及一種估算四倍體物種基因組的二倍化程度的方法。
背景技術：
四倍體基因組是正常二倍體基因組在自然情況或人工操作過程中，經自我加倍產生；或者兩個親緣關系相對較近物種通過自然或人工雜交并染色體加倍產生。前者往往稱之為同源四倍體，后者稱之為異源四倍體。四倍體物種的細胞核內包含組染色體，可以組成兩套二倍體基因組。在進化過程中，四倍體基因組中的兩套二倍體基因組會逐漸趨于形成更多的差異，向兩個方向發展，使得同源四倍體異源化最終成為二倍體基因組。這個過程可以稱之為四倍體基因組的二倍化。
當前基因組學研究中，并沒有專門針對四倍體基因組二倍化程度進行直接分析的方法。在對四倍體基因組進行組裝分析時，將四倍體基因組中兩套基因組分別或同時組裝出來，通過組裝出來的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％。
以上所述僅為本發明的較佳實施例，并不用以限制本發明，凡在本發明的精神和原則之內，所作的任何修改、等同替換、改進等，均應包含在本發明的保護范圍之內。