CN103761453B - 一种基于簇图结构的并行基因拼接方法 - Google Patents
一种基于簇图结构的并行基因拼接方法 Download PDFInfo
- Publication number
- CN103761453B CN103761453B CN201310666751.2A CN201310666751A CN103761453B CN 103761453 B CN103761453 B CN 103761453B CN 201310666751 A CN201310666751 A CN 201310666751A CN 103761453 B CN103761453 B CN 103761453B
- Authority
- CN
- China
- Prior art keywords
- scaffold
- cluster
- read
- sequence
- gene
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供一种基于簇图结构的并行基因拼接算法。本发明以多个其他基因拼接算法的拼接得到的长序列（scaffold）和双端测序仪生成的短读长基因序列（read‑pair）为输入，通过构建索引、映射read‑pair、scaffold聚簇、构建簇图、搜索路径等步骤将具有互补的scaffold拼接成更长的序列。构建索引和读长映射两个步骤旨在通过读长找到不同算法获得的长序列scaffold之间的相关性和匹配度，然后通过相关性和匹配度进行聚簇，簇内的所有scaffold具有互补性，是潜在的可拼接的序列。最后构建簇图，求解图的全局最长路径，得到拼接的长基因序列。
Description
技术领域
本发明属于生物信息学技术领域，具体涉及一种新的基于簇图结构的并行基因拼接算法。
背景技术
自从2006年5月18日《自然》杂志报道称，科学家已对含有2.23亿个碱基对占人类基因组中碱基对总量的8%左右的人类第一号染色体完成测序宣告持续16年的人类基因组计划全部完成。作为人类自然科学史上重要的里程碑，“人类基因组”的研究已从“结构基因组”阶段进入“功能基因组”阶段。在人类基因组计划后相继推出的水稻基因组计划、马铃薯基因组计划、草鱼基因组计划等和快速增长的微生物基因测序“海量”的基因信息的积累催生了“功能基因组”时代的来临。针对充分利用“海量”基因组信息的生物信息学不仅应运而生而且为以注释、阐明基因功和利用基因生物学功能的“后基因组时代”的研究发挥了重大作用。
基因组测序的目的就是要确定DNA分子的碱基序列，而DNA序列拼接则是基因组测序的关键技术之一。DNA序列拼接的定义可概括为：从DNA片段集合F中重构该DNA序列S，其中F为该DNA序列S的子序列。无模版拼接算法，是指在没有参考基因作为模板的情况下,根据F重构S。这些重构的DNA源序列可以被进一步的评估和分析，可以成为解决生物问题的线索，如寻找致病病毒、进行药物设计、研究如何将纤维物质转化为生物燃料、揭示生物遗传和变异的。另外，对进行基因诊断、基因治疗、药物设计都有巨大的作用。
基因组拼接的挑战在于将碎片状的读长进行重建得到原始的基因组。基于化学的第一代测序技术--桑格测序方法得到的读长的长度范围从大约500至1000个碱基。像Illumina，Complete Genomics公司、Helicos、454 Life Sciences、SOLID、Ion Torrent公司等这些新一代的技术是以牺牲读长的长度为代价获得高产量。这种海量的、短小的、包含错误的读长数据导致了拼接的高难度。
随着新一代基因组测序技术的推广使用，全基因组Shotgun拼接算法和软件得到了广泛的研究。当前的基因测序技术获得的DNA序列数据相对于第一代测序方法--Sanger测序表现为：高通量、高覆盖率、低成本，与此同时还具有短读长、更多类型的错误等特点，而且普通高等生物的基因组碱基数目巨大，如人类基因组总长约30亿bp。另外，高等生物的基因还具有非常复杂的重复结构，因而基因组的无模板拼接具有很大难度。自从2005年以后，出现了多种基于下一代测序平台基因序列的从头拼接算法软件包，包括:
1. Telescoper(http:// sourceforge.net/p/telescoper)
2. Velvet（http://www.ebi.ac.uk/～zebino/velvet/）
3. ABySS(http://www.bcgsc.ca/platform/bioinfo/software/abyss)
4. AllPath2
5. SOAPdenovo(http://soap.genomics.org.cn/soapdenovo.html)
6. EULER-USR
7. Cortex。
然而，2011年ALKAN等人发表在的Nature报告中指出，使用短读长进行人类基因组无模板拼接的结果比使用长读长得到的拼接结果还短16%。因此，很容易认识到设计出更好的算法进行基因拼接还有很大的发展空间。
本发明提供一种基于簇图结构的并行基因拼接算法进行全基因组拼接。本发明创建了一种适于并行的、快速的簇图结构，可运行于多种操作系统平台（Linux、Windows、Mac），可将大多数基因拼接算法得到的长基因序列进行聚簇并进行再拼接，从而得到更长的基因序列。以多个其他基因拼接算法的拼接得到的长序列和双端测序仪生成的短读长基因序列为输入，通过构建索引、映射短读长、scaffold聚簇、构建簇图、搜索路径等步骤将具有互补的scaffold拼接成更长的序列。构建索引和读长映射两个步骤旨在通过read-pair找到不同算法获得的长序列scaffold之间的相关性和匹配度，然后通过这个相关性和匹配度进行聚簇，簇内的所有scaffold具有互补性，是潜在的可拼接的序列。最后构建簇图，求解图的全局最长路径，得到拼接的结果。
发明内容
针对上述现存的基因拼接方法在真实数据集上仍然难以得到生物基因的完整（100%）序列，本发明提供一种基于簇图结构的并行基因拼接算法，算法相对于上述现存的基因拼接方法在多项指标上都有较大的提高。特别的，在大肠杆菌基因数据集上测试，得到的最长序列长度提高的百分比超过了50%。另外，本发明设计的并行计算框架使得当输入大数据集时，有较短的运行时间。对以上多个方法的结果进行详细的分析，我们提出了一种基于簇图的并行基因拼接算法，在使用原始短读长作为输入的同时，运用多个其他拼接方法产生的scaffold，通过构建索引、读长映射、scaffold聚簇、构建簇图等步骤将scaffold拼接成更长的称为基因序列。构建索引和读长映射两个步骤旨在通过读长（read）找到不同算法获得的scaffold之间的相关性，然后通过这个相关性进行聚簇，簇内的所有scaffold具有互补性，是潜在的可拼接的序列。最后构建簇图，求解图的全局最长路径，得到拼接的结果。
实验结果表明，算法得到最长的scaffold序列的长度和scaffold N50等两项指标，相对于目前拼接效果最好的算法Velvet、ABySS、SOAPdenovo等增长的比例高达50%。当更多的基算法结果加入到我们的算法中，结果将会有更大的提高。本文提出的方法大幅度的提高了scaffold的序列的长度，将为进一步对基因评估和分析降低了难度，为解决生物问题提供了更好的线索，将迅速加快生物基因组研究的步伐。
为实现上述目的，本发明公开了如下的技术方案：
一种基于簇图结构的并行基因拼接算法，其特征在于所述基因拼接算法包含创建簇图和搭建并行框架；
其中创建簇图指的是：根据原始基因数据（read-pair）与其他算法生成结果长序列（scaffold）之间的映射结果对scaffold进行相似性和匹配度计算，然后进行聚簇，簇中的两个匹配的scaffold构成scaffold对（scaffold-pair）,所有scaffold-pair中具有多个匹配的区域，以这些区域作为节点，他们之间的连接构成边，创建簇图；
搭建并行框架指的是：贯穿在整个基因拼接算法的各个步骤中，包括读写文件、构建索引、短读长映射、scaffold聚簇、构建簇图、搜索路径等步骤；采用的并行框架对每个步骤中的任务进行分割、执行、合并，执行过程中节省了大量的时间；包括以下步骤：
（1）数据准备：准备本方法所有的输入数据，包括两种数据，一是原始的双端读长（read-pair）数据，这个可以在NCBI上获得；二是来自其他拼接算法的结果数据scaffold；这两类数据分别要进行预处理；
（2）构建索引： 构建索引就是要将来自其他拼接算法的结果数据scaffold所包含的序列建立一个索引结构，这个索引结构为下一步读长映射提供基础。索引构建完毕，将得到每个算法的索引文件；
（3）读长映射：利用索引将原始数据库中的全部的基因序列映射到scaffold上。首先将上一步中生成的索引文件读入到内存，接下来就是对读长进行映射了，映射的方式并没有采用读长序列中所有的碱基，而是只使用了读长对的内侧的一部分（L=3*k-mer），所谓内侧是指left read的右端和right read的左端。规定只有这部分映射成功之后，整个读长对就可映射成功，映射结果表现为一个scaffold的不同的位置上有多个read与之映射；
（4）Scaffold聚簇：为了下一步进行拼接生成簇图，我们首先对所有的scaffold进行聚簇。每个scaffold的特征由上一步映射结果得到的read的集合 来体现，根据计算不同scaffold对应的read集合之间的相关性和scaffold之间的匹配程度，我们找到互补的、潜在的、可拼接的scaffold对（scaffold-pair），并将他们聚到同一个簇中，对于每一个簇将会通过构建簇图并寻找最优路径的方式得到算法的产物scaffold；
（5）构建簇图：构建簇图的过程包括生成子图和合并子图两个步骤，即对于簇中第个scaffold生成子图，然后将簇中所有合并成能表示一个簇的最终图，最后求解簇图的最长路径。最长路径所包含的碱基序列即为我们算法拼接之后的结果；
（6）生成拼接结果：得到簇图之后，通过计算簇图的最长路径，便可以输出算法结果，即由多个Contig组装成的Scaffold。
其中所述的数据准备步骤要下载的基因序列原始数据要求是来自Illumina测序平台生成的双端短序列，文件格式要是fasta或fastq格式，其他格式的文件需要先进行转换；对于从NCBI官方网站上下载的数据一般是SRA格式，需要使用SRA Toolkit工具包将下载得到的*.sra文件转化成要求的fastq和fasta格式的数据文件。运行命令：
$ fastq-dump --split-files
转换得到两个文件，是读长对（read-pair）分别存储的左读长（left reads）和右读长（right reads）的fastq文件；
数据准备步骤要准备的第二类数据是长序列scaffold文件；该文件是其他基因拼接算法的结果文件，所以需要配置并运行这些算法，并得到最终结果，这些拼接算法可以是Velvet、AbySS、SOAPdenovo、Ray。
数据预处理要求处理未知碱基，DNA序列中碱基只有四种，即A、C、G、T；然而由于测序过程中的一些技术限制或错误导致了未能准确区别两种碱基，从而生成了不确定的非A、C、G、T碱基，需要采用一定的方法将这些不确定的碱基进行确定话。
其中所述的索引结构是首先共享的索引可供多个线程访问，其次索引结构是采用的是预分配空间直接存取的方式，这种方式节省了映射过程中查找序列的时间。
其中所述的短读长（read）高通量基因测序平台产生的序列，一次测序中仪器读取的核苷酸序列，该序列是原始DNA序列经过随机打断生成的碎片序列，基因序列的无模板拼接（de novo assembly）就是要将这些碎片序列拼接成更长的序列。高通量测序中read一般会成对出现，也就是以read-pair形式存在。
其中所述的Scaffold，是由其他拼接方法产生的更长的序列，在实际情况中，scaffold和scaffold之间并不能直接连接起来，很多情况下是它们之间只有通过它们内部的一些小的read之间的某些距离信息或者mate信息进行连接，它借助其他reads之间的关系信息，把contig直接的缝隙进行填充。
其中所述的实现的算法软件包可以运行在64位或32位Linux、Mac、Windows 等多类型的操作系统中，推荐使用64-bit，系统需要的软件包依赖是Java、R、rJava 包；其中Java支持32-bit版本，推荐使用64-bit。版本选用JDK1.6版以上（包括1.6）,下载网址：http://java.com/en/download/manual.jsp。运行软件包时可以修改相关的配置文件，以软件包分配合适的运行时内存。Linux系统下可安装OpenJDK1.6版以上（包括1.6）。
其中所述的运行算法的主机内存要求110G以上，所需内存大小主要是由基因数据集的测试深度和物种的基因组序列的长度决定的，实验中用到的是测试深度约为500、物种的基因组序列的长度大约为数据集(ERR022075)大约消耗内存110Gb。
其中所述的运行算法的主机处理器要求是多核的。本发明中设计的并行框架将合理的调用多个不同的处理器，并为其分配各自的任务，最终将完成任务之后的结果进行汇总。对于大型的数据集，本发明设计的并行框架节省了大量的运行时间，相对于其他算法也具有很明显的优势。
其中所述的软件包依赖，其中R软件包含2.5.X版本以及以上版本，下载网址:http://www.r-project.org/。
其中所述的包依赖，其中rJava包：在R中安装rJava软件包，命令：install.packages("rJava")。R和rJava包是用于绘图，提供了用于绘制簇图和相关性能分析和评价的可视化接口。
本发明更加详细的方法如下：
（1）数据准备。准备本方法所有的输入数据，包括两种数据，一是原始的双端读长（read-pair）数据，这个可以在NCBI上获得；二是来自其他拼接算法的结果数据scaffold；这两类数据分别要进行预处理。
（2）构建索引。 构建索引就是要将来自其他拼接算法的结果数据scaffold所包含的序列建立一个索引结构，这个索引结构为下一步读长映射提供基础。索引构建完毕，将得到每个算法的索引文件。
（3）读长映射。利用索引将原始数据库中的全部的基因序列映射到scaffold上。首先将上一步中生成的索引文件读入到内存，接下来就是对读长进行映射了，映射的方式并没有采用读长序列中所有的碱基，而是只使用了读长对的内侧的一部分（L=3*k-mer），所谓内侧是指left read的右端和right read的左端。规定只有这部分映射成功之后，整个读长对就可映射成功。映射结果表现为一个scaffold的不同的位置上有多个read与之映射。
（4）Scaffold聚簇。为了下一步进行拼接生成簇图，我们首先对所有的scaffold进行聚簇。每个scaffold的特征由上一步映射结果得到的read的集合来体现，根据计算两个scaffold对应的read集合之间的相关性（具体计算方法见公式（1））和scaffold之间的匹配程度，我们找到互补的、潜在的、可拼接的scaffold对（scaffold-pair），并将他们聚到同一个簇中，对于每一个簇将会通过构建簇图并寻找最优路径的方式得到算法的结果序列。表示和匹配罚分公式如下：
公式（1）
我们发明了公式（1）用于定义scaffold-pair中两个scaffold之间的相关性，其中和分别scaffold-pair中的两个scaffold映射的read集合。相关性定义为两个集合的公共read-pair的个数与scaffold-pair中短序列的长度之比。
公式（2）是公认的求解最长公共子序列问题过程中的递归公式。用于对相关性达到阈值的scaffold-pair进行求解最长公共最长子序列以及其匹配的位置。
（5）构建簇图。构建簇图的过程包括生成子图和合并子图两个步骤，即对于簇中第个scaffold生成子图，然后将簇中所有合并成能表示一个簇的最终图，最后求解簇图的最长路径。最长路径所包含的碱基序列即为我们算法拼接之后的结果。
（6）生成拼接结果。得到簇图之后，通过计算簇图的最长路径，便可以输出算法结果，即由多个scaffold拼接成的更长的序列。
基因序列拼接一直以来是一个未解决的问题，其难度较大且十分有意义，其中无模版的序列拼接难度远远大于重测序，也更具挑战性。在深入分析当前基因测序软件的特性以及拼接数据高通量、短序列的特点之后，结合当前序列拼接问题的研究现状，针对当前序列拼接的结果序列中的contig和scaffold之间的互补性等方面的问题，对现有的拼接结果contig和scaffold进行了重新拼接，提出提出通过构建簇图的方式解决当前基因序列拼接的问题。
本发明公开的基于簇图结构的并行基因拼接算法相对于现有技术的主要特征包括以下几个方面:
1．通过将基本拼接算法拼接得到的scaffold经过处理之后，以k-mer大小的序列进行构建索引，并生成索引文件，索引文件提供下一步映射使用。将索引由内存中写入硬盘文件中，可减少大量的内存消耗。而且方便的实现了索引的可重用性，下次需要用时直接载入就行，不用每次都要重建索引了；
2．读长映射过程，使用mate-pair两端各只用30bp(原始读长的长度是100bp)进行映射。映射时容许一定的错误率，要求不能30bp中不超过2bp，如此一来我们就可以通过将30bp分成3段，只要其中有一段能够100%匹配，就说明该读长可以映射成功。这种方式不仅减少了序列比对的时间，而且省去了错误率统计的操作；
3．通过使用多进程并行读取基因数据文件的方式，实现了基因数据的读写的动态拆分和合并，解决了单个计算机使用大规模基因数据内存不足的问题,实现了基因数据资源的负载均衡分布；
4．针对寻找scaffold之间的重复区域的问题，我们没有采用直接的序列比对方式，而是通过统计两条scaffold被映射上的相同的短读长的数目实现的。当两个scaffold被映射上的相同的读长超过一定的数量，就初步认为scaffold是相关的，这样的两个scaffold就是一个scaffold-pair，下一步将通过比对序列的方式计算这个scaffold-pair的最长公共子序列；
5．找到所有的scaffold-pair之后，采用对所有的scaffold-pair进行聚簇的方式而不是直接进行图的构建。通过scaffold-pair聚簇，能够将不相关的、不可拼接的scaffold-pair过滤，从而大大降低了下一步构建簇图的复杂性，也减小了回溯法求解最长路径的解空间。
6．最后，通过寻找簇图中的无前驱的节点，以这些节点为根节点，利用回溯法搜索解空间，求出簇图的最长路径。再根据簇图的最长路径中的scaffold信息得到拼接之后的长序列。
通过与其他拼接算法的比较，本发明在所向指标上都有显著的提高，尤其是最长的scaffold和scaffold N50两个指标在多个拼接算法的基础上提高了高达50%的效果。本发明是基于其他不同算法结果之间的互补性，因此当两种互补性较高的算法的结果加进来时，本发明的优势就更加显著。
附图说明
图1是短序列read-pair示意图；
图2是本发明中的提出的算法框架；
图3是本发明与其他三个方法进行相比提高的比例柱状图；
图4是本发明中的提出的算法构建的子图示意图；
图5是本发明中的提出的算法构建的簇图示意图。
具体实施方式
为了简单和清楚的目的，下文恰当的省略了公知技术的描述，以免那些不必要的细节影响对本技术方案的描述。以下结合较佳实施例，对本发明做进一步的描述。
实施例1
一种基于簇图结构的并行基因拼接算法，法包含创建簇图和搭建并行框架；
其中创建簇图指的是：根据原始基因数据（短读长）与其他算法生成结果长序列（scaffold）之间的映射结果对scaffold进行相似性和匹配度计算，然后进行聚簇，簇中的两个匹配的scaffold构成scaffold对（scaffold-pair）,所有scaffold-pair中具有多个匹配的区域，以这些区域作为节点，他们之间的连接构成边，创建簇图；
搭建并行框架指的是：贯穿在整个基因拼接算法的各个步骤中，包括读写文件、构建索引、短读长映射、scaffold聚簇、构建簇图、搜索路径等步骤；采用的并行框架对每个步骤中的任务进行分割、执行、合并，执行过程中节省了大量的时间；包括以下步骤：
（1）数据准备：准备本方法所有的输入数据，包括两种数据，一是原始的双端读长（read-pair）数据，这个可以在NCBI上获得；二是来自其他拼接算法的结果数据scaffold；这两类数据分别要进行预处理；
（2）构建索引： 构建索引就是要将来自其他拼接算法的结果数据scaffold所包含的序列建立一个索引结构，这个索引结构为下一步读长映射提供基础。索引构建完毕，将得到每个算法的索引文件；
（3）读长映射：利用索引将原始数据库中的全部的基因序列映射到scaffold上。首先将上一步中生成的索引文件读入到内存，接下来就是对读长进行映射了，映射的方式并没有采用读长序列中所有的碱基，而是只使用了读长对的内侧的一部分（L=3*k-mer），所谓内侧是指left read的右端和right read的左端。规定只有这部分映射成功之后，整个读长对就可映射成功，映射结果表现为一个scaffold的不同的位置上有多个read与之映射；
（4）Scaffold聚簇：为了下一步进行拼接生成簇图，我们首先对所有的scaffold进行聚簇。每个scaffold的特征由上一步映射结果得到的read的集合来体现，根据计算不同scaffold对应的read集合之间的相关性和scaffold之间的匹配程度，我们找到互补的、潜在的、可拼接的scaffold对（scaffold-pair），并将他们聚到同一个簇中，对于每一个簇将会通过构建簇图并寻找最优路径的方式得到算法的产物scaffold；
（5）构建簇图：构建簇图的过程包括生成子图和合并子图两个步骤，即对于簇中第个scaffold生成子图，然后将簇中所有合并成能表示一个簇的最终图，最后求解簇图的最长路径。最长路径所包含的碱基序列即为我们算法拼接之后的结果；
（6）生成拼接结果：得到簇图之后，通过计算簇图的最长路径，便可以输出算法结果，即由多个Contig组装成的Scaffold。
所述的数据准备步骤要下载的基因序列原始数据要求是来自Illumina测序平台生成的双端短序列，文件格式要是fasta或fastq格式，其他格式的文件需要先进行转换；对于从NCBI官方网站上下载的数据一般是SRA格式，需要使用SRA Toolkit工具包将下载得到的*.sra文件转化成要求的fastq和fasta格式的数据文件。运行命令：
$ fastq-dump --split-files
转换得到两个文件，是读长对（read-pair）分别存储的左读长（left reads）和右读长（right reads）的fastq文件；
数据准备步骤要准备的第二类数据是长序列scaffold文件；该文件是其他基因拼接算法的结果文件，所以需要配置并运行这些算法，并得到最终结果，这些拼接算法可以是Velvet、AbySS、SOAPdenovo、Ray。
数据预处理要求处理未知碱基，DNA序列中碱基只有四种，即A、C、G、T；然而由于测序过程中的一些技术限制或错误导致了未能准确区别两种碱基，从而生成了不确定的非A、C、G、T碱基，需要采用一定的方法将这些不确定的碱基进行确定话。
所述的索引结构是首先共享的索引可供多个线程访问，其次索引结构是采用的是预分配空间直接存取的方式，这种方式节省了映射过程中查找序列的时间。
所述的短读长（read）高通量基因测序平台产生的序列，一次测序中仪器读取的核苷酸序列，该序列是原始DNA序列经过随机打断生成的碎片序列，基因序列的无模板拼接（de novo assembly）就是要将这些碎片序列拼接成更长的序列。高通量测序中read一般会成对出现，也就是以read-pair形式存在。
所述的Scaffold，是由其他拼接方法产生的更长的序列，在实际情况中，scaffold和scaffold之间并不能直接连接起来，很多情况下是它们之间只有通过它们内部的一些小的read之间的某些距离信息或者mate信息进行连接，它借助其他reads之间的关系信息，把contig直接的缝隙进行填充。
所述的方法实现的算法软件包可以运行在64位或32位Linux/Mac/Windows 等多类型的操作系统中，推荐使用64-bit，系统需要的软件包依赖是Java、R、rJava 包；其中Java支持32-bit版本，推荐使用64-bit。版本选用JDK1.6版以上（包括1.6）。运行软件包时可以修改相关的配置文件，以软件包分配合适的运行时内存。Linux系统下可安装OpenJDK1.6版以上（包括1.6）。其内存要求110G以上，所需内存大小主要是由基因数据集的测试深度和物种的基因组序列的长度决定的，实验中用到的是测试深度约为500、物种的基因组序列的长度大约为数据集(ERR022075)大约消耗内存110Gb。
对于处理器应该是多核的，核数的多少直接影响拼接执行的时间。软件包依赖，其中R包含2.5.X版本以及以上版本，下载网址。软件包依赖，其中rJava包：在R中安装rJava软件包，命令：install.packages("rJava")。R和rJava包是用于绘图，提供了用于绘制簇图和相关性能分析和评价的可视化接口。
实施例2
本发明提出的一种基于簇图结构的并行基因拼接算法，可以在多类型操作系统（Linux、Mac、Windows）上运行，运行的方式很简单。所述方案的具体运行方式包括以下步骤：
（1）在操作系统上安装权利要求中的所有软件包依赖；
（2）准备两类数据，数据一是原始双端测序基因短序列，二是以数据一作为多个其他基因拼接算法的输入得到的输出（长序列）；
（3）修改config.cfg文件中的路径和参数；
#------------input----------------
#########Mapping reads#####
Kmer_Size=30
Available_Processor_Num=20
Read_1=/home/ub/genome/realdata/SRR034959/fasta/SRR034959_1.fasta
Read_2=/home/ub/genome/realdata/SRR034959/fasta/SRR034959_2.fasta
Contig_File=/home/ub/genome/realdata/SRR034959/abyss/k64/SRR034959-scaffolds.fa
Map_Output_Dir=/home/ub/genome/realdata/SRR034959/xk/map/abyss_k64
LogFile_Path=/home/ub/genome/realdata/SRR034959/xk/map/abyss_k64/log.txt
#######Related contigs#####
#Test option
Test=false
#two contigs in algrithm:0-23-1-227
Test_Param=1-0-0-369
Gap_Size=30
ComSeq_Min_Count=100
ComReadpairs_Count=1000
Readpairs_Count=17404920
#contig or scaffolds file from read-pairs map file
MapFile_Path0=D:\share\data\SRR034959\ray\map\ray_k50\all.l-r.res
MapFile_Path1=D:\share\data\SRR034959\velvet\map\velvet_k50\all.l-r.res
#------------output----------------
RCtg_Output_Dir=D:\share\data\SRR034959\xk\Rctg\Rctg_1000
#LogFile_Path=D:\share\data\SRR034959\xk\Rctg\Rctg_1000\log.log
（4）运行创建索引和read映射程序（例如Linux版）
$ ./Mapread.sh
（4）运行scaffold聚簇程序；
（5）运行簇图构建程序
（6）得到结果文件。
实施例3
下表是本发明的方法与现有的三个常规基因拼接算法(ABySS, Velvet,SOAPdenove)
在大肠杆菌K-12 MG1655 (NCBI SRA accession
ERR022075, http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgiview=run_browser&run=ERR022075) 数据集的结果的对比，其中#Aby表示算法ABySS的实验编号、#Vel表示算法Velvet的实验编号、#Soa表示算法SOAPdenove的实验编号，#Cob表示本发明的实验编号，显然本发明的优势非常明显。
结论：
（1）本发明的方法大幅度的提高了scaffold的序列的长度。经测试在大肠杆菌基因数据集上得到的最长序列长度提高的百分比超过了其他常规算法50%。
（2）本发明的方法进行基因拼接得到的长序列为进一步对基因评估和分析降低了难度，为解决生物问题提供了更好的线索，将迅速加快生物基因组研究的步伐。
（3）本发明设计的方法将杂乱无章的序列进行聚簇，将不相关的、不可拼接的序列分离，并将基因序列拼接问题转化为构建多个簇图结构和寻找路径的问题，从而简化了常规算法中使用到的复杂的de Bruijn图结构，降低了解决问题的复杂度。
（4）本发明设计的并行计算框架在多个步骤中进行任务分配和合并。例如读取基因数据文件的方式，实现了基因数据的读写的动态拆分和合并。解决了单个计算机使用大规模基因数据内存不足的问题，实现了基因数据资源的负载均衡分布。
上述技术方案只是本发明的一种运行方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选地，而并不具有限制性的意义。
Claims (8)
1.一种基于簇图结构的并行基因拼接方法，其特征在于所述基因拼接方法包含创建簇图和搭建并行框架；
其中创建簇图指的是：根据原始基因数据与Velvet、ABySS、SOAPdenovo、Ray这些算法
生成结果长序列之间的映射结果对scaffold进行相似性和匹配度计算，然后进行聚簇，簇中的两个匹配的scaffold构成scaffold对，所有scaffold对中具有多个匹配的区域，以这些区域作为节点，它们之间的连接构成边，创建簇图；搭建并行框架指的是：贯穿在整个基因拼接方法的各个步骤中，包括读写文件、构建索引、短读长映射、scaffold聚簇、构建簇图、搜索路径步骤；采用的并行框架对每个步骤中的任务进行分割、执行、合并，执行过程中节省了大量的时间，包括以下步骤：
（1）数据准备：准备所有的输入数据，包括两种数据，一是原始的双端读长数据，这个在NCBI上获得；二是Velvet、ABySS、SOAPdenovo、Ray这些算法的结果数据scaffold；这两类数据分别要进行预处理；
（2）构建索引： 构建索引就是要将结果数据scaffold所包含的序列建立一个索引结构，这个索引结构为下一步读长映射提供基础；索引构建完毕，将得到每个算法的scaffold的索引文件；
（3）读长映射：利用索引将read-pair映射到scaffold上，首先将上一步中生成的索引文件读入到内存，接下来就是对读长进行映射了，映射的方式并没有采用读长序列中所有的碱基，而是只使用了读长对的内侧的一部分，L=3*k-mer，所谓内侧是指left read的右端和right read的左端；规定只有这部分映射成功之后，整个读长对就可映射成功，映射结果表现为一个scaffold的不同的位置上有多个read与之映射；
（4）scaffold聚簇：为了下一步进行拼接生成簇图，首先对所有的scaffold进行聚簇；每个scaffold的特征由上一步映射结果得到的read的集合来体现，根据计算不同scaffold对应的read集合之间的相关性和scaffold之间的匹配程度，找到互补的、潜在的、可拼接的scaffold对，并将它们聚到同一个簇中，对于每一个簇将会通过构建簇图并寻找最长路径的方式得到长序列；
（5）构建簇图：构建簇图的过程包括生成子图和合并子图两个步骤，即对于簇中第个contig生成子图，然后将簇中所有合并成能表示一个簇的最终图，最后求解簇图的最长路径；最长路径所包含的碱基序列即为方法拼接之后的结果；
（6）生成拼接结果：得到簇图之后，通过计算簇图的最长路径，根据路径信息得到拼接成的基因序列。
2.根据权利要求1所述的基于簇图结构的并行基因拼接方法，其特征在于数据准备步骤要下载的基因序列原始数据要求是来自Illumina测序平台生成的双端短序列，文件格式要是fasta或fastq格式，其它格式的文件需要先进行转换；对于从NCBI官方网站上下载的数据一般是SRA格式，需要使用SRA Toolkit工具包将下载得到的*.sra文件转化成要求的fastq和fasta格式的数据文件；
运行命令：
转换得到两个文件，是读长对分别存储的左读长和右读长的fastq文件；
数据准备步骤要准备的第二类数据是长序列scaffold文件；该scaffold文件是其他基因拼接算法的结果文件，所以需要配置并运行这些算法，并得到最终结果，这些拼接算法是Velvet、ABySS、SOAPdenovo、Ray。
3.根据权利要求1所述的基于簇图结构的并行基因拼接方法，其特征在于，所述索引结构是首先共享的索引可供多个线程访问，其次索引结构采用的是预分配空间直接存取的方式，这种方式节省了映射过程中查找序列的时间。
4.根据权利要求1所述的基于簇图结构的并行基因拼接方法，其特征在于作为方法输入数据的DNA序列是经过随机打断生成的碎片序列，基因序列的无模板拼接就是要将这些碎片序列拼接成更长的序列；高通量测序中read一般会成对出现，也就是以read-pair形式存在。
5.根据权利要求1所述的基于簇图结构的并行基因拼接方法，其中所述scaffold是由Velvet、ABySS、SOAPdenovo、Ray算法产生的更长的序列，在实际情况中，scaffold和scaffold之间并不能直接连接起来，很多情况下是它们之间只有通过它们内部的一些小的read之间的某些距离信息或者mate信息进行连接，它借助其他reads之间的关系信息，把contig之间的缝隙进行填充。
6.根据权利要求1所述的基于簇图结构的并行基因拼接方法，其特征在于所述方法实现的方法软件包可以运行在64位或32位Linux/Mac/Windows 多类型的操作系统中，系统需要的软件包依赖是Java、R、rJava 包。
7.根据权利要求1所述的基于簇图结构的并行基因拼接方法，其特征在于其内存要求110G以上，所需内存大小主要是由基因数据集的测试深度和物种的基因组序列的长度决定的，实验中用到的是测试深度为500、物种的基因组序列的长度为数据集，ERR022075消耗内存110Gb。
8.根据权利要求1所述的基于簇图结构的并行基因拼接方法，其特征在于方法的索引结构，映射过程，scaffold之间的拼接这些部分均采用并行计算的方式。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201310666751.2A CN103761453B (zh) | 2013-12-09 | 2013-12-09 | 一种基于簇图结构的并行基因拼接方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201310666751.2A CN103761453B (zh) | 2013-12-09 | 2013-12-09 | 一种基于簇图结构的并行基因拼接方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN103761453A CN103761453A (zh) | 2014-04-30 |
| CN103761453B true CN103761453B (zh) | 2017-10-27 |
Family
ID=50528689
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201310666751.2A Expired - Fee Related CN103761453B (zh) | 2013-12-09 | 2013-12-09 | 一种基于簇图结构的并行基因拼接方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN103761453B (zh) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104200133B (zh) * | 2014-09-19 | 2017-03-29 | 中南大学 | 一种基于读数和距离分布的基因组De novo序列拼接方法 |
| CA2963868A1 (en) * | 2014-10-10 | 2016-04-14 | Invitae Corporation | Methods, systems and processes of de novo assembly of sequencing reads |
| US20160246921A1 (en) * | 2015-02-25 | 2016-08-25 | Spiral Genetics, Inc. | Multi-sample differential variation detection |
| CN104965999B (zh) * | 2015-06-05 | 2016-08-17 | 西安交通大学 | 一种中短基因片段测序的分析拼接方法及设备 |
| CN109817280B (zh) * | 2016-04-06 | 2023-04-14 | 晶能生物技术（上海）有限公司 | 一种测序数据组装方法 |
| CN107858408A (zh) * | 2016-09-19 | 2018-03-30 | 深圳华大基因科技服务有限公司 | 一种基因组二代序列组装方法和系统 |
| CN107944221B (zh) * | 2017-11-21 | 2020-12-29 | 南京溯远基因科技有限公司 | 一种并行分离核酸片断的拼接算法及其应用 |
| CN110317856B (zh) * | 2018-03-28 | 2023-08-11 | 中国科学院分子植物科学卓越创新中心 | 基于表观组信息低成本组装解析生物核心基因组信息 |
| CN109710314B (zh) * | 2018-12-20 | 2019-11-12 | 四川新网银行股份有限公司 | 一种基于图结构分布式并行模式构建图的方法 |
| CN111028897B (zh) * | 2019-12-13 | 2023-06-20 | 内蒙古农业大学 | 一种基于Hadoop的基因组索引构建的分布式并行计算方法 |
| CN112599195B (zh) * | 2020-11-30 | 2024-04-19 | 中国科学院深圳先进技术研究院 | 一种基因序列拼接方法及应用 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102867134A (zh) * | 2012-08-16 | 2013-01-09 | 盛司潼 | 一种对基因序列片段进行拼接的系统和方法 |
| CN103093121A (zh) * | 2012-12-28 | 2013-05-08 | 深圳先进技术研究院 | 双向多步deBruijn图的压缩存储和构造方法 |
| CN103258145A (zh) * | 2012-12-22 | 2013-08-21 | 中国科学院深圳先进技术研究院 | 一种基于De Bruijn图的并行基因拼接方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4716651B2 (ja) * | 2003-10-29 | 2011-07-06 | 株式会社メイズ | スプライシングバリアントの特定方法 |
| JP2005176730A (ja) * | 2003-12-19 | 2005-07-07 | Hitachi Ltd | ｃＤＮＡ配列をゲノム配列にマッピングする方法 |
-
2013
- 2013-12-09 CN CN201310666751.2A patent/CN103761453B/zh not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102867134A (zh) * | 2012-08-16 | 2013-01-09 | 盛司潼 | 一种对基因序列片段进行拼接的系统和方法 |
| CN103258145A (zh) * | 2012-12-22 | 2013-08-21 | 中国科学院深圳先进技术研究院 | 一种基于De Bruijn图的并行基因拼接方法 |
| CN103093121A (zh) * | 2012-12-28 | 2013-05-08 | 深圳先进技术研究院 | 双向多步deBruijn图的压缩存储和构造方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN103761453A (zh) | 2014-04-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN103761453B (zh) | 一种基于簇图结构的并行基因拼接方法 | |
| Raghavan et al. | A simple guide to de novo transcriptome assembly and annotation | |
| Schbath et al. | Mapping reads on a genomic sequence: an algorithmic overview and a practical comparative analysis | |
| Schmidt et al. | Next-generation sequencing: big data meets high performance computing | |
| Chou et al. | A comparative study of SVDquartets and other coalescent-based species tree estimation methods | |
| Canzar et al. | Short read mapping: an algorithmic tour | |
| AU2014340461B2 (en) | Systems and methods for using paired-end data in directed acyclic structure | |
| Narzisi et al. | Comparing de novo genome assembly: the long and short of it | |
| Bao et al. | Evaluation of next-generation sequencing software in mapping and assembly | |
| Shi et al. | MSOAR 2.0: Incorporating tandem duplications into ortholog assignment based on genome rearrangement | |
| AU2014340461A1 (en) | Systems and methods for using paired-end data in directed acyclic structure | |
| Aparicio et al. | Extending the applicability of graphlets to directed networks | |
| Wei et al. | DBH: a de Bruijn graph-based heuristic method for clustering large-scale 16S rRNA sequences into OTUs | |
| Wei et al. | smsMap: mapping single molecule sequencing reads by locating the alignment starting positions | |
| Henry et al. | WGDTree: a phylogenetic software tool to examine conditional probabilities of retention following whole genome duplication events | |
| Firtina et al. | BLEND: A fast, memory-efficient, and accurate mechanism to find fuzzy seed matches | |
| Runge et al. | Rnabench: A comprehensive library for in silico rna modelling | |
| Vasimuddin et al. | Identification of significant computational building blocks through comprehensive investigation of NGS secondary analysis methods | |
| Saeed et al. | A high performance multiple sequence alignment system for pyrosequencing reads from multiple reference genomes | |
| Saggese et al. | STAble: a novel approach to de novo assembly of RNA-seq data and its application in a metabolic model network based metatranscriptomic workflow | |
| Cascitti et al. | RNACache: A scalable approach to rapid transcriptomic read mapping using locality sensitive hashing | |
| Mokhov et al. | Language and hardware acceleration backend for graph processing | |
| Chen et al. | Constructing consensus genetic maps in comparative analysis | |
| Expósito et al. | BigDEC: A multi-algorithm Big Data tool based on the k-mer spectrum method for scalable short-read error correction | |
| Köster et al. | Massively parallel read mapping on GPUs with the q-group index and PEANUT |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20171027 Termination date: 20191209 |