CN102789553B - 利用长转录组测序结果装配基因组的方法及装置 - Google Patents

Abstract

本发明涉及一种利用长转录组测序结果装配基因组的方法及装置，该方法是将同一物种的转录组测序读段与基因组片段进行比对，去除仅比对到1个基因组片段的转录组测序读段，对保留下来的转录组测序读段上的查询区段以规定条件进行筛选，然后按照规定条件获得与保留下来的查询区段相关的区段连接，基于该区段连接对基因组片段进行连接，从而完成对基因组序列的组装。根据本发明的利用长转录组测序结果装配基因组的方法，能够利用长片段测序数据，包括现有公开的大量Sanger数据，从而能够利用长的转录组测序读段来组装基因组序列。

Description

利用长转录组测序结果装配基因组的方法及装置

技术领域

本发明涉及一种利用长转录组测序结果装配基因组的方法及装置。

背景技术

目前，在基因组和转录组研究中，至少采用4种测序技术，其中包括传统的第一代Sanger测序技术以及作为第二代高通量测序技术的罗氏454测序技术、Illumina测序技术和AB公司的SOLiD技术。Sanger测序技术和Roche454测序技术产生的测序读长平均在300bp以上，而Illumina测序技术和SOLiD技术产生的测序读长低于150bp。

在基因组序列组装过程中，需要构建不同长度的基因组测序文库，来连接两个基因组片段。使用基因组片段比较短的测序文库，能够更准确地连接两个基因组片段，使用基因组片段比较长的测序文库，能够使组装后的基因组长度更长。因此，目前，在基因组组装过程中，通常先利用短测序文库组装基因组，再使用长测序文库组装基因组，使得组装后的基因组片段长度越来越大。目前第二代测序技术能够准确地构建短于20kb的基因组测序文库，但难以构建长于20kb的测序文库，所以，必须利用酵母人工染色体和细菌人工染色体，产生更长的测序文库，例如30kb和150kb的测序文库。但利用酵母人工染色体和细菌人工染色体的技术存在产生的数据量少、成本高、周期长的缺陷。为了加快基因组组装进度，人们尝试采用新方法替代上述构建长片段测序文库的方法。但是迄今尚未开发出低廉的长片段测序文库的构建方法。

成熟的转录本是通过以连续性基因组为模板进行转录，然后剪切去除内含子，拼接剩余的外显子而形成。如果基因组片段未组装起来，即基因组不完整，则可能将一条成熟的转录本分成2个或者更多个转录本片段。利用这些转录本片段，能够重新将对应的基因组片段串联起来，形成更长的基因组序列。而两个转录本片段之间的距离是被剪切去除的内含子。据报道内含子的平均长度为90kb，最长可达到200kb，因此利用转录组组装基因组的方法与利用长测序文库组装基因组相似。

第二代测序技术已经广泛应用于转录组测序，例如，当利用Illumina公司的Solexa GAIIx仪器获得1G的测序结果时，能产生千万级的转录组读段，并且仅需1周左右的时间，从而降低成本。因此与基于长测序文库的测序相比，使用转录组数据的优势在于容量高、成本低和周期短。

基于这种思想，研究人员尝试利用转录组来拼接基因组。目前ERANGE软件中的RNA-PATH模块可利用Illumina的转录组双端测序结果来组装基因组，具体组装过程为，首先直接将转录组的双端读段比对到基因组片段中，然后去除仅比对到同一基因组片段中的读段，利用比对到不同基因组片段的双端读段，将两个基因组片段连接在一起。

但是，上述软件存在以下问题：

（1）目前RNA-PATH模块仅能支持illumina测序技术产生的双端数据，而无法利用其他测序技术产生的单向测序数据，特别是无法利用现有公开的大量单向的Sanger数据和Roche454测序数据，导致使用范围受到限制。

（2）RNA-PATH模块根据读段长度分别使用Eland、Bowtie和BLAT三种比对程序将转录组读段比对到预先装配的基因组上。但是，上述三种比对程序可比对的转录组数据有限，导致实际有效的转录组数据少于其他比对程序。

因此，本领域期待一种能够利用长的转录组测序读段组装基因组序列的方法。

发明内容

本发明提供一种能够利用长的转录组测序读段组装基因组序列的方法。

为了对本发明作出清楚的说明，首先针对本说明书中使用的技术术语如下进行定义。

查询序列（query sequence），也称为输入序列，在本说明书中是指转录组测序读段。

匹配序列（target sequence），是与上述查询序列相似的序列，在本说明书中是指与转录组测序读段相似的基因组片段。

比对区域（alignment region），在本说明书中是指查询序列与匹配序列相似或者一致的区域。由于基因组装配尚不完整，所以，一条查询序列可能被分割成多个比对区域，这些比对区域分别位于多个基因组片段上。

查询区段（query block），在本说明书中是指比对区域在转录组测序读段上的序列。

匹配区段（target block），在本说明书中是指比对区域在基因组片段上的序列。

查询序列的长度（qSize），是指查询序列的所有碱基数。在本说明书中，具体为转录组测序读段的所有碱基的总和。

匹配序列的长度（tSize），是指匹配序列的所有碱基数。在本说明书中，具体为基因组片段的所有碱基的总和。

比对区域的相对起始位置（qStart）：是指比对区域相对于整条查询序列的起始位置。在本说明书中，具体为比对区域在转录组测序读段中的起始位置。以下简称为相对起始位置。

比对区域的相对终止位置（qEnd），是指比对区域相对于整条查询序列的终止位置。在本说明书中，具体为比对区域在转录组测序读段中的终止位置。以下简称为相对终止位置。

