CN103761453A - 一种基于簇图结构的并行基因拼接算法 - Google Patents

一种基于簇图结构的并行基因拼接算法 Download PDF

Info

Publication number
CN103761453A
CN103761453A CN201310666751.2A CN201310666751A CN103761453A CN 103761453 A CN103761453 A CN 103761453A CN 201310666751 A CN201310666751 A CN 201310666751A CN 103761453 A CN103761453 A CN 103761453A
Authority
CN
China
Prior art keywords
scaffold
read
sequence
bunch
data
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.)
Granted
Application number
CN201310666751.2A
Other languages
English (en)
Other versions
CN103761453B (zh
Inventor
陈科
徐魁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tianjin Polytechnic University
Original Assignee
Tianjin Polytechnic University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Tianjin Polytechnic University filed Critical Tianjin Polytechnic University
Priority to CN201310666751.2A priority Critical patent/CN103761453B/zh
Publication of CN103761453A publication Critical patent/CN103761453A/zh
Application granted granted Critical
Publication of CN103761453B publication Critical patent/CN103761453B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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的集合                                                
Figure 664927DEST_PATH_IMAGE001
来体现,根据计算不同scaffold对应的read集合之间的相关性和scaffold之间的匹配程度,我们找到互补的、潜在的、可拼接的scaffold对(scaffold-pair),并将他们聚到同一个簇中,对于每一个簇将会通过构建簇图并寻找最优路径的方式得到算法的产物scaffold;
(5)构建簇图:构建簇图的过程包括生成子图和合并子图两个步骤,即对于簇
Figure 833445DEST_PATH_IMAGE003
中第
Figure 929577DEST_PATH_IMAGE004
个scaffold生成子图,然后将簇中所有合并成能表示一个簇的最终图
Figure 769860DEST_PATH_IMAGE006
,最后求解簇图的最长路径。最长路径所包含的碱基序列即为我们算法拼接之后的结果;
(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的集合
Figure 353288DEST_PATH_IMAGE007
来体现,根据计算两个scaffold对应的read集合
Figure 665321DEST_PATH_IMAGE007
之间的相关性
Figure 650594DEST_PATH_IMAGE008
(具体计算方法见公式(1))和scaffold之间的匹配程度,我们找到互补的、潜在的、可拼接的scaffold对(scaffold-pair),并将他们聚到同一个簇中,对于每一个簇将会通过构建簇图并寻找最优路径的方式得到算法的结果序列。表示和匹配罚分公式如下:
       公式(1)
我们发明了公式(1)用于定义scaffold-pair中两个scaffold之间的相关性,其中
Figure 271435DEST_PATH_IMAGE007
Figure 48898DEST_PATH_IMAGE010
分别scaffold-pair中的两个scaffold映射的read集合。相关性
Figure 529558DEST_PATH_IMAGE008
定义为两个集合的公共read-pair的个数与scaffold-pair中短序列的长度之比。
Figure 821999DEST_PATH_IMAGE011
          
公式(2)是公认的求解最长公共子序列问题过程中的递归公式。用于对相关性达到阈值
Figure 863117DEST_PATH_IMAGE012
的scaffold-pair进行求解最长公共最长子序列以及其匹配的位置。
(5)构建簇图。构建簇图的过程包括生成子图和合并子图两个步骤,即对于簇
Figure 557404DEST_PATH_IMAGE003
中第
Figure 146648DEST_PATH_IMAGE004
个scaffold生成子图
Figure 191964DEST_PATH_IMAGE005
,然后将簇中所有合并成能表示一个簇的最终图
Figure 322917DEST_PATH_IMAGE006
,最后求解簇图的最长路径。最长路径所包含的碱基序列即为我们算法拼接之后的结果。
(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的集合
Figure 410959DEST_PATH_IMAGE001
来体现,根据计算不同scaffold对应的read集合
Figure 615676DEST_PATH_IMAGE001
之间的相关性和scaffold之间的匹配程度,我们找到互补的、潜在的、可拼接的scaffold对(scaffold-pair),并将他们聚到同一个簇中,对于每一个簇将会通过构建簇图并寻找最优路径的方式得到算法的产物scaffold;
(5)构建簇图:构建簇图的过程包括生成子图和合并子图两个步骤,即对于簇
Figure 814576DEST_PATH_IMAGE003
中第
Figure 217875DEST_PATH_IMAGE004
个scaffold生成子图
Figure 663769DEST_PATH_IMAGE005
,然后将簇中所有合并成能表示一个簇的最终图
Figure 92793DEST_PATH_IMAGE006
,最后求解簇图的最长路径。最长路径所包含的碱基序列即为我们算法拼接之后的结果;
(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表示本发明的实验编号,显然本发明的优势非常明显。
Figure 616178DEST_PATH_IMAGE013
结论:
(1)本发明的方法大幅度的提高了scaffold的序列的长度。经测试在大肠杆菌基因数据集上得到的最长序列长度提高的百分比超过了其他常规算法50%。
(2)本发明的方法进行基因拼接得到的长序列为进一步对基因评估和分析降低了难度,为解决生物问题提供了更好的线索,将迅速加快生物基因组研究的步伐。
(3)本发明设计的方法将杂乱无章的序列进行聚簇,将不相关的、不可拼接的序列分离,并将基因序列拼接问题转化为构建多个簇图结构和寻找路径的问题,从而简化了常规算法中使用到的复杂的de Bruijn图结构,降低了解决问题的复杂度。
(4)本发明设计的并行计算框架在多个步骤中进行任务分配和合并。例如读取基因数据文件的方式,实现了基因数据的读写的动态拆分和合并。解决了单个计算机使用大规模基因数据内存不足的问题,实现了基因数据资源的负载均衡分布。
上述技术方案只是本发明的一种运行方式,对于本领域内的技术人员而言,在本发明公开了应用方法和原理的基础上,很容易做出各种类型的改进或变形,而不仅限于本发明上述具体实施方式所描述的方法,因此前面描述的方式只是优选地,而并不具有限制性的意义。

Claims (10)

1.一种基于簇图结构的并行基因拼接算法,其特征在于所述基因拼接算法包含创建簇图和搭建并行框架;
其中创建簇图指的是:根据原始基因数据(read-pair)与其他算法生成结果长序列(scaffold)之间的映射结果对scaffold进行相似性和匹配度计算,然后进行聚簇,簇中的两个匹配的scaffold构成scaffold对(scaffold-pair),所有scaffold-pair中具有多个匹配的区域,以这些区域作为节点,他们之间的连接构成边,创建簇图;
搭建并行框架指的是:贯穿在整个基因拼接算法的各个步骤中,包括读写文件、构建索引、短读长映射、scaffold聚簇、构建簇图、搜索路径等步骤;采用的并行框架对每个步骤中的任务进行分割、执行、合并,执行过程中节省了大量的时间;
包括以下步骤:
(1)数据准备:准备本方法所有的输入数据,包括两种数据,一是原始的双端读长(read-pair)数据,这个可以在NCBI上获得;二是来自其他拼接算法的结果数据scaffold;这两类数据分别要进行预处理;
(2)构建索引: 构建索引就是要将来自其他拼接算法的结果数据scaffold所包含的序列建立一个索引结构,这个索引结构为下一步读长映射提供基础;
索引构建完毕,将得到每个算法的scaffold的索引文件;
(3)读长映射:利用索引将read-pair映射到scaffold上;
首先将上一步中生成的索引文件读入到内存,接下来就是对读长进行映射了,映射的方式并没有采用读长序列中所有的碱基,而是只使用了读长对的内侧的一部分(L=3*k-mer),所谓内侧是指left read的右端和right read的左端;
规定只有这部分映射成功之后,整个读长对就可映射成功,映射结果表现为一个scaffold的不同的位置上有多个read与之映射;
(4)Scaffold聚簇:为了下一步进行拼接生成簇图,首先对所有的scaffold进行聚簇;
每个scaffold的特征由上一步映射结果得到的read的集合                                               
Figure 525976DEST_PATH_IMAGE002
来体现,根据计算不同scaffold对应的read集合
Figure DEST_PATH_IMAGE004A
之间的相关性和scaffold之间的匹配程度,我们找到互补的、潜在的、可拼接的scaffold对(scaffold-pair),并将他们聚到同一个簇中,对于每一个簇将会通过构建簇图并寻找最长路径的方式得到的长序列;
(5)构建簇图:构建簇图的过程包括生成子图和合并子图两个步骤,即对于簇
Figure 976680DEST_PATH_IMAGE005
中第
Figure 250666DEST_PATH_IMAGE006
个contig生成子图
Figure DEST_PATH_IMAGE008AA
,然后将簇中所有
Figure DEST_PATH_IMAGE009
合并成能表示一个簇的最终图
Figure DEST_PATH_IMAGE011
,最后求解簇图的最长路径;最长路径所包含的碱基序列即为我们算法拼接之后的结果;
(6)生成拼接结果:得到簇图之后,通过计算簇图的最长路径,根据路径信息得到拼接成的基因序列。
2.根据权利要求1所述的方法,其特征在于数据准备步骤要下载的基因序列原始数据要求是来自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碱基,需要采用一定的方法将这些不确定的碱基进行确定话。
3.根据权利要求1所述的方法,其特征在于,所述索引结构是首先共享的索引可供多个线程访问,其次索引结构是采用的是预分配空间直接存取的方式,这种方式节省了映射过程中查找序列的时间。
4.根据权利要求1所述的方法,其特征在于,所述短读长(read)高通量基因测序平台产生的序列,一次测序中仪器读取的核苷酸序列,该序列是原始DNA序列经过随机打断生成的碎片序列,基因序列的无模板拼接(de novo assembly)就是要将这些碎片序列拼接成更长的序列;高通量测序中read一般会成对出现,也就是以read-pair形式存在。
5.根据权利要求1所述的方法,其中所述Scaffold,是由其他拼接方法产生的更长的序列,在实际情况中,scaffold和scaffold之间并不能直接连接起来,很多情况下是它们之间只有通过它们内部的一些小的read之间的某些距离信息或者mate信息进行连接,它借助其他reads之间的关系信息,把contig直接的缝隙进行填充。
6.根据权利要求1所述的方法,其特征在于所述方法实现的算法软件包可以运行在64位或32位Linux/Mac/Windows 等多类型的操作系统中,推荐使用64-bit,系统需要的软件包依赖是Java、R、rJava 包;其中Java支持32-bit版本,推荐使用64-bit;版本选用JDK1.6版以上(包括1.6);运行软件包时可以修改相关的配置文件,以软件包分配合适的运行时内存;
Linux系统下可安装OpenJDK1.6版以上(包括1.6)。
7.根据权利要求1所述所述的方法,其特征在于其内存要求110G以上,所需内存大小主要是由基因数据集的测试深度和物种的基因组序列的长度决定的,实验中用到的是测试深度约为500、物种的基因组序列的长度大约为数据集(ERR022075)大约消耗内存110Gb。
8.根据权利要求1所述所述的方法,其特征在于其处理器是多核的,核数的多少直接影响拼接执行的时间。
9.根据权利要求1所述的方法,其特征在于所述软件包依赖,其中R包含2.5.X版本以及以上版本,下载网址。
10.根据权利要求1所述的方法,其特征在于所述软件包依赖,其中rJava包:在R中安装rJava软件包,命令:install.packages("rJava");R和rJava包是用于绘图,提供了用于绘制簇图和相关性能分析和评价的可视化接口。
CN201310666751.2A 2013-12-09 2013-12-09 一种基于簇图结构的并行基因拼接方法 Expired - Fee Related CN103761453B (zh)

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 true CN103761453A (zh) 2014-04-30
CN103761453B 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)

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104200133A (zh) * 2014-09-19 2014-12-10 中南大学 一种基于读数和距离分布的基因组De novo序列拼接方法
CN104965999A (zh) * 2015-06-05 2015-10-07 西安交通大学 一种中短基因片段测序的分析拼接方法及设备
CN106795568A (zh) * 2014-10-10 2017-05-31 因维蒂公司 测序读段的de novo组装的方法、系统和过程
CN107858408A (zh) * 2016-09-19 2018-03-30 深圳华大基因科技服务有限公司 一种基因组二代序列组装方法和系统
CN107944221A (zh) * 2017-11-21 2018-04-20 南京溯远基因科技有限公司 一种并行分离核酸片断的拼接算法及其应用
CN108140070A (zh) * 2015-02-25 2018-06-08 螺旋遗传学公司 多样品差分变异检测
CN109710314A (zh) * 2018-12-20 2019-05-03 四川新网银行股份有限公司 一种基于图结构分布式并行模式构建图的方法
CN109817280A (zh) * 2016-04-06 2019-05-28 晶能生物技术(上海)有限公司 一种测序数据组装方法
CN110317856A (zh) * 2018-03-28 2019-10-11 中国科学院上海生命科学研究院 基于表观组信息低成本组装解析生物核心基因组信息
CN111028897A (zh) * 2019-12-13 2020-04-17 内蒙古农业大学 一种基于Hadoop的基因组索引构建的分布式并行计算方法
CN112599195A (zh) * 2020-11-30 2021-04-02 中国科学院深圳先进技术研究院 一种基因序列拼接方法及应用

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005135053A (ja) * 2003-10-29 2005-05-26 Maze:Kk スプライシングバリアントの特定方法
US20050159898A1 (en) * 2003-12-19 2005-07-21 Hitachi, Ltd. Method that aligns cDNA sequences to genome sequences
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图的并行基因拼接方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005135053A (ja) * 2003-10-29 2005-05-26 Maze:Kk スプライシングバリアントの特定方法
US20050159898A1 (en) * 2003-12-19 2005-07-21 Hitachi, Ltd. Method that aligns cDNA sequences to genome sequences
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图的压缩存储和构造方法

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104200133A (zh) * 2014-09-19 2014-12-10 中南大学 一种基于读数和距离分布的基因组De novo序列拼接方法
CN104200133B (zh) * 2014-09-19 2017-03-29 中南大学 一种基于读数和距离分布的基因组De novo序列拼接方法
CN106795568A (zh) * 2014-10-10 2017-05-31 因维蒂公司 测序读段的de novo组装的方法、系统和过程
CN108140070A (zh) * 2015-02-25 2018-06-08 螺旋遗传学公司 多样品差分变异检测
CN104965999A (zh) * 2015-06-05 2015-10-07 西安交通大学 一种中短基因片段测序的分析拼接方法及设备
CN109817280A (zh) * 2016-04-06 2019-05-28 晶能生物技术(上海)有限公司 一种测序数据组装方法
CN107858408A (zh) * 2016-09-19 2018-03-30 深圳华大基因科技服务有限公司 一种基因组二代序列组装方法和系统
CN107944221A (zh) * 2017-11-21 2018-04-20 南京溯远基因科技有限公司 一种并行分离核酸片断的拼接算法及其应用
CN110317856A (zh) * 2018-03-28 2019-10-11 中国科学院上海生命科学研究院 基于表观组信息低成本组装解析生物核心基因组信息
CN110317856B (zh) * 2018-03-28 2023-08-11 中国科学院分子植物科学卓越创新中心 基于表观组信息低成本组装解析生物核心基因组信息
CN109710314A (zh) * 2018-12-20 2019-05-03 四川新网银行股份有限公司 一种基于图结构分布式并行模式构建图的方法
CN109710314B (zh) * 2018-12-20 2019-11-12 四川新网银行股份有限公司 一种基于图结构分布式并行模式构建图的方法
CN111028897A (zh) * 2019-12-13 2020-04-17 内蒙古农业大学 一种基于Hadoop的基因组索引构建的分布式并行计算方法
CN112599195A (zh) * 2020-11-30 2021-04-02 中国科学院深圳先进技术研究院 一种基因序列拼接方法及应用
CN112599195B (zh) * 2020-11-30 2024-04-19 中国科学院深圳先进技术研究院 一种基因序列拼接方法及应用

Also Published As

Publication number Publication date
CN103761453B (zh) 2017-10-27

Similar Documents

Publication Publication Date Title
CN103761453A (zh) 一种基于簇图结构的并行基因拼接算法
Minkin et al. Scalable multiple whole-genome alignment and locally collinear block construction with SibeliaZ
Lee et al. Bioinformatics tools and databases for analysis of next-generation sequence data
Chaisson et al. Short read fragment assembly of bacterial genomes
Gu et al. Using SOAPaligner for short reads alignment
Page et al. BamBam: genome sequence analysis tools for biologists
D'Antonio et al. WEP: a high-performance analysis pipeline for whole-exome data
ES2966889T3 (es) Sistemas de procesamiento paralelo y métodos para el análisis altamente escalable de datos de secuencia biológica
WO2016025818A1 (en) Systems and methods for genetic analysis
Guzzi et al. coresnp: Parallel processing of microarray data
Zhao et al. Cloud computing for next-generation sequencing data analysis
Xu et al. An efficient algorithm for DNA fragment assembly in MapReduce
Expósito et al. SMusket: Spark-based DNA error correction on distributed-memory systems
Vasimuddin et al. Identification of significant computational building blocks through comprehensive investigation of NGS secondary analysis methods
Firtina et al. BLEND: A fast, memory-efficient, and accurate mechanism to find fuzzy seed matches
Outten et al. Methods and developments in graphical pangenomics
Ghosh et al. Pakman: a scalable algorithm for generating genomic contigs on distributed memory machines
Zhong et al. GRASP2: fast and memory-efficient gene-centric assembly and homolog search for metagenomic sequencing data
Wang et al. SaAlign: Multiple DNA/RNA sequence alignment and phylogenetic tree construction tool for ultra-large datasets and ultra-long sequences based on suffix array
US10629292B2 (en) Generation and use of simulated genomic data
Lee et al. BulkAligner: A novel sequence alignment algorithm based on graph theory and Trinity
Ahmed et al. A comparative analysis of parallel computing approaches for genome assembly
Kaye Approaches to genome analysis through the application of graph theory
Clemente et al. PROJECTION algorithm for motif finding on gPUs
Sachdeva Efficient parallel algorithms for error correction and transcriptome assembly of biological sequences

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