CN106055925B - 基于转录组双端测序数据组装基因组序列的方法和装置 - Google Patents

Abstract

本发明提供一种基于转录组双端测序数据组装基因组序列的方法，所述方法包括将转录组双端测序序列比对到基因组上，保留双端测序序列分别仅能比对唯一的不同基因组序列以及基于最多转录组双端测序序列连接证据的基因组序列拼接筛选、形成新的基因组序列等步骤。本发明还提供实现上述方法的装置。利用本发明方法，通过将转录组双端测序序列比对到基因组上，获得基于最多双端测序比对结果的连接证据，从而进行基因组序列的拼接组装，以提升基因组的完整性。转录组双端测序数据既可以是公共数据库中该物种的转录组双端测序序列，也可以是实验产生的转录组双端测序数据。

Description

基于转录组双端测序数据组装基因组序列的方法和装置

技术领域

本发明涉及基因组学、遗传学和生物信息学领域，特别涉及一种基于转录组双端测序数据组装基因组序列的方法和装置。

背景技术

目前物种的全基因组装配主要依赖于鸟枪法策略。在构建多个插入片段长度不一的文库后，先利用插入片段短的文库组装基因组，逐步利用插入片段更长的文库组装基因组，使得基因组长度逐步增长。但是采用鸟枪法策略构建的基因组无法完整覆盖全部的基因。

DNA转录过程以连续性基因组为模板进行转录，形成转录本。如果基因组不完整，可能导致转录本被分割到不同的基因组contig上。目前利用Illumina测序技术进行转录组测序是较为常见的测序方法，主要包括以下步骤：(1)提取样品总RNA，利用oligo-dT逆转录富集有polyA尾巴的转录本；或者利用ribo-zero法，去除rRNA后，利用随机引物逆转录富集除rRNA之外的所有转录本；(2)将逆转录后获得的cDNA打断成特定大小的片段，构建Illumina pair-end文库；(3)采用双端测序策略，分别从一条转录本的两端开始测序，获得双端测序序列。

如果基因组不完整，那么来自同一条转录本的双端序列将被分别比对到两个不同的基因组序列上。利用这些区域及其在转录组双端测序的位置，能够重新将对应的基因组序列串联起来，形成更长的基因组序列。因此，开发利用转录组双端测序数据组装基因组序列的方法及装置具有可行性，该方法及装置的开发将能够提高基因组的完整性。

发明内容

本发明旨在解决目前全基因组序列组装中存在的拼接不完整、基因区域不完整等问题，提供一种基于转录组双端测序数据组装基因组序列的方法和装置。

为对本发明做出清楚说明，对发明中涉及的技术术语进行如下定义：

contig：预先拼接的基因组序列；

转录组双端测序序列包括左端序列和右端序列；

基因组连接：由两个contig按照前后序列排列在一起，排列在前的序列，称为起点序列；排列在后的序列，称为终止序列；

比对区域：是指转录组双端测序序列与基因组序列相似或一致的区域；由于基因组装配尚不完整，因此一条转录本的两端序列可能被分割到两个不同的contig上；

转录组双端测序序列总长度：是指两端测序序列所含的总碱基数；

比对区域的绝对位置：是指比对区域相对于基因组序列的位置；

比对区域之间的间隔：是指前后两个比对区域(i和j)在对应的基因组序列(A和B)的距离差，即

比对区域之间的间隔＝基因组序列A的长度-比对区域i在A序列的位置+比对区域j在B序列的位置；

序列覆盖度：比对区域长度与转录组双端测序序列总长度的比值。

基因组序列拼接：两个contig按照在转录组双端测序序列中的位置，先后排序而成。

为了实现本发明目的，本发明提供的一种基于转录组双端测序数据组装基因组序列的方法，所述方法包括将转录组双端测序序列比对到基因组上，保留双端测序序列分别仅能比对唯一的不同基因组序列以及基于最多转录组双端测序序列连接证据的基因组序列拼接筛选、形成新的基因组序列等步骤。

具体地，所述方法包括如下步骤：

(1)转录组双端测序序列的清洗

用SolexaQA软件中的dynamictrim和lengthsort模块分别对转录组双端测序结果fastq文件进行清洗，去除低质量序列以及短片段序列；

(2)第一轮转录组双端测序序列的比对

将转录组双端测序序列与contig进行第一轮序列比对，获得双端测序序列在所有可比对上的contig上的位置信息，保留双端测序序列的任一端都比对到唯一且不同的contig上的结果；