比对区域的相对长度（query block length），等于相对终止位置与相对起始位置的差，表示比对区域在转录组测序读段中的长度。以下简称为相对长度。

比对区域的绝对起始位置（tStart），是指匹配区域相对于整条匹配序列的起始位置。在本说明书中，具体为比对区域在基因组片段中的起始位置。以下简称为绝对起始位置。

比对区域的绝对终止位置（tEnd），是指比对区域相对于整条匹配序列的终止位置。在本说明书中，具体为比对区域在基因组片段中的终止位置。以下简称为绝对终止位置。

比对区域的绝对长度（target block length），等于绝对终止位置与绝对起始位置的差，表示比对区域在基因组片段中的长度。以下简称为绝对长度。

一致碱基数（matches），是指比对区域中查询区段和匹配区段中碱基相同的数量。

错配碱基数（misMatches），是指比对区域中查询序列与匹配序列不一致的碱基数。

查询序列比对插入的次数（qNumInsert），是指查询序列比对匹配序列的过程中，查询序列中错配碱基的插入次数，需要说明的是，当2个以上碱基连续发生错配时，该2个以上连续错配碱基作为一个整体，视为插入一次。在本说明书中，具体为将转录组测序读段比对到基因组片段上的过程中，转录组测序读段中错配碱基的插入次数。

匹配序列比对插入的次数（tNumInsert），是指查询序列比对匹配序列的过程中，匹配序列中错配碱基的插入次数，需要说明的是，当2个以上碱基连续发生错配时，该2个以上连续错配碱基作为一个整体，视为插入一次。在本说明书中，具体为将转录组测序读段比对到基因组片段上的过程中，基因组片段中错配碱基的插入次数。

区段连接：两个或者多个查询区段在整条查询序列中的先后顺序。

本发明的利用长转录组测序结果装配基因组的方法的具体方案如下。

一种利用长转录组测序结果装配基因组的方法，包括以下步骤：

（1）将同一物种的转录组测序读段与基因组片段进行比对；

（2）去除仅比对到1个基因组片段的转录组测序读段；

（3）针对（2）中保留下来的各转录组测序读段，去除相似程度低的转录组测序读段，优选去除相似程度低于90%的转录组测序读段；

（4）针对（3）中保留下来的各转录组测序读段，分别在每一个转录组测序读段上，按照各查询区段的相对起始位置从小到大依次排列，并根据该排列顺序对各查询区段由小到大进行编号，然后以编号相对小的查询区段（优选编号最小的查询区段）作为参考区段，分别将后续的查询区段（也称为后续查询区段）与所述参考区段进行比较，保留符合下述两个条件中任意一个条件的后续的查询区段，去除下述两个条件都不符合的后续的查询区段，

ⅰ后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且后续的查询区段的相对终止位置与参考区段的相对终止位置的差的绝对值小于10，

ⅱ后续的查询区段的相对终止位置与所述参考区段的相对终止位置之差大于等于10，

将保留下来的后续查询区段（优选该转录组测序读段上保留下来的第一个后续查询区段）作为新参考区段，继续将该新参考区段的后续查询区段（优选该新参考区段的所有后续查询区段）与该新参考区段进行比较，保留符合上述条件之一的查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段，

对保留下来的各查询区段按照上述编号方法重新进行偏号，然后针对每一转录组测序读段，以编号相对小的查询区段（优选编号最小的查询区段）作为参考区段，分别将后续查询区段（优选所有后续查询区段）与所述参考区段进行比较，如果后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且该后续的查询区段的相对终止位置和参考区段的相对终止位置的差的绝对值小于10，则去除该后续查询区段，而保留不满足上述条件的后续查询区段，完成所有比较后，如果存在满足上述条件的后续查询区段，则将该参考区段也去除，

接下来，将保留的后续查询区段（优选该转录组测序读段上保留下来的第一个后续查询区段）作为新参考区段，继续将该新参考区段的后续查询区段（优选该新参考区段的所有后续查询区段）与该新参考区段进行比较，以相同条件判断去除或保留该新参考区段和后续查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段，

需要说明的是，在上述步骤（4）中，当后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且该后续的查询区段的相对终止位置和参考区段的相对终止位置的差的绝对值小于10时，则认为该后续查询区段与参考区段相似，基于该参考区段一直比较至该转录组测序读段上最后一个查询区段，在比较过程中，去除相似的后续查询区段，保留不满足上述条件的后续查询区段，并且，如果存在与参考区段相似的查询区段，则在比较结束时去除该参考区段，然后将保留的后续查询区段作为新参考区段，继续将该新参考区段的后续查询区段与该新参考区段进行比较，以相同条件判断去除还是保留该新参考区段和后续查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段；

（5）针对（4）中所保留下来的所有查询区段按照与（4）相同的编号方法进行编号，然后将每一转录组测序读段中编号相对小的查询区段（优选编号最小的查询区段）作为起点区段，将其后面的所有查询区段分别与其进行比较，如果在后续查询区段中存在相对起始位置与该起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则将该后续的查询区段和所述起点区段的组合作为一个区段连接保留下来，并且计算出这两个区段之间的距离，即该后续查询区段的相对起始位置与所述起点区段的相对终止位置的差值，如果在后续的查询区段中不存在相对起始位置与起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则去除所述的起点区段，然后以该起点区段后续的查询区段（优选最靠近该起点区段的后续查询区段）作为新的起点区段，以同样条件进行与上述相同的比较，以确定与所述新的起点区段相关的区段连接，并且如上所述地计算出这两个区段之间的距离，如果不存在与所述新的起点区段相关的区段连接，则去除所述新的起点区段，如此反复地进行比较，直至该转录组测序读段上最后一个查询区段作为起点区段，

