本發明涉及分子生物學技術領域，尤其涉及一種基于基因拷貝數變異位點的基因型分型方法，更具體涉及一種基于基因拷貝數變異位點對林木自然群體的基因型的分型方法。

背景技術：

拷貝數變異(Copy Number Variations,CNVs)是指不同個體基因組中長度范圍從1kb至數Mb的DNA區段，與參考基因組比較而發生拷貝數變異的亞微觀染色體結構變異。在人類中，CNV已被證實與許多復雜疾病相關。而在植物基因組中，由于CNV的存在影響了植株的表型性狀、代謝過程、生理過程和適應性進化等。因此，對群體中基因拷貝數變異進行全面研究顯得尤為重要。盡管在過去的幾年里，高通量測序平臺的問世，使得大量檢測CNV的方法和技術也陸續出現和不斷完善，但對于精確的發現和檢測CNV，尤其是對CNV位點基因型檢測的準確計算仍具有強烈的挑戰。

目前，利用高通量測序結果對CNV位點基因型進行分型大多是基于CNV位點附近的SNP位點的基因型進行的分型，該方法操作復雜、結果分析困難，不能精確的檢測CNV位點的基因型，使其應用受到很大的局限。此外，在植物領域，尤其是在林木自然群體中對CNV位點基因型進行檢測的方法研究更是空白。因此，研究植物的CNV分型方法對培育高產、優質、抗病性強等優良品種具有深遠的意義。

技術實現要素：

鑒于此，本發明的目的是提供一種適用于不同群體大小的林木(或植物)CNV位點的基因型分型方法，操作簡便易行，并且較精確地發現和檢測CNV位點的基因型。

為了實現上述發明目的，本發明提供以下技術方案：

本發明提供了基于基因拷貝數變異位點對林木群體基因型的分型方法，包括以下步驟：

1)將林木物種中已經公布的基因組序列作為參考序列，操作界面的每個窗口上顯示連續的、非重疊的、大小相等區域，每個窗口能空間上左右調整；

2)將所述林木物種的不同個體測序，將所述測序得到的reads比對到所述步驟1)的參考序列上，將每個窗口內比對到所述窗口上的reads數目作為讀取深度信號；

3)根據每個窗口內reads的GC含量和偏差，空間上左右調整得到所述窗口的讀取深度信號值，根據調整后的讀取深度信號值得到不同個體在每個窗口內的讀取深度信號值形成的數據集，再根據所述數據集計算得到中值；

4)根據步驟3)得到的中值對所述步驟2)得到的讀取深度信號值進行數據標準化處理，得到校正后的讀取深度信號值；

5)根據步驟4)得到的校正后的讀取深度信號值計算不同個體在所述窗口中拷貝數，以拷貝數為2的拷貝數設為正常拷貝，正常拷貝的基因型為(1；1)；

6)當所述步驟5)中得到的拷貝數發生變異，且所述窗口內讀取深度信號與無拷貝數變異的窗口相比有顯著地增高時，屬于基因拷貝數發生重復的結構變異；對重復拷貝數CN為3或4的基因位點進行分型；

7)所述步驟6)的分型方法具體為：將標準化校正后的讀取深度信號值在1.25<讀取深度信號值<1.75范圍內的基因拷貝數變異位點設為雜合重復，雜合重復的基因型為(1；2)；將標準化校正后的讀取深度信號值在1.75<讀取深度信號值<2.25范圍內的基因拷貝數變異位點設為純合重復，純合重復的基因型為(2；2)；

8)當所述步驟5)中得到的拷貝數發生變異，且所述窗口內讀取深度信號與無拷貝數變異的窗口相比有顯著地降低時，屬于基因拷貝數發生缺失的結構變異；

9)將所述步驟8)中缺失的結構變異進行分型，所述分型的方法具體為：將標準化校正后的讀取深度信號值<0.10的基因拷貝數變異位點設為純合缺失，設定純合缺失基因型為(0；0)；將標準化校正后的讀取深度信號值為0.10<讀取深度信號值<0.75的位點設為雜合缺失，設定雜合缺失的基因型為(0；1)。

