CN105046105A - 染色体跨度的单体型图及其构建方法 - Google Patents

Abstract

本发明提供了一种染色体跨度的单体型图及其构建方法。构建方法包括对包含待测样本基因组交联位点的测序文库进行测序，得到PE？reads；将PE？reads分别与参考基因组上的序列进行第一次比对，得到第一比对结果；根据第一比对结果构建PE？reads的一致性序列，获得高质量的SNP位点；根据所获取的高质量的SNP位点对每一条reads进行筛选，并从中提取至少含有2个杂合SNP位点的reads；并利用杂合SNP位点构建染色体跨度的单体型图。这种基于全基因组的交联位点的测序数据更全面，因而得到的染色体跨度的单体型图的准确性和分辨率都很高，且具有物种适应范围广的优势。

Description

染色体跨度的单体型图及其构建方法

技术领域

本发明涉及单体型图谱构建领域，具体而言，涉及一种染色体跨度的单体型图及其构建方法。

背景技术

SNP是指在染色体基因组水平上单个核苷酸的变异引起的DNA序列多态性。位于一条染色体上的或某一区域的一组SNP等位位点被称为单体型(haplotype)。单体型是单倍体型的简称，是指在同一染色体上进行共同遗传的多个基因座上等位基因的组合。按照某一指定基因座上基因重组发生的数量，单体型甚至可以指至少两个基因座或整个染色体。高质量的单体型在物种进化历史、表型性状的关联分析等方面有重要的价值。

目前，国内在单体型的构建方面，采用的方法主要是基于对长片段序列(长片段序列是指从头测到尾的几Kb～几十Kb的片段)、大片段文库(大片段文库是指2K、5K、10K或20K等长度的mate-pair文库，一般进行双末端测序，例如，用PE100的测序策略，那么测出的数据是首端100bp和尾端的100bp)、fosmid文库等的测序数据组装构建单体型，其准确性和分辨率非常高，但是所构建的单体型的长度仅有几百Kb～1Mb。

此外，也有染色体跨度的单体型构建的报道，但这种染色体跨度的单体型的构建方法是基于染色体显微切割或者单精子等单倍体测序来构建的，其构建方法具体如下：1)首先，通过显微切割获得微量的单倍体DNA遗传物质；2)采用MALBAC或者多重置换扩增(MDA)等全基因组扩增技术提高DNA的起始量；3)构建小片段文库，进行二代测序获得短序列数据(PEreads)；4)短序列数据(PEreads)与已经发表的参考基因组比对，检测单核苷酸位点多态性(SNPs)，构建单倍体的单体型。

上述染色体跨度的单体型构建方法，对于实验技术和实验仪器有着很高的要求，在取材方面非常受限，只能用于一些符合严格减数分裂的物种。最重要的是，为了达到建库的最低起始量，需要对DNA进行扩增。在扩增过程中会引入一些错误，并且随着不断的扩增会导致误差的进一步放大。而且，扩增过程中的碱基偏好性都会影响SNP检测的准确性和全面性，从而对最终单体型构建的结果造成一定的误差。此外，上述方法整个实验过程时间较长，通常需要3个月以上。

因此，仍需要对现有技术进行改进，以提供一种适用范围广、准确性高、染色体跨度大的个体单体型的构建方法。

发明内容

本发明的主要目的在于提供一种染色体跨度的单体型图及其构建方法，以提供一种适用范围广、准确性高以及染色体跨度大的个体单体型图的构建方法。

为了实现上述目的，根据本发明的一个方面，提供了一种染色体跨度的单体型图的构建方法，该构建方法包括：步骤S1，对包含待测样本基因组交联位点的测序文库进行测序，得到PEreads；步骤S2，将PEreads分别与参考基因组上的序列进行第一次比对，得到第一比对结果；步骤S3，根据第一比对结果构建来源于待测样本的PEreads的一致性序列，获得高质量的SNP位点；以及步骤S4，根据所获取的高质量的SNP位点对每一条reads进行筛选，并提取至少含有2个杂合SNP位点的reads；并利用杂合SNP位点构建染色体跨度的单体型图。

进一步地，在步骤S2之前，构建方法还包括判断测序文库是否合格的步骤；判断测序文库是否合格的步骤包括：步骤a，抽取部分PEreads作为待质控文库；步骤b，将待质控文库与参考基因组上的序列进行第二次比对，得到第二比对结果；步骤c，利用第二比对结果计算待质控文库的插入片段在参考基因组上的跨度；步骤d，若跨度大于10kb，且跨度大于10kb的插入片段在待质检文库中的比例为50％以上，则判断测序文库合格。

进一步地，在步骤S1之后，以及在步骤S2之前，构建方法还包括对PEreads进行质控的步骤；质控的步骤包括：检测PEreads中是否存在外源样本污染，和/或对PEreads进行低质量数据过滤。

进一步地，在步骤S2之前，还包括对参考基因组上的序列进行前处理的步骤；前处理的步骤包括：基于BWA和SAMtools对参考基因组上的序列进行索引文件的构建；对参考基因组上的序列的长度、碱基含量和空缺的比例进行统计；以及获取构建测序文库时所使用的限制性内切酶在参考基因组的序列上的分布位置和数目。

进一步地，步骤S2包括：利用BWA软件的mem模块，将PEreads通过第一次比对分别比对到建好索引文件的参考基因组的序列上，得到第一比对结果。