对于双端测序序列中任一端，如果(i)比对到多个contig上，或(ii)没有比对到contig上，或(iii)左端、右端序列比对到相同的contig上，则去除该双端测序序列，不做后续分析，本步骤仅保留双端测序序列任一端仅比对到唯一contig，且两个contig不同的结果；

(3)第二轮转录组双端测序序列的比对

将步骤(2)获得的转录组双端测序序列与contig进行第二轮序列比对，过滤掉左端序列和右端序列比对到contig上的序列覆盖度至少90％，且双端序列比对到相同contig或者任一端比对到多个contig上的结果；

(4)最可靠基因组序列连接的筛选

将经过上述两轮比对筛选后，获得的转录组双端测序序列及其比对位置作为contig拼接的连接证据；

在所有的基因组连接中，contig可以有三种角色：(i)仅作为起点序列；(ii)仅作为终止序列；(iii)既可以作为起点序列，也可以作为终止序列；有后续contig与之连接的序列为起点序列；之前有contig与之连接的序列为终止序列；

对于所有作为起点序列的contig，为每个contig选择有最多连接证据的contig；

对于所有作为终止序列的contig，为每个contig选择有最多连接证据的contig；

对于上述保留的所有连接，判断1个连接中的两个contig是否为彼此最多连接证据的序列，如果是，则保留该连接；

保留具有最多连接证据的所述起点序列和所述终止序列；将后续连接新contig，而之前没有连接新contig的基因组序列作为起始点，将之前有连接新contig，而后续没有连接新contig的基因组序列作为终结点，将之前既连接新contig，后续又连接新contig的基因组序列作为中间点；

(5)新基因组序列的组装

对于步骤(4)最终保留的contig只能被分配到如下三种数据集合之一：(i)仅作为起点序列的contig；(ii)仅作为终止序列的contig；(iii)既可以作为起点序列，也可以作为终止序列的contig；

从(i)集合中依次挑选contig作为起始点，从(ii)和(iii)集合中选择后续的连接contig，为该contig进一步选择新的连接contig，直至最后在(ii)集合中找到连接contig为止，至此构建一条完整的组装通路，最后形成的组装通路数量等于(i)集合中contig数目；

即，根据步骤(4)最终保留的序列，将每个只能作为起始点的基因组片段，分别将其作为起始点，选择后续的中间点，为这个中间点进一步选择新的中间点，直至找到终结点为止，从而形成一条完整的组装通路；根据上述组装通路，将各contig串联组装成更长的基因组序列。

前述的方法，步骤(1)用SolexaQA软件中的dynamictrim模块过滤低质量转录组双端测序序列，默认保留测序质量p值＜0.05的序列；用SolexaQA软件中的lengthsort模块过滤掉长度小于25个碱基的reads。

前述的方法，步骤(2)进行第一轮序列比对采用的软件为hisat2。

前述的方法，步骤(3)进行第二轮序列比对采用的软件为blat。

本发明中使用的转录组双端测序序列来源于该物种已公开的转录组双端测序序列，或通过实验方法获得的该物种转录组双端测序序列。例如，所述转录组双端测序序列为：①基于oligo-dT逆转录获得的RNA-seq双端测序序列，②基于ribo-zero方法构建的RNA-seq双端测序序列。

具体来说，例如，将一对转录组双端测序序列(假定为左端序列为a，右端序列为b)比对到基因组上。获得双端序列对应的contig(假定为a对应A，b对应B)及其在contig上的绝对位置。按照本发明步骤(2)提供的筛选方法，过滤后的转录组双端序列为唯一且不同的比对，其特征为，转录组双端序列分别比对到唯一contig，且双端序列比对到的contig不同。通过上述筛选后保留下来的转录组双端测序序列，可作为后续contig拼接的连接证据。

由于第一轮转录组双端序列可能错误比对到contig上，或第一轮比对未穷尽所有可比对的contig，因此将保留下来的转录组双端序列与contig再进行比对。如果比对区域的序列覆盖度大于90％且双端序列对应同一contig，或比对区域的序列覆盖度大于90％且任一端转录组序列对应contig，则之前的比对区域拼接认为是不可靠的，对应的转录组双端序列予以去除。对于上述比对区域，一对转录组双端测序a、b将作为contig拼接A-＞B的连接证据。

