本發(fā)明涉及分子生物學技術(shù)領(lǐng)域，尤其涉及一種基因組胞嘧啶位點表觀基因型分型方法。

背景技術(shù)：

在高等真核生物中，DNA甲基化僅發(fā)生在Cp G二核苷酸G5′側(cè)的C上。當位于基因啟動子區(qū)內(nèi)富含Cp G序列(Cp G島)時，這種修飾作用則對基因的表達有著重要的調(diào)控作用。此外，它還同基因組印記、女性X染色體的基因滅活、細胞增殖、分化發(fā)育、腫瘤的發(fā)生和發(fā)展以及遺傳的不穩(wěn)定性等密切有關(guān)。

近年來甲基化測序技術(shù)逐漸發(fā)展完善，包括重硫酸鹽處理基因組甲基化測序，利用重硫酸鹽處理基因組DNA，使未發(fā)生甲基化的胞嘧啶脫氨基變成尿嘧啶，而發(fā)生了甲基化的胞嘧啶則不會發(fā)生變化。比對經(jīng)重硫酸鹽處理和未經(jīng)處理的樣本，可以檢測到甲基化位點。進一步結(jié)合高通量測序技術(shù)，能夠從全基因組水平和單堿基精度來分析5＇甲基胞嘧啶，由此能夠發(fā)現(xiàn)很多傳統(tǒng)的基因組學研究所不能檢測到的甲基化位點。

目前，甲基化測序數(shù)據(jù)分析涉及到基因組甲基化水平，甲基化分布類型以及分布傾向等基本特征的分析，但是目前還無法實現(xiàn)表觀基因型分型。

技術(shù)實現(xiàn)要素：

有鑒于此，本發(fā)明的目的在于提供一種基因組胞嘧啶位點表觀基因型分型方法。

為了實現(xiàn)上述發(fā)明目的，本發(fā)明提供以下技術(shù)方案：

本發(fā)明提供了一種基因組胞嘧啶位點表觀基因型分型方法，包括以下步驟：1)對待測樣品父母本和子代樣本進行重亞硫酸鹽全基因組甲基化測序，獲得父母本和子代樣本基因組序列；2)將父母本和子代樣本基因組序列與參考基因組比對獲得比對結(jié)果，確定待測胞嘧啶位點；3)利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，將比對結(jié)果進行染色體坐標排序、reads去重復處理，再通過GATK2-V3.2對已知胞嘧啶位點上下游5～10bp的序列進行Call SNPs，從而區(qū)分子代等位基因序列；4)將已獲得的被區(qū)分過子代等位基因序列與父母本基因組序列比對，完成胞嘧啶表觀基因型分型。

優(yōu)選的，步驟1)中所述的待測樣品為有參考基因組的物種。

優(yōu)選的，所述參考基因組為待測樣品本物種已測序基因組或待測樣品近緣物種已測序基因組。

優(yōu)選的，所述步驟1)具體包括以下步驟：1.1)用CTAB法提取待測樣品父母本和子代樣本基因組DNA；1.2)對提取得到的基因組DNA樣品進行質(zhì)量、純度和濃度檢測篩選獲得合格的父母本和子代樣本基因組DNA樣品；1.3)重亞硫酸鹽法構(gòu)建上述合格的父母本和子代樣本基因組DNA樣品測序文庫；1.4)質(zhì)檢篩選合格基因組DNA樣品測序文庫，保證所述文庫插入片段為320-520bp，文庫有效濃度>2nM；1.5)對合格的DNA樣品文庫進行雙末端Hiseq測序，獲得父母本和子代樣本基因組序列。

優(yōu)選的，步驟1.3)構(gòu)建父母本和子代樣本基因組DNA樣品測序文庫時加入比例為建庫DNA起始量的1/1000的陰性對照lambdaDNA。

優(yōu)選的，步驟1.3)構(gòu)建父母本和子代樣本基因組DNA樣品測序文庫包括以下步驟：隨機打斷基因組DNA至200～300bp獲得DNA片段；對DNA片段進行平末端修復后加尾巴A堿基獲得帶尾巴A的DNA片段；在帶尾巴A的DNA片段上連接測序接頭后進行Bisulfite處理；最后進行PCR擴增獲得基因組DNA樣品測序文庫。

優(yōu)選的，步驟3)中所述染色體坐標排序采用picard-tools工具中的SortSam完成按照染色體坐標順序從小到大排序。

