CN106845152B - 一种基因组胞嘧啶位点表观基因型分型方法 - Google Patents

Abstract

本发明提供了一种基因组胞嘧啶位点表观基因型分型方法，包括以下步骤：1)对待测样品父母本和子代样本进行重亚硫酸盐全基因组甲基化测序，获得父母本和子代样本基因组序列；2)将父母本和子代样本基因组序列与参考基因组比对获得比对结果，确定待测胞嘧啶位点；3)将比对结果进行染色体坐标排序、reads去重复处理，再通过GATK2‑V3.2对已知胞嘧啶位点上下游5～10bp的序列进行Call SNPs，从而区分子代等位基因序列；4)将已获得的被区分过子代等位基因序列与父母本基因组序列比对，完成胞嘧啶表观基因型分型。本发明首次完成了基因组胞嘧啶位点表观基因型分型，技术成熟，成本低，易于操作与推广应用。

Description

一种基因组胞嘧啶位点表观基因型分型方法

技术领域

本发明涉及分子生物学技术领域，尤其涉及一种基因组胞嘧啶位点表观基因型分型方法。

背景技术

在高等真核生物中，DNA甲基化仅发生在Cp G二核苷酸G5′侧的C上。当位于基因启动子区内富含Cp G序列(Cp G岛)时，这种修饰作用则对基因的表达有着重要的调控作用。此外，它还同基因组印记、女性X染色体的基因灭活、细胞增殖、分化发育、肿瘤的发生和发展以及遗传的不稳定性等密切有关。

近年来甲基化测序技术逐渐发展完善，包括重硫酸盐处理基因组甲基化测序，利用重硫酸盐处理基因组DNA，使未发生甲基化的胞嘧啶脱氨基变成尿嘧啶，而发生了甲基化的胞嘧啶则不会发生变化。比对经重硫酸盐处理和未经处理的样本，可以检测到甲基化位点。进一步结合高通量测序技术，能够从全基因组水平和单碱基精度来分析5＇甲基胞嘧啶，由此能够发现很多传统的基因组学研究所不能检测到的甲基化位点。

目前，甲基化测序数据分析涉及到基因组甲基化水平，甲基化分布类型以及分布倾向等基本特征的分析，但是目前还无法实现表观基因型分型。

发明内容

有鉴于此，本发明的目的在于提供一种基因组胞嘧啶位点表观基因型分型方法。

为了实现上述发明目的，本发明提供以下技术方案：

本发明提供了一种基因组胞嘧啶位点表观基因型分型方法，包括以下步骤：1)对待测样品父母本和子代样本进行重亚硫酸盐全基因组甲基化测序，获得父母本和子代样本基因组序列；2)将父母本和子代样本基因组序列与参考基因组比对获得比对结果，确定待测胞嘧啶位点；3)利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，将比对结果进行染色体坐标排序、reads去重复处理，再通过GATK2-V3.2对已知胞嘧啶位点上下游5～10bp的序列进行Call SNPs，从而区分子代等位基因序列；4)将已获得的被区分过子代等位基因序列与父母本基因组序列比对，完成胞嘧啶表观基因型分型。

优选的，步骤1)中所述的待测样品为有参考基因组的物种。

优选的，所述参考基因组为待测样品本物种已测序基因组或待测样品近缘物种已测序基因组。

优选的，所述步骤1)具体包括以下步骤：1.1)用CTAB法提取待测样品父母本和子代样本基因组DNA；1.2)对提取得到的基因组DNA样品进行质量、纯度和浓度检测筛选获得合格的父母本和子代样本基因组DNA样品；1.3)重亚硫酸盐法构建上述合格的父母本和子代样本基因组DNA样品测序文库；1.4)质检筛选合格基因组DNA样品测序文库，保证所述文库插入片段为320-520bp，文库有效浓度>2nM；1.5)对合格的DNA样品文库进行双末端Hiseq测序，获得父母本和子代样本基因组序列。

优选的，步骤1.3)构建父母本和子代样本基因组DNA样品测序文库时加入比例为建库DNA起始量的1/1000的阴性对照lambdaDNA。