接下来，针对每个查询区段，根据计算的两个查询区段之间的距离，保留与该查询区段相关且两个区段之间的距离最小的区段连接，去除其余与该查询区段相关的区段连接；

（6）将（5）中保留下来的每一个区段连接作为与其对应的两个基因组片段连接的支持证据；

（7）将（6）中所保留下来的每个基因组片段分别作为起始基因组片段，并在与其连接的所有基因组片段中选择支持证据最多的基因组片段作为终止片段，形成一个基因组片段连接关系，

针对（6）中所保留下来的每个基因组片段都进行上述选择终止基因组片段并形成一个基因组片段连接关系的操作，需要说明的是，这里只是形成一种连接关系，而不是将两个基因组片段实际连接在一起，对于各基因组片段而言，存在三种情况：（i）只能连接在其他基因组片段之前，作为起始基因组片段，（ii）只能连接在其他基因组片段之后，作为终止基因组片段，（iii）既能够连接在其他基因组片之前作为起始基因组片段，又能够连接在其他基因组片段之后作为终止基因组片段，

另外，在步骤（7）中，如果起始基因组片段有两个或两个以上支持证据数量相同且最多的连接基因组片段，则去除该起始基因组片段

（8）针对（7）中每个只能作为起始基因组片段的基因组片段（即，步骤（7）中的（i）中描述的基因组片段），分别将其作为起始点，从只能连接在其他基因组片段之后作为终止基因组片段的基因组片段（即，步骤（7）中的（ii）中描述的基因组片段），以及既能够连接在其他基因组片之前作为起始基因组片段，又能够连接在其他基因组片段之后作为终止基因组片段的基因组片段（即，步骤（7）中的（iii）中描述的基因组片段）中，寻找可连接的基因组片段，形成基因组片段连接，将该基因组片段连接作为新的起始点，进一步如上所述那样寻找可连接的基因组片段，直至没有可连接的基因组片段为止，根据上述各基因组片段连接的前后顺序将各基因组片段连接组装成更长的基因组片段。

优选在上述步骤（2）和步骤（3）之间或者在步骤（3）和步骤（4）之间，还包括去除序列覆盖度高于90%的转录组测序读段的步骤。

在本发明的利用长转录组测序结果装配基因组的方法中，所述转录组序列读段可以由Sanger测序技术或Roche 454测序技术获得，并且可以采用序列比对程序Blat的单机版模式，以程序默认参数进行比对。并且，所述转录组序列读段与所述基因组片段优选来源于斑马鱼或罗非鱼。

本发明还提供一种组装基因组序列的装置，包括以下模块：

比对模块，将同一物种的转录组测序读段与基因组片段进行比对；

去除模块1，去除仅比对到1个基因组片段的转录组测序读段；

去除模块2，针对去除模块1中保留下来的各转录组测序读段，去除相似程度低的转录组测序读段；

查询区段选择模块，针对去除模块2中保留下来的各转录组测序读段，分别在每一个转录组测序读段上，按照各查询区段的相对起始位置从小到大依次排列，并根据该排列顺序对各查询区段由小到大进行编号，然后以编号相对小的查询区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较，保留符合下述两个条件中任意一个条件的后续的查询区段，去除下述两个条件都不符合的后续查询区段，

ⅰ后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且后续的查询区段的相对终止位置与参考区段的相对终止位置的差的绝对值小于10，

ⅱ后续的查询区段的相对终止位置与所述参考区段的相对终止位置之差大于等于10，

将保留下来的后续查询区段作为新参考区段，继续将该新参考区段后续的查询区段与该新参考区段进行比较，保留符合上述条件之一的后续查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段，

对保留下来的各查询区段按照上述编号方法重新进行偏号，然后针对每一转录组测序读段，以编号相对小的查询区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较，如果后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且该后续的查询区段的相对终止位置和参考区段的相对终止位置的差的绝对值小于10，则去除该后续的查询区段，而保留不满足上述条件的后续查询区段，完成所有比较后，如果存在满足上述条件的后续的查询区段，则将该参考区段也去除，

接下来，将保留的后续查询区段作为新参考区段，继续将该新参考区段后续的查询区段与该新参考区段进行比较，以相同条件判断去除或保留该新参考区段和后续的查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段；

区段连接选择模块，针对查询区段选择模块中所保留下来的所有查询区段按照与查询区段选择模块中相同的编号方法进行编号，然后将每一转录组测序读段中编号相对小的查询区段作为起点区段，将其后续的所有查询区段分别与其进行比较，如果在后续的查询区段中存在相对起始位置与该起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则将该后续的查询区段和所述起点区段的组合作为一个区段连接保留下来，并且计算出这两个区段之间的距离，即该后续的查询区段的相对起始位置与所述起点区段的相对终止位置的差值，如果在后续的查询区段中不存在相对起始位置与起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则去除所述起点区段，然后以该起点区段后续的区段作为新的起点区段，以同样条件进行与上述相同的比较，以确定与所述新的起点区段相关的区段连接，并且如上所述地计算出这两个区段之间的距离，如果不存在与所述新的起点区段相关的区段连接，则去除所述新的起点区段，如此反复地进行比较，直至该转录组测序读段上最后一个查询区段作为起点区段，

接下来，针对每个查询区段，根据计算的两个查询区段之间的距离，保留与该查询区段相关且两个区段之间的距离最小的区段连接，去除其余与该查询区段相关的区段连接；