優(yōu)選的，步驟3)中所述reads去重復處理采用picard-tools完成。

優(yōu)選的，步驟3)中所述Call SNPs采用UnifiedGenotyper工具完成。

優(yōu)選的，步驟3)中所述Call SNPs完成后還包括對所述Call SNPs結(jié)果進行過濾。

本發(fā)明的有益效果：本發(fā)明提供的基因組胞嘧啶位點表觀基因型分型方法將樣本基因組經(jīng)重亞硫酸鹽處理，全基因組甲基化測序，篩選測定后的胞嘧啶位點，經(jīng)序列比對，SNP Calling后，最終將子代與父母本序列比對完成胞嘧啶表觀基因型分型。本發(fā)明首次完成了基因組胞嘧啶位點表觀基因型分型，技術(shù)成熟，成本低，易于操作與推廣應(yīng)用。

附圖說明

圖1為無降解，無RNA污染檢測合格的DNA樣品電泳條帶圖；

圖2為嚴重降解檢測不合格的DNA樣品電泳條帶圖；

圖3為父母本胞嘧啶甲基化表觀基因型分型示意圖；

圖4為子代胞嘧啶分型純合甲基化表觀基因型示意圖；

圖5為子代胞嘧啶分型雜合甲基化表觀基因型示意圖；

圖6為子代胞嘧啶分型純合非甲基化表觀基因型示意圖。

具體實施方式

本發(fā)明提供了一種基因組胞嘧啶位點表觀基因型分型方法，包括以下步驟：1)對待測樣品父母本和子代樣本進行重亞硫酸鹽全基因組甲基化測序，獲得父母本和子代樣本基因組序列；2)將父母本和子代樣本基因組序列與參考基因組比對獲得比對結(jié)果，確定待測胞嘧啶位點；3)利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，將比對結(jié)果進行染色體坐標排序、reads去重復處理，再通過GATK2-V3.2對已知胞嘧啶位點上下游5～10bp的序列進行Call SNPs，從而區(qū)分子代等位基因序列；4)將已獲得的被區(qū)分過子代等位基因序列與父母本基因組序列比對，完成胞嘧啶表觀基因型分型。

本發(fā)明中所述待測樣品優(yōu)選的為有參考基因組的物種，所述參考基因組為待測樣品本物種已測序基因組或待測樣品近緣物種已測序基因組。所述待測樣品優(yōu)選的為人類，斑馬魚，毛果樣，胡楊，小葉楊或擬南芥等有參考基因組的物種。

本發(fā)明在確定待測樣品后，對待測樣品父母本和子代樣本進行重亞硫酸鹽全基因組甲基化測序，獲得父母本和子代樣本基因組序列。本發(fā)明中優(yōu)選的先進行待測樣品父母本和子代樣本基因組DNA的提取；所述基因組DNA的提取優(yōu)選的采用CTAB法，本發(fā)明中CTAB法提取基因組DNA的具體步驟參見本領(lǐng)域常規(guī)的CTAB法提取基因組DNA的步驟，優(yōu)選的采用文獻(Chang et al.,1993)中所記載的方法。

本發(fā)明在得到待測樣品父母本和子代樣本基因組DNA后，對所述的基因組DNA進行質(zhì)量、純度和濃度檢測，篩選獲得合格的父母本和子代樣本基因組DNA樣品。在本發(fā)明中所述基因組DNA質(zhì)量檢測優(yōu)選的使用瓊脂糖凝膠電泳進行，根據(jù)電泳條帶分析基因組DNA降解程度以及是否有RNA污染。具體的結(jié)果如圖1和圖2，圖1中樣品為檢測合格的無降解，無RNA污染樣品電泳條帶，圖2中樣品4號為嚴重降解；5號有重度RNA污染。本發(fā)明篩選質(zhì)量合格的如圖1所示的基因組DNA進行下一步操作，若基因組DNA樣品質(zhì)量不合格需要重新提取基因組DNA。本發(fā)明在檢測基因組DNA質(zhì)量合格后，優(yōu)選的進一步檢測其純度是否合格。在本發(fā)明中，優(yōu)選的采用Nanodrop檢測基因組DNA的純度，當OD260/280為1.6～1.8之間，說明基因組DNA的純度合格，若所述基因組DNA的OD260/280不在上述范圍內(nèi)，則所述基因組DNA純度不合格，需要重新提取。本發(fā)明在篩選獲得質(zhì)量和純度合格的基因組DNA后，優(yōu)選的采用Qubit對基因組DNA濃度進行精確定量，為后續(xù)操作提供依據(jù)。