優選的，所述步驟1)中窗口的大小固定，所述窗口的大小為500bp。

優選的，所述步驟2)中讀取深度信號值由CNVnator軟件計算得到。

優選的，所述步驟3)調整所述窗口的讀取深度信號值的具體方法是：調整讀取深度信號值使GC含量達到48～52％和偏差小于5％。

優選的，所述步驟3)中值的計算方法為：將每個個體得到的讀取深度信號值的數據集按從小到大的順序排列，位于該數列中間位置的數值為該個體的讀取深度信號值的中值。

優選的，所述步驟4)中數據標準化處理按照式I計算；所述式I為x′＝x_i-Me/Std(x)，其中x′為得到的新數據，x_i為原始數據，Me為該列數據集中的中值，Std(x)為標準差。

優選的，所述步驟5)中計算不同個體在窗口中拷貝數的方法為：以2個拷貝為中心，將標準化后的讀取深度信號值經過四舍五入后，取最接近的整數作為該個體在該窗口中的拷貝數。

優選的，所述步驟5)，步驟7)和步驟9)中整個步驟是基于自然群體中的平衡選擇定律。

本發明提供的基于基因拷貝數變異位點對林木群體基因型的分型方法，是一種適用于不同群體大小的林木CNV位點的基因型分型方法，利用高通量測序結果的讀取深度信號確定CNV位點的基因型，其算法復雜度較低，操作簡便易行，可以較精確地發現和檢測CNV位點的基因型。

另外，本發明提供的分型方法，利用讀取深度信號值數據集的中值對讀深信號值進行標準化處理，可以減少測序過程中產生的測序深度、技術和人員操作等方面的誤差，降低假陽性，使整體數值更為均一化，分型結果更為精確。

說明書附圖

圖1為實施例1中利用讀深信號值的方法在4個樣本中檢測到的缺失示意圖；

圖2為實施例1中利用讀深信號值的方法在4個樣本中檢測到的重復示意圖。

具體實施方式

本發明提供了基于基因拷貝數變異位點對林木群體基因型的分型方法，包括以下步驟：

1)將林木物種中已經公布的基因組序列作為參考序列，操作界面的每個窗口上顯示連續的、非重疊的、大小相等區域，每個窗口能空間上左右調整；

2)將所述林木物種的不同個體測序，將所述測序得到的reads比對到所述步驟1)的參考序列上，將每個窗口內比對到所述窗口上的reads數目作為讀取深度信號；

3)根據每個窗口內reads的GC含量和偏差，空間上左右調整得到所述窗口的讀取深度信號值，根據調整后的讀取深度信號值得到不同個體在每個窗口內的讀取深度信號值形成的數據集，再根據所述數據集計算得到中值；

4)根據步驟3)得到的中值對所述步驟2)得到的讀取深度信號值進行數據標準化處理，得到校正后的讀取深度信號值；

5)根據步驟4)得到的校正后的讀取深度信號值計算不同個體在所述窗口中拷貝數，以拷貝數為2的拷貝數設為正常拷貝，即無拷貝數變異；所述無拷貝數變異的基因型為(1；1)；

6)當所述步驟5)中得到的拷貝數發生變異，且所述窗口內讀取深度信號與無拷貝數變異的窗口相比有顯著地增高時，屬于基因拷貝數發生重復的結構變異；對重復拷貝數CN為3或4的基因位點進行分型；

7)所述步驟6)的分型方法具體為：將標準化校正后的1.25<讀取深度信號值<1.75范圍內的基因拷貝數變異位點設為雜合重復，其基因型為(1；2)；將標準化校正后的1.75<讀取深度信號值<2.25的基因拷貝數變異位點設為純合重復，其基因型為(2；2)；

8)當所述步驟5)中得到的拷貝數發生變異，且所述窗口內讀取深度信號與無拷貝數變異的窗口相比有顯著地降低時，屬于基因拷貝數發生缺失的結構變異；

