一种基于基因数值化表达的二代从头组装方法及系统

1.本发明属于生物信息学技术领域，尤其涉及一种基于基因数值化表达的二代从头组装方法及系统。


背景技术：

2.自2005年罗氏推出第一款二代测序仪罗氏454后，生命科学正式开始进入高通量测序时代。illumina（因美纳）系列测序平台的推出，极大地降低了二代测序的价格，使得高通量测序在生命科学各个研究领域得到了广泛的普及。至今，第二代短读长测序技术在全球测序市场上仍然占有绝对的优势地位。
3.在生物信息学领域，基因组组装是基因组分析的核心问题。第二代测序技术虽然大大提高了测序的通量，但其获得的单条序列长度很短，往往只在50~300bp，并且因为测序的覆盖范围更深，导致基因组项目中需要处理的序列数量迅速增加，需要将这些较短的单条序列拼接为较长的完整序列用于进一步分析，提高物种注释分析的准确性。
4.目前由于组装技术的限制和实际情况的复杂性，最终组装得到的序列与真实基因组序列之间仍可能存在差异，甚至只能得到若干条无法进一步连接起来的序列。目前主流的组装算法是dbg（de brujin grapph，德布莱英图）组装算法，二代测序得到的reads（测序片段）长度较短，数量较多，不易发现reads之间的重叠关系，因此采用将reads转化为定长k-mer（长度为k的核苷酸序列）从而寻找重叠关系，构建de bruijn graph，把短序列拼接问题转化为de bruijn graph中的欧拉路径问题。但基于k-mer的dbg组装算法仍然存在较大缺陷，如kmer序列较短，易受测序错误影响组装结果，寻找碱基之间的overlap（重叠）关系无法实现快排等。


技术实现要素：