本發(fā)明在獲得用合格的父母本和子代樣本基因組DNA樣品后，采用重亞硫酸鹽法構(gòu)建上述合格的父母本和子代樣本基因組DNA樣品測序文庫；在本發(fā)明中，構(gòu)建父母本和子代樣本基因組DNA樣品測序文庫時加入比例為建庫DNA起始量的1/1000的陰性對照lambda DNA。所述基因組DNA樣品測序文庫的構(gòu)建具體的包括以下步驟：首先隨機打斷基因組DNA至200～300bp獲得DNA片段，所述隨機打斷基因組DNA優(yōu)選的使用Covaris S220進行；然后對DNA片段進行平末端修復，所述平末端修復優(yōu)選的使用內(nèi)切酶和外切酶將黏性末端修復成為平末端；隨后在平末端DNA片段后加尾巴A堿基獲得帶尾巴A的DNA片段，所述加尾巴A堿基步驟采用本領(lǐng)域常規(guī)的手段即可，無其他特殊要求；本發(fā)明在帶尾巴A的DNA片段上連接測序接頭后進行Bisulfite處理；所述的Bisulfite處理優(yōu)選的采用Zymo Research的EZ DNA Methylation GoldKit進行，本發(fā)明中所述Bisulfite處理的作用是將未發(fā)生甲基化的C變成U，在后續(xù)的PCR擴增過程后變?yōu)門，而甲基化的C則保持不變。本發(fā)明在Bisulfite處理后進行PCR擴增獲得基因組DNA樣品測序文庫。

本發(fā)明在獲得基因組DNA樣品測序文庫后，對文庫進行質(zhì)檢篩選獲得合格基因組DNA樣品測序文庫。本發(fā)明中在質(zhì)檢文庫之前優(yōu)選的使用Qubit2.0對文庫DNA濃度進行初步定量，然后根據(jù)文庫DNA的濃度將其進行稀釋，文庫稀釋后的濃度優(yōu)選的為1ng/ul；在文庫稀釋后，優(yōu)選的使用Agilent2100對文庫的插入片段長度進行檢測，篩選后的文庫插入片段在320～520bp為合格；使用Q-PCR方法對文庫的有效濃度進行準確定量，所述文庫有效濃度>2nM為合格。

本發(fā)明在獲得合格基因組DNA樣品測序文庫后，對合格的DNA樣品文庫進行雙末端Hiseq測序，獲得父母本和子代樣本基因組序列。所述的雙末端Hiseq測序優(yōu)選的委托測序公司完成，在本發(fā)明具體實施例中具體的委托北京諾禾致源生物信息科技有限公司通過全基因組甲基化測序完成。

本發(fā)明在獲得父母本和子代樣本基因組序列后，將所獲子代基因組序列與參考基因組比對，本發(fā)明中所述參考基因組優(yōu)選的為所測樣品本物種基因組；如本物種未測序，優(yōu)選的選取已測序近緣物種基因組做參考基因組，確定待測胞嘧啶位點。

本發(fā)明在獲得待測胞嘧啶位點后，利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，將比對結(jié)果進行染色體坐標排序、reads去重復處理，再通過GATK2-V3.2對已知胞嘧啶位點上下游5～10bp的序列進行Call SNPs，從而區(qū)分子代等位基因序列，由于序列經(jīng)由重亞硫酸鹽甲基化測序，所以胞嘧啶上的甲基化修飾狀況已被記錄。

本發(fā)明中上述操作優(yōu)選的委托北京諾禾致源生物信息科技有限公司完成。具體的操作如下：1)對原始下機fastq文件進行過濾和比對(mapping)，對Illumina下機數(shù)據(jù)使用bwa進行mapping。