9)將所述步驟8)中缺失的結構變異進行分型，所述分型的方法具體為：將標準化校正后的讀取深度信號值<0.10的基因拷貝數變異位點設為純合缺失，設定純合缺失的基因型為(0；0)；將標準化校正后的讀取深度信號值為0.10<讀取深度信號值<0.75的位點設為雜合缺失，設定雜合缺失的基因型為(0；1)。

本發明將林木物種中已經公布的基因組序列作為參考序列，操作界面的每個窗口上顯示連續的、非重疊的、大小相等區域，每個窗口能空間上左右調整。

本發明中，林木物種中已經公布的基因組序列優選為從NCBI網站中下載得到。

本發明中，提供所述操作界面的軟件優選為Linux系統的Shell窗口。所述Linux系統的Shell窗口是基于大型計算機服務器。

本發明中，所述操作界面的每個窗口上顯示的連續的、非重疊的、大小相等區域優選利用CNVnator算法劃分。所述CNVnator算法優選利用均值漂移技術(mean-shift technique)將讀取深度信號分成有潛在CNV的小片段中。所述窗口的大小固定，所述窗口的大小優選為500bp。

得到參考序列后，本發明將所述林木物種的不同個體測序，將所述測序得到的reads比對到所述參考序列上，將每個窗口內比對到所述窗口上的reads數目作為讀取深度信號。所述林木物種優選為楊樹，更優選為毛白楊。

本發明中，利用了林木基因雜合度高、DNA序列多態性豐富的特點，結合Illumina的pair-end和454的mate-pair兩種測序方法對物種的不同個體進行測序，這種測序方案很好地結合了短序列插入片段和長序列插入片段的各自突出優勢，更好地提高了測序結果的準確性，為后續的精確分型提供了保障。

本發明中，所述林木物種的不同個體測序具體是測定林木物種中不同個體的基因組。所述林木物種個體的個數優選為400～500個，更優選為435個。

本發明中，所述比對優選包括以下步驟：

Ⅰ.利用軟件BWA-0.7.8中的aln算法將不同個體得到的大量reads分別比對到參考基因組序列上；

Ⅱ.以picard軟件包中的Markduplicate工具標記可能的PCR重復；

Ⅲ.利用軟件Samtools文件包中的flagstat工具統計得到所有樣品的比對信息文件，并以bam格式進行保存。

比對結束后，本發明將同一個個體中比對到參考序列上gap位置的reads刪除。

本發明中，所述窗口的大小優選為固定值，所述窗口的大小優選為500bp。

本發明中，所述讀取深度信號值優選由CNVnator軟件劃分的連續的、非重疊的、大小相等區域中映射的reads數目計算得到。利用讀取深度信號值分型的思路為：假定讀取深度是一個泊松分布，然后利用讀深信號的隨機分布檢測目標樣本中的重復和缺失類型的CNV，并進行分型。

得到每個窗口的讀取深度信號值后，本發明根據每個窗口內的GC含量和偏差，空間上左右調整得到所述窗口的讀取深度信號值，根據調整后的讀取深度信號值得到不同個體在每個窗口內的讀取深度信號值形成的數據集，再根據所述數據集計算得到中值。

本發明中，所述調整所述窗口的讀取深度信號值的具體方法優選是：調整讀取深度值使GC含量達到48～52％或偏差小于5％。

本發明中，所述步驟3)中值的計算方法優選為：將每個個體得到的讀取深度信號值的數據集按從小到大的順序排列，位于該數列中間位置的數值為該個體的讀取深度信號值的中值，也稱中位數。

得到中值后，本發明根據中值對所述讀取深度信號值進行數據標準化處理，得到校正后的讀取深度信號值。

本發明中，所述數據標準化處理按照式I計算；

x′＝x_i-Me/Std(x) 式I

其中為x′得到的新數據；

x_i為原始數據；