支持证据收集模块，将区段连接选择模块中保留下来的每一个区段连接作为与其对应的两个基因组片段连接的支持证据；

基因组片段连接关系选择模块，将支持证据收集模块中所保留下来的每个基因组片段分别作为起始基因组片段，并在与其连接的所有基因组片段中选择支持证据最多的基因组片段作为终止片段，形成一个基因组片段连接关系；

基因组片段连接确定模块，针对基因组片段连接关系选择模块中每个只能作为起始基因组片段的基因组片段，分别将其作为起始点，从只能连接在其他基因组片段之后作为终止基因组片段的基因组片段，以及既能够连接在其他基因组片之前作为起始基因组片段，又能够连接在其他基因组片段之后作为终止基因组片段的基因组片段中，寻找可连接的基因组片段，形成基因组片段连接，将该基因组片段连接作为新的起始点，进一步如上所述那样寻找可连接的基因组片段，直至没有可连接的基因组片段为止，根据上述各基因组片段连接的前后顺序将各基因组片段连接组装成更长的基因组片段。

优选在所述去除模块2中，去除相似程度低于90%的转录组测序读段。

另外，优选在所述去除模块1和所述去除模块2之间或者在所述去除模块2和所述查询区段选择模块之间，还包括去除模块3，用于去除序列覆盖度高于90%的转录组测序读段。

在本发明的组装基因组序列的装置中，优选所述查询区段选择模块针对去除模块2或去除模块3中保留下来的各转录组测序读段，分别在每一个转录组测序读段上，按照各查询区段的相对起始位置从小到大依次排列，对于相对起始位置相同的查询区段，进一步按照该查询区段的相对终止位置从大到小依次排列，并根据该排列顺序对各查询区段由小到大进行编号，然后以编号最小的查询区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较，保留符合下述两个条件中任意一个条件的后续的查询区段，去除下述两个条件都不符合的后续查询区段，

ⅰ后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且后续的查询区段的相对终止位置与参考区段的相对终止位置的差的绝对值小于10，

ⅱ后续的查询区段的相对终止位置与所述参考区段的相对终止位置之差大于等于10，

将保留下来的第一个后续查询区段作为新参考区段，继续将该新参考区段后续的所有查询区段与该新参考区段进行比较，保留符合上述条件之一的查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段，

对保留下来的各查询区段按照上述编号方法重新进行偏号，然后针对每一转录组测序读段，以编号最小的查询区段作为参考区段，分别将后续的所有查询区段与所述参考区段进行比较，如果后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且该后续的查询区段的相对终止位置和参考区段的相对终止位置的差的绝对值小于10，则去除该后续的查询区段，而保留不满足上述条件的后续查询区段，完成所有比较后，如果存在满足上述条件的后续查询区段，则将该参考区段也去除，

接下来，将该转录组测序读段上保留的第一个后续查询区段作为新参考区段，继续将该新参考区段后续的所有查询区段与该新参考区段进行比较，以相同条件判断去除或保留该新参考区段和后续的查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段。

在本发明的组装基因组序列的装置，还优选所述区段连接选择模块针对查询区段选择模块中所保留下来的所有查询区段按照与查询区段选择模块中相同的编号方法进行编号，然后将每一转录组测序读段中编号最小的查询区段作为起点区段，将其后续的所有查询区段分别与其进行比较，如果在后续的查询区段中存在相对起始位置与该起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则将该后续的查询区段和所述起点区段的组合作为一个区段连接保留下来，并且计算出这两个区段之间的距离，如果在后续的查询区段中不存在相对起始位置与起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则去除所述起点区段，然后以该起点区段的第一个后续区段作为新的起点区段，以同样条件进行与上述相同的比较，以确定与所述新的起点区段相关的区段连接，并且如上所述地计算出这两个区段之间的距离，如果不存在与所述新的起点区段相关的区段连接，则去除所述新的起点区段，如此反复地进行比较，直至该转录组测序读段上最后一个查询区段作为起点区段，

接下来，针对每个查询区段，根据计算的两个查询区段之间的距离，保留与该查询区段相关且两个区段之间的距离最小的区段连接，去除其余与该查询区段相关的区段连接。

在所述基因组片段连接关系选择模块中，如果起始基因组片段有两个或两个以上支持证据数量相同且最多的连接基因组片段，则去除该起始基因组片段。

在上述比对模块中，优选所述转录组序列读段由Sanger测序技术或Roche454测序技术获得，采用序列比对程序Blat，并选择单机版模式，以程序默认参数进行比对。

在使用本发明的装置组装基因组序列时，所述转录组序列读段与所述基因组片段优选来源于斑马鱼或罗非鱼。

根据本发明的利用长转录组测序结果装配基因组的方法，能够利用长片段测序数据，包括Sanger测序技术和Roche454测序技术产生的转录组测序读段，从而能够利用读长较长的转录组测序读段来组装基因组片段，进而获得长度更长的基因组片段。

具体实施方式

下面通过实施例具体说明本发明的利用长转录组测序结果装配基因组的方法，但本发明并不限于所列举的实施例。

实施例1

利用斑马鱼转录组Sanger测序读段来组装斑马鱼的基因组序列

材料：从美国生物技术信息中心（NCBI，National Center for BiotechonlogyInformation）网站（http://www.ncbi.nlm.nih.gov/dbEST/index.html）下载1546467条FASTA格式的斑马鱼转录组Sanger测序读段。从英国桑格研究院（Sangerinstitute,http://www.sanger.ac.uk/Projects/D_rerio/wgs.shtml）下载37298条预先装配的斑马鱼基因组片段，这些基因组片段的平均长度为143274bp。