所述Bwa比對步驟如下：Ⅰ.對參考基因組構(gòu)建索引：Ⅱ.尋找輸入reads文件的SA坐標。Ⅲ.生成sam格式的比對文件。如果一條read比對到多個位置，會隨機選擇一種。2)對sam文件進行進行重新排序(reorder)，由BWA生成的sam文件時按字典式排序法進行的排序(lexicographically)進行排序的(chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chrX，chrY)，但是GATK在進行callsnp的時候是按照染色體組型(karyotypic)進行的(chrM，chr1，chr2…chr22，chrX，chrY)，因此要對原始sam文件進行reorder?？梢允褂胮icard-tools中的ReorderSam完成。

3)將sam文件轉(zhuǎn)換成bam文件(bam是二進制文件，運算速度快)，這一步可使用samtools view完成。4)對bam文件進行sort排序處理，這一步是將sam文件中同一染色體對應(yīng)的條目按照坐標順序從小到大進行排序?？梢允褂胮icard-tools中SortSam完成。5)對bam文件進行加頭(head)處理GATK2.0以上版本不支持無頭文件的變異檢測。加頭這一步可以在BWA比對的時候進行，通過-r參數(shù)的選擇可以完成。如果在BWA比對期間沒有選擇-r參數(shù)，可以增加這一步驟?？墒褂胮icard-tools中AddOrReplaceReadGroups完成。6)Merge：如果一個樣本分為多個lane進行測序，那么在進行下一步之前可以將每個lane的bam文件合并。7)Duplicates Marking：在制備文庫的過程中，由于PCR擴增過程中會存在一些偏差，也就是說有的序列會被過量擴增。這樣，在比對的時候，這些過量擴增出來的完全相同的序列就會比對到基因組的相同位置。而這些過量擴增的reads并不是基因組自身固有序列，不能作為變異檢測的證據(jù)，因此，要盡量去除這些由PCR擴增所形成的duplicates，這一步可以使用picard-tools來完成。去重復的過程是給這些序列設(shè)置一個flag以標志它們，方便GATK的識別。還可以設(shè)置REMOVE_DUPLICATES＝true來丟棄duplicated序列。對于是否選擇標記或者刪除，對結(jié)果應(yīng)該沒有什么影響，GATK官方流程里面給出的例子是僅做標記不刪除。這里定義的重復序列是這樣的：如果兩條reads具有相同的長度而且比對到了基因組的同一位置，那么就認為這樣的reads是由PCR擴增而來，就會被GATK標記。8)對上一步得到的結(jié)果生成索引文件：可以用samtools完成，生成的索引后綴是bai。9)Local realignment around indels：這一步的目的就是將比對到indel附近的reads進行局部重新比對，將比對的錯誤率降到最低。一般來說，絕大部分需要進行重新比對的基因組區(qū)域，都是因為插入/缺失的存在，因為在indel附近的比對會出現(xiàn)大量的堿基錯配，這些堿基的錯配很容易被誤認為SNP。還有，在比對過程中，比對算法對于每一條read的處理都是獨立的，不可能同時把多條reads與參考基因組比對來排錯。因此，即使有一些reads能夠正確的比對到indel，但那些恰恰比對到indel開始或者結(jié)束位置的read也會有很高的比對錯誤率，這都是需要重新比對的。Local realignment就是將由indel導致錯配的區(qū)域進行重新比對，將indel附近的比對錯誤率降到最低。主要分為兩步：第一步，通過運行RealignerTargetCreator來確定要進行重新比對的區(qū)域。第二步，通過運行IndelRealigner在這些區(qū)域內(nèi)進行重新比對。10)Base quality score recalibration：這一步是對bam文件里reads的堿基質(zhì)量值進行重新校正，使最后輸出的bam文件中reads中堿基的質(zhì)量值能夠更加接近真實的與參考基因組之間錯配的概率。這一步適用于多種數(shù)據(jù)類型，包括illunima、solid、454、CG等數(shù)據(jù)格式。在GATK2.0以上版本中還可以對indel的質(zhì)量值進行校正，這一步對indel calling非常有幫助BQSR主要有三步：第一步：利用工具BaseRecalibrator，根據(jù)一些known sites，生成一個校正質(zhì)量值所需要的數(shù)據(jù)文件，GATK網(wǎng)站以“.grp”為后綴命名。第二步：利用第一步生成的ChrALL.100.sam.recal.08-1.grp來生成校正后的數(shù)據(jù)文件，也是以“.grp”命名，這一步主要是為了與校正之前的數(shù)據(jù)進行比較，最后生成堿基質(zhì)量值校正前后的比較圖，如果不想生成最后BQSR比較圖，這一步可以省略。第三步：利用工具PrintReads將經(jīng)過質(zhì)量值校正的數(shù)據(jù)輸出到新的bam文件中，用于后續(xù)的變異檢測。11)分析和評估BQSR結(jié)果：這一步會生成評估前后堿基質(zhì)量值的比較結(jié)果，可以選擇使用圖片和表格的形式展示。13)Reduce bam file：這一步是使用ReduceReads這個工具將bam文件進行壓縮，生成新的bam文件，新的bam文件仍然保持bam文件的格式和所有進行變異檢測所需要的信息。這樣不僅能夠節(jié)省存儲空間，也方便后續(xù)變異檢測過程中對數(shù)據(jù)的處理。14)Variant Calling：GATK在這一步里面提供了兩個工具進行變異檢測——UnifiedGenotyper和HaplotypeCaller。我們使用UnifiedGenotyper。UnifiedGenotyper是集合多種變異檢測方法而成的一種Variants Caller，既可以用于單個樣本的變異檢測，也可以用于群體的變異檢測。UnifiedGenotyper使用貝葉斯最大似然模型，同時估計基因型和基因頻率，最后對每一個樣本的每一個變異位點和基因型都會給出一個精確的后驗概率。15)對原始變異檢測結(jié)果進行過濾(hard filter andVQSR)：這一步的目的就是對上一步call出來的變異位點進行過濾，去掉不可信的位點。這一步可以有兩種方法，一種是通過GATK的VariantFiltration，另一種是通過GATK的VQSR(變異位點質(zhì)量值重新校正)進行過濾。通過GATK網(wǎng)站上提供的最佳方案可以看出，GATK是推薦使用VASR的，但使用VQSR數(shù)據(jù)量一定要達到要求，數(shù)據(jù)量太小無法使用高斯模型。還有，在使用VAQR時，indel和snp要分別進行。Ⅰ.VariantRecalibrator，VariantRecalibrator：通過大量的高質(zhì)量的已知變異集合的各個注釋(包括很多種，后面介紹)的值來創(chuàng)建一個高斯混合模型，然后用于評估所有的變異位點。這個文件最后將生成一個recalibration文件。Ⅱ.ApplyRecalibration，ApplyRecalibration：這一步將模型的各個參數(shù)應(yīng)用于原始vcf文件中的每一個變異位點，這時，每一個變異位點的注釋信息列中都會出現(xiàn)一個VQSLOD值，然后模型會根據(jù)這個值對變異位點進行過濾，過濾后的信息會寫在vcf文件的filter一列中。