接下来，根据本发明，每个contig在序列拼接中有三种属性：起点序列，终止序列，既可以作为起点序列，也可以作为终止序列。例如，两个contig的拼接A-＞B中，A为起点序列，而B为终止序列。对于每个contig，作为起点序列，可能有多个contig与之拼接。根据本发明，仅保留有最多连接证据的基因组拼接。例如，对于contig序列A，作为起点序列，可能有多种拼接方式，例如A-＞B、A-＞K，和A-＞M。每种拼接的双端测序序列连接证据为5、3和2，即A-＞B连接有5对双端测序序列支持，A-＞K和A-＞M连接分别有3对和2对双端测序序列支持。则应该保留A-＞B。同理，对于每个contig，作为终止序列，也采取同上的操作步骤。例如，对于基因组序列B，作为终点序列，可能有多种拼接方式，例如A-＞B、F-＞B，和G-＞B。每种拼接方式的蛋白连接证据为5、3和2，应保留A-＞B。

最后，将保留后的基因组拼接串联起，形成新的基因组序列。针对上一步中每个只能作为起始点的基因组片段，分别将其作为起始点，从保留的基因组拼接中，选择后续的中间点；为这个中间点进一步选择新的中间点，直至找到终结点为止。根据上述各基因组序列连接的前后顺序将各基因组片段连接组装成更长的基因组片段。例如，保留下的基因组拼接A-＞B和B-＞D，则串联后形成的基因组顺序为A-＞B-＞D。

本发明还提供实现上述方法的装置，所述装置包括如下单元：

1)转录组双端测序序列清洗单元

用SolexaQA软件中的dynamictrim和lengthsort模块分别对转录组双端测序结果fastq文件进行清洗和配对；

其中，用SolexaQA软件中的dynamictrim模块过滤低质量转录组双端测序序列，默认保留测序质量p值＜0.05的序列；用SolexaQA软件中的lengthsort模块过滤掉长度小于25个碱基的reads。

2)第一轮转录组双端测序序列比对结果的保留单元

用hisat2软件进行第一轮序列比对，获得转录组双端测序序列上的两端序列能比对到contig的数量，及其在不同contig上的绝对位置。

所述保留单元包括：(a1)左端序列比对到唯一contig上的筛选模块；(a2)右端序列比对到唯一contig上的筛选模块；(a3)区分两个不同contig的模块。

满足上述保留单元的双端序列，要保留。

3)第二轮转录组双端测序序列比对结果的筛选单元

用blat软件进行第二轮序列比对，获得序列覆盖度至少90％的，且转录组双端测序序列上的两端序列能比对到contig的数量，及其在不同contig上的绝对位置。

所述筛选单元包括：(b1)左端序列和右端序列比对到contig的序列覆盖度至少90％的筛选模块；(b2)从(b1)筛选得到的序列中，判断比对到相同contig上的筛选模块；(b3)从(b1)筛选得到的序列中，判断任一端序列比对到多个contig上的筛选模块。

满足(b2)或者(b3)的双端序列，要剔出。将通过第二轮转录组双端测序序列比对结果的筛选单元保留的转录组双端测序序列及其比对结果作为contig拼接的连接证据。

4)最可靠基因组序列连接的筛选单元

所述筛选单元包括以下三个模块：

(c1)起点序列最可靠连接contig的筛选模块：对于每条作为起点序列的contig，从3)保留的比对结果中选择与其有最多连接证据的终止序列，并保留对应的连接；

(c2)终止序列最可靠连接contig的筛选模块：对于每条作为终止序列的contig，从(c1)保留的连接中，选择与其有最多连接证据的起点序列，并保留对应的连接；

(c3)双向最可靠连接contig的筛选模块：在(c2)保留的连接中，对于每条起始序列最可靠连接的终止序列，如果该终止序列的最可靠连接恰好也是该起始序列，则保留该连接。

5)新基因组序列的组装单元

根据4)保留的连接，将每个只能作为起始点的基因组片段，分别将其作为起始点，选择后续的中间点，为这个中间点进一步选择新的中间点，直至找到终结点为止，从而形成一条完整的组装通路。

根据上述组装通路，将各contig串联组装成更长的基因组序列。

利用本发明提供的基于转录组双端测序数据组装基因组序列的方法及装置，通过将转录组双端测序序列比对到基因组上，获得基于最多双端测序比对结果的连接证据，从而进行基因组序列的拼接组装，以提升基因组的完整性。转录组双端测序数据既可以是公共数据库中该物种的转录组双端测序序列，也可以是实验产生的转录组双端测序数据。

具体实施方式

以下实施例用于说明本发明，但不用来限制本发明的范围。若未特别指明，实施例中所用的技术手段为本领域技术人员所熟知的常规手段，所用原料均为市售商品。