优选的，步骤1.3)构建父母本和子代样本基因组DNA样品测序文库包括以下步骤：随机打断基因组DNA至200～300bp获得DNA片段；对DNA片段进行平末端修复后加尾巴A碱基获得带尾巴A的DNA片段；在带尾巴A的DNA片段上连接测序接头后进行Bisulfite处理；最后进行PCR扩增获得基因组DNA样品测序文库。

优选的，步骤3)中所述染色体坐标排序采用picard-tools工具中的SortSam完成按照染色体坐标顺序从小到大排序。

优选的，步骤3)中所述reads去重复处理采用picard-tools完成。

优选的，步骤3)中所述Call SNPs采用UnifiedGenotyper工具完成。

优选的，步骤3)中所述Call SNPs完成后还包括对所述Call SNPs结果进行过滤。

本发明的有益效果：本发明提供的基因组胞嘧啶位点表观基因型分型方法将样本基因组经重亚硫酸盐处理，全基因组甲基化测序，筛选测定后的胞嘧啶位点，经序列比对，SNP Calling后，最终将子代与父母本序列比对完成胞嘧啶表观基因型分型。本发明首次完成了基因组胞嘧啶位点表观基因型分型，技术成熟，成本低，易于操作与推广应用。

附图说明

图1为无降解，无RNA污染检测合格的DNA样品电泳条带图；

图2为严重降解检测不合格的DNA样品电泳条带图；

图3为父母本胞嘧啶甲基化表观基因型分型示意图；

图4为子代胞嘧啶分型纯合甲基化表观基因型示意图；

图5为子代胞嘧啶分型杂合甲基化表观基因型示意图；

图6为子代胞嘧啶分型纯合非甲基化表观基因型示意图。

具体实施方式

本发明提供了一种基因组胞嘧啶位点表观基因型分型方法，包括以下步骤：1)对待测样品父母本和子代样本进行重亚硫酸盐全基因组甲基化测序，获得父母本和子代样本基因组序列；2)将父母本和子代样本基因组序列与参考基因组比对获得比对结果，确定待测胞嘧啶位点；3)利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，将比对结果进行染色体坐标排序、reads去重复处理，再通过GATK2-V3.2对已知胞嘧啶位点上下游5～10bp的序列进行Call SNPs，从而区分子代等位基因序列；4)将已获得的被区分过子代等位基因序列与父母本基因组序列比对，完成胞嘧啶表观基因型分型。

本发明中所述待测样品优选的为有参考基因组的物种，所述参考基因组为待测样品本物种已测序基因组或待测样品近缘物种已测序基因组。所述待测样品优选的为人类，斑马鱼，毛果样，胡杨，小叶杨或拟南芥等有参考基因组的物种。

本发明在确定待测样品后，对待测样品父母本和子代样本进行重亚硫酸盐全基因组甲基化测序，获得父母本和子代样本基因组序列。本发明中优选的先进行待测样品父母本和子代样本基因组DNA的提取；所述基因组DNA的提取优选的采用CTAB法，本发明中CTAB法提取基因组DNA的具体步骤参见本领域常规的CTAB法提取基因组DNA的步骤，优选的采用文献(Chang et al.,1993)中所记载的方法。

本发明在得到待测样品父母本和子代样本基因组DNA后，对所述的基因组DNA进行质量、纯度和浓度检测，筛选获得合格的父母本和子代样本基因组DNA样品。在本发明中所述基因组DNA质量检测优选的使用琼脂糖凝胶电泳进行，根据电泳条带分析基因组DNA降解程度以及是否有RNA污染。具体的结果如图1和图2，图1中样品为检测合格的无降解，无RNA污染样品电泳条带，图2中样品4号为严重降解；5号有重度RNA污染。本发明筛选质量合格的如图1所示的基因组DNA进行下一步操作，若基因组DNA样品质量不合格需要重新提取基因组DNA。本发明在检测基因组DNA质量合格后，优选的进一步检测其纯度是否合格。在本发明中，优选的采用Nanodrop检测基因组DNA的纯度，当OD260/280为1.6～1.8之间，说明基因组DNA的纯度合格，若所述基因组DNA的OD260/280不在上述范围内，则所述基因组DNA纯度不合格，需要重新提取。本发明在筛选获得质量和纯度合格的基因组DNA后，优选的采用Qubit对基因组DNA浓度进行精确定量，为后续操作提供依据。