本發(fā)明在獲得區(qū)分子代等位基因序列后，將已獲得的被區(qū)分過子代等位基因序列與父母本基因組序列比對，完成胞嘧啶表觀基因型分型。結(jié)果如圖3～6所示：圖3是父母本樣品胞嘧啶甲基化表觀基因型分型圖示，X代表胞嘧啶分型片段側(cè)翼序列中的任意堿基(可以為ATCG中任意一個)，m代表修飾在胞嘧啶上的甲基基團，子代胞嘧啶分型片段經(jīng)過與親本比對最終確定的甲基化表觀基因型可分為純合甲基化表觀基因型，雜合甲基化表觀基因型和純合非甲基化型三種。

下面結(jié)合實施例對本發(fā)明提供的基因組胞嘧啶位點表觀基因型分型方法進行詳細的說明，但是不能把它們理解為對本發(fā)明保護范圍的限定。

實施例1

以人類為研究對象，將包含父母本及其子代在內(nèi)的樣本進行重硫酸鹽全基因組甲基化測序，獲得基因組序列；

1)樣本基因組DNA提取，使用CTAB法；2)DNA樣品檢測，使用瓊脂糖凝膠電泳分析DNA降解程度以及是否有RNA污染，檢測參考圖1和圖2，圖1中樣品為檢測合格的樣品電泳條帶如下：無降解，無RNA污染，圖2中樣品4號為嚴重降解；5號有重度RNA污染。再通過Nanodrop檢測DNA純度(OD260/280)，最后使用Qubit對DNA濃度進行精確定量；3)文庫構(gòu)建，樣品檢測為A或B時認為合格，合格后加入比例為建庫起始量的1/1000的陰性對照(lambda DNA)，首先使用Covaris S220將基因組DNA隨機打斷至200-300bp。對打斷后的DNA片段進行末端修復、加尾巴A，并連接上所有胞嘧啶均經(jīng)過甲基化修飾的測序接頭。隨后進行Bisulfite處理(采用EZ DNAMethylation GoldKit，Zymo Research)，經(jīng)過處理，未發(fā)生甲基化的C變成U(PCR擴增后變?yōu)門)，而甲基化的C保持不變，最后進行PCR擴增，得到最終的DNA文庫；4)文庫質(zhì)檢，文庫構(gòu)建完成后，先使用Qubit2.0進行初步定量，稀釋文庫至1ng/ul，隨后使用Agilent2100對文庫的插入片段長度進行檢測，符合預(yù)期后，使用Q-PCR方法對文庫的有效濃度進行準確定量(文庫有效濃度>2nM)，以保證文庫質(zhì)量；5)上機測序，庫檢合格后，把不同文庫按照有效濃度及目標下機數(shù)據(jù)量的需求，pooling后進行Hiseq測序，測序策略為雙末端測序。