步骤01：从美国加州大学圣克鲁斯分校（http://hgdownload.cse.ucsc.edu/admin/exe/）下载BLAT（BlAST-like alignment tool）程序，选择单机版模式，以各转录组测序读段作为查询序列，以基因组片段作为匹配序列，按照该比对程序的默认参数，将1546467条Sanger测序读段与37298条基因组片段进行比对。结果显示有1420506条转录组Sanger测序读段比对到基因组片段上。

步骤02：根据如下公式分别计算比对区域的序列覆盖度（length coverage）和相似程度（identity）。

如下计算序列覆盖度：

$\mathrm{lengthcoverage}=\frac{\mathrm{qEnd}-\mathrm{qStart}}{\mathrm{qSize}}*100\%$

如下计算相似程度：

（1）当（（qEnd-qStart）－（tEnd-tStart））<0时：

$\mathrm{identity}=\frac{100-\frac{1000*(\mathrm{misMatch}+\mathrm{qNumInsert}+3*\log \left(1\right))}{\mathrm{matches}+\mathrm{misMatches}}*0.1}{100}$

（2）当（（qEnd-qStart）－（tEnd-tStart））≥0时：

$\mathrm{identity}=$

$\frac{100-\frac{1000*(\mathrm{misMatches}+\mathrm{qNumInsert}+3*\log (1+((\mathrm{qEnd}-\mathrm{qStart})-(\mathrm{tEnd}-\mathrm{tStart}))))}{\mathrm{matches}+\mathrm{misMatches}}*0.1}{100}$

以相似程度90%对步骤01中保留下来的1420506条转录组测序读段进行筛选。保留相似程度大于等于90%的转录组测序读段。经过该步骤，保留了1419653条转录组测序读段。对该保留的转录组测序读段进行筛选，保留序列覆盖度（length coverage）小于90%的转录组测序读段。

步骤03：对于步骤02中保留下来的转录组测序读段，去除仅比对到1个基因组片段的转录组测序读段，保留比对到2个以上基因组片段的转录组测序读段。经步骤03后，有397079条转录组测序读段被保留下来，对应有27933071个查询区段。

步骤04：对于步骤03中保留下来的查询区段，在同一转录组测序读段上，按照各查询区段的相对起始位置从小到大依次排列并对其依次编号为第一区段、第二区段……，对于相对起始位置相同的查询区段，进一步按照该查询区段的相对终止位置从大到小依次排列。例如，在同一转录组测序读段上包括3个查询区段a、b、c，查询区段a在转录组测序读段上的相对起始位置小于查询区段b、c，查询区段b、c的相对起始位置相同，但查询区段b的末端碱基的位点数大于查询区段c，此时查询区段a、b、c在该转录组测序读段上依次编号为第一区段、第二区段、第三区段。

然后以编号为第一区段的查询区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较，保留符合下述两个条件中任意一个条件的后续的查询区段，去除下述两个条件都不符合的后续查询区段。

ⅰ后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10,并且后续的查询区段的相对终止位置和参考区段的相对终止位置的差的绝对值小于10。

ⅱ后续的查询区段的相对终止位置与所述参考区段的相对终止位置之差大于等于10。

将保留下来的第一个后续查询区段作为新参考区段，继续将该新参考区段后续的查询区段与该新参考区段进行比较，保留符合上述条件之一的查询区段。反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段。

针对步骤03中保留下来的所有转录组测序读段分别进行上述比较。该比较结束后剩余376875个查询区段。

接下来，对于保留下来的376875个查询区段，在同一转录组测序读段中按照与上述相同的编号方法对查询区段重新进行编号。然后以第一区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较。如果后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且该后续的查询区段的相对终止位置和参考区段的相对终止位置的差的绝对值小于10，则认为该后续查询区段与参考区段相似，基于该参考区段一直比较至该转录组测序读段上最后一个查询区段，在比较过程中，去除相似的后续查询区段，保留不满足上述条件的后续查询区段，并且，如果存在与参考区段相似的查询区段，则在比较结束时去除该参考区段，然后将保留下来的第一个后续查询区段作为新参考区段，继续将该新参考区段的后续查询区段与该新参考区段进行比较，以相同条件判断去除还是保留该新参考区段和后续查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段；

针对其他转录组测序读段上的查询区段也相同地进行编号，并分别进行上述比较。最终剩余293660个查询区段。

步骤05：针对步骤04中所保留下来的所有查询区段按照与步骤04相同的编号方法进行编号。然后将同一转录组测序读段中的第一区段作为起点区段，将其后面的所有查询区段分别与其进行比较，如果在后续的查询区段中存在相对起始位置与该起点区段的相对终止位置之差的绝对值小于30、并且后续查询区段的绝对起始位置与起点区段的绝对终止位置之间的长度小于等于200kb，则将后续所述起点区段和该后续的查询区段的组合作为一个区段连接保留下来，并且计算出这两个区段之间的距离（即该后续的查询区段的相对起始位置与所述起点区段的相对终止位置的差值）。如果在后续的查询区段中不存在相对起始位置与起点区段的相对终止位置之差的绝对值小于30、并且绝对起始位置与起点区段的绝对终止位置之间的长度小于等于200kb的查询区段，则去除所述的起点区段。然后以该转录组测序读段中的第二区段作为新的起点区段，进行与上述第一区段相同的比较，以确定与第二区段相关的区段连接，并且如上所述地计算出这两个区段之间的距离，同样，如果不存在与第二区段相关的区段连接，则去除所述第二区段。如此反复进行与上述第一区段相同的比较，直至该转录组测序读段上最后一个查询区段作为起点区段。