Me為該列數據集中的中值；

Std(x)為標準差。

得到校正后的讀取深度信號值x’后，本發明根據校正后的讀取深度信號值計算不同個體在所述窗口中的拷貝數，以CN為2的拷貝數設為正常拷貝，即無拷貝數變異，所述無拷貝數變異的基因型為(1；1)。

本發明中，所述計算不同個體在窗口中拷貝數的方法優選為：以2個拷貝(即二倍體)為中心，將標準化后的讀取深度信號值經過四舍五入后，取最接近的整數作為該個體在該窗口中的拷貝數。

本發明中，當個體某個基因存在結構變異時，個體樣本所測得到的reads映射到窗口中的讀取深度信號較無拷貝數變異的區域會有顯著地增高或降低，說明在該區域發生了基因拷貝數的重復(duplication)或缺失(deletion)。

當所述窗口內讀取深度信號與無拷貝數變異的窗口相比有顯著地增高時，屬于基因拷貝數發生重復的結構變異；本發明僅對重復拷貝數CN為3或4的基因位點進行分型。所述分型方法具體為：將標準化校正后的1.25<讀取深度信號值<1.75范圍內的基因拷貝數變異位點設為雜合重復，雜合重復的基因型為(1；2)；將標準化校正后的1.75<RD-value<2.25的CNV位點設為純合重復，純合重復的基因型為(2；2)。本發明中，所述顯著地增高是指發生結構變異的窗口中映射匹配的reads數目顯著的比參考基因組上對應的該窗口中的reads數目較多。

當所述窗口內讀取深度信號與無拷貝數變異的窗口相比有顯著地降低時，屬于基因拷貝數發生缺失的結構變異；將所述缺失的結構變異進行分型，所述分型的方法具體為：將標準化校正后的讀取深度信號值<0.10的基因拷貝數變異位點設為純合缺失，設定純合缺失的基因型為(0；0)；將標準化校正后的讀取深度信號值為0.10<讀取深度信號值<0.75的位點設為雜合缺失，設定雜合缺失的基因型為(0；1)。本發明中，所述顯著地降低是指發生結構變異的窗口中映射匹配的reads數目顯著的比參考基因組上對應的該窗口中的reads數目較少。

本發明中，所述顯著地增高或降低具體是將讀取深度信號值采用數學上的顯著性統計算法進行計算。所述算法為單樣本t-檢驗(one-sample t-test)(p<0.05)。

本發明中，所述檢測標準及分型方法均基于自然群體中的平衡選擇定律。

下面結合實施例對本發明提供的一種基于基因拷貝數變異位點對林木群體基因型分型方法進行詳細的說明，但是不能把它們理解為對本發明保護范圍的限定。

實施例1

(1)原材料的獲得：從毛白楊自然分布區收集了435株個體作為研究對象。利用CTAB法提取每株個體的基因組DNA后，送至上海伯豪生物技術有限公司進行測序。測序選用Illumina的pair-end和454的mate-pair兩種測序方法對毛白楊個體進行測序。

(2)比對：利用上述比對的工具、軟件和算法，將每株個體得到的測序片段即reads與參考基因組序列進行比對，去除PCR重復、冗余和測序過程中引入的接頭序列。

(3)統計：利用CNVnator軟件和其算法統計每個具有潛在CNV位點區域的相關信息，如CNV的起始-結束坐標、CNV的長度、類型(deletion或duplication)、讀取深度信號值(RD-value)、可進行數學上顯著性統計的P值、確定候選CNV假陽性的q₀值等。為了提高分型結果的精確性，避免假陽性，本發明中選擇p<0.01和q₀<0.5的CNV位點進行后續分析。

結合上一步驟得到的讀取深度信號值，根據當前窗口內的GC含量和偏差，基于不同個體在該窗口內的讀取深度信號值的數據集的中值，對該窗口的讀深信號值進行數據標準化處理。

(4)校正：利用不同個體在某個窗口內的讀深信號值數據集的中值對讀深信號進行標準化校正，減少誤差，降低假陽性。