本发明在获得用合格的父母本和子代样本基因组DNA样品后，采用重亚硫酸盐法构建上述合格的父母本和子代样本基因组DNA样品测序文库；在本发明中，构建父母本和子代样本基因组DNA样品测序文库时加入比例为建库DNA起始量的1/1000的阴性对照lambdaDNA。所述基因组DNA样品测序文库的构建具体的包括以下步骤：首先随机打断基因组DNA至200～300bp获得DNA片段，所述随机打断基因组DNA优选的使用Covaris S220进行；然后对DNA片段进行平末端修复，所述平末端修复优选的使用内切酶和外切酶将黏性末端修复成为平末端；随后在平末端DNA片段后加尾巴A碱基获得带尾巴A的DNA片段，所述加尾巴A碱基步骤采用本领域常规的手段即可，无其他特殊要求；本发明在带尾巴A的DNA片段上连接测序接头后进行Bisulfite处理；所述的Bisulfite处理优选的采用Zymo Research的EZ DNAMethylation GoldKit进行，本发明中所述Bisulfite处理的作用是将未发生甲基化的C变成U，在后续的PCR扩增过程后变为T，而甲基化的C则保持不变。本发明在Bisulfite处理后进行PCR扩增获得基因组DNA样品测序文库。

本发明在获得基因组DNA样品测序文库后，对文库进行质检筛选获得合格基因组DNA样品测序文库。本发明中在质检文库之前优选的使用Qubit2.0对文库DNA浓度进行初步定量，然后根据文库DNA的浓度将其进行稀释，文库稀释后的浓度优选的为1ng/ul；在文库稀释后，优选的使用Agilent2100对文库的插入片段长度进行检测，筛选后的文库插入片段在320～520bp为合格；使用Q-PCR方法对文库的有效浓度进行准确定量，所述文库有效浓度>2nM为合格。

本发明在获得合格基因组DNA样品测序文库后，对合格的DNA样品文库进行双末端Hiseq测序，获得父母本和子代样本基因组序列。所述的双末端Hiseq测序优选的委托测序公司完成，在本发明具体实施例中具体的委托北京诺禾致源生物信息科技有限公司通过全基因组甲基化测序完成。

本发明在获得父母本和子代样本基因组序列后，将所获子代基因组序列与参考基因组比对，本发明中所述参考基因组优选的为所测样品本物种基因组；如本物种未测序，优选的选取已测序近缘物种基因组做参考基因组，确定待测胞嘧啶位点。

本发明在获得待测胞嘧啶位点后，利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，将比对结果进行染色体坐标排序、reads去重复处理，再通过GATK2-V3.2对已知胞嘧啶位点上下游5～10bp的序列进行Call SNPs，从而区分子代等位基因序列，由于序列经由重亚硫酸盐甲基化测序，所以胞嘧啶上的甲基化修饰状况已被记录。

本发明中上述操作优选的委托北京诺禾致源生物信息科技有限公司完成。具体的操作如下：1)对原始下机fastq文件进行过滤和比对(mapping)，对Illumina下机数据使用bwa进行mapping。

所述Bwa比对步骤如下：Ⅰ.对参考基因组构建索引：Ⅱ.寻找输入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。可以使用picard-tools中的ReorderSam完成。