此时，针对每个查询区段，根据之前所计算的两个查询区段之间的距离，保留与该查询区段相关且两个区段之间的距离最小的区段连接，去除其余与该查询区段相关的区段连接。例如对于某一转录组测序读段中的第n查询区段（n≥1），存在多个与该第n查询区段相关的区段连接，在这些区段连接中该第n查询区段与其后续的查询区段之间的距离分别为10、12、18……，此时保留距离为10的区段连接，而去除其他与该第n查询区段相关的区段连接。如果在与该第一查询区段相关的多个区段连接中存在两个以上区段之间距离相同且最小的区段连接，那么这些区段连接都不保留。针对所有查询区段进行相同处理，即保留与该查询区段相关的两个区段间距离最小的区段连接，而去除与该查询区段相关的其余区段连接。

步骤06：将经步骤05所保留下来的每一区段连接作为对应的基因组片段的连接证据。

步骤07：将步骤06中保留下来的每个基因组片段分别作为起始基因组片段，在与其连接的所有基因组片段中选择支持证据最多的基因组片段作为终止片段，形成一个基因组片段连接关系。这些基因组片段可以分为三类：（i）只能连接在其他基因组片段之前，作为起始基因组片段，（ii）只能连接在其他基因组片段之后，作为终止基因组片段，（iii）既能够连接在其他基因组片之前作为起始基因组片段，又能够连接在其他基因组片段之后作为终止基因组片段。如果起始基因组片段有两个或两个以上支持证据数量相同且最多的连接基因组片段，则去除该起始基因组片段。本步骤结束后产生了4419个基因组片段连接关系。

步骤08：针对步骤07中属于（i）类的每个基因组片段（即，只能作为起始基因组片段的基因组片段），分别将其作为起始点，从属于（ii）类和（iii）类的基因组片段（即，只能连接在其他基因组片段之后作为终止基因组片段的基因组片段，以及既能够连接在其他基因组片之前作为起始基因组片段，又能够连接在其他基因组片段之后作为终止基因组片段的基因组片段）中，寻找可连接的基因组片段，形成基因组片段连接，将该基因组片段连接作为新的起始点，进一步如上所述那样寻找可连接的基因组片段，直至没有可连接的基因组片段为止，根据上述各基因组片段连接的前后顺序将各基因组片段连接组装成更长的基因组片段，从而完成基因组组装过程。本步骤结束后产生3473个基因组片段。

结果：组装后的斑马鱼基因组序列为33625条，较原来减少了9.85%；平均长度为165154bp，增长了15.27%。

实例2

利用罗非鱼转录组454测序读段来组装罗非鱼的基因组片段

材料：从美国Broad机构（http://bouillabase.org/）下载5900条罗非鱼全基因组序列片段，其平均长度为2.8M。从美国生物技术信息中心网站下载fastq格式的罗非鱼454测序读段（http://www.ncbi.nlm.nih.gov/sra/SRX078333和http://www.ncbi.nlm.nih.gov/sra/SRX078329）。

方法：首先利用Solexa QA软件包（solexaqa.sourceforge.net），以默认参数过滤低质量的转录组测序读段和短的转录组测序读段。然后使用srtoolbox软件包（http://brianknaus.com/software/srtoolbox/）中的fastq2fasta.pl将fastq格式转换为fasta格式。

接下来，按照实施例1中的步骤01~步骤07进行罗非鱼基因组片段的组装。

结果：预先装配前的罗非鱼基因组片段的平均长度为2.8M，数量为5900条。使用转录组454测序读段组装后，罗非鱼基因组片段的平均长度为3.1M，组装长度增加了11%，数量减少至5597条。

上述实施例中是以斑马鱼和罗非鱼为例进行了说明，当然，本发明的利用长转录组测序结果装配基因组的方法还可以用于其他生物的基因组序列的组装。另外，上述实施例中采用的所述转录组序列读段由Sanger测序技术或Roche454测序技术获得，也可以采用其他测序技术获得的读长较长的转录组测序读段，还可以采用预先装配Illumina短测序读段而得到的更长读段，来组装基因组片段。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的构思和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (10)

1.一种利用长转录组测序结果装配基因组的方法，其特征在于，包括以下步骤：

(1)将同一物种的转录组测序读段与基因组片段进行比对；

(2)去除仅比对到1个基因组片段的转录组测序读段；

(3)针对(2)中保留下来的各转录组测序读段，去除相似程度低的转录组测序读段；

(4)针对(3)中保留下来的各转录组测序读段，分别在每一个转录组测序读段上，按照各查询区段的相对起始位置从小到大依次排列，并根据该排列顺序对各查询区段由小到大进行编号，对于相对起始位置相同的查询区段，进一步按照该查询区段的相对终止位置从大到小依次排列，然后以编号相对小的查询区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较，保留符合下述两个条件中任意一个条件的后续的查询区段，去除下述两个条件都不符合的后续查询区段，

ⅰ后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且后续的查询区段的相对终止位置与参考区段的相对终止位置的差的绝对值小于10，

ⅱ后续的查询区段的相对终止位置与所述参考区段的相对终止位置之差大于等于10，

将保留下来的后续查询区段作为新参考区段，继续将该新参考区段后续的查询区段与该新参考区段进行比较，保留符合上述条件之一的后续查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段，