进一步地，步骤S2在得到第一比对结果后，还包括对PEreads进行覆盖均匀性进行检测的步骤，覆盖均匀性检测的步骤包括：检测在不同测序深度下，PEreads对参考基因组的覆盖程度和覆盖深度，覆盖程度是指PEreads覆盖参考基因组上的序列的长度与参考基因组上的序列的总长度的比值；覆盖深度是指参考基因组相应位置被PEreads覆盖到的次数。

进一步地，步骤S3包括：根据第一比对结果，利用SAMtools软件构建PEreads的一致性序列，以及获取一致性序列中存在的高质量的SNP位点。

进一步地，步骤S4包括：从高质量的SNP位点中提取至少含有2个杂合SNP位点的reads；以杂合SNP位点为点，以两个杂合SNP之间的片段长度为边，构建单体型块；基于交联位点之间的线性距离小于30Mb的两侧片段上的杂合SNP位点，将多个单体型块进行连接，得到候选单体型图；对候选单体型图的完整性和分辨率进行评估，并对评估结果依次按照完整性和分辨率的高低进行排序，挑取在完整性最高条件下分辨率最高的候选单体型图作为染色体跨度的单体型图；其中，挑选交联位点之间的线性距离小于30Mb的两侧片段的步骤包括：根据染色体的位置，对第一比对结果进行排序，得到PEreads中的每一条reads在参考基因组的序列上的位置信息；根据PEreads中的每一条reads在参考基因组的序列上的位置信息，将构成交联位点两端的reads进行配对，得到交联位点在参考基因组上的线性距离；根据线性距离，挑选得到交联位点之间的线性距离小于30Mb的两侧片段。

进一步地，在步骤S1之前，构建方法还包括：利用空间构象捕获的方法构建得到包含待测样本基因组交联位点的测序文库；利用空间构象捕获的方法构建测序文库的步骤包括：对待测样本的DNA进行交联固定；利用限制性内切酶对DNA进行酶切，产生具有酶切缺口的DNA片段；利用生物素标记的寡核苷酸将具有酶切缺口的DNA片段上的酶切缺口进行补平；利用核酸连接酶将DNA片段进行连接，得到连接交联DNA片段；对连接交联DNA片段进行解交联，得到解交联DNA；以及对解交联DNA进行片段化文库构建，得到测序文库。

根据本发明的另一方面，提供了一种染色体跨度的单体型图，该单体型图利用上述任一种构建方法构建而成。

应用本发明的技术方案，通过以待测样本全基因组上带有交联位点的测序数据为处理对象，能够获取全基因组上存在的重组交换位点，利用这些重组交换位点中的SNP位点构建而成的单体型图在基因组上的跨度较大，能够达到染色体跨度水平，且这种基于全基因组的交联位点的测序数据更全面，因而得到的染色体跨度的单体型图的准确性和分辨率都很高。相比现有技术的构建方法，该方法不局限于某些特定的物种，具有物种适应范围广的优势；且实验过程避免了扩增带来的偏好性及误差，相对提高了本发明所构建的单体型图的准确性和分辨率。

附图说明

构成本申请的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1和图2示出了根据本发明的优选实施例中将测序数据与参考基因组上的序列进行比对并构建成一致性序列的示意图；其中，

图1示出了具体测序序列与参考基因组上的序列的详细序列比对结果，其中第一行表示参考基因组上的序列的碱基位置，第二行表示参考基因组上的序列的碱基信息；第三行表示一致性序列；第四行及其以下行表示测序所得序列；

图2示出了存在SNP位点的序列比对信息，其中，左侧框中表示第一个单体型块；右侧框中表示另一个单体型块；

图3示出了在本发明一种优选的实施例中单体型块及候选单体型图的构建流程示意图；

图4示出了多个候选单体型图的各种参数值；

图5示出了本发明的一种优选实施例中的构建方法所构建的单体型图在完整性方面与现有技术的比对结果图；

图6示出了本发明的一种优选实施例中的构建方法所构建的单体型图在准确度方面与现有技术的比对结果图；以及

图7示出了本发明的一种优选实施例中的构建方法所构建的单体型图在分辨率方面与现有技术的比对结果图。

具体实施方式

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将结合实施例来详细说明本发明。

名称解释：

单体型图：单体型图是指在同一染色体上进行共同遗传的多个基因座上等位位点的基因型的直线排列图。

Read:指高通量测序时，测序仪器每读取一次所能读取到的待测序列的长度，后把每次读取得到的一个短片段的序列称为一个read。

双端reads(PEreads)：双末端测序序列，是指对所构建的测序文库中的插入片段从两端进行测序，得到包括插入片段两端序列信息的测序序列。

单端reads:是指双端测序序列中某一端的测序序列。

一致性序列：Consensussequence，通过序列比较发现相似但不一定完全相同的核苷酸序列或氨基酸序列。在本发明中，即个体测序序列在与参考序列相同位置上的基因型所组成的序列。

一致性序列的构建：根据参考基因组上的位置及相应位置出现频率最大的测序序列基因型，将位置与参考序列一致的测序个体的碱基序列构建成一致性序列。