3)将sam文件转换成bam文件(bam是二进制文件，运算速度快)，这一步可使用samtools view完成。4)对bam文件进行sort排序处理，这一步是将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序。可以使用picard-tools中SortSam完成。5)对bam文件进行加头(head)处理GATK2.0以上版本不支持无头文件的变异检测。加头这一步可以在BWA比对的时候进行，通过-r参数的选择可以完成。如果在BWA比对期间没有选择-r参数，可以增加这一步骤。可使用picard-tools中AddOrReplaceReadGroups完成。6)Merge：如果一个样本分为多个lane进行测序，那么在进行下一步之前可以将每个lane的bam文件合并。7)Duplicates Marking：在制备文库的过程中，由于PCR扩增过程中会存在一些偏差，也就是说有的序列会被过量扩增。这样，在比对的时候，这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。而这些过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由PCR扩增所形成的duplicates，这一步可以使用picard-tools来完成。去重复的过程是给这些序列设置一个flag以标志它们，方便GATK的识别。还可以设置REMOVE_DUPLICATES＝true来丢弃duplicated序列。对于是否选择标记或者删除，对结果应该没有什么影响，GATK官方流程里面给出的例子是仅做标记不删除。这里定义的重复序列是这样的：如果两条reads具有相同的长度而且比对到了基因组的同一位置，那么就认为这样的reads是由PCR扩增而来，就会被GATK标记。8)对上一步得到的结果生成索引文件：可以用samtools完成，生成的索引后缀是bai。9)Local realignmentaround indels：这一步的目的就是将比对到indel附近的reads进行局部重新比对，将比对的错误率降到最低。一般来说，绝大部分需要进行重新比对的基因组区域，都是因为插入/缺失的存在，因为在indel附近的比对会出现大量的碱基错配，这些碱基的错配很容易被误认为SNP。还有，在比对过程中，比对算法对于每一条read的处理都是独立的，不可能同时把多条reads与参考基因组比对来排错。因此，即使有一些reads能够正确的比对到indel，但那些恰恰比对到indel开始或者结束位置的read也会有很高的比对错误率，这都是需要重新比对的。Local realignment就是将由indel导致错配的区域进行重新比对，将indel附近的比对错误率降到最低。主要分为两步：第一步，通过运行RealignerTargetCreator来确定要进行重新比对的区域。第二步，通过运行IndelRealigner在这些区域内进行重新比对。10)Base quality score recalibration：这一步是对bam文件里reads的碱基质量值进行重新校正，使最后输出的bam文件中reads中碱基的质量值能够更加接近真实的与参考基因组之间错配的概率。这一步适用于多种数据类型，包括illunima、solid、454、CG等数据格式。在GATK2.0以上版本中还可以对indel的质量值进行校正，这一步对indel calling非常有帮助BQSR主要有三步：第一步：利用工具BaseRecalibrator，根据一些known sites，生成一个校正质量值所需要的数据文件，GATK网站以“.grp”为后缀命名。第二步：利用第一步生成的ChrALL.100.sam.recal.08-1.grp来生成校正后的数据文件，也是以“.grp”命名，这一步主要是为了与校正之前的数据进行比较，最后生成碱基质量值校正前后的比较图，如果不想生成最后BQSR比较图，这一步可以省略。第三步：利用工具PrintReads将经过质量值校正的数据输出到新的bam文件中，用于后续的变异检测。11)分析和评估BQSR结果：这一步会生成评估前后碱基质量值的比较结果，可以选择使用图片和表格的形式展示。13)Reducebam file：这一步是使用ReduceReads这个工具将bam文件进行压缩，生成新的bam文件，新的bam文件仍然保持bam文件的格式和所有进行变异检测所需要的信息。这样不仅能够节省存储空间，也方便后续变异检测过程中对数据的处理。14)Variant Calling：GATK在这一步里面提供了两个工具进行变异检测——UnifiedGenotyper和HaplotypeCaller。我们使用UnifiedGenotyper。UnifiedGenotyper是集合多种变异检测方法而成的一种VariantsCaller，既可以用于单个样本的变异检测，也可以用于群体的变异检测。UnifiedGenotyper使用贝叶斯最大似然模型，同时估计基因型和基因频率，最后对每一个样本的每一个变异位点和基因型都会给出一个精确的后验概率。15)对原始变异检测结果进行过滤(hardfilter andVQSR)：这一步的目的就是对上一步call出来的变异位点进行过滤，去掉不可信的位点。这一步可以有两种方法，一种是通过GATK的VariantFiltration，另一种是通过GATK的VQSR(变异位点质量值重新校正)进行过滤。通过GATK网站上提供的最佳方案可以看出，GATK是推荐使用VASR的，但使用VQSR数据量一定要达到要求，数据量太小无法使用高斯模型。还有，在使用VAQR时，indel和snp要分别进行。Ⅰ.VariantRecalibrator，VariantRecalibrator：通过大量的高质量的已知变异集合的各个注释(包括很多种，后面介绍)的值来创建一个高斯混合模型，然后用于评估所有的变异位点。这个文件最后将生成一个recalibration文件。Ⅱ.ApplyRecalibration，ApplyRecalibration：这一步将模型的各个参数应用于原始vcf文件中的每一个变异位点，这时，每一个变异位点的注释信息列中都会出现一个VQSLOD值，然后模型会根据这个值对变异位点进行过滤，过滤后的信息会写在vcf文件的filter一列中。