对保留下来的各查询区段按照上述编号方法重新进行偏号，然后针对每一转录组测序读段，以编号相对小的查询区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较，如果后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且该后续的查询区段的相对终止位置和参考区段的相对终止位置的差的绝对值小于10，则去除该后续的查询区段，而保留不满足上述条件的后续查询区段，完成所有比较后，如果存在满足上述条件的后续的查询区段，则将该参考区段也去除，

接下来，将保留的后续查询区段作为新参考区段，继续将该新参考区段后续的查询区段与该新参考区段进行比较，以相同条件判断去除或保留该新参考区段和后续的查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段；

(5)针对(4)中所保留下来的所有查询区段按照与(4)相同的编号方法进行编号，然后将每一转录组测序读段中编号相对小的查询区段作为起点区段，将其后续的所有查询区段分别与其进行比较，如果在后续的查询区段中存在相对起始位置与该起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则将该后续的查询区段和所述起点区段的组合作为一个区段连接保留下来，并且计算出这两个区段之间的距离，即该后续的查询区段的相对起始位置与所述起点区段的相对终止位置的差值，如果在后续的查询区段中不存在相对起始位置与起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则去除所述起点区段，然后以该起点区段后续的区段作为新的起点区段，以同样条件进行与上述相同的比较，以确定与所述新的起点区段相关的区段连接，并且如上所述地计算出这两个区段之间的距离，如果不存在与所述新的起点区段相关的区段连接，则去除所述新的起点区段，如此反复地进行比较，直至该转录组测序读段上最后一个查询区段作为起点区段，

接下来，针对每个查询区段，根据计算的两个查询区段之间的距离，保留与该查询区段相关且两个区段之间的距离最小的区段连接，去除其余与该查询区段相关的区段连接；

(6)将(5)中保留下来的每一个区段连接作为与其对应的两个基因组片段连接的支持证据；

(7)将(6)中所保留下来的每个基因组片段分别作为起始基因组片段，并在与其连接的所有基因组片段中选择支持证据最多的基因组片段作为终止片段，形成一个基因组片段连接关系；

(8)针对(7)中每个只能作为起始基因组片段的基因组片段，分别将其作为起始点，从只能连接在其他基因组片段之后作为终止基因组片段的基因组片段，以及既能够连接在其他基因组片之前作为起始基因组片段，又能够连接在其他基因组片段之后作为终止基因组片段的基因组片段中，寻找可连接的基因组片段，形成基因组片段连接，将该基因组片段连接作为新的起始点，进一步如上所述那样寻找可连接的基因组片段，直至没有可连接的基因组片段为止，根据上述各基因组片段连接的前后顺序将各基因组片段连接组装成更长的基因组片段。

2.根据权利要求1所述的方法，其特征在于，在步骤(3)中，去除相似程度低于90％的转录组测序读段。

3.根据权利要求2所述的方法，其特征在于，在步骤(2)和步骤(3)之间或者在步骤(3)和步骤(4)之间，还包括去除序列覆盖度高于90％的转录组测序读段的步骤。

4.根据权利要求1至3中任一项所述的方法，其特征在于，如下所述进行步骤(4)，针对(3)中保留下来的各转录组测序读段，或者针对去除序列覆盖度高于90％的转录组测序读段后保留下来的各转录组测序读段，分别在每一个转录组测序读段上，按照各查询区段的相对起始位置从小到大依次排列，对于相对起始位置相同的查询区段，进一步按照该查询区段的相对终止位置从大到小依次排列，并根据该排列顺序对各查询区段由小到大进行编号，然后以编号最小的查询区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较，保留符合下述两个条件中任意一个条件的后续的查询区段，去除下述两个条件都不符合的后续查询区段，

ⅰ后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且后续的查询区段的相对终止位置与参考区段的相对终止位置的差的绝对值小于10，

ⅱ后续的查询区段的相对终止位置与所述参考区段的相对终止位置之差大于等于10，

将保留下来的第一个后续查询区段作为新参考区段，继续将该新参考区段后续的所有查询区段与该新参考区段进行比较，保留符合上述条件之一的查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段，

对保留下来的各查询区段按照上述编号方法重新进行偏号，然后针对每一转录组测序读段，以编号最小的查询区段作为参考区段，分别将后续的所有查询区段与所述参考区段进行比较，如果后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且该后续的查询区段的相对终止位置和参考区段的相对终止位置的差的绝对值小于10，则去除该后续的查询区段，而保留不满足上述条件的后续查询区段，完成所有比较后，如果存在满足上述条件的后续查询区段，则将该参考区段也去除，

接下来，将该转录组测序读段上保留的第一个后续查询区段作为新参考区段，继续将该新参考区段后续的所有查询区段与该新参考区段进行比较，以相同条件判断去除或保留该新参考区段和后续的查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段。

5.根据权利要求1至3中任一项所述的方法，其特征在于，如下所述进行步骤(5)，针对(4)中所保留下来的所有查询区段按照与(4)相同的编号方法进行编号，然后将每一转录组测序读段中编号最小的查询区段作为起点区段，将其后续的所有查询区段分别与其进行比较，如果在后续的查询区段中存在相对起始位置与该起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则将该后续的查询区段和所述起点区段的组合作为一个区段连接保留下来，并且计算出这两个区段之间的距离，如果在后续的查询区段中不存在相对起始位置与起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则去除所述起点区段，然后以该起点区段的第一个后续区段作为新的起点区段，以同样条件进行与上述相同的比较，以确定与所述新的起点区段相关的区段连接，并且如上所述地计算出这两个区段之间的距离，如果不存在与所述新的起点区段相关的区段连接，则去除所述新的起点区段，如此反复地进行比较，直至该转录组测序读段上最后一个查询区段作为起点区段，