杂合SNP位点：Singlenucleotidepolymorphism，即单核苷酸多态性，是指测序个体基因型和参考序列不一致的位点，对于二倍体，是指来源于父母本的两个碱基型不一致的位点。

Block：即单体型块，由于有些杂合SNP位点的深度不足或者不同杂合SNP位点之间没有捕获到交联序列，导致单体型断成多个块。

定型位点数：即phase的位点数目，在相应的单体型block上可以被明确区分到两个单体型上的杂合SNP位点的数目。

覆盖深度：指参考基因组序列相应位置被测序序列覆盖到的次数。

覆盖程度：指被覆盖的基因组序列的长度与参考基因组总长度的比值。

低质量数据：包括(1)带接头(adapter)的readspair；(2)当单端的read中碱基类型为N的数目占该条read总碱基数目的比例超过10％时，包含该单端read的pairedreads都属于低质量数据；(3)当单端测序read中含有的低质量(Q<＝5)碱基数超过该条read长度比例的50％时，包含该单端read的pairedreads都属于低质量数据。

高质量的SNP位点：高质量的SNP位点是指覆盖深度≧4，比对质量≧20的SNP位点。

比对质量：高通量测序中，每测一个碱基会给出一个相应的质量值，这个质量值是衡量测序准确度的。如果测序错误率用e表示，illuminaHiSeqTM2500/MiSeq的碱基质量值用Qphred表示，则碱基比对质量值与测序错误率之间存在下列关系：Qphred＝-10log10(e)；碱基的测序质量值为13时，对应的测序错误率为5％；碱基的测序质量值为20时，对应的测序错误率为1％；碱基的测序质量值为30时，对应的测序错误率为0.1％。Q20与Q30则表示质量值大于等于20或30的碱基所占百分比。比如一共测了1G的数据量，其中有0.9G的碱基质量值大于或等于20，那么Q20则为90％。

由于现有技术中染色体跨度的单体型构建方法，在实验技术上具有很大的难度，且广谱适用性差，不适合大多数的物种。为了改善上述技术问题，在本发明一种典型的实施方式中，提供了一种构建染色体跨度的单体型图的方法，该构建方法包括：步骤S1，对包含待测样本基因组交联位点的测序文库进行测序，得到PEreads；步骤S2，将PEreads分别与参考基因组上的序列进行第一次比对，得到第一比对结果；步骤S3，根据第一比对结果，构建PEreads的一致性序列，获得高质量的SNP位点；步骤S4，根据所获取的高质量的SNP位点对每一条reads进行筛选，并从中获取至少含有2个杂合SNP位点的reads；并利用杂合SNP位点构建染色体跨度的单体型图。

本发明的上述方法通过以待测样本全基因组上带有交联位点的测序数据为处理对象，能够获取全基因组上存在的重组交换位点，利用这些重组交换位点中的SNP位点构建而成的单体型图在基因组上的跨度较大，能够达到染色体跨度水平，且这种基于全基因组的交联位点的测序数据更全面，因而得到的染色体跨度的单体型图的准确性和分辨率都很高。相比现有技术的构建方法，该方法不局限于某些特定的物种，具有物种适应范围广的优势；且实验过程避免了扩增带来的偏好性及误差，相对提高了本发明所构建的单体型图的准确性和分辨率。

上述第一次比对只要是能够将来源于待测样本的测序数据---PEreads与参考基因组上的序列进行比对即可，即可为单端比对亦可为双端比对。

上述构建方法中，在对包含待测样本基因组交联位点的测序文库进行测序之前，在构建文件的步骤中会对各步骤进行监控，以增加文库构建成功的概率。同样，为了进一步确保文库所产出的数据能够用于单体型构建，在本发明一种优选的实施例中，在上述步骤S2之前，该构建方法还包括判断测序文库是否合格的步骤；其中，判断测序文库是否合格的步骤包括：步骤a，抽取部分PEreads作为待质检文库；步骤b，将待质检文库与参考基因组上的序列进行第二次比对，得到第二比对结果；步骤c，利用第二比对结果计算待质检文库的插入片段在参考基因组上的跨度；步骤d，若跨度大于10kb，且跨度大于10kb的插入片段在待质检文库中插入片段总数的比例为50％以上，则判断测序文库合格。

上述优选实施例中，通过对测序数据中的一部分进行抽查，并根据所构建的文库的插入片段在基因组上的染色体跨度大小是否大于10kb来判断所构建的文库是否合格。由于基因上存在交联位点的片段之间的物理距离相对较大，而小于10kb的插入片段可能不是交联位点两端的序列。因而通过该步骤对测序文库的质控检测，在确定所用的文库的数据相对可靠的基础上进行后续分析，有助于提高最终所构建的单体型的准确度。上述第二次比对为双端比对。

在上述构建方法中，从插入片段在染色体上的跨度大小来监测文库已经能够提到所用测序数据的可靠性。为了进一步提高后续步骤所用数据的可靠性和准确性，在本发明另一种优选的实施例中，在上述步骤S1之后，以及在步骤S2之前，该构建方法还包括对PEreads进行质控的步骤；且质控的步骤包括：检测PEreads中是否存在外源样本污染，和/或对PEreads进行低质量数据过滤。