將所獲基因組序列包含父母本與子代基因組與人類基因組比對；

利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，將比對結(jié)果進行染色體坐標排序、reads去重復處理，再通過GATK2(V3.2)對已知胞嘧啶位點附近10bp的序列進行Call SNPs，從而區(qū)分等位基因序列(此時，由于序列經(jīng)由重硫酸鹽甲基化測序，所以胞嘧啶上的甲基化修飾狀況已被記錄)；

將已獲得的被區(qū)分過的子代等位基因序列與父母本基因組序列比對，進而完成胞嘧啶表觀基因型分型如圖3～6所示。圖3是父母本樣品胞嘧啶甲基化表觀基因型分型圖示，X代表胞嘧啶分型片段側(cè)翼序列中的任意堿基(可以為ATCG中任意一個)，m代表修飾在胞嘧啶上的甲基基團，子代胞嘧啶分型片段經(jīng)過與親本比對最終確定的甲基化表觀基因型可分為純合甲基化表觀基因型，雜合甲基化表觀基因型和純合非甲基化型三種。

實施例2

以胡楊為研究對象，將包含父母本及其子代在內(nèi)的樣本進行重硫酸鹽全基因組甲基化測序，獲得基因組序列；

1)樣本基因組DNA提取，使用CTAB法；2)DNA樣品檢測，使用瓊脂糖凝膠電泳分析DNA降解程度以及是否有RNA污染，檢測參考圖1和圖2，圖1中樣品為檢測合格的樣品電泳條帶如下：無降解，無RNA污染，圖2中樣品4號為嚴重降解；5號有重度RNA污染。再通過Nanodrop檢測DNA純度(OD260/280)，最后使用Qubit對DNA濃度進行精確定量；3)文庫構(gòu)建，樣品檢測為A或B時認為合格，合格后加入比例為建庫起始量的1/1000的陰性對照(lambda DNA)，首先使用Covaris S220將基因組DNA隨機打斷至200-300bp。對打斷后的DNA片段進行末端修復、加尾巴A，并連接上所有胞嘧啶均經(jīng)過甲基化修飾的測序接頭。隨后進行Bisulfite處理(采用EZ DNAMethylation GoldKit，Zymo Research)，經(jīng)過處理，未發(fā)生甲基化的C變成U(PCR擴增后變?yōu)門)，而甲基化的C保持不變，最后進行PCR擴增，得到最終的DNA文庫；4)文庫質(zhì)檢，文庫構(gòu)建完成后，先使用Qubit2.0進行初步定量，稀釋文庫至1ng/ul，隨后使用Agilent2100對文庫的插入片段長度進行檢測，符合預(yù)期后，使用Q-PCR方法對文庫的有效濃度進行準確定量(文庫有效濃度>2nM)，以保證文庫質(zhì)量；5)上機測序，庫檢合格后，把不同文庫按照有效濃度及目標下機數(shù)據(jù)量的需求，pooling后進行Hiseq測序，測序策略為雙末端測序。

將所獲基因組序列包含父母本與子代基因組與胡楊基因組比對；

利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，將比對結(jié)果進行染色體坐標排序、reads去重復處理，再通過GATK2(V3.2)對已知胞嘧啶位點附近10bp的序列進行Call SNPs，從而區(qū)分等位基因序列(此時，由于序列經(jīng)由重硫酸鹽甲基化測序，所以胞嘧啶上的甲基化修飾狀況已被記錄)；