本发明在获得区分子代等位基因序列后，将已获得的被区分过子代等位基因序列与父母本基因组序列比对，完成胞嘧啶表观基因型分型。结果如图3～6所示：图3是父母本样品胞嘧啶甲基化表观基因型分型图示，X代表胞嘧啶分型片段侧翼序列中的任意碱基(可以为ATCG中任意一个)，m代表修饰在胞嘧啶上的甲基基团，子代胞嘧啶分型片段经过与亲本比对最终确定的甲基化表观基因型可分为纯合甲基化表观基因型，杂合甲基化表观基因型和纯合非甲基化型三种。

下面结合实施例对本发明提供的基因组胞嘧啶位点表观基因型分型方法进行详细的说明，但是不能把它们理解为对本发明保护范围的限定。

实施例1

以人类为研究对象，将包含父母本及其子代在内的样本进行重硫酸盐全基因组甲基化测序，获得基因组序列；

1)样本基因组DNA提取，使用CTAB法；2)DNA样品检测，使用琼脂糖凝胶电泳分析DNA降解程度以及是否有RNA污染，检测参考图1和图2，图1中样品为检测合格的样品电泳条带如下：无降解，无RNA污染，图2中样品4号为严重降解；5号有重度RNA污染。再通过Nanodrop检测DNA纯度(OD260/280)，最后使用Qubit对DNA浓度进行精确定量；3)文库构建，样品检测为A或B时认为合格，合格后加入比例为建库起始量的1/1000的阴性对照(lambdaDNA)，首先使用Covaris S220将基因组DNA随机打断至200-300bp。对打断后的DNA片段进行末端修复、加尾巴A，并连接上所有胞嘧啶均经过甲基化修饰的测序接头。随后进行Bisulfite处理(采用EZ DNAMethylation GoldKit，Zymo Research)，经过处理，未发生甲基化的C变成U(PCR扩增后变为T)，而甲基化的C保持不变，最后进行PCR扩增，得到最终的DNA文库；4)文库质检，文库构建完成后，先使用Qubit2.0进行初步定量，稀释文库至1ng/ul，随后使用Agilent2100对文库的插入片段长度进行检测，符合预期后，使用Q-PCR方法对文库的有效浓度进行准确定量(文库有效浓度>2nM)，以保证文库质量；5)上机测序，库检合格后，把不同文库按照有效浓度及目标下机数据量的需求，pooling后进行Hiseq测序，测序策略为双末端测序。

将所获基因组序列包含父母本与子代基因组与人类基因组比对；

利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，将比对结果进行染色体坐标排序、reads去重复处理，再通过GATK2(V3.2)对已知胞嘧啶位点附近10bp的序列进行Call SNPs，从而区分等位基因序列(此时，由于序列经由重硫酸盐甲基化测序，所以胞嘧啶上的甲基化修饰状况已被记录)；

将已获得的被区分过的子代等位基因序列与父母本基因组序列比对，进而完成胞嘧啶表观基因型分型如图3～6所示。图3是父母本样品胞嘧啶甲基化表观基因型分型图示，X代表胞嘧啶分型片段侧翼序列中的任意碱基(可以为ATCG中任意一个)，m代表修饰在胞嘧啶上的甲基基团，子代胞嘧啶分型片段经过与亲本比对最终确定的甲基化表观基因型可分为纯合甲基化表观基因型，杂合甲基化表观基因型和纯合非甲基化型三种。

实施例2

以胡杨为研究对象，将包含父母本及其子代在内的样本进行重硫酸盐全基因组甲基化测序，获得基因组序列；