由于上机测序之前的文库混合步骤的误差、测序过程中测序仪器的稳定性或者其他人为因素均可能导致测序数据质量的变化，因而上述优选实施例中，通过检测测序数据中是否存在其他外源样本的序列信息来判断待测样本是否存在污染，和/或去除低质量的测序数据来进一步提高可用数据的质量。其中，所需要去除的低质量数据包括：构建文库时所添加的接头序列以及低质量的PEreads；属于低质量PEreads的情形包括：其中一个单端read中碱基类型为N(表示不确定)的碱基数占该单端read碱基数的比例超过10％时；或者，单端read中含有的低质量(Q<＝5)碱基数占该单端read碱基数的比例超过50％时。

在上述构建方法中，将待测样本的测序数据与参考基因组上的序列进行比对的步骤可以采用现有技术的比对方法，但为了使比对操作更方便，在本发明又一种优选的实施例中，在步骤S2之前，还包括对参考基因组上的序列进行前处理的步骤；前处理的步骤包括：基于BWA和SAMtools对参考基因组上的序列进行索引文件的构建；对参考基因组上的序列的长度、碱基含量和空缺(gap)的比例进行统计；以及获取构建测序文库时所使用的限制性内切酶在参考基因组的序列上的分布位置和数目。

上述优选的实施例中，前处理的步骤中，建立索引文件有助于对基因组各碱基的位置和类型进行明确标记，以便后续很容易读取测序序列的位置和碱基类型；而对参考基因组上的序列的长度、碱基含量和空缺(gap)的比例进行统计便于掌握基因组的组装程度，用于后续对覆盖度进行统计；获取构建测序文库时所使用的限制性内切酶在参考基因组的序列上的分布位置和数目的作用是用于在测序文库质控时判断含有的限制性内切酶的酶切位点的序列是否是交联形成的。该前处理的步骤可以在得到测序数据之前或之后，只要在进行比对之前进行即可。

上述步骤S2中，采用现有可用的比对软件能够实现对比的功能即可，比如可用基于Novoalign软件进行测序数据的比对。Novoalign软件在比对时，对未能对比到参考基因组上的测序片段所采取的策略是，将该测序片段尾端剪切掉5bp，然后再进行比对，如此操作，直到能够与参考基因组比对上，或者剪切后的片段的长度小于25bp时，停止比对。基于该软件的上述比对方式完全能够实现比对的功能，但为了进一步提高比对的质量和比对数据的有效性，在本发明一种优选的实施例中，上述步骤S2包括：利用BWA软件的mem模块，将PEreads通过第一次比对分别比对到建好索引文件的参考基因组的序列上，得到第一比对结果。

上述优选实施例中通过采用BWA软件进行比对，最小的BWA软件中的mem模块能够处理含有交联位点的跳跃序列，即同一条reads的前后两部分来源于基因组上物理距离较远的两个位置，并可以保留所有比对的结果，一方面提高相应比对上的位点的reads支持度，另外可以获得更多的杂合模块片段，有效提高数据的可利用率。

在上述优选的实施例中，根据参考基因组上的位置及相应位置上出现频率最大的read的基因型，构建出位置与参考序列一致的待测样本的碱基序列，这样便建立了待测样本与参考基因组的一致性序列。

在上述优选实施例中，在得到第一比对结果之后，还包括对PEreads进行覆盖均匀性检测的步骤，优选覆盖均匀性检测的步骤包括：检测在不同测序深度下，PEreads对参考基因组的覆盖程度和覆盖深度，覆盖程度是指PEreads覆盖参考基因组上的序列的长度与参考基因组上的序列的总长度的比值；覆盖深度是指参考基因组相应位置被PEreads覆盖到的次数。

通过在如1×、4×、10×以及20×等不同的测序深度下，对PEreads的覆盖均匀性进行检测，有助于进一步确定所得到的测序数据的可靠性和准确性，进而有助于提高后续所构建的单体型图的准确度和分辨率。

上述构建方法中，构建一致性序列的步骤采用常规的构建步骤进行构建即可。在本发明一种优选的实施例中，上述步骤S3包括：根据第一比对结果，利用SAMtools软件构建来源待测样本的PEreads的一致性序列，以及获取一致性序列中存在的高质量的SNP位点。

上述优选实施例中，在利用BWA软件进行对比得到比对结果的基础上，能够用来构建一致性序列的reads的数量相对较多，且未经过末端剪切，且各read的长度也相对较长，因而通过利用SAMtools软件所构建的一致性序列的长度也相对较长。同时，也大大增加了SNP的检测几率。另外，通过对一致性序列中的SNP位点的质量进行检测，获取得到覆盖次数相对较高的SNP用于后续的单体型构建，能够提高单体型图的准确度和精密度。选择覆盖深度至少4次以上的SNP位点，是从根据经验和概率的角度考虑，4次以下的SNP位点的可靠性和准确性较低，利用这样的SNP位点容易降低单体型图的准确度。