將已獲得的被區(qū)分過的子代等位基因序列與父母本基因組序列比對，進而完成胞嘧啶表觀基因型分型如圖3～6所示。圖3是父母本樣品胞嘧啶甲基化表觀基因型分型圖示，X代表胞嘧啶分型片段側(cè)翼序列中的任意堿基(可以為ATCG中任意一個)，m代表修飾在胞嘧啶上的甲基基團，子代胞嘧啶分型片段經(jīng)過與親本比對最終確定的甲基化表觀基因型可分為純合甲基化表觀基因型，雜合甲基化表觀基因型和純合非甲基化型三種。

實施例3

以小葉楊為研究對象，將包含父母本及其子代在內(nèi)的樣本進行重硫酸鹽全基因組甲基化測序，獲得基因組序列；

1)樣本基因組DNA提取，使用CTAB法；2)DNA樣品檢測，使用瓊脂糖凝膠電泳分析DNA降解程度以及是否有RNA污染，檢測參考圖1和圖2，圖1中樣品為檢測合格的樣品電泳條帶如下：無降解，無RNA污染，圖2中樣品4號為嚴重降解；5號有重度RNA污染。再通過Nanodrop檢測DNA純度(OD260/280)，最后使用Qubit對DNA濃度進行精確定量；3)文庫構(gòu)建，樣品檢測為A或B時認為合格，合格后加入比例為建庫起始量的1/1000的陰性對照(lambda DNA)，首先使用Covaris S220將基因組DNA隨機打斷至200-300bp。對打斷后的DNA片段進行末端修復、加尾巴A，并連接上所有胞嘧啶均經(jīng)過甲基化修飾的測序接頭。隨后進行Bisulfite處理(采用EZ DNAMethylation GoldKit，Zymo Research)，經(jīng)過處理，未發(fā)生甲基化的C變成U(PCR擴增后變?yōu)門)，而甲基化的C保持不變，最后進行PCR擴增，得到最終的DNA文庫；4)文庫質(zhì)檢，文庫構(gòu)建完成后，先使用Qubit2.0進行初步定量，稀釋文庫至1ng/ul，隨后使用Agilent2100對文庫的插入片段長度進行檢測，符合預(yù)期后，使用Q-PCR方法對文庫的有效濃度進行準確定量(文庫有效濃度>2nM)，以保證文庫質(zhì)量；5)上機測序，庫檢合格后，把不同文庫按照有效濃度及目標下機數(shù)據(jù)量的需求，pooling后進行Hiseq測序，測序策略為雙末端測序。

因小葉楊基因組尚未被測序，所以將所獲基因組序列包含父母本與子代基因組與小葉楊近緣物種毛果楊的基因組比對；

利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，將比對結(jié)果進行染色體坐標排序、reads去重復處理，再通過GATK2(V3.2)對已知胞嘧啶位點附近10bp的序列進行Call SNPs，從而區(qū)分等位基因序列(此時，由于序列經(jīng)由重硫酸鹽甲基化測序，所以胞嘧啶上的甲基化修飾狀況已被記錄)；

將已獲得的被區(qū)分過的子代等位基因序列與父母本基因組序列比對，進而完成胞嘧啶表觀基因型分型如圖3～6所示。圖3是父母本樣品胞嘧啶甲基化表觀基因型分型圖示，X代表胞嘧啶分型片段側(cè)翼序列中的任意堿基(可以為ATCG中任意一個)，m代表修飾在胞嘧啶上的甲基基團，子代胞嘧啶分型片段經(jīng)過與親本比對最終確定的甲基化表觀基因型可分為純合甲基化表觀基因型，雜合甲基化表觀基因型和純合非甲基化型三種。

由以上實施例可知，本發(fā)明提供的基因組胞嘧啶位點表觀基因型分型方法將樣本基因組經(jīng)重亞硫酸鹽處理，全基因組甲基化測序，篩選測定后的胞嘧啶位點，經(jīng)序列比對，SNP Calling后，最終將子代與父母本序列比對完成胞嘧啶表觀基因型分型。本發(fā)明首次完成了基因組胞嘧啶位點表觀基因型分型，技術(shù)成熟，成本低，易于操作與推廣應(yīng)用。

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