1)样本基因组DNA提取，使用CTAB法；2)DNA样品检测，使用琼脂糖凝胶电泳分析DNA降解程度以及是否有RNA污染，检测参考图1和图2，图1中样品为检测合格的样品电泳条带如下：无降解，无RNA污染，图2中样品4号为严重降解；5号有重度RNA污染。再通过Nanodrop检测DNA纯度(OD260/280)，最后使用Qubit对DNA浓度进行精确定量；3)文库构建，样品检测为A或B时认为合格，合格后加入比例为建库起始量的1/1000的阴性对照(lambdaDNA)，首先使用Covaris S220将基因组DNA随机打断至200-300bp。对打断后的DNA片段进行末端修复、加尾巴A，并连接上所有胞嘧啶均经过甲基化修饰的测序接头。随后进行Bisulfite处理(采用EZ DNAMethylation GoldKit，Zymo Research)，经过处理，未发生甲基化的C变成U(PCR扩增后变为T)，而甲基化的C保持不变，最后进行PCR扩增，得到最终的DNA文库；4)文库质检，文库构建完成后，先使用Qubit2.0进行初步定量，稀释文库至1ng/ul，随后使用Agilent2100对文库的插入片段长度进行检测，符合预期后，使用Q-PCR方法对文库的有效浓度进行准确定量(文库有效浓度>2nM)，以保证文库质量；5)上机测序，库检合格后，把不同文库按照有效浓度及目标下机数据量的需求，pooling后进行Hiseq测序，测序策略为双末端测序。

将所获基因组序列包含父母本与子代基因组与胡杨基因组比对；

利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，将比对结果进行染色体坐标排序、reads去重复处理，再通过GATK2(V3.2)对已知胞嘧啶位点附近10bp的序列进行Call SNPs，从而区分等位基因序列(此时，由于序列经由重硫酸盐甲基化测序，所以胞嘧啶上的甲基化修饰状况已被记录)；

将已获得的被区分过的子代等位基因序列与父母本基因组序列比对，进而完成胞嘧啶表观基因型分型如图3～6所示。图3是父母本样品胞嘧啶甲基化表观基因型分型图示，X代表胞嘧啶分型片段侧翼序列中的任意碱基(可以为ATCG中任意一个)，m代表修饰在胞嘧啶上的甲基基团，子代胞嘧啶分型片段经过与亲本比对最终确定的甲基化表观基因型可分为纯合甲基化表观基因型，杂合甲基化表观基因型和纯合非甲基化型三种。

实施例3

以小叶杨为研究对象，将包含父母本及其子代在内的样本进行重硫酸盐全基因组甲基化测序，获得基因组序列；

1)样本基因组DNA提取，使用CTAB法；2)DNA样品检测，使用琼脂糖凝胶电泳分析DNA降解程度以及是否有RNA污染，检测参考图1和图2，图1中样品为检测合格的样品电泳条带如下：无降解，无RNA污染，图2中样品4号为严重降解；5号有重度RNA污染。再通过Nanodrop检测DNA纯度(OD260/280)，最后使用Qubit对DNA浓度进行精确定量；3)文库构建，样品检测为A或B时认为合格，合格后加入比例为建库起始量的1/1000的阴性对照(lambdaDNA)，首先使用Covaris S220将基因组DNA随机打断至200-300bp。对打断后的DNA片段进行末端修复、加尾巴A，并连接上所有胞嘧啶均经过甲基化修饰的测序接头。随后进行Bisulfite处理(采用EZ DNAMethylation GoldKit，Zymo Research)，经过处理，未发生甲基化的C变成U(PCR扩增后变为T)，而甲基化的C保持不变，最后进行PCR扩增，得到最终的DNA文库；4)文库质检，文库构建完成后，先使用Qubit2.0进行初步定量，稀释文库至1ng/ul，随后使用Agilent2100对文库的插入片段长度进行检测，符合预期后，使用Q-PCR方法对文库的有效浓度进行准确定量(文库有效浓度>2nM)，以保证文库质量；5)上机测序，库检合格后，把不同文库按照有效浓度及目标下机数据量的需求，pooling后进行Hiseq测序，测序策略为双末端测序。

因小叶杨基因组尚未被测序，所以将所获基因组序列包含父母本与子代基因组与小叶杨近缘物种毛果杨的基因组比对；