在本发明的上述构建方法中，在得到高质量的SNP位点之后，从高质量的SNP位点中提取至少含有2个杂合SNP位点的reads；利用这些reads中的杂合位点构建单体型图的步骤采用本领域的常规构建步骤即可。在本发明一种优选的实施例中，利用至少含有2个杂合SNP位点的reads的杂合位点构建单体型图的步骤包括：以杂合SNP位点为点，以两个杂合SNP之间的片段长度为边，构建单体型块；基于交联位点之间的线性距离小于30Mb的两侧片段上的杂合SNP位点，将多个单体型块进行连接，得到候选单体型图；对候选单体型图的完整性和分辨率进行评估，并对评估结果依次按照完整性和分辨率的高低进行排序，挑取在完整性最高条件下分辨率最高的候选单体型图作为染色体跨度的单体型图；其中，挑选交联位点之间的线性距离小于30Mb的两侧片段的步骤包括：根据染色体的位置，对第一比对结果进行排序，得到每一条reads在参考基因组的序列上的位置信息；根据每一条reads在参考基因组的序列上的位置信息，将构成交联位点两端的reads进行配对，得到交联位点在参考基因组上的线性距离；根据线性距离，挑选得到交联位点之间的线性距离小于30Mb的两侧片段。

由于本发明的单体型图的构建主要是在存在交联位点的测序片段之间建立联系。若某个测序片段与其余测序片段之间没有交联位点，则该测序片段上所含有的SNP位点无法与其他测序片段上的SNP位点联系起来。因而在构建单体型图的过程中被无交联位点的片段割裂成一个个的单体型块。为了得到完整性和分辨率更高的单体型图，在本发明的上述优选实施例中，通过对测序数据中交联位点之间线性距离小于30Mb的片段进行判断，在将短的单体型块连接成更长的候选单体型图的步骤中，选择那些存在交联位点且交联位点之间的线性距离小于30Mb的片段，对短的单体型块进行连接进而形成多个候选单体型图。最后，通过对多个候选单体型图的完整性和分辨率进行评估排序，并选择完整性最高条件下分辨率相对最高的候选单体型图来作为本发明所称的染色体跨度的单体型图。这种利用未剪切过末端的测序数据直接比对得到的SNP位点构建而成的单体型图，具有数据利用率高，所形成的单体型图的完整性好，分辨率更高的优势。

上述优选实施例中，选择线性距离小于30Mb的两侧片段的原因在于，发明人根据以前的研究得知，该线性距离的长度越长，姐妹染色单体之间发生交联(h-trans)的概率越大，H-trans会影响单体型构建的准确度。而该长度在30Mb以下的时候，h-trans的比例≤5％，有利于提高单体型图构建的准确度。

在本发明的上述构建方法中，得到包含待测样本基因组交联位点的测序文库的方法有很多种，具体制备方法可根据具体研究目的的不同选择合适的方法制备得到。为了进一步降低操作难度、提高物种适用范围，提高单体型图的准确性，在本发明一种优选的实施例中，在步骤S1之前，上述构建方法还包括：利用空间构象捕获的方法构建包含待测样本基因组交联位点的测序文库。

通过采用空间构象捕获的方法所构建的测序文库能够覆盖整个待测样本的全基因组上所有交联位点的序列信息，因而这种方法所得到的测序文库的测序数据所涵盖的序列更全面，构建单体型图所可用的SNP位点的数量和质量也相对较高，有助于提高单体型图的准确度。

上述利用空间构象捕获的方法构建测序文库的步骤采用现有的制备步骤进行制备即可。在本发明另一种优选的实施例中，上述利用空间构象捕获的方法构建测序文库的步骤包括：A，对待测样本的DNA进行交联固定；B，利用限制性内切酶对DNA进行酶切，产生具有酶切缺口的DNA片段；C，利用生物素标记的寡核苷酸将具有酶切缺口的DNA片段上的酶切缺口进行补平；D，利用核酸连接酶将DNA片段进行连接，得到连接交联DNA片段；E，对连接交联DNA片段进行解交联，得到解交联DNA；以及F，对解交联DNA进行片段化文库构建，得到测序文库。

上述空间构象捕获的方法(Hi-C)所构建的测序文库能够获得在基因组上线性距离分布广泛的PEreads。通过利用上述构建方法进行单体型图的构建，不仅有利于快速、高效地获得单体型图，而且所得到的单体型图的准确性和分辨率都较高，且跨度能够达到染色体水平。

在本发明的上述构建方法中，利用至少含有2个杂合SNP位点的单端reads中的杂合SNP位点构建染色体跨度的单体型图的步骤可以采用现有的方法进行构建。在本发明一种优选的实施例中，上述构建染色体跨度的单体型图的步骤包括：以杂合SNP为点，两个杂合SNP之间的片段为边构建候选单体型图，运用greedymax-cutheuristic算法和MEC打分函数对候选单体型图中的每部分的block进行确定及评估，得到多个候选单体型的block；基于交联位点距离小于30M的两侧片段上的杂合SNPs将block进行连接，构建得到染色体跨度的候选单体型图；然后根据基因组长度、block物理位置跨度、SNPs数目、定型位点数目以及MEC分值对每个block的完整性和分辨率进行比对，选择比对结果最优的单体型图作为最终的染色体跨度的单体型图。

在本发明另一种典型的实施方式中，提供了一种染色体跨度的单体型图，该单体型图利用上述任一种构建方法构建而成。由上述方法构建得到的单体型图不仅可以达到染色体的跨度，而且具有更高的准确度和分辨率。

下面以拟南芥为待测样本，来进一步说明本发明的染色体跨度的单体型图的详细构建流程。以下如无详细说明，使用试剂或药品均来自于NEB公司。

一、文库构建及测序步骤：