(5)合并：由于所采用的材料毛白楊本身具有基因雜合度較高，DNA序列多態性豐富的特點，再加上測序過程中不可避免的因為測序技術、實驗人員的操作等造成一定的試驗誤差。因此，即使同一個區域中的CNV位點也會因上述存在的問題而出現每株個體與每株個體之間的CNV起始和結束坐標并非完全一致。為了解決這一問題，本發明采取小于等于5個bin窗口(每個bin＝500bp，5個bin的長度＝500*5＝2,500bp)的算法進行合并，合并后再利用上述步驟5)、6)和7)中的檢測標準和分型方法進行CNV位點基因型的檢測。該步驟具體如下(以本發明中的實際例子說明)：

如，第1號染色體上發現了一個deletion類型的CNV位點，發生該CNV位點的毛白楊株數為144株，經過上述所有步驟檢測到該CNV位點發生在染色體上的位置為(1260001-1265000，即1260001為起始坐標，1265000為結束坐標)，其中有30株個體以1260001為起始坐標，3株個體以1260501為起始坐標，111株個體以1261001為起始坐標。基于上述小于等于5個bin窗口的算法對該deletion類型的CNV位點的染色體位置進行合并如下：chr01：1261001-1265000。再如：同樣在第1號染色體上發現了一個duplication類型的CNV，發生該CNV位點的毛白楊株數為114株，經過上述所有步驟和算法檢測到該CNV發生在染色體上的位置為(1292001-1327500)，其中12株個體的結束坐標為1325000，6株個體的結束坐標為1326000，15株個體的結束坐標為1326500，66株個體的結束坐標為1327000，15株個體的結束坐標為1327500。基于上述小于等于5個bin窗口的算法對該duplication類型的CNV在染色體上的位置進行合并如下，chr01：1292001-1325000。即，如果某CNV位點存在不一致的起始坐標或結束坐標時，其起始坐標總是以發生該CNV位點最大的起始坐標作為該CNV位點的起始坐標，其結束坐標總是以發生該CNV位點的最小結束坐標作為該CNV位點的結束坐標。基于此原則，本發明對毛白楊基因組中的19條染色體上發生的CNV位點進行合并，然后進行分型。

(6)分型：利用標準化后的讀深信號值對完成合并的每個CNV位點進行基因型分型。其具體操作如下：

i.對于deletion類型的CNV，將標準化校正后的RD-value<0.10的CNV位點設為純合缺失，其基因型為(0；0)；將標準化校正后的0.10<RD-value<0.75的位點設為雜合缺失，其基因型為(0；1)；

ii.對于沒有發生結構變異的基因位點，將其設為正常拷貝(CN＝2)，其基因型為(1；1)；

iii.對于duplication類型的CNV，將標準化校正后的1.25<RD-value<1.75的CNV位點設為雜合重復，其基因型為(1；2)；將標準化校正后的1.75<RD-value<2.25的CNV位點設為純合重復，其基因型為(2；2)。

統計結果顯示：本發明在毛白楊基因組的前4條染色體上共成功分型到1,628個CNV位點，其中第1號染色體上成功分型了743個CNV位點，包括523個deletion和220個duplication；第2號染色體成功分型了321個CNV位點，包括199個deletion和122個duplication；第3號成功分型了284個CNV位點，包括175個deletion和109個duplication；第4號染色體上分別成功分型了280個CNV位點，包括195個deletion和85個duplication。

由以上實施例可知，本發明提供的基于基因拷貝數變異位點對林木群體基因型的分型方法，利用讀深信號進行林木自然群體中的CNV位點的基因型檢測，對CNV與數量性狀的關聯分析具有重要意義。

以上所述僅是本發明的優選實施方式，應當指出，對于本技術領域的普通技術人員來說，在不脫離本發明原理的前提下，還可以做出若干改進和潤飾，這些改進和潤飾也應視為本發明的保護范圍。