接下来，针对每个查询区段，根据计算的两个查询区段之间的距离，保留与该查询区段相关且两个区段之间的距离最小的区段连接，去除其余与该查询区段相关的区段连接。

6.根据权利要求1至3中任一项所述的方法，其特征在于，在步骤(7)中，如果起始基因组片段有两个或两个以上支持证据数量相同且最多的连接基因组片段，则去除该起始基因组片段。

7.根据权利要求1至3中任一项所述的方法，其特征在于，在步骤(1)中，所述转录组序列读段由Sanger测序技术或Roche 454测序技术获得，采用序列比对程序Blat，并选择单机版模式，以程序默认参数进行比对。

8.根据权利要求1至3中任一项所述的方法，其特征在于，所述转录组序列读段与所述基因组片段来源于斑马鱼或罗非鱼。

9.一种利用长转录组测序结果装配基因组的装置，其特征在于，包括以下模块：

比对模块，将同一物种的转录组测序读段与基因组片段进行比对；

去除模块1，去除仅比对到1个基因组片段的转录组测序读段；

去除模块2，针对去除模块1中保留下来的各转录组测序读段，去除相似程度低的转录组测序读段；

查询区段选择模块，针对去除模块2中保留下来的各转录组测序读段，分别在每一个转录组测序读段上，按照各查询区段的相对起始位置从小到大依次排列，并根据该排列顺序对各查询区段由小到大进行编号，对于相对起始位置相同的查询区段，进一步按照该查询区段的相对终止位置从大到小依次排列，然后以编号相对小的查询区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较，保留符合下述两个条件中任意一个条件的后续的查询区段，去除下述两个条件都不符合的后续查询区段，

ⅰ后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且后续的查询区段的相对终止位置与参考区段的相对终止位置的差的绝对值小于10，

ⅱ后续的查询区段的相对终止位置与所述参考区段的相对终止位置之差大于等于10，

将保留下来的后续查询区段作为新参考区段，继续将该新参考区段后续的查询区段与该新参考区段进行比较，保留符合上述条件之一的后续查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段，

对保留下来的各查询区段按照上述编号方法重新进行偏号，然后针对每一转录组测序读段，以编号相对小的查询区段作为参考区段，分别将后续的查询区段与所述参考区段进行比较，如果后续的查询区段的相对起始位置与参考区段的相对起始位置之差小于等于10，并且该后续的查询区段的相对终止位置和参考区段的相对终止位置的差的绝对值小于10，则去除该后续的查询区段，而保留不满足上述条件的后续查询区段，完成所有比较后，如果存在满足上述条件的后续的查询区段，则将该参考区段也去除，

接下来，将保留的后续查询区段作为新参考区段，继续将该新参考区段后续的查询区段与该新参考区段进行比较，以相同条件判断去除或保留该新参考区段和后续的查询区段，反复进行上述比较，直至该转录组测序读段上最后一个查询区段作为参考区段；

区段连接选择模块，针对查询区段选择模块中所保留下来的所有查询区段按照与查询区段选择模块中相同的编号方法进行编号，然后将每一转录组测序读段中编号相对小的查询区段作为起点区段，将其后续的所有查询区段分别与其进行比较，如果在后续的查询区段中存在相对起始位置与该起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则将该后续的查询区段和所述起点区段的组合作为一个区段连接保留下来，并且计算出这两个区段之间的距离，即该后续的查询区段的相对起始位置与所述起点区段的相对终止位置的差值，如果在后续的查询区段中不存在相对起始位置与起点区段的相对终止位置之差的绝对值小于30、且绝对起始位置与该起点区段的绝对终止位置之差小于200Kb的查询区段，则去除所述起点区段，然后以该起点区段后续的区段作为新的起点区段，以同样条件进行与上述相同的比较，以确定与所述新的起点区段相关的区段连接，并且如上所述地计算出这两个区段之间的距离，如果不存在与所述新的起点区段相关的区段连接，则去除所述新的起点区段，如此反复地进行比较，直至该转录组测序读段上最后一个查询区段作为起点区段，

接下来，针对每个查询区段，根据计算的两个查询区段之间的距离，保留与该查询区段相关且两个区段之间的距离最小的区段连接，去除其余与该查询区段相关的区段连接；

支持证据收集模块，将区段连接选择模块中保留下来的每一个区段连接作为与其对应的两个基因组片段连接的支持证据；

基因组片段连接关系选择模块，将支持证据收集模块中所保留下来的每个基因组片段分别作为起始基因组片段，并在与其连接的所有基因组片段中选择支持证据最多的基因组片段作为终止片段，形成一个基因组片段连接关系；

基因组片段连接确定模块，针对基因组片段连接关系选择模块中每个只能作为起始基因组片段的基因组片段，分别将其作为起始点，从只能连接在其他基因组片段之后作为终止基因组片段的基因组片段，以及既能够连接在其他基因组片之前作为起始基因组片段，又能够连接在其他基因组片段之后作为终止基因组片段的基因组片段中，寻找可连接的基因组片段，形成基因组片段连接，将该基因组片段连接作为新的起始点，进一步如上所述那样寻找可连接的基因组片段，直至没有可连接的基因组片段为止，根据上述各基因组片段连接的前后顺序将各基因组片段连接组装成更长的基因组片段。

10.根据权利要求9所述的装置，其特征在于，

在所述去除模块2中，去除相似程度低于90％的转录组测序读段。