1、交联：用甲醛等细胞交联剂处理培养到一定数目的细胞，使全基因组范围内的具有物理上相互作用的DNA和其结合蛋白交联而得到固定；

2、限制性内切酶酶切：选择HindIII限制性核酸内切酶处理上述交联的DNA，使交联点两侧产生缺口；

3、生物素标记：用生物素标记的寡核苷酸将连接点两侧酶切缺口补平；

4、邻近DNA片段的连接：采用T4核酸连接酶将邻近的DNA片段连接；

5、解交联、DNA片段的随机打断及生物素标记片段的回收：选取蛋白酶K将蛋白质和DNA解除交联状态，然后用超声波将DNA随机打断，采用链霉亲和素将含有生物素标记的DNA片段回收；

6、小片段文库构建：将回收的DNA片段构建成可在高通量测序仪上测序的小片段文库，并进行文库检测及定量；

7、测序：基于新一代的测序平台Hiseq对经库检合格的DNA文库进行双末端测序，获得可用于后续分析的PEreads数据；

二、信息分析步骤：

1、参考序列处理：

a)索引文件构建：基于Burrows-WheelerAlignment(BWA)软件和SAMtools对参考序列进行索引文件的构建；

b)基因组特征统计：对于基因组的长度、碱基含量、空缺(gap)比例等基本特征进行统计，统计结果见下表1；

表1：

附：N50和N90是评价基因组拼接质量好坏的一个指标，通常该值越大，表示拼接质量越好。其中，N50是指：将测序得到的若干条reads进行拼接，然后把拼接得到的片段从长到短进行排列，然后相加，当恰好加到该基因组总长度的50％时所对应的那条片段的长度即为N50的长度。同样，当恰好加到该基因组总长度的90％时所对应的那条片段的长度即为N90的长度。

c)酶切位点分析：基于Nhel限制性内切酶酶切位点的序列特征分析其在基因组上分布的位置、数目等；

2、质控

a)样本：抽取一部分PEreads与核酸数据库(nt)进行blast比对，检测是否有外源微生物等污染，判断样本是否合格；

b)测序：对PEreads进行Q20、Q30、N含量等分析，过滤低质量数据(高通量测序中,每测一个碱基会给出一个相应的质量值，这个质量值是衡量测序准确度的。碱基的质量值13,错误率为5％；20的错误率为1％；30的错误率为0.1％。Q20与Q30则表示质量值大于等于20或30的碱基所占百分比。比如，一共测了1G的数据量，其中有0.9G的碱基质量值大于或等于20，那么Q20则为90％。)，获得高质量的数据用于后续的分析；

c)文库：抽取一部分PEreads与参考基因组比对，对插入片段在基因组上的跨度进行分析，大于10K的比例在50％以上，则上述Hi-C文库构建合格；

3、比对：如图1所示，基于Burrows-WheelerAlignment(BWA)软件将PEreads比对到参考基因组上，根据染色体位置对比对结果进行排序，并对低质量的比对进行过滤，从而得到PEreads在参考基因组上的位置信息以及染色体的覆盖深度；

4、SNPs检测：基于比对的结果和SAMtools进行一致性序列的构建和SNPs检测，通过覆盖深度、SNPs质量等对结果进行过滤，获得高质量的SNPs。

如图1和图2所示，其中，第一行表示参考基因组序列的位置信息；第二行表示参考基因组序列的碱基信息；第三行表示一致性序列，其中，点表示与参考基因组上的序列的碱基类型一致；字母表示和参考基因组组上是序列的间接类型不一致，即SNP的位点，用简并碱基来表示；以下其余行表示对比上的测序序列(reads)的信息。

图2中，第三中，字母表示的简并碱基及点表示的与参考基因组序列上的碱基类型一致的碱基构成一致性序列，简并碱基构成一个单体型块。根据第三行以下的序列信息可以看出，由于待测样本为二倍体，来源于同一个亲本的简并位置处的碱基类型构成了GGACGCCT；而简并位置处的碱基类型构成TAGGATTC的序列与参考基因组上的序列相同，因而来自于另一亲本。

5、杂合模块片段提取：基于比对结果和SNPs信息，提取至少含有两个杂合SNPs的reads。

6、候选单体型图构建：以杂合SNP为点，两个杂合SNP之间的片段为边构建候选单体型图，运用greedymax-cutheuristic算法和MEC打分函数对候选单体型图中的每个block进行确定及评估，得到多个候选单体型图的block。

根据图2的一致性序列可以看出，哪些测序序列来自于第一亲本，哪些测序序列来自与第二亲本。在候选单体型图构建的步骤中，如图3所示(说明：B，为Block的简写，表示单体型块。碱基，表示基因组上杂合SNP位点；点，表示基因组上非杂合SNP的位点，问号，表示非定型的位点，上一行表示来自同一亲本；下一行表示来自另一个亲本)，首先，利用来源于同一亲本的含有2或2个以上杂合SNP的同一测序序列中的杂合SNP位点构建短的单体型块。其次，基于含有交联位点或者来源于交联问点两侧的测序序列上的SNP将短的单体型块不断地连接成更大的单体型块，最后形成含有单体型块1、2和4的候选单体型图(如图2中，左框中的单体型块与右框中的单体型块基于共同的位点Y，可以将两个单体型块连接成更大的单体型块)。由于单体型块B3上的杂合SNP位点与其他的单体型块之间没有捕获到交联序列，导致单体型块B3没有连接到更大的单体型块上，即在更大的单体型块中属于无定型的位点，因而用问号来表示。