利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，将比对结果进行染色体坐标排序、reads去重复处理，再通过GATK2(V3.2)对已知胞嘧啶位点附近10bp的序列进行Call SNPs，从而区分等位基因序列(此时，由于序列经由重硫酸盐甲基化测序，所以胞嘧啶上的甲基化修饰状况已被记录)；

将已获得的被区分过的子代等位基因序列与父母本基因组序列比对，进而完成胞嘧啶表观基因型分型如图3～6所示。图3是父母本样品胞嘧啶甲基化表观基因型分型图示，X代表胞嘧啶分型片段侧翼序列中的任意碱基(可以为ATCG中任意一个)，m代表修饰在胞嘧啶上的甲基基团，子代胞嘧啶分型片段经过与亲本比对最终确定的甲基化表观基因型可分为纯合甲基化表观基因型，杂合甲基化表观基因型和纯合非甲基化型三种。

由以上实施例可知，本发明提供的基因组胞嘧啶位点表观基因型分型方法将样本基因组经重亚硫酸盐处理，全基因组甲基化测序，筛选测定后的胞嘧啶位点，经序列比对，SNP Calling后，最终将子代与父母本序列比对完成胞嘧啶表观基因型分型。本发明首次完成了基因组胞嘧啶位点表观基因型分型，技术成熟，成本低，易于操作与推广应用。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (8)

1.一种基因组胞嘧啶位点表观基因型分型方法，其特征在于，包括以下步骤：

1)对待测样品父母本和子代样本进行重亚硫酸盐全基因组甲基化测序，获得父母本和子代样本基因组序列；

2)将所述父母本和子代样本基因组序列与待测样品已知的参考基因组比对，确定待测胞嘧啶位点；

3)利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96，将所述待测胞嘧啶位点的序列数据进行染色体坐标排序后进行reads去重复处理，获得校正后的序列数据，再通过GATK2-V3.2对所述校正后的序列数据中已知胞嘧啶位点上下游5～10bp的序列进行CallSNPs，区分出子代等位基因序列；

4)将所述子代等位基因序列与父母本基因组序列比对，得到胞嘧啶表观基因型分型结果；

所述步骤1)具体包括以下步骤：

1.1)用CTAB法提取待测样品父母本和子代样本基因组DNA；

1.2)对所述提取到的基因组DNA进行质量、纯度和浓度检测，筛选出合格的父母本和子代样本基因组DNA样品；

1.3)采用重亚硫酸盐法构建所述合格的父母本和子代样本基因组DNA样品测序文库；

1.4)质检筛选合格基因组DNA样品测序文库，所述合格基因组DNA样品测序文库的插入片段为320～520bp，有效浓度＞2nM；

1.5)对所述合格基因组DNA样品测序文库进行双末端Hiseq测序，获得父母本和子代样本基因组序列。

2.根据权利要求1所述的方法，其特征在于，所述参考基因组为待测样品本物种已测序基因组或待测样品近缘物种已测序基因组。

3.根据权利要求1所述的方法，其特征在于，步骤1.3)构建父母本和子代样本基因组DNA样品测序文库时加入建库DNA起始量1/1000的阴性对照lambdaDNA。

4.根据权利要求1或3所述的方法，其特征在于，步骤1.3)构建测序文库包括以下步骤：

随机打断基因组DNA至200～300bp，获得DNA片段；

对所述DNA片段进行平末端修复后加尾巴A碱基，获得带尾巴A的DNA片段；

在所述带尾巴A的DNA片段上连接测序接头后进行Bisulfite处理；

将所述Bisulfite处理后的DNA片段进行PCR扩增，获得基因组DNA样品测序文库。

5.根据权利要求1所述的方法，其特征在于，步骤3)中所述染色体坐标排序为：采用picard-tools工具中的SortSam按照染色体坐标顺序从小到大排序。

6.根据权利要求1所述的方法，其特征在于，步骤3)中所述reads去重复处理采用picard-tools完成。

7.根据权利要求1所述的方法，其特征在于，步骤3)中所述CallSNPs采用UnifiedGenotyper工具完成。

8.根据权利要求1所述的方法，其特征在于，步骤3)中所述CallSNPs后还包括：对所述CallSNPs结果进行过滤；所述过滤采用GATK工具的变异位点质量值重新校正-VQSR进行。