实施例1 利用人的转录组双端测序数据组装人基因组序列

从美国国立生物技术信息中心(NCBI，http://www.ncbi.nlm.nih.gov/)网站SRA数据库下载人转录组双端测序序列(Accession：ERR420387，共27318482对双端测序)和36437条FASTA格式的人基因组contig序列(N50：148715bp)。

1、清洗转录组双端测序序列

从http://solexaqa.sourceforge.net/网站下载SolexaQA程序，用dynamictrim模块对转录组双端序列进行清洗，去除低质量序列，默认保留测序质量p值＜0.05的序列。然后用lengthsort模块去除长度小于25个碱基的测序序列，最后保留26247926对高质量的转录组双端测序序列。

2、第一轮转录组双端测序序列的比对

(1)从http://www.ccb.jhu.edu/下载hisat2程序，用hisat2-build给参考基因组建立索引。

(2)用hisat2将转录组双端测序序列与基因组contig序列进行比对。获得转录组双端序列比对上的所有contig及其在所述contig上的绝对位置。

(3)筛选出双端序列分别仅比对到一条基因组序列且双端所对比上的基因组序列不同的转录组双端测序序列，得到唯一不同比对的转录组双端测序比对结果。经过该步处理后，有150114对转录组双端序列满足筛选要求，并进行第二轮序列比对。

3、第二轮转录组双端测序序列的比对

(4)从http://hgdownload.cse.ucsc.edu/admin/exe/下载BLAT程序，选择单机版模式，以所述唯一不同比对的转录组双端测序序列作为查询序列，以基因组片段作为匹配序列，参数设置为-noHead。每一条序列的比对覆盖度大于90％。

(5)去除任一端比对到多个contig，或两端比对同一contig的转录组序列。经过该步处理后，有90992对转录组双端序列满足筛选要求，用于后续contig拼接。

4、基因组contig拼接筛选

对步骤3中保留下来的每个contig进行连接，按照本发明提供的方法，分别为其选择连接证据最多的起点序列和终止序列。本步骤结束后产生了4873个可靠的基因组序列拼接关系。

将这些基因组序列分为(i)仅作为起点序列的contig；(ii)仅作为终止序列的contig；(iii)既可以作为起点序列，也可以作为终止序列的contig三类。

5、形成新的基因组序列

针对步骤4中属于(i)类的每个contig，分别将其作为起始点，从(ii)类和(iii)类的contig中，寻找可拼接的contig，形成基因组序列连接；将找到的contig作为新的起始点，进一步如上所述操作，寻找可连接的contig，直至没有可连接的contig为止。根据上述各基因组序列连接的前后顺序拼接组装成更长的基因组序列，从而完成基因组组装过程。本步骤结束后产生3131个新的基因组序列。

结果：组装后的人基因组序列为31564条，较原来减少了13.37％；N50长度为169805bp，增长了14.18％。

虽然，上文中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (6)

1.一种基于转录组双端测序数据组装基因组序列的方法，其特征在于，所述方法包括将转录组双端测序序列比对到基因组上，保留双端测序序列分别仅能比对唯一的不同基因组序列以及基于最多转录组双端测序序列连接证据的基因组序列拼接筛选、形成新的基因组序列的步骤；

所述方法包括如下步骤：

（1）转录组双端测序序列的清洗

用SolexaQA软件中的dynamictrim和lengthsort模块分别对转录组双端测序结果fastq文件进行清洗，去除低质量序列以及短片段序列；

（2）第一轮转录组双端测序序列的比对

将转录组双端测序序列与预先拼接的基因组序列contig进行第一轮序列比对，获得双端测序序列在所有可比对上的contig上的位置信息，保留双端测序序列的任一端都比对到唯一且不同的contig上的结果；

（3）第二轮转录组双端测序序列的比对

将步骤（2）获得的转录组双端测序序列与预先拼接的基因组序列contig进行第二轮序列比对，过滤掉左端序列和右端序列比对到contig上的序列覆盖度至少90%的结果，且过滤掉双端序列比对到相同contig或者任一端比对到多个contig上的结果；

（4）最可靠基因组序列连接的筛选

将经过上述两轮比对筛选后，获得的转录组双端测序序列及其比对位置作为contig拼接的连接证据；

有后续contig与之连接的序列为起点序列；之前有contig与之连接的序列为终止序列；

对于所有作为起点序列的contig，为每个contig选择有最多连接证据的contig，作为终止序列；

对于所有作为终止序列的contig，为每个contig选择有最多连接证据的contig，作为起点序列；

保留具有最多连接证据的所述起点序列和所述终止序列；将后续连接新contig，而之前没有连接新contig的基因组序列作为起始点，将之前有连接新contig，而后续没有连接新contig的基因组序列作为终结点，将之前既连接新contig，后续又连接新contig的基因组序列作为中间点；

（5）新基因组序列的组装

根据步骤（4）最终保留的序列，将每个只能作为起始点的基因组片段，分别将其作为起始点，选择后续的中间点，为这个中间点进一步选择新的中间点，直至找到终结点为止，从而形成一条完整的组装通路；

根据上述组装通路，将各contig串联组装成更长的基因组序列；

所述转录组双端测序序列来源于某物种已公开的转录组双端测序序列，或通过实验方法获得的该物种转录组双端测序序列；所述转录组双端测序序列为：①基于oligo-dT逆转录获得的RNA-seq双端测序序列，②基于ribo-zero方法构建的RNA-seq双端测序序列。

2.如权利要求1所述的方法，其特征在于，步骤（1）用SolexaQA软件中的dynamictrim模块过滤低质量转录组双端测序序列，默认保留测序质量p值< 0.05的序列；用SolexaQA软件中的lengthsort模块过滤掉长度小于25个碱基的reads。

3.如权利要求1所述的方法，其特征在于，步骤（2）进行第一轮序列比对采用的软件为hisat2。

4.如权利要求1-3任一项所述的方法，其特征在于，步骤（3）进行第二轮序列比对采用的软件为blat。

5.一种基于转录组双端测序数据组装基因组序列的装置，其特征在于，所述装置包括如下单元：

1）转录组双端测序序列清洗单元

用SolexaQA软件中的dynamictrim和lengthsort模块分别对转录组双端测序结果fastq文件进行清洗和配对；

2）第一轮转录组双端测序序列比对结果的保留单元

用hisat2软件进行第一轮序列比对，获得转录组双端测序序列上的两端序列能比对到contig的数量，及其在不同contig上的绝对位置；

所述保留单元包括：（a 1）左端序列比对到唯一contig上的筛选模块；（a2）右端序列比对到唯一contig上的筛选模块；（a3）区分两个不同contig的模块；

满足上述保留单元的双端序列，要保留；

3）第二轮转录组双端测序序列比对结果的筛选单元

用blat软件进行第二轮序列比对，获得双端序列的覆盖度都超过90%的，且转录组双端测序序列上的两端序列能比对到contig的数量，及其在不同contig上的绝对位置；

所述第二轮转录组双端测序序列比对结果的筛选单元包括：（b1）左端序列和右端序列比对到contig的序列覆盖度至少90%的筛选模块；（b2）从（b1）筛选得到的序列中，判断比对到相同contig上的筛选模块；（b3）从（b1）筛选得到的序列中，判断任一端序列比对到多个contig上的筛选模块；

满足（b2）或者(b3)的双端序列，要剔除；

将通过第二轮转录组双端测序序列比对结果的筛选单元保留的转录组双端测序序列及其比对结果作为contig拼接的连接证据；

4）最可靠基因组序列连接的筛选单元

所述最可靠基因组序列连接的筛选单元包括以下三个模块：

（c1）起点序列最可靠连接contig的筛选模块：对于每条作为起点序列的contig，从3）保留的比对结果中选择与其有最多连接证据的终止序列，并保留对应的连接；

（c2）终止序列最可靠连接contig的筛选模块：对于每条作为终止序列的contig，从（c1）保留的连接中，选择与其有最多连接证据的起点序列，并保留对应的连接；

（c3）双向最可靠连接contig的筛选模块：在（c2）保留的连接中，对于每条起始序列最可靠连接的终止序列，如果该终止序列的最可靠连接恰好也是该起始序列，则保留该连接；

5）新基因组序列的组装单元

根据4）保留的连接，将每个只能作为起始点的基因组片段，分别将其作为起始点，选择后续的中间点，为这个中间点进一步选择新的中间点，直至找到终结点为止，从而形成一条完整的组装通路；

根据上述组装通路，将各contig串联组装成更长的基因组序列。

6.如权利要求5所述的装置，其特征在于，对于所述转录组双端测序序列清洗单元，用SolexaQA软件中的dynamictrim模块过滤低质量转录组双端测序序列，默认保留测序质量p值< 0.05的序列；用SolexaQA软件中的lengthsort模块过滤掉长度小于25个碱基的reads。