7、结果统计：根据基因组长度、单体型块的物理位置跨度、SNPs数目、定型位点数目、MEC打分等评估每个单体型块的完整性和分辨率，选择最佳的单体型作为染色体跨度的单体型图。

如图4所示，Block，表示单体型块；offset表示第一个变异位点在VCF文件中的行号；len表示单体型块在基因组相应区域内所有的杂合SNP数目；phased：表示定型位点数目；SPAN：表示候选单体型在基因组上的跨度；MECscore表示minimunerrorcorrection，用来进行候选单体型组装的计算模型；fragments：表示支持该单体型块的片段数。其中，各列的具体含义如下：

第一列，表示variantID，即变异位点所在的行号(在cvf文件中)。

第二列，表示haplotype_1，即单体型1的基因型(0代表与参考基因组上的序列想基因型一致，1代表与变异基因型一致)。

第三列，表示haplotype_2，即单体型2的基因型(0代表与参考基因组上的序列想基因型一致，1代表与变异基因型一致)。

第四列，表示chromosome，即染色体编号。

第五列，表示position，即在染色体中的位置。

第六列，表示refallele，即在参考基因组上的序列的碱基类型。

第七列，表示variantallele，即变异碱基的类型。

第八列，表示genotype，即基因型(0/1代表杂合)。

第九列，表示allelecounts，即此位点的reads支持数；genotypelikelihoods，即基于比对上的reads计算三种二倍体基因型(0/0,0/1,1/1)可能性的似然估计值；delta，即GL(0/0)-GL(0/1),GL(1/1)-GL(0/1)中最大的值，代表此位点有可能不是杂合位点的概率，这个值越大代表这个位点不是杂合位点的可能性越大；MECvariant；即移除此位点后block的MEC分数降低的值，这个值越大代表这个位点不是变异位点或者杂合位点的可能性越大。

对图4中的来源于1号染色体的两个候选单体型图进行完整性和分辨率的评估，评估结果如下表2：

表2：

/	完整性	分辨率
			候选单体型图A	221/30427671	2/2＝1
候选单体型图B	63992/30427671	63992/70523＝0.91

由上表2可以看出，候选单体型图2即为该样本的染色体跨度的单体型图。

从以上的描述中可知，本发明上述的实施例实现了如下技术效果：上述实施例通过Hi-C的方法获得在基因组上线性距离分布广泛的PEreads，并通过上述信息分析方法准确提取出同一条染色体内部的双末端reads，根据reads上的杂合位点将单体型模块连接起来形成了具有染色体跨度的单体型图。本发明是基于长片段测序和单体型构建的思想进行改进的，提供了个体染色体跨度的单体型构建的方法。

另外，发明人还利用现有技术中基于Novoalign软件进行测序数据的比对分析的方法对上述拟南芥的测序数据进行了单体型图构建。并对两种构建方法进行了比对，具体比对结果如图5至图7。

从图5可以看出，单体型图的完整性随着测序覆盖深度(简称，测序深度)的升高也会相应的增加，但是测序深度达到一定的程度，完整性就会达到饱和；本发明和现有技术的方法所构建的单体型图的完整性可达到饱和的最低覆盖深度分别为：28X和41X。可见，本发明的构建方法能够在相对较低的覆盖深度下达到与现有技术相同的完整性。

从图6可以看出，准确度随着覆盖深度的升高也会相应的增加，但是测序深度达到一定的程度，准确度就会趋于平稳，随着深度的增加提升不显著；在20X以上的同样的覆盖深度条件下，本发明的方法所构建的单体型图的准确度高于现有技术的准确度，因此在进行单体型构建的分析时，比对检测SNP可以采用将测序reads直接进行BWA比对，无需对测序的reads进行剪切(tream)。

从图7可以看出，在同样覆盖深度下，本发明的构建方法所构建的单体型图的分辨率都高于现有技术所构建的单体型图。而且，当测序覆盖深度达到55x，本发明的完整性可达到约95％、准确度可达到96％的效果，分辨率为33％左右，而现有技术最高的分辨率小于30％。

在本发明的上述构建方法中，针对于基于杂合SNP进行候选单体型构建的部分也可以采用其他的现有方法进行分析。而且，基于上述比较结果可以看出，本发明提供的构建方法具有以下优点：

1.物种适应范围广：可用于一切已经发表参考基因组的物种，取材简便，只需要进行目标个体的细胞培养即可。

2.实验操作误差小：基于Hi-C的实验，可用足够的细胞量，避免扩增带来的偏好性及误差。

3.分析结果准确性高：基于Hi-C的实验方法可得到线性距离分布广泛的双末端reads。利用reads内部的杂合位点构建出短的block，并根据双末端reads上杂合位点之间的相互连接关系将这些block串联到一起，最终得到染色体跨度，并且准确性和分辨率都很高的个体单体型结果。

4.分析操作简便性好：一键化实现个体染色体跨度的单体型重建。

5.全面性：结果统计部分可全面性评估单体型构建的效果。

6.时间短：本发明可在30天内有效的完成实验及分析的工作，得到高质量的单体型结果。在质量和效率上都有了可观的提升。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种染色体跨度的单体型图的构建方法，其特征在于，所述构建方法包括：