5.本发明的目的在于提供一种基于基因数值化表达的二代从头组装方法及系统，以解决在基因组组装中，测序片段reads排列速度较慢、生物信息学组装的内存过大、组装错误风险较高的技术问题。
6.本发明的目的是采用以下技术方案实现的：一种基于基因数值化表达的二代从头组装方法，包括如下步骤：s1：通过软磁盘阵列raid实现碱基测序数据管理；s2：基因解析碱基序列并自定义编号；s3：计算每条测序片段reads两组特定长度的前后质心值；s4：通过每条测序片段reads的碱基序列和编号read id构建b 树索引；s5：应用人工智能算法获得前后质心值与碱基序列对齐关系的匹配模型，或使用数值快速排列算法将质心值顺次排序，据此快排结果生成id对接表；s6：根据id对接表多线程组装小型重叠群small contig，输出线程标识表、小型重叠群small contig片段及其编号和单核苷酸变异snp信息表；
s7：根据线程标识构建德布莱英图de brujin grapph并寻路，将小型重叠群small contig组装为重叠群contig；s8：基于重叠群contig组装框架序列scaffold。
7.进一步的，所述步骤s1具体包括：通过内置raid卡或nfs协议，将磁盘中存储的碱基测序数据，以数据流形式陆续调入由高iops的硬盘组成的伪虚拟内存中，并且分散存储，所述碱基测序数据包括质心数据和基因文件。
8.进一步的，所述步骤s2具体包括：通过cpu从高随机读写能力固态硬盘中读取基因文件，得到多条基因测序片段reads，并对每条基因测序片段reads进行编号，得到随机编号id，然后，通过cuda流将基因测序片段reads和随机编号id传输至gpu。
9.进一步的，所述步骤s3具体包括：通过gpu将碱基测序数据分别进行宏层标识，将每条测序片段reads前后分别去m和n个碱基，其中，m和n的数值不同，得到两组共四条拆分后的前后层片段，采取无监督机器学习分别将两组前层与后层序列表达为硬聚类计算的质心值。
10.进一步的，所述步骤s4具体包括：在两组特定长度的前后质心值，将每条测序片段reads的碱基序列及其对应的编号read id构建起b 树索引，以实现碱基序列的随机提取。
11.进一步的，所述步骤s5具体包括：在得到所有测序序列两组不同片段的前后质心值后，使用人工智能算法，将两条序列的两组前后层质心和序列长度作为输入向量，两条序列是否能匹配作为输出向量，应用深度学习神经网络（dnn）获得前后质心值与碱基序列对齐关系的匹配模型，据此生成高准确的id对接表，或使用数值快排算法，将质心顺次排序，据此生成id对接表。
12.进一步的，所述人工智能算法包括以下子步骤：s51：将两条序列的两组前后质心值及其序列长度变换至[-1,1]以内，作为数据集的输入向量x；s52：将两条序列使用kmp算法进行匹配对齐，若能对齐匹配，则返回true，y=1，不能对齐匹配，返回false，y=0，y作为数据集的输出向量；s53：训练前，首先设置一个均值为0的随机分布初始化网络权重，训练时先将输入向量x送至dnn网络的输入层，进行前向传播，经过若干隐层后得到前向输出y'；s54：计算y'与真实值y之间的误差，再从输出层开始将误差反向传播至第一个隐层，按照减小误差的方向调整网络的各个连接权值；s55：不断循环以上调整权值的过程，直到训练误差收敛到极小值，即得到符合要求的dnn网络权重。
[0013]
进一步的，所述步骤s6具体包括：多线程读取id对接表中的编号read id，使用b 树索引分别从硬盘或虚拟内存中提取对应的编号read id，并依次与基准测序片段reads进行匹配，将相似的测序片段reads与基准测序片段reads分别匹配并延伸，组成小型重叠群small contig。
[0014]
进一步的，所述步骤s7具体包括：生成的线程标识表中前后存在线程重叠，其中包含了小型重叠群small contig的重叠信息，通过图论方法构建线程标识的前后关系，构建成一个多级多节的德布莱英图de brujin grapph线程路径，并按照该线程路径进行序列拼接，将小型重叠群small contig组装为重叠群contig。
[0015]
进一步的，所述步骤s8包括以下子步骤：s81：提取每条重叠群contig碱基片段的特征；s82：对重叠群contig碱基片段的特征进行比较，将特征相同的对应的重叠群contig进行部位嵌合，遍历纠错；s83：输出最终组装框架序列scaffold。
[0016]
一种基于基因数值化表达的二代从头组装系统，包括内存管理模块、测序片段reads对接模块、小型重叠群small contig组装模块、重叠群contig组装模块和框架序列scaffold组装模块，其中，所述内存管理模块用以实现碱基测序数据管理；所述测序片段reads对接模块用以基因解析碱基序列并自定义编号、计算每条测序片段reads两组特定长度的前后质心、通过每条测序片段reads的碱基序列和编号read id构建b 树索引，应用人工智能算法获得前后质心值与碱基序列对齐关系的匹配模型，或使用数值快速排列算法将质心值顺次排序，据此生成id对接表；所述小型重叠群small contig组装模块用以根据id对接表多线程组装小型重叠群small contig，输出线程标识表、小型重叠群small contig片段及其编号和单核苷酸变异snp信息表；所述重叠群contig组装模块用以根据线程标识构建德布莱英图de brujin grapph并寻路，将小型重叠群small contig组装为重叠群contig；所述框架序列scaffold组装模块用以基于重叠群contig组装框架序列scaffold。
[0017]
本发明的有益效果在于：本发明创造性地提出了一种基于基因数值化表达的二代从头组装方法及系统，完成二代测序的低成本高性能组装。相较于以往的组装算法，本方法没有将reads拆分为多段kmer，通过序列间overlap构建de bruijn graph有向图，而是通过无监督机器学习等方法将测序碱基数值化，将reads之间overlap的排列顺序转换为数值对应顺序，不仅提高了reads排列的速度，还一定程度上避免了组装过程中可能产生kmer间gap位置的产生，减少反复使用reads回帖比对的次数，在保证准确率的基础上大大地提高了运算速度。在组装scaffold过程中，没有采用read重新回帖比对增加算力成本和io开销的方法，而是采用酶切位点标识辅助组装的方法，提取长片段序列的酶切位点特征，将基因文件进行嵌合，组装为更长的片段。在组装过程中，同时完成了snp变异信息的检测。使用cuda异构编程，采取gpu加速的策略，协调cpu、gpu共同运算，提高并行速度。总而言之，本方法是一种集高精准度、低算力需求、高效率于一体的二代组装方案，大幅提高了基因组装分析的生物信息学分析流程的速度，有助于生信分析人员以更快的速度、更短的时间获得准确的分析结果。
[0018]
依托于基因测序技术以及生物信息学等科学领域的快速发展，本方法为低成本、快速高效的基因组组装分析提供了有力保障，在真正意义上实现了以计算机科学技术为工具，对生物信息进行储存、检索与分析。同时，生物信息学是21世纪自然科学的核心领域之一，本方法作为生物信息学中dna-seq分析的关键技术，在基因组学中具有极高技术优势和市场价值，适合于基因分析的实际应用与技术上的进一步推广。
[0019]
本发明的核心突破点至少包括：1、彻底破除了生物信息学组装的内存过大难题。2、实现了完全遍历组装，极大程度降低了组装错误风险。3、运用异构异步计算实现了计算效率的极大提升。
附图说明
[0020]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
[0021]
图1为发明流程图；图2为de bruijn graph路径泡状结构图；图3为de bruijn graph路径单出入度结构图；图4为de bruijn graph路径环状结构图；图5为头部对齐示意图；图6为尾部对齐或参差对齐示意图；图7为小型重叠群small contig组装流程图；图8为重叠群contig组装流程图；图9为框架序列scaffold组装流程图；图10为本发明系统框图。
具体实施方式
[0022]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
[0023]
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
[0024]
下面结合附图，对本发明的一些实施方式作详细说明。在不冲突的情况下，下述的实施例及实施例中的特征可以相互组合。
[0025]
实施例1：参阅图1，一种基于基因数值化表达的二代从头组装方法，包括如下步骤：s1：通过软磁盘阵列raid实现碱基测序数据管理；s2：基因解析碱基序列并自定义编号；s3：计算每条测序片段reads两组特定长度的前后质心值；s4：通过每条测序片段reads的碱基序列和编号read id构建b 树索引；s5：应用人工智能算法获得前后质心值与碱基序列对齐关系的匹配模型，或使用数值快速排列算法将质心值顺次排序，据此生成id对接表；s6：根据id对接表多线程组装小型重叠群small contig，输出线程标识表、小型重叠群small contig片段及其编号和单核苷酸变异snp信息表；s7：根据线程标识构建德布莱英图de brujin grapph并寻路，将小型重叠群small contig组装为重叠群contig；s8：基于重叠群contig组装框架序列scaffold。
[0026]
在本实施例当中，所述步骤s1具体包括：通过内置raid卡或nfs协议，将磁盘中存
储的碱基测序数据，以数据流形式陆续调入由高iops的硬盘组成的伪虚拟内存中，并且分散存储，所述碱基测序数据包括质心数据和基因文件。所述步骤s2具体包括：通过cpu从高随机读写能力固态硬盘中读取基因文件，得到多条基因测序片段reads，并对每条基因测序片段reads进行编号，得到随机编号id，然后，通过cuda流将基因测序片段reads和随机编号id传输至gpu。所述步骤s3具体包括：通过gpu将碱基测序数据分别进行宏层标识，将每条测序片段reads前后各去n个碱基（n值一大一小），得到两组共四条拆分后的前后层片段，采取无监督机器学习分别将两组前层与后层序列表达为硬聚类计算的质心值，如kmeans聚类。所述步骤s4具体包括：两组特定长度的前后质心值，将每条测序片段reads的碱基序列及其对应的编号read id构建起b 树索引，以实现碱基序列的随机提取。所述步骤s5具体包括：在得到所有测序序列两组不同片段的前后质心值后，使用人工智能算法，将两条序列的两组前后层质心和序列长度作为输入向量，两条序列是否能匹配作为输出向量，应用深度学习神经网络（dnn）获得前后质心值与碱基序列对齐关系的匹配模型，据此生成高准确的id对接表，或使用数值快排算法，将质心顺次排序，据此生成id对接表。
[0027]
进一步的，步骤s1~s5，可以看作测序片段reads的对接，reads对接大体上可以分为三步，目的是要得到根据碱基重叠部分初步排列出来的reads顺序。第一步实现碱基测序数据管理，通过内置raid卡或nfs协议等将磁盘中储存的数据调入虚拟内存中，或者使用windows系统自带的虚拟内存。同时调用gpu启动无监督机器学习核函数，将每条测序片段reads分别用前后层质心进行标识，或者计算每条测序片段reads两组特定长度的前后质心值。第二步使用每条测序片段reads的碱基序列和编号read id构建b 树索引，并对质心值顺次进行快排。第三步通过设定阈值或计算数据相似度等方法将相似度高的前后层质心值对应的随机编号id对接到一起，或将两条序列的两组前后层质心和序列长度作为输入向量，两条序列是否能匹配作为输出向量，应用深度学习神经网络（dnn）获得前后质心值与碱基序列对齐关系的匹配模型，得到更为准确的id对接表，至此，得到了全部测序数据测序片段reads的大致排列顺序，即olc（overlap-layout-consensus，重叠-布局-共识）碱基序列。对于小型基因数据，可以首先采用多线程方案对每个对接表进行一次全序列io，删除无法匹配的编号read id，将id对接表进行优化删减，得到一个完全正确的clean id对接表。
[0028]
测序片段reads的对接需要使用到基因数值化表达的方法，该方法主要分为两个层次：第一层次为宏层标识层，下机后的测序fastq文件（基因文件），经过质控去除接头、低质量碱基等预处理后得到高质量的clean data，针对预处理后的fastq文件（即clean data），摒弃传统dbg组装算法中将测序片段reads打断为固定长度k-mer的方法，将r1文件和r2（基因文件）文件分别进行宏层标识，将一条测序片段reads前后分别去m和n个碱基，其中，m和n的数值不同，得到两组共四条拆分后的前后层片段，采取无监督机器学习分别将两组前层与后层序列表达为kmeans计算的质心值，使用前层质心与后层质心共同标识一条测序片段reads。例如一段预处理后的clean read为atttcgggtcgaaa，长度为14bp，首尾各裁去3个碱基，拆为atttcgggtcg和tcgggtcgaaa，分别计算其前后层质心，则前层质心和后层质心能够共同标识一条测序片段reads。具体地，在宏层质心计算中，设计a、t为垂直坐标的正负轴，g、c为水平坐标的正负轴。考虑碱基的位置，将拆分后的前后序列中每一个碱基表达为坐标系中的点，使得碱基的坐标点散布在笛卡尔坐标系的坐标轴上。例如一个前层序
列：tgctt，t为第一个碱基，其坐标为（0，-1），g为第二个碱基，其坐标为（2，0），以此类推，可将第一段前层序列转换为{（0，-1），（2，0），（-3，0），（0，-4），（0，-5）}的坐标数组。随后运用kmeans聚类，以样本点间欧几里得距离的均方差作为准则函数，计算出这组坐标点的质心，其可以作为该段前层序列的特征值，同理计算后层序列特征值。计算过程中通过调整权重，保证前层质心对前20bp左右碱基的变化更为敏感，其余部分少数几个碱基变化不影响最终生成的前层质心值，后层质心对最后20bp左右碱基的变化更为敏感，其余部分少数几个碱基变化不影响最终生成的后层质心值，用相同质心表达尽量多的相似测序片段reads。
[0029]
本发明提取前后层序列特征标识测序片段reads的方法不限于k均值聚类算法，还包括其他聚类算法，如基于高斯混合模型的最大期望聚类算法、均值漂移聚类算法等，以及傅里叶系数提取特征等，本发明的优选实施例为kmeans聚类分析算法进行质心计算。
[0030]
第二层次为质心快排层，宏层质心计算之后得到的文件包括测序片段reads的碱基序列和前后层质心，对其进行自定义编号，作为编号read id。并且将编号read id对应的前后层质心按照大小进行数值快排，形成快排表。通过内置raid卡（磁盘阵列卡）或者nfs协议（network file system，网络文件系统），将磁盘中存储的质心数据及基因文件以数据流形式陆续调入由高iops（input/output operating per second，每秒进行读写（i/o）操作的次数）的硬盘组成的伪虚拟内存中，并且分散存储，构建软raid。将fastq文件的reads id和软raid中随机储存的id对应的碱基序列构建b 树索引，方便后续能够根据read id从硬盘中快速提取碱基序列。
[0031]
对于大型基因组，基于软raid和特种硬盘进行提取，对于小型基因组组装，可以将所有中间文件以及测序数据置于windows系统的虚拟内存中，方便随时调取。
[0032]
参阅图7，所述步骤s6具体包括：多线程读取id对接表中的编号read id，使用b 树索引分别从硬盘或虚拟内存中提取对应的编号read id，并依次与1号测序片段reads进行匹配，将相似的测序片段reads与1号测序片段reads分别匹配并延伸，组成小型重叠群small contig。小型重叠群small contig组装的目的是要得到高准确性的small contig片段和snp信息表，首先多线程读取id对接表中的read id，使用b 树索引分别从硬盘或虚拟内存中提取对应的碱基序列，并依次进行匹配，将极相似的测序片段reads分别匹配并延伸，组成small contig并输出small contig的唯一线程标识、组装好的small contig和snp信息表。
[0033]
参阅图8，所述步骤s7具体包括：生成的线程标识表中前后存在线程重叠，其中包含了小型重叠群small contig的重叠信息，通过图论方法构建线程标识的前后关系，构建成一个多级多节的德布莱英图de brujin grapph线程路径，并按照该线程路径进行序列拼接，将小型重叠群small contig组装为重叠群contig。进一步的，多线程基因组重叠群contig组装大体可以分为三步，目的是要得到高准确性的小型重叠群small contig片段和snp信息表，首先读取完整打上线程标识的id对接表以及线程标识中零散的编号reads id，并将其替换成对接表中的线程标识。接下来按照线程标识的匹配信息构建de bruijn graph有向图，其中每一条路径都包含了contig组装的顺序。最后按照de bruijn graph中的路径提取线程标识对应的小型重叠群small contig，将其组装为重叠群contig，并修改snp对接表中的变异碱基位点以及变异碱基类型等信息。
[0034]
进一步的，步骤s6~s7需要使用到多线程构建基因重叠群contig的方法以及snp变
异检测的方法，其中，多线程构建基因重叠群contig的方法总体可以分为三个层次：第一个层次为id对接层：在将测序片段reads的质心值初步排列完成后，多线程读取1-n号的前层质心和后层质心（假设1号read为0.000001 、0.000001；2.359876 、5.164616），在生成的快排结果中找到与1号后层质心值相似误差均在一定范围内（以
±
0.002为例）的前层质心，并提取该部分的编号reads id，生成id对接表。即1-3、9、18、66、77、105。
[0035]
本发明中在进行id对接时，除了直接设定指定阈值范围（即误差范围）外，还可以采取相似性度量算法，通过计算质心值之间的相似度，将相似度高的质心值对应的id对接到一起。这里的相似度计算方法包括计算前后层质心点的欧几里得距离、余弦相似度、曼哈顿距离等。
[0036]
为了优化id对接表，使reads之间的前后重叠匹配关系更加精确而不冗余，使用深度学习神经网络对数据集进行训练，将质心值的对应规律与碱基序列的相似性相关联。具体地，对每条read裁剪前后各n个碱基计算前后质心值，再提取每条read前60bp及后60bp（不足60bp的按照碱基序列整体计算）并计算前后质心值，由此得到四组质心值。使用kmp将所有碱基序列进行两两匹配，若能对齐匹配（相似）则返回true，不能对齐匹配（不相似）返回false。得到两两进行匹配的序列各自的质心坐标值（4组共8个质心数据）以及两条序列的长度后，将该十个分量（第一条序列去掉头部m个碱基的质心坐标（x1,y1），第一条序列去掉尾部m个碱基的质心值坐标（x2,y2），第一条序列去掉头部n个碱基的质心坐标（x3,y3），第一条序列去掉尾部n个碱基的质心值坐标（x4,y4）,第二条序列去掉头部m个碱基的质心坐标（x5,y5），第二条序列去掉尾部m个碱基的质心值坐标（x6,y6），第二条序列去掉头部n个碱基的质心坐标（x7,y7），第二条序列去掉尾部n个碱基的质心值坐标（x8,y48），第一条序列的长度length1，低二条序列的长度length2）变换至[-1,1]以内，作为数据集的输入向量x。根据匹配计算的结果确定两条序列相似或不相似，相似则y=1，不相似则y=0，y作为数据集的输出向量。在dnn训练前，首先设置一个均值为0的随机分布初始化网络权重。训练时先将输入向量x送至dnn网络的输入层，进行前向传播，经过若干隐层后得到前向输出y'。随后计算y'与真实值y之间的误差，再从输出层开始将误差反向传播至第一个隐层，按照减小误差的方向调整网络的各个连接权值。不断循环往复整个过程，直到训练误差收敛到极小值，即得到符合要求的dnn网络权重。
[0037]
当dnn模型训练完毕后，对目标质心文件进行判别，从而确定不同序列之间的一对多衔接关系，输出碱基序列的id对接表。
[0038]
第二个层次为序列匹配层：多线程读取id对接表中的编号reads id，使用b 树索引分别从硬盘或虚拟内存中提取对应的编号reads id（3、9、18、66、77、105），并依次与1号测序片段reads进行匹配，将极相似的测序片段reads与1号分别匹配并延伸，组成小型重叠群small contig。
[0039]
具体的多线程实施方案为：线程开启计算时，即在纳入线程范围的id对接表中的编号reads id前打上标签，表示进行1-n号碱基序列的逐一匹配。在第一个线程进行计算，匹配1号与其他对接碱基序列时，则在id对接表的第一行（1-3、9、18、66、77、105）前打上1-1标签作为线程标识数，表示第一个线程的第二次计算，并按照编号reads id根据b 树索引提取相应序列进行逐一匹配。当3号、66号测序片段reads的碱基与1号大部分碱基相同，只存在少量碱基的snp变异时，即可认为1号和3号、66号序列存在前后overlap区域，可以合并
为一条。9号、18号、105号序列与1号序列的overlap区域存在大量不同碱基，则舍弃这三条序列。当3号能和1号匹配时，继续搜索3号序列的id对接表时，发现在id对接表前已有第三个线程开启计算的3-1线程标识数，则停止计算，输出本次计算的线程标识数{1-1：3-1，66，77}以及该标识对应的1号碱基序列，该条小型重叠群samll contig序列标识为contig1-1。
[0040]
1号测序片段reads计算完毕后，第一个线程开启第二次计算，寻找id对接表中没有标识的测序片段reads，例如9-101、102，并在其前面打上线程标识1-2。若101号、102号能和9号匹配，且均未被其他线程开启标识，则同时提取101号、102号对接的碱基序列，如101-1007、1008，102-1009、1010进行拼接，直到检索到某一条序列前已有标识时终止计算，若1010号已有其他线程8-15开启过计算，则本次计算结束，并输出标识数{1-2：8-15，1007，1008，1009}以及9号和101、102拼接后的碱基序列，该条小型重叠群small contig序列标识为contig1-2。依次类推。
[0041]
当所有id对接表均计算完毕后，检索全部标记上线程标识的id对接表，将输出的线程标识表进行优化，例如分别有4-7、8-9线程开启过read 66、read77的计算，则将标识表中的{1-1：3-1，66，77}修正为{1-1：3-1，4-7，8-9}。其中线程的标识同时也代表了输出的小型重叠群small contig的编号。
[0042]
第三个层次为图论拼接层：生成的标识表中前后存在线程重叠，这其中包含了小型重叠群small contig的重叠信息。利用图论方法构建其线程标识的前后关系，构成一个多级多节点的de bruijn graph线程路径，并按照该路径进行更进一步的序列拼接，形成更长的重叠群contig片段。
[0043]
具体地，得到一系列的线程标识表之后，构建的de bruijn graph路径图如图2~图4所示。其中，图2为泡状结构、图3为单出入度结构、图4为环状结构，对于泡状结构的图，按照不同路径分别组装，若最后多条路径的组装结果相同，则合并去重之后统一输出，若不同，则分别输出。对于单出入度的图。按照不同路径分别输出组装结果。对于环状结构，则说明可能存在重复序列，按照不同路径组装后单独输出一次，做为重复序列文件，在后期组装时重点关注。在按照id对接表进行匹配组装时，质心值计算匹配的前后碱基过长，可能存在前后10~20bp能够匹配，但是质心值代表的范围过大，从而在id对接表中没有对接上的问题。针对这类重叠群contig片段的组装，需要生成不同长度的质心值，使用dnn（深度神经网络）计算两组代表碱基长短不一致的质心值的匹配关系，从而优化生成更为准确的id对接表。
[0044]
对于本发明中的多线程组装，上述方法适用于基因组数据大的生物组装重叠群contig，若基因数据较小，对于生成的id对接表，可以首先采用多线程方案对每个对接表进行一次全序列io，删除无法匹配的编号reads id，将id对接表进行优化删减，得到一个完全正确的clean id对接表，再根据id对接表整体进行de bruijn graph构建，最后按照路径提取碱基序列进行拼接组装，生成重叠群contig片段。
[0045]
进一步的，snp变异检测的方法：在本发明中，除了完成基因序列的从头组装外，还能检测基因的snp（单核苷酸）变异。在多线程构建重叠群contig时，通过提取测序片段reads进行逐一匹配，则能输出snp的位点信息：在多线程组装小型重叠群small contig时第一次输出snp位点信息，包括该条小型重叠群small contig的标识、全长、snp变异碱基位点、输出重叠群contig在该变异位点的碱基、发生变异后替换的碱基等信息。在图论构建完
成之后，提取小型重叠群small contig组装重叠群contig时，snp变异信息要做相应调整。
[0046]
具体的，当根据路径提取小型重叠群small contig时，其对齐情况有图5和图6所示两种情况，其中，图5为头部对齐，图6为尾部对齐或参差对齐，图5中k1、k2、k3分别为原本snp信息表中的变异位点，如果对齐情况为头部对齐，则将原来的snp变异信息直接输出，如果基准序列和其余序列的变异位点有重合，则合并替换的碱基信息后输出。图6中k1、k2、k3、k4、k5分别为原本snp信息表中的变异位点，l为基准序列的全长，o为比对定位后上下两条碱基的重合度。如果对齐情况为尾部对齐（参差对齐），则基准序列的变异位点不变，其余序列的变异位点依次调整为l-o k3、l-o k4、l-o k5。当基准序列和和其余序列的变异位点有重合，则合并替换的碱基信息后输出。
[0047]
在整体发明中，针对小基因组装，提供cpu计算和gpu计算的环境检测功能，为无gpu的个人电脑提供cpu适配组装方法。整体算法可使用cpu的avx-512深度学习指令集代替短片段进行更复杂的分支计算。
[0048]
参阅图9，所述步骤s8包括以下子步骤：s81：提取每条重叠群contig碱基片段的特征；s82：对重叠群contig碱基片段的特征进行比较，将对应的重叠群contig进行部位嵌合，遍历纠错；s83：输出最终组装框架序列scaffold。
[0049]
进一步的，基于重叠群contig组装框架序列scaffold，在分别得到r1和r2（基因文件）组装好的重叠群contig序列后，将r2的碱基进行互补倒置，转换为5＇－3＇端的序列，并通过指定的酶切识别例如5＇gggcc^c 3＇、5＇g^gatcc 3＇等位点，可采用单酶切标记或双酶切标记，对所有重叠群contig进行位点标记，并输出每条重叠群contig的位点标记间隔表。最后根据间隔数目的信息匹配将位点信息间隔相同的重叠群contig进行进一步组装，形成更长的框架序列scaffold片段，并一步修正snp信息表。除利用酶切位点作为重叠群contig的特征外，还可以使用双向循环神经网络、trasnformer多头注意机制等深度学习方法生成每条重叠群contig的循环特征矩阵，并且将相似特征的重叠群contig提取进行框架序列scaffold拼接。
[0050]
本发明步骤s8采用了酶切位点识别辅助组装scaffold（框架序列）的方法，该方法针对组装完毕后的重叠群contig，摒弃所有测序片段reads重新回帖到重叠群contig上后利用pair-end信息组装框架序列scaffold的方法，使用酶切位点识别辅助组装的方法，通过识别不同的酶切位点，输出标记位点间隔数，从而将重叠群contig组装为框架序列scaffold。
[0051]
具体地，在双端测序中，r1文件通常为5＇－3＇端的测序结果，r2文件通常为3＇－5＇端的测序结果。在分别得到r1和r2组装好的重叠群contig序列后，将r2的碱基进行互补倒置，转换为5＇－3＇端的序列，并通过指定的酶切识别位点，对所有重叠群contig进行位点标记，并输出每条重叠群contig的位点标记表。例如：apaⅰ酶剪切的碱基位点为5＇gggcc^c 3＇，通过识别每条重叠群contig中5＇gggcc^c 3＇并输出位点间隔的数目，以此作为每条重叠群contig的特征，在进行组装时，将间隔数目相同的重叠群contig提取组装为框架序列scaffold。除单酶切外，还可以同时识别两种酶切位点，将两种酶切位点的间隔数目分别用正负数表达，从而进行比较。
[0052]
本发明中组装框架序列scaffold的方法主要依赖重叠群contig序列的特征，除了使用酶切位点作为特征进行组装外，还可以使用双向循环神经网络（bidirectional rnn, bi-rnn）或transformer多头注意机制对重叠群contig序列向前向后提取特征，组成循环特征网络，根据相似特征进行重叠群contig组装。
[0053]
参阅图10，基于同一发明构思，本发明还提出了一种基于基因数值化表达的二代从头组装系统，用以实现上述于基因数值化表达的二代从头组装方法，该系统包括内存管理模块、测序片段reads对接模块、小型重叠群small contig组装模块、重叠群contig组装模块和框架序列scaffold组装模块，其中，所述内存管理模块用以实现碱基测序数据管理；所述测序片段reads对接模块用以基因解析碱基序列并自定义编号、计算每条测序片段reads两组特定长度的前后质心、通过每条测序片段reads的碱基序列和编号read id构建b 树索引，应用人工智能算法获得前后质心值与碱基序列对齐关系的匹配模型，或使用数值快速排列算法将质心值顺次排序，据此生成id对接表；所述小型重叠群small contig组装模块用以根据id对接表多线程组装小型重叠群small contig，输出线程标识表、小型重叠群small contig片段及其编号和单核苷酸变异snp信息表；所述重叠群contig组装模块用以根据线程标识构建德布莱英图de brujin grapph并寻路，、将小型重叠群small contig组装为重叠群contig；所述框架序列scaffold组装模块用以基于重叠群contig组装框架序列scaffold。
[0054]
可以想到的，本发明的目的：由于基因组组装是生物信息学领域的核心问题，是进行基因功能分析、开展下游群体遗传和结构差异等一系列研究的起点。本发明基于坐标偏移算法、无监督机器学习硬聚类算法，使用双层质心标识reads序列；使用线程标签代替测序文件拆分的kmer序列，利用图论方法构建汉密尔顿路径或欧拉路径，寻找reads间的顺序；采用cuda数据流，配合异步操作，避免gpu（graphics processing unit，图形处理器）浪费及内存开销，提高计算速度；采用cuda异构编程，保证cpu（central processing unit/processor，中央处理器）、gpu间高速通讯，协同计算，实现高准确度、高效率的基因测序数据从头组装，缩短基因组分析的时间成本；使用cpu的avx-512指令集，提高多分支复杂计算的效率。
[0055]
本发明至少具有以下技术效果：本发明创造性地提出了一种基于基因数值化表达的二代从头组装方法及系统，完成二代测序的低成本高性能组装。相较于以往的组装算法，本方法没有将reads拆分为多段kmer，通过序列间overlap构建de bruijn graph有向图，而是通过无监督机器学习等方法将测序碱基数值化，将reads之间overlap的排列顺序转换为数值快排顺序，不仅提高了reads排列的速度，还一定程度上避免了组装过程中可能产生kmer间gap位置的产生，减少反复使用reads回帖比对的次数，在保证准确率的基础上大大地提高了运算速度。在组装scaffold过程中，没有采用read重新回帖比对增加算力成本和io开销的方法，而是采用酶切位点辅助组装的方法，提取长片段序列的酶切位点特征，将基因文件进行嵌合，组装为更长的片段。在组装过程中，同时完成了snp变异信息的检测。使用cuda异构编程，采取gpu加速的策略，协调cpu、gpu共同运算，提高并行速度。总而言之，本方法是一种集高精准度、低算力需求、高效率于一体的二代组装方案，大幅提高了基因组装分析的生物信息学分析流程的速度，有助于生信分析人员以更快的速度、更短的时间获得准确的分析结果。
[0056]
依托于基因测序技术以及生物信息学等科学领域的快速发展，本方法为低成本、快速高效的基因组组装分析提供了有力保障，在真正意义上实现了以计算机科学技术为工具，对生物信息进行储存、检索与分析。同时，生物信息学是21世纪自然科学的核心领域之一，本方法作为生物信息学中dna-seq分析的关键技术，在基因组学中具有极高技术优势和市场价值，适合于基因分析的实际应用与技术上的进一步推广。
[0057]
本发明的核心突破点至少包括：1、彻底破除了生物信息学组装的内存过大难题。2、实现了完全遍历组装，极大程度降低了组装错误风险。3、运用异构计算实现了计算效率的极大提升。
[0058]
需要说明的是，对于前述的实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本技术并不受所描述的动作顺序的限制，因为依据本技术，某一些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例属于优选实施例，所涉及的动作并不一定是本技术所必须的。
[0059]
上述实施例中，描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。