步骤S1，对包含待测样本基因组交联位点的测序文库进行测序，得到PEreads；

步骤S2，将所述PEreads分别与参考基因组上的序列进行第一次比对，得到第一比对结果；

步骤S3，根据所述第一比对结果构建所述PEreads的一致性序列，获得高质量的SNP位点；以及

步骤S4，根据获取的所述高质量的SNP位点对所述PEreads中的每一条reads进行筛选，获得至少含有2个杂合SNP位点的reads；利用所述杂合SNP位点构建所述染色体跨度的单体型图。

2.根据权利要求1所述的构建方法，其特征在于，在所述步骤S2之前，所述构建方法还包括判断所述测序文库是否合格的步骤，所述判断测序文库是否合格的步骤包括：

步骤a，抽取部分所述PEreads作为待质控文库；

步骤b，将所述待质控文库与所述参考基因组上的序列进行第二次比对，得到第二比对结果；

步骤c，利用所述第二比对结果计算所述待质控文库的插入片段在所述参考基因组上的跨度；以及

步骤d，若所述跨度大于10kb，且所述跨度大于10kb的插入片段在所述待质控文库中所有插入片段总数的比例为50％以上，则判断所述测序文库合格。

3.根据权利要求1或2所述的构建方法，其特征在于，在所述步骤S1之后，以及在所述步骤S2之前，所述构建方法还包括对PEreads进行质控的步骤；所述质控的步骤包括：

检测所述PEreads中是否存在外源样本污染，和/或

对所述PEreads进行低质量数据过滤。

4.根据权利要求1或2所述的构建方法，其特征在于，在所述步骤S2之前，还包括对所述参考基因组上的序列进行前处理的步骤；所述前处理的步骤包括：

基于BWA和SAMtools对所述参考基因组上的序列进行索引文件的构建；

对所述参考基因组上的序列的长度、碱基含量和空缺的比例进行统计；以及

获取构建所述测序文库时所使用的限制性内切酶在所述参考基因组序列上的分布位置和数目。

5.根据权利要求4所述的构建方法，其特征在于，所述步骤S2包括：

利用BWA软件的mem模块，将所述PEreads通过第一次比对分别比对到建好所述索引文件的所述参考基因组的序列上，得到所述第一比对结果。

6.根据权利要求5所述的构建方法，其特征在于，所述步骤S2在得到所述第一比对结果后，还包括对所述PEreads进行覆盖均匀性检测的步骤，所述覆盖均匀性检测的步骤包括：

检测在不同测序深度下，所述PEreads对所述参考基因组的覆盖程度和覆盖深度，所述覆盖程度是指所述PEreads覆盖所述参考基因组上的序列的长度与所述参考基因组上的序列的总长度的比值；所述覆盖深度是指所述参考基因组相应位置被所述PEreads覆盖到的次数。

7.根据权利要求1所述的构建方法，其特征在于，所述步骤S3包括：

根据所述第一比对结果，利用SAMtools软件构建所述PEpeads的一致性序列，以及

获取所述一致性序列中存在的所述高质量的SNP位点。

8.根据权利要求5至7中任一项所述的构建方法，其特征在于，所述步骤S4包括：

从所述高质量的SNP位点中获取至少含有2个杂合SNP位点的reads；

以所述杂合SNP位点为点，以两个所述杂合SNP之间的片段长度为边，构建单体型块；

基于所述交联位点之间的线性距离小于30Mb的两侧片段上的所述杂合SNP位点，将多个所述单体型块进行连接，得到候选单体型图；

对所述候选单体型图的完整性和分辨率进行评估，并对评估结果依次按照完整性和分辨率的高低进行排序，挑取在完整性最高条件下分辨率最高的所述候选单体型图作为所述染色体跨度的单体型图；

其中，挑选所述交联位点之间的线性距离小于30Mb的两侧片段的步骤包括：

根据染色体的位置，对所述第一比对结果进行排序，得到所述PEreads中的每一条reads在所述参考基因组的序列上的位置信息；

根据所述PEreads中的每一条reads在所述参考基因组的序列上的所述位置信息，将构成所述交联位点两端的reads进行配对，得到所述交联位点在所述参考基因组上的所述线性距离；

根据所述线性距离，挑选得到所述交联位点之间的所述线性距离小于30Mb的所述两侧片段。

9.根据权利要求5至7中任一项所述的构建方法，其特征在于，在所述步骤S1之前，所述构建方法还包括：利用空间构象捕获的方法构建所述包含待测样本基因组交联位点的测序文库；所述利用空间构象捕获的方法构建所述测序文库的步骤包括：

对所述待测样本的DNA进行交联固定；

利用限制性内切酶对所述DNA进行酶切，产生具有酶切缺口的DNA片段；

利用生物素标记的寡核苷酸将所述具有酶切缺口的DNA片段上的所述酶切缺口进行补平；

利用核酸连接酶将所述DNA片段进行连接，得到连接交联DNA片段；

对所述连接交联DNA片段进行解交联，得到解交联DNA；以及

对所述解交联DNA进行片段化文库构建，得到所述测序文库。

10.一种染色体跨度的单体型图，其特征在于，所述单体型图利用权利要求1至9中任一项所述的构建方法构建而成。