CN115440302A - 基因组叠阵、基因组架构、基因组序列组装方法及系统 - Google Patents

基因组叠阵、基因组架构、基因组序列组装方法及系统 Download PDF

Info

Publication number
CN115440302A
CN115440302A CN202110620095.7A CN202110620095A CN115440302A CN 115440302 A CN115440302 A CN 115440302A CN 202110620095 A CN202110620095 A CN 202110620095A CN 115440302 A CN115440302 A CN 115440302A
Authority
CN
China
Prior art keywords
sequence
sequencing
array
genome
stacked
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.)
Pending
Application number
CN202110620095.7A
Other languages
English (en)
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.)
Academy of Mathematics and Systems Science of CAS
Original Assignee
Academy of Mathematics and Systems Science of CAS
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 Academy of Mathematics and Systems Science of CAS filed Critical Academy of Mathematics and Systems Science of CAS
Priority to CN202110620095.7A priority Critical patent/CN115440302A/zh
Publication of CN115440302A publication Critical patent/CN115440302A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/20—Sequence assembly
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/30—Detection of binding sites or motifs

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Biotechnology (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Theoretical Computer Science (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Molecular Biology (AREA)
  • Genetics & Genomics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明提供一种基因组叠阵组装方法及系统、一种基因组架构组装方法及系统以及一种基因组序列组装方法及系统，属于生物信息技术领域。本发明通过稳健回归技术优化一个全局损失函数来确定测序序列之间的叠落关系，消除假阳性比对的干扰，从而得到更准确的组装叠阵集，避免组装基因组上拷贝数的缺失或者错误拼接；基于回归矩阵计算进行叠阵排列的估计，对异常值的容忍度更高，除非歧义的样本信息数量多于真实样本的数量，不会出现“崩溃”的情形；通过重抽样技术对测序数据进行多次独立地采样、组装、改进、评估，然后整合成最终的组装基因组，可以减少由测序数据中的噪音给组装结果带来的不确定性，降低组装结果对组装参数选择的敏感性。

Description

基因组叠阵、基因组架构、基因组序列组装方法及系统
技术领域
本发明涉及生物信息技术领域，具体地涉及一种基因组叠阵组装方法、一种基因组架构组装方法、一种基因组序列组装方法、一种基因组叠阵组装系统、一种基因组架构组装系统以及一种基因组序列的组装系统。
背景技术
基因组测序是开展分子生物学研究的重要技术。研究人员可以通过基因组测序来获得一个生命体的基因组碱基序列，它作为这个生命体的遗传序列模板，为基因的转录、调控、修饰等层面的研究提供了重要参照，从而帮助解释生命现象背后的分子机制。通过将被测物种或者个体之间的基因组进行比较，研究人员可以发现它们在基因组水平上的差异，从而为揭示物种演化历史，改良培育作物、诊断遗传疾病、优化药物治疗等提供指导。
目前应用比较广泛的测序技术是第二代和第三代测序技术。第二代测序技术碱基辨识质量高，价格便宜，所测到的序列长度(也就是碱基对数目)通常在100-300bp之间。它的一个重要特点就是可以从两端对一个很长的片段进行测序，得到这个长片段两端的碱基序列，也就是双末端测序序列。第三代测序技术能够测量的序列长度大幅增加(通常在10-100kbp之间)，但同时伴随着测序成本和碱基辨识错误率的增加。
将测序序列组装成基因组是计算生物学领域的基本问题。因为测序仪所能测量的序列长度远小于基因组长度，所以需要在测序后对序列进行组装，推断出它们的相对位置，从而还原出被测基因组。基因组组装面临的一个关键挑战是基因组上有很多相似度很高、或是重复出现的区段，它们的存在给推测测序序列的相对位置增加了很大的不确定性。
现有的基因组组装方法在原理上主要分为两类。一类是基于De Bruijn图的方法，该方法的主要操作是：对于每一个测序序列，每隔一个碱基切割出一个特定长度的子序列(通常称为k-mer，k代表子序列的长度)；利用所有被切割出来的子序列构造De Bruijn图；进行一定的纠错操作后，在图上寻找路径，每条路径被推断为被测基因组上的片段。另一类方法是基于测序序列叠落关系的，该方法对每两条序列进行比对，然后根据比对结果推断序列的叠落关系，进而得到组装叠阵集，其中每个叠阵对应被测基因组上的一个片段。
得到若干个叠阵(或者片段)以后，需要将他们进行定向和排列，组装成架构。一般的方法是将测序序列与所有的叠阵一致序列进行比对，利用双末端测序序列的映射和库长信息，或者单分子长读长测序序列映射到多个叠阵一致序列的坐标信息，来确定不同叠阵之间的顺序、方向关系和距离范围，从而得到组装基因组架构。
基于De Bruijn图的方法对重复度低的基因组组装效果会比较好，而对于重复度高的基因组不是很理想，因为所切割出的子序列长度明显短于测序序列，所以特异性会降低，导致在图上寻找路径时出现错误。
基于测序序列叠落关系的方法以测序序列为单位，而不是k-mer，所以保留了更多基因组上的位置信息，更有利于重复序列的还原。但是目前的方法主要基于贪婪策略，通过依次选择局部最优的比对来确定序列之间的叠落关系，这类方法容易受到假阳性比对的干扰，得到错误的叠落关系，从而导致组装基因组上拷贝数的缺失，或者错误拼接。但是如果通过使用更严格的序列比对准则来降低比对的假阳性，就会造成比对假阴性的增加，从而降低组装基因组的连续性。
由于双末端测序的库长有时会出现较大偏差，同时，基因组上的重复序列、杂合区域会使叠阵的排列信息产生歧义，现有的基因组架构方法在确定叠阵之间的方向、相对位置时难以得到兼顾长度和准确度的基因组架构。
由于样品制备和碱基辨识过程中产生的偏差，测序序列中可能会存在与被测基因组不一致的错误，从而导致组装结果不准确。对于同一个基因组组装方法和同一个被测基因组，输入不同的测序数据集往往会输出不同的组装结果，而现有的基因组组装系统通常将所有测序序列作为输入，无法消除由测序数据给组装结果带来的不确定性，并且组装结果往往也依赖于组装参数的选择，而参数的调试通常是基于经验的。
发明内容
本发明实施方式的目的是提供一种基因组叠阵组装方法和系统，以消除假阳性比对对于确定测序序列叠落关系的误导，从而增加组装基因组的准确性和连续性。
为了实现上述目的，本发明提供一种基因组叠阵组装方法，该方法首先从样品测序序列中随机抽取一定深度的第一测序序列集；然后通过预定的序列比对算法将第一测序序列集中的测序序列进行两两比对，得到测序序列的两两比对信息；根据测序序列的两两比对信息将所有满足叠落条件的比对共同表示为第一线性回归模型；通过迭代重加权最小二乘算法对所述第一线性回归模型求解，得到全局损失函数最小的解和第一测序序列集中的测序序列在被测基因组上的坐标的稳健估计；根据所述全局损失函数最小的解构造的无向图将所述测序序列划分为多个连通分量，按照各个连通分量内的测序序列的坐标的稳健估计对测序序列进行排列，生成叠阵；根据测序序列的两两比对信息对所述叠阵进行合并，得到初步组装叠阵集；获取所述初步组装叠阵集中各个叠阵的一致序列，将所述测序序列比对映射到所述一致序列，得到所述测序序列向一致序列的映射结果，所述映射结果为新的组装叠阵集。
本发明还提供一种基因组叠阵组装系统，该系统应用前述的基因组叠阵组装方法进行基因组叠阵组装。
本发明实施方式的另一个目的是提供一种基因组架构组装方法及系统，以消除叠阵排列的歧义信息对架构组装的干扰，从而兼顾组装架构的长度和准确度。
为了实现上述目的，本发明提供一种基因组架构组装方法，该方法首先从样品测序序列中随机抽取一定深度的第二测序序列集；然后根据所述第二测序序列集组装基因组叠阵，得到组装叠阵集
Figure BDA0003099519500000031
将所述第二测序序列集与组装叠阵集
Figure BDA0003099519500000032
的一致序列集进行比对，根据比对结果筛选出单映射测序序列叠阵，得到唯一映射的测序序列信息和筛选后的一致序列集；基于唯一映射的测序序列信息和筛选后的一致序列集构建无向有权叠阵图，将所述无向有权叠阵图划分为多个连通分量；采用组合优化模型确定各个连通分量上的每个顶点所对应的叠阵的方向；建立回归模型，估计叠阵在被测基因组上的排列坐标，根据各个顶点对应的叠阵的方向和所述排列坐标对同一连通分量中的所有叠阵进行排序，得到该连通分量的准架构，根据多个连通分量的准架构得到架构集；对所述架构集中每个架构的叠阵排列坐标进行梳理，消除坐标重叠或者冗余的叠阵，得到组装基因组架构。
本发明还提供一种基因组架构组装系统，该系统应用前述的基因组架构组装方法进行基因组架构组装。
本发明实施方式的再一个目的是提供一种基因组序列组装方法及系统，以减少由测序数据中的噪音给组装结果带来的不确定性，并且降低组装结果对组装参数选择的敏感性。
为了实现上述目的，本发明提供一种基因组序列组装方法，该方法首先从样品测序序列中随机抽取一定深度的第三测序序列集和另一深度的第四测序序列集；采用所述的基因组叠阵组装方法对第三测序序列集进行处理得到组装叠阵集；采用所述的基因组架构组装方法对组装叠阵集进行处理得到组装基因组架构；采用第四测序序列集更新改进所述组装基因组架构，并对改进后的组装基因组进行质量评估，若质量达到标准，则将改进后的组装基因组储存为一条候选组装基因组；重复前述步骤，直到得到预定数量的候选组装基因组，将所有候选组装基因组进行多重序列比对，将比对的每一列中频率最大的碱基连接起来得到最终的组装基因组。
本发明还提供一种基因组序列组装系统，该系统应用前述的基因组序列组装方法进行基因组序列组装。
通过上述技术方案，能够消除假阳性比对对于确定测序序列叠落关系的误导，从而增加组装基因组的准确性和连续性，同时消除叠阵排列的歧义信息对架构组装的干扰，从而兼顾组装架构的长度和准确度；减少由测序数据中的噪音给组装结果带来的不确定性，并且降低组装结果对组装参数选择的敏感性。
本发明实施方式的其它特征和优点将在随后的具体实施方式部分予以详细说明。
附图说明
附图是用来提供对本发明实施方式的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明实施方式，但并不构成对本发明实施方式的限制。在附图中：
图1是根据本发明实施例的利用测序序列两两比对信息建立线性回归模型的示意图；
图2是根据本发明实施例的基于测序序列两两比对信息将回归结果导出的连通分量进行必要合并的示意图；
图3是根据本发明实施例的将测序序列重新映射到一致序列的示意图；
图4是根据本发明实施例的基因组架构的组装方法的总体流程图；
图5是根据本发明实施例的从双末端测序序列映射数据中提取叠阵一致序列两两之间方向关系和起点距离信息的示意图；
图6是根据本发明实施例的从单分子长读长测序序列映射数据中提取叠阵一致序列两两之间方向关系和起点距离信息的示意图；
图7是根据本发明实施例的在叠阵图的一个连通分量上用组合优化模型优化各顶点方向指派的算法流程图；
图8a-图8b是根据本发明实施例的计算准架构两两之间方向关系和距离信息的八种情况的示意图；
图9是根据本发明实施例的利用回归估计坐标、连接信息以及重叠区域的局部比对梳理最终的架构排列的示意图；
图10是根据本发明实施例的基于重抽样的基因组序列组装系统的框图。
具体实施方式
以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。
为了便于理解本发明实施例的基因组序列的组装方法，下面先对本发明实施例的相关技术名词作出如下定义和解释：
基因组：包含一个生物体所有遗传信息的遗传序列，由四种碱基排列而成，通常表示为含有A，C，G，T四种字符的序列，每一种字符代表一种碱基。
被测基因组：需要确定碱基序列的基因组。
测序序列：由测序平台输出的某一生物体基因组子序列的测量结果，由若干个代表四种碱基的字符(A，C，G，T)组成，代表基因组的一个子序列。测序序列的一端为5’端，另一端为3’端；通常在测序平台给出的文件中，测序序列的左端为5’端，右端为3’端。
bp：英语basepair的缩写。100bp指一个碱基序列有100个碱基。1kbp指一个碱基序列有1000个碱基。
深度：测序序列的总碱基数与被测基因组大小的比值，也就是被测基因组上单个碱基被读取的平均次数。比如某样品的测序深度为30X，就是说该样品的基因组上每个单碱基平均被读取了30次。
第二代测序技术：也称作新一代测序技术，英文名称为Next GenerationSequencing，英文缩写为NGS；其特点是能并行地产生大量的双末端测序序列数据，每一条测序序列的长度大多为100-300bp。
双末端测序序列：第二代测序技术在进行操作时，将被测的序列切割成若干片段，片段的长度可小可大；之后在片段的两端进行测序，在一端得到一个测序序列，在另一端得到一个测序序列；如果片段过长，中间部分一般无法测到；所得到的两个测序序列对应于同一个片段，称为双末端测序序列。
片段库、库长：在测序平台给出的同一个双末端测序序列文件中，所有双末端测序序列所属的片段长度相近，测序平台会估计出一个平均值；在下文中，称测序平台给出的每一个双末端测序序列文件为片段库，称其估计出的片段长度平均值为库长。
短库、长库：在双末端测序过程中，如果被测的序列片段先被环化然后打断测序，得到的片段库称为长库(英文名称为mate-paired library)，库长通常在2-10kbp，且测序方向为从片段内部向两头；如果被测的序列片段未经过环化，得到的片段库称为短库(英文名称为paired-end library)，库长通常在200-500bp，且测序方向为从片段两头向内部。
第三代测序技术：指的是单分子测序技术，英文名称为Third GenerationSequencing或Single Molecule Sequencing或Long-read Sequencing；目前主要包括两大技术方向：单分子荧光测序和纳米孔测序；其特点是能产生长读长的测序序列数据，每一条测序序列的长度大约在10-100kbp之间，称为单分子长读长测序序列，但测序错误率较高，一般在10％左右。
质量值：在测序过程中，测序序列的每一个碱基都有被误读的可能性，质量值是这一可能性的度量。测序序列的每一个碱基都对应一个这样的质量值，质量值越高，被误读的可能性越低。
插入：指相对于被测基因组，测序序列的某两个相邻的碱基之间额外插入了一段碱基序列。
删失：指相对于被测基因组，测序序列丢失一段或若干段碱基序列。
INDEL：一个插入或者删失。
碱基替换：测序序列被映射到参考基因组上后某些碱基和参考基因组上对应的碱基不相同。
测序序列映射：对于一个测序序列，在参考基因组上寻找一个与其长度大体一致的子序列，该子序列与测序序列完全一样，或者存在差异，但差异在预先设定的标准之内。通常以碱基替换个数，INDEL长度来衡量二者的差异。
映射结果：对于一个测序序列，如果参考基因组上存在子序列，且子序列与测序序列的差异在设定的标准之内，则该测序序列映射成功。映射结果包括：(1)参考基因组上的子序列最右端碱基的坐标，该值作为测序序列在参考基因组上的映射坐标；(2)测序序列的映射方向，可以按5’至3’的方向映射，也可以按3’至5’的方向映射；(3)测序序列与参考基因组上子序列的比对信息，包括两者之间存在的碱基替换，INDEL情况。当参考基因组上存在多个与测序序列差异在设定标准之内的子序列时，该测序序列具有多个映射结果，每一个映射结果对应于参考基因组上的一个子序列。
基因组组装：利用测序序列还原被测基因组，组装后得到的基因组称为组装基因组。
测序序列叠阵(可以简称为叠阵)：英文术语为contig或layout；表示由若干测序序列形成的集合，在该集合当中，每两条测序序列之间的叠落关系或者相对距离是确定的。测序序列叠阵可以用含有A，C，G，T，Ф五种字符的矩阵表示；矩阵的某个元素为Ф表示矩阵的该位置上不存在碱基；矩阵的每一列对应于叠阵的每一个位点且每一列至少含有一个非Ф字符；矩阵的列数为叠阵的长度；矩阵的每一行代表一个测序序列，每行最后一个非Ф字符所在列的列标为该行所代表的测序序列在叠阵中的坐标。
叠阵集：由叠阵组成的集合。
一致序列(或叠阵一致序列)：英文术语为consensus；对于叠阵的每一列，根据其包含的所有碱基信息推断出一个碱基，并将所推断出的碱基连接后得到的序列；推断碱基的一个可行的方法为，对叠阵的每一列，在其包含的所有碱基中取频率最大者。
一致序列集(叠阵一致序列集)：由一致序列组成的集合。
叠阵图：以叠阵一致序列为顶点，基于测序序列向叠阵一致序列映射结果建立连接顶点的边，所构造的图。
架构(或组装基因组架构)：英文术语为scaffold；它由一系列按照一定顺序排列后的叠阵组成，其中所有叠阵的方向和相对位置关系已知，作为组装基因组的一个框架。
架构集：由架构组成的集合。
间隙：英文术语为Gap；架构中确定先后顺序的相邻叠阵之间未组装出来的空白区域，在组装结果中一般以N填充，N的数目代表了间隙长度的估计。
本发明以双末端测序序列或者单分子长读长测序序列作为输入，其中，双末端测序序列和单分子长读长测序序列来自被测基因组，双末端测序序列可以包含多个具有不同库长的片段库。
本发明一方面提供一种基因组叠阵组装方法，所述方法包括：
S01：从输入的样品测序序列中抽取一定深度的第一测序序列集，深度为d1，第一测序序列集中的序列记为
Figure BDA0003099519500000081
S02：通过序列比对算法对所述第一测序序列集中的测序序列进行两两比对，得到测序序列的两两比对信息，比对前可以根据测序序列的长度和质量分布情况，以及对比对灵敏度、特异度的要求设置比对所用的参数和判断比对得分是否显著的标准。
S03：根据测序序列的两两比对信息将所有满足叠落条件的比对共同表示为第一线性回归模型。具体的：
S031：叠落条件设定为：比对得分达到显著性标准，并且叠落区域两端未比对上的子序列长度小于给定阈值；
S032：分别沿着被测样品DNA的两条单链从5’端到3’端建立两条坐标轴；第一测序序列集中的任意一条测序序列Ri在其所在的DNA单链上的坐标定义为Ri的最右端碱基的坐标，记为βi；同样可以定义Ri的反向互补序列(记为
Figure BDA0003099519500000082
)在互补单链上的坐标
Figure BDA0003099519500000083
Figure BDA0003099519500000084
为待估参数；对于任意两条满足叠落条件的测序序列Ri和Rj，假定测序序列Ri的最右端碱基对齐到测序序列Rj内部，Rj右侧悬垂部分含有yi,j个碱基，由此得到一个观测样本：yi,j＝βj-βi+εi,j；其中，εi,j表示观测样本的观测误差，βj是Rj在其所在的DNA单链上的坐标；
S032：假定第二测序序列Rj的最左端碱基对齐到第一测序序列Ri内部，Ri左侧悬垂部分含有
Figure BDA0003099519500000085
个碱基，由此得到一个对偶观测样本：
Figure BDA0003099519500000086
其中，
Figure BDA0003099519500000087
表示对偶观测样本的观测误差，
Figure BDA0003099519500000088
是Rj的反向互补序列在互补单链上的坐标；
S033：根据测序序列的两两比对信息将所有满足叠落条件的比对的观测样本和对偶观测样本整合为一个矩阵方程：Y＝Xβ+ε；其中，每一行为一个观测样本，向量
Figure BDA0003099519500000089
为待估参数，X是一个大型稀疏矩阵，X的每一行只有两个非零元素1和-1，X储存为稀疏矩阵格式。
需要说明的是，在本实施例中，观测误差εi,j具体是指Rj右侧悬垂部分相对于被测基因组插入的碱基数减去删失的碱基数；同理，
Figure BDA0003099519500000091
表示对偶观测样本的观测误差，具体是指Ri左侧悬垂部分相对于被测基因组插入的碱基数减去删失的碱基数；矩阵方程中的ε表示向量(ε1,...,εm)T，其中m表示观测样本和对偶观测样本的总数。
在其他一些实施例中，可以将测序序列的坐标定义为该测序序列最左端碱基的坐标，则观测样本表示为：yi,j＝βj-βi+εi,j；式中yi,j表示Ri左侧悬垂部分的碱基数目，εi,j表示Ri左侧悬垂部分相对于被测基因组插入的碱基数减去删失的碱基数。
图1为根据测序序列比对信息建立线性回归模型的示意图，如图1所示，每两条满足叠落条件的测序序列给出模型的一个观测样本，同时这两条序列的反向互补序列给出模型的另一个对偶观测样本。
S04：通过迭代重加权最小二乘算法对所述第一线性回归模型求解，得到全局损失函数最小的解和所述第一测序序列集中的测序序列在被测基因组上的坐标的稳健估计，包括：
S041：利用大型稀疏线性方程组的加速求解算法(例如LGMRES算法)求解关于β的稀疏线性方程组：XTX·β＝XTY，得到的解作为
Figure BDA0003099519500000092
的初始值；
S042：在每一轮迭代中，根据当前的
Figure BDA0003099519500000093
计算残差向量
Figure BDA0003099519500000094
以及
Figure BDA0003099519500000095
各分量绝对值的最大值
Figure BDA0003099519500000096
S043：如果
Figure BDA0003099519500000097
则计算权重矩阵
Figure BDA0003099519500000098
其中，α为预设的收敛阈值，测序序列的碱基辨识错误率越高，给α设定的值越大，反之则越小；权重函数为
Figure BDA0003099519500000099
S044：根据计算得到的W，取互为对偶的两个样本的权重值中的较小值赋予这两个样本，即：假定互为对偶的两个样本的权重分别为a1和a2，其中a1＞a2，则将a2作为这两个样本的权重；
S045：利用大型稀疏线性方程组的加速求解算法(例如LGMRES算法)求解关于β的稀疏线性方程组：XTWX·β＝XTWY；
S046：用得到的解更新
Figure BDA0003099519500000101
删除矩阵X和Y中对应权重为0的行，返回进行下一轮迭代；
S047：直到
Figure BDA0003099519500000102
停止迭代，记录最后一轮迭代的矩阵X；
S048：对最后一轮
Figure BDA0003099519500000103
的每个分量取整后输出
Figure BDA0003099519500000104
即为所述第一测序序列集中的测序序列在被测基因组上的坐标的稳健估计。
其中，函数ρ可以取为任意一个同时满足如下3个条件的函数，且在实际应用中取得良好效果：
(1)ρ(x)≥0对所有χ成立，且ρ(x)在χ＝0处取到最小值；
(2)ρ(x)关于Y轴对称；
(3)χ≥0时，ρ(x)关于χ单调非降，且χ增大时ρ(x)不会变得太大。
在一些实施例中，函数ρ采用双权损失函数ρB：
Figure BDA0003099519500000105
对应的权重函数为：
Figure BDA0003099519500000106
在本发明中，参数c应当是关于
Figure BDA0003099519500000107
的单调递增函数，且
Figure BDA0003099519500000108
在一些实施例中，经测试效果理想的取法的函数为如下的阶梯函数：
Figure BDA0003099519500000109
S05：根据所述全局损失函数最小的解构造的无向图将所述测序序列划分为多个连通分量，按照各个连通分量内的测序序列的坐标的稳健估计对测序序列进行排列，生成叠阵，包括：
S051：以第一测序序列集中的测序序列及所述测序序列对应的反向互补序列作为顶点，共包含2n1个顶点，分别表示测序序列
Figure BDA0003099519500000111
根据最后一轮迭代的矩阵X中每一行的两个非零元素确定两条测序序列，将这两条测序序列对应的顶点之间建立一条边，构造得到无向图；
S052：遍历所述无向图，将所述无向图划分为多个连通分量，得到2k个连通分量，记为C1，...，Ck，
Figure BDA0003099519500000112
其中Ci与
Figure BDA0003099519500000113
所包含的测序序列一一对应，互为反向互补序列；如果某个连通分量中同时包括某个Ri和
Figure BDA0003099519500000114
则将这个连通分量去除；
S053：将各个连通分量中的测序序列根据其在被测基因组上的坐标的稳健估计进行排列，得到叠阵集：L1，...，Lk，
Figure BDA0003099519500000115
叠阵Li和
Figure BDA0003099519500000116
分别对应于DNA双链上互为反向互补的两条子序列。
需要说明的是，一般情况下测序序列是由左侧开始按照从小到大的顺序进行排列，在其他一些实施例中，测序序列可以从右侧开始，从大到小进行排列。以连通分量C1为例，将其中的测序序列按照坐标估计从大到小的顺序依次排开，对于任意两个相邻的测序序列，他们右端相隔的碱基数目等于他们坐标估计的差值，得到的排列即为叠阵L1。
S06：根据测序序列的两两比对信息对所述叠阵进行合并，得到初步组装叠阵集，包括：
S061：以所述叠阵集中的叠阵作为顶点共包括2k个顶点，分别表示叠阵L1，...，Lk，
Figure BDA0003099519500000117
S062：对于任意两个叠阵，不失一般性假设为Li和Lj，若Li3’端s个测序序列分别和Lj5’端s个测序序列具有大于t个满足叠落条件的比对，则建立一条Li对应顶点指向Lj对应顶点的有向边，构建有向图；其中，边的权重设定为满足叠落条件的比对数，s和t的取值正相关于抽样深度d1；
S063：从所述有向图中包含最多测序序列的叠阵Lx开始，确定从叠阵Lx出发的权重最大的边所指向的叠阵Ly；
S064：从Lx的3’端选取一条测序序列Ri，和从Ly的5’端选取一条测序序列Rj，Ri和Rj之间的比对满足叠落条件；
S065：假定Ri的最右端碱基对齐到Rj内部，Rj右侧悬垂部分含有Δ个碱基，则将Ly中所有测序序列的坐标的稳健估计都加上
Figure BDA0003099519500000121
假定Ly中的任意一条测序序列Ru，其坐标估计由
Figure BDA0003099519500000122
变为
Figure BDA0003099519500000123
S067：将Lx和Ly中的测序序列按照其坐标的稳健估计从小到大进行排列合并成一个叠阵，同时将Lx和
Figure BDA0003099519500000124
在有向图中对应的顶点去除；
S068：从Ly开始重复步骤S063-S067，直到没有可以合并的叠阵为止，再从有向图中剩下的顶点中包含测序序列最多的叠阵开始重复步骤S063-S067，直到遍历有向图中所有顶点为止，得到的合并后的叠阵作为初步组装叠阵集，记为D1,...,Dm。
图2为基于测序序列比对信息将叠阵进行合并的示意图，如图2所示，三个叠阵通过5’端和3’端序列之间的比对所提供的连接信息合并成一个叠阵。
S07：获取所述初步组装叠阵集中各个叠阵的一致序列，将所述测序序列比对映射到所述一致序列，得到所述测序序列向一致序列的映射结果，所述映射结果为新的组装叠阵集，包括：
S071：对于所述初步组装叠阵集中的任意一个叠阵Di，取每一列中频率最大的碱基生成叠阵Di对应的一致序列，记为Ei；
S072：对于叠阵Di中的任意一条测序序列Rj的最左端碱基在Di中的坐标记为a1，Rj的最右端碱基在Di中的坐标记为a2，将Ei中第a1-τ个碱基到第a2+τ个碱基之间的子序列取出与Rj进行局部序列比对；
S073：将所述局部序列比对一致的部分作为Rj向Ei的映射结果，Di中的所有测序序列向Ei的映射结果构成了一个新的叠阵
Figure BDA0003099519500000125
Figure BDA0003099519500000126
为更新后的组装叠阵集。
图3为将测序序列重新映射到一致序列的示意图，如图3所示，对于叠阵中任意一条测序序列，找到其对应在一致序列中的位置，将一致序列在所述位置处的子序列连同其两边预定长度的碱基序列取出与该测序序列进行比对。
本发明第二方面提供一种基因组叠阵组装系统，该系统应用所述的基因组叠阵组装方法进行基因组叠阵组装。
本实施例提供的基因组叠阵组装方法通过稳健回归技术优化一个全局损失函数来确定测序序列之间的叠落关系，可以消除假阳性比对的干扰，从而得到更准确的组装叠阵集，避免了组装基因组上拷贝数的缺失，或者错误拼接；优化算法基于求解大型稀疏线性方程组，有充足的理论研究来支持对于计算时间和存储空间的节约利用，并且可以通过GPU或者分布式计算实现加速；基于回归矩阵计算而不是图的计算进行叠阵排列的估计，对异常值的容忍度更高，除非歧义的样本信息数量多于真实样本的数量，不会出现“崩溃”的情形：即估计排列值与真实排列相差甚远的情形。
本发明还提供一种基因组架构组装方法，如图4所示，所述方法包括：
S11：从输入的样品测序序列中抽取一定深度的第二测序序列集，深度为d1，第二测序序列集中的序列记为
Figure BDA0003099519500000131
S12：根据所述第二测序序列集组装基因组叠阵，得到组装叠阵集
Figure BDA0003099519500000132
S13：将所述第二测序序列集R1,R2,...,Rn1与组装叠阵集
Figure BDA0003099519500000133
的一致序列集进行比对，根据比对结果筛选出单映射测序序列叠阵，根据设定的条件对所述单映射的测序序列及其映射上的叠阵进行筛选，得到唯一映射的测序序列信息和筛选后的一致序列集。具体的，
若输入的测序序列为双末端测序序列，则：
S1301：从比对结果中挑选出双末端序列的两端映射到同一个叠阵一致序列上的映射结果，以估计双末端测序片段库的库长，所述映射结果包含每一端的序列比对到叠阵一致序列上的起始点坐标，分别记为s1,t1和s2,t2，库长的观测值记为：in＝max(s1,t1,s2,t2)-min(s1,t1,s2,t2)；对每个片段库分别预定数量的观测值，取观测值的中位数M作为库长的观测，并计算其标准差σ，在本发明的一个实施例中，每个片段库分别取30000个观测值，再取这30000个观测值的中位数作为库长的观测；
S1302：将分别映射到两个不同的叠阵一致序列上的双末端测序序列的映射坐标仍分别记为s1,t1和s2,t2；
S1303：根据映射坐标判断所述双末端测序序列与其映射上的叠阵一致序列的方向是否一致：由二代测序原理决定短库的测序方向为从两头向内部，而长库的测序方向是由内部向两头，我们根据两端序列起始点映射坐标的大小与上述方向是否一致可以知道两端序列与其映射上的叠阵一致序列的方向是否一致(短库左端序列起点坐标小于终点坐标，右端序列起点坐标大于终点坐标，则说明方向一致，长库则正好相反)；根据双末端测序序列与其映射上的叠阵一致序列的方向判断结果判断两个不同的叠阵一致序列的方向是否一致；根据两个不同的叠阵一致序列的方向判断结果和所述库长计算得到两个叠阵起始点之间的距离；同时筛除其中一端映射到叠阵上的位置超出预设的插入长度范围的连接，得到测序序列到两两叠阵一致序列的连接信息，图5为从双末端序列映射数据中提取叠阵一致序列两两之间方向关系和起点距离信息的示意图，如图5所示，具体包括：
定义所述双末端测序序列的左端序列映射上的叠阵方向为正，且位置在前，假定左右端映射上的叠阵起点分别为F1,F2，则：
a)若s1＜t1且s2＜t2，则两个叠阵一致序列的方向各自与双端测序序列一致，即两叠阵的方向相反，则距离为：F2-F1＝in+s1+s2，筛除len1-t1＞in+3σ或者len2-t2＞in+3σ的连接；
b)若s1＜t1且s2＞t2，则两个叠阵一致序列的方向各自与双端测序序列相反，即两叠阵的方向一致，则距离为：F2-F1＝in+s1-s2，筛除len1-t1＞in+3σ或者s2＞in+3σ的连接；
c)若s1＞t1且s2＞t2，则两个叠阵一致序列的方向各自与双端测序序列一致，即两叠阵的方向相反，则距离为：F2-F1＝in-s1-s2，筛除s1＞in+3σ或者s2＞in+3σ的连接；
d)若s1＞t1且s2＜t2，则两个叠阵一致序列的方向各自与双端测序序列相反，即两叠阵的方向一致，则距离为：F2-F1＝in-s1+s2，筛除s1＞in+3σ或者len2-t2＞in+3σ的连接。
若输入的测序序列为单分子长读长测序序列，则考虑映射到多个叠阵一致序列的测序序列，首先假定比对到多个叠阵一致序列上的比对区间在测序序列上的坐标按起点排序依次为[qs1,qe1],[qs2,qe2],...,[qsk,qek]，其对应的叠阵一致序列上的坐标依次为[cs1,ce1],[cs2,ce2],...,[csk,cek]；根据上述坐标得到任意两个比对到同一测序序列上的叠阵一致序列的方向关系和起点距离作为两两叠阵一致序列的连接信息；图6为从单分子长读长测序序列映射数据中提取叠阵一致序列两两之间方向关系和起点距离信息的示意图。
其中，若叠阵与测序序列方向一致，则：Fi＝qsi-csi；若叠阵与测序序列方向相反，则：Fi＝qsi-csi+leni；距离为：Fj-Fi，其中，Fi和Fj测序序列映射上的叠阵序列的起点；最后筛除同一测序序列上相互交叠的连接。
S14：基于唯一映射的测序序列信息和筛选后的一致序列集构建无向有权叠阵图，将所述无向有权叠阵图划分为多个连通分量，包括：
S1401：定义叠阵之间形成边的连接数阈值参数m，阈值参数的设置主要取决于测序数据的深度，可设计为测序数据总深度的一半，最低不应低于3；
S1402：以组装叠阵集的一致序列集作为顶点，单映射测序序列叠阵之间的连接数大于阈值参数m的形成一条边，构建无向有权叠阵图；这条对应着n个连接的边包含了这两个顶点之间的方向信息和距离信息，有一部分连接指示两个顶点方向相同，另一部分则指示其方向相反；
S1403：将无向有权叠阵图中顶点i和顶点j之间的同向连接数和异向连接数分别记为aij和bij，对应的距离信息分别录入同向和异向的叠阵距离列表；
S1404：将所述无向有权叠阵图划分为多个连通分量。
S15：采用组合优化模型确定各个连通分量上的每个顶点所对应的叠阵的方向，包括：
S1501：以aij和bij中的极大值给连通分量上的每条边赋予权重并标记该边的方向，若aij≥bij，则权重wij＝aij，Dij＝1；否则权重wij＝bij，Dij＝-1；遍历所有顶点得到连通分量的每个顶点的初始方向；
S1502：迭代优化当前连通分量上的矛盾方向连接总数，每轮迭代中分别计算改变每个顶点，导致目标函数的单点变化量，若存在某顶点的单点变化量为负，就选取单点变化量与总连接数比值最大的顶点改变其方向，更新方向后继续下一轮迭代，直到改变任一顶点都无法使得总矛盾数更小，停止迭代，得到每个顶点所对应的叠阵的方向。图7示出了在叠阵图的一个连通分量上用组合优化模型优化各顶点方向指派的算法流程。
在一些实施例中，采用权重递降的深度优先遍历方法遍历所有顶点得到连通分量的每个顶点的初始方向：从零顶点(方向设为1)出发，总是在访问当前顶点(访问的含义：从顶点i访问顶点j，令Dj＝DijDj)后选择访问当前顶点权重最大的邻居，依次遍历得到所有顶点的初始方向。
在其他一些实施例中，采用最大生成树方法遍历所有顶点得到连通分量的每个顶点的初始方向：先找到当前连通分量的一个最大权重生成树，以生成树的边来遍历连通分量，得到所有顶点的初始方向。
S16：建立回归模型估计叠阵在被测基因组上的排列坐标，根据各个顶点对应的叠阵的方向和所述排列坐标对同一连通分量中的所有叠阵进行排序，得到该连通分量的准架构，根据多个连通分量的准架构得到架构集，包括：
S1601：建立回归模型；
S1602：根据所述两两叠阵一致序列的连接信息以及每个顶点所对应的叠阵的方向，确定每条边的类型，类型包括同向边和异向边；
S1603：根据每条边的类型，取出每条边对应的叠阵距离列表的中位数作为回归模型的观测值；
S1604：设置加权截断最小二乘估计的残差允许阈值max Error，为了保证稳健回归的准确性，我们将此阈值尽量设小，但考虑到库长方差和迭代效率，设置为100；
S1605：以库长的方差估计的倒数与当前边所包含的符合当前顶点方向的连接数的乘积作为初始权重，即
Figure BDA0003099519500000161
计算叠阵一致序列起点坐标的加权最小二乘估计
Figure BDA0003099519500000162
其中系数矩阵X除第一行为(1,0,...,0)，其余行表示对应顶点i与j之间的连边，故第i项为-1，第j项为1；
S1606：根据上一轮的估计计算所有样本残差，当残差最大值小于残差允许阈值max Error，停止迭代，得到可行估计
Figure BDA0003099519500000163
否则选取样本中残差小于(100-iter)％分位数的样本，其中，iter表示当前迭代次数；
判断选取的样本对应的系数矩阵是否列满秩：若不满秩，根据系数矩阵将当前子图进一步划分成若干个连通分量，在子图的每个连通分量上从步骤S1501重新开始执行；若满秩，则在选取的样本子集上再次进行加权最小二乘估计，得到
Figure BDA0003099519500000164
再重复步骤S1606，直到得到可行估计
Figure BDA0003099519500000171
S1607：根据每个顶点所对应的叠阵的方向将当前连通分量中的所有叠阵按照其对应的可行估计的大小进行排序，整理到同一个坐标系中，得到一条准架构，根据多个连通分量的准架构得到架构集。
S17：对所述架构集中每个架构的叠阵排列坐标进行梳理，消除坐标重叠或者冗余的叠阵，得到组装基因组架构，包括：
S1701：将架构中的所有叠阵一致序列按每个顶点所对应的叠阵的方向进行调整，将反向的叠阵的起点与终点坐标对换，使得架构排列中的所有叠阵起点坐标均在终点坐标之前，将叠阵按起点坐标由小到大排序；
S1702：从排序后的第一个叠阵开始，依次向后，对坐标上被当前选入的最后一个叠阵包含或交叉的叠阵进行局部序列比对，比对成功，则去除被包含的叠阵，或将被交叉的叠阵并入当前叠阵；选择与当前叠阵终点坐标紧邻又有连接信息的叠阵加入当前架构，留置与当前叠阵终点坐标紧邻但无连接信息的叠阵；其中，
1)若两个叠阵一致序列i与j的坐标满足：
0＜startj-endi≤2kbp，且endj＞endi；则两个叠阵紧邻；
2)若两个叠阵一致序列i与j的坐标满足：
starti＜startj＜endj＜endi；则i包含j；
3)若两个叠阵一致序列i与j的坐标满足：
starti＜startj＜endi＜endj；则两个叠阵交叉；
其中，start,end分别表示叠阵一致序列在准架构中的起点和终点坐标；
S1703：更新当前叠阵为最新加入当前架构的叠阵，重复步骤S1702，直至到达该架构上的最后一个叠阵，结束当前架构，输出到最终的组装架构文件；
S1704：在被留置的叠阵中重复步骤S1702-S1703，确定下一个架构排列并输出至最终的组装架构文件；重复此步骤直至所有叠阵被安排进某一个最终的架构之中，得到组装基因组架构、叠阵在架构中的位置和方向信息、紧邻叠阵之间间隙长度的估计及估计方差、由坐标交叉且比对成功而合并得到的新的叠阵一致序列。
图9为利用回归估计坐标、连接信息以及重叠区域的局部比对梳理最终的架构排列的示意图。
在其他一些实施例中，为了得到更长的架构，将所述准架构组成的准架构集代叠阵一致序列集，重复步骤S14-S16，得到架构集。在迭代过程中需要注意：
A.叠阵图中的连接直接来源于测序序列比对到两个不同叠阵一致序列上的信息，准架构图中，则需选取比对到两个不同准架构上的测序序列信息，即挑选出两个叠阵一致序列不在同一个准架构中的映射信息；
B.由A中连接的定义更新，两个准架构之间的连接指示其同向还是异向需补充考虑两端映射的叠阵一致序列在各自准架构中的方向，若两端映射的叠阵一致序列同向且在各自准架构中的方向相同，则两准架构的方向相同；若两端映射的叠阵一致序列异向且在各自准架构中的方向相反，则两准架构的方向也相同；否则两准架构的方向相反；由A中连接定义的不同，计算两个准架构序列起点距离观测值的公式也需分情况替换，共分为如下8种情况，其中start,end分别表示叠阵一致序列在准架构中的起点和终点坐标，diff为步骤S1303中计算得到的叠阵一致序列起点之间的距离：
a)若两端映射的叠阵一致序列异向且在各自准架构中的方向相反，前者正后者负，则：
Figure BDA0003099519500000181
b)若两端映射的叠阵一致序列异向且在各自准架构中的方向相同，均为正方向，则：
Figure BDA0003099519500000182
c)若两端映射的叠阵一致序列同向且在各自准架构中的方向相同，均为正方向，则：
Figure BDA0003099519500000183
d)若两端映射的叠阵一致序列同向且在各自准架构中的方向相反，前者正后者负，则：
Figure BDA0003099519500000184
e)若两端映射的叠阵一致序列同向且在各自准架构中的方向相同，均为负方向，则：
Figure BDA0003099519500000185
f)若两端映射的叠阵一致序列同向且在各自准架构中的方向相反，前者负后者正，则：
Figure BDA0003099519500000186
g)若两端映射的叠阵一致序列异向且在各自准架构中的方向相反，前者负后者正，则：
Figure BDA0003099519500000191
h)若两端映射的叠阵一致序列异向且在各自准架构中的方向相同，均为负方向，则：
Figure BDA0003099519500000192
图8a-图8b是计算准架构两两之间方向关系和距离信息的八种情况的示意图，如图8a-图8b所示分别考虑了双末端测序序列两端映射到的叠阵一致序列在各自准架构中的方向。
C.在参数选择上也要适当变化。设计准架构顶点之间形成边的连接数阈值参数m’为m的1.5倍，设计加权截断最小二乘估计的残差允许阈值max Error′为max Error的1.5倍。
本发明第四方面提供一种基因组架构组装系统，该系统应用所述的基因组架构组装方法进行基因组架构组装。
本发明提供的基因组架构组装方法采用分层迭代算法充分利用所有有价值的测序信息，并采取先严后松的策略，保证底层准架构的准确性基础上进行迭代延长；基于回归方法的叠阵排列计算还能给出相邻叠阵之间的间隙长度估计的误差，即给出的是一个区间估计，为后续间隙填补提供一个可靠的参考区间；能同时利用二代双末端测序和三代单分子长读长测序，也能只用其中一种测序数据，对平台的适应性强。
本发明还提供一种基因组序列组装方法，如图10所示，所述方法包括：
S21：抽取一定深度的第三测序序列集，用于根据前述的基因组叠阵组装方法生成组装叠阵集，所述组装叠阵集用于根据前述的基因组架构组装方法生成组装基因组架构。
S22：抽取另一深度的第四测序序列集更新改进所述组装基因组架构，并对改进后的组装基因组进行质量评估，若质量达到标准，则将改进后的组装基因组储存为一条候选组装基因组，包括：
S2201：将第四测序序列
Figure BDA0003099519500000193
映射到组装基因组，去除其中同时映射到多个位置的测序序列，形成单映测序序列叠阵；对于叠阵中的每一列，如果其被多于d0条单映序列覆盖，则取该列中频率最大的碱基更新组装基因组的相应位置的碱基，得到改进后的组装基因组；
S2202：将第四测序序列
Figure BDA0003099519500000201
映射到改进后的组装基因组，根据映射结果对其质量进行评估，若质量达到标准，则将其储存为一条候选组装基因组。
S23：重复步骤S21和S22，直到得到预定数量的候选组装基因组，将所有候选组装基因组进行多重序列比对，将比对的每一列中频率最大的碱基连接起来得到最终的组装基因组。
最终的组装基因组被存储为标准格式的文件。
本发明还提供了一种基因组序列组装系统，该系统应用前述的基因组序列组装方法进行基因组序列组装。
本发明基于回归思想实现基因组叠阵组装，同样基于回归思想实现基因组架构组装，最后通过重抽样技术对测序数据进行多次独立地采样、组装、改进、评估，然后整合成最终的组装基因组，可以减少由测序数据中的噪音给组装结果带来的不确定性，并且降低组装结果对组装参数选择的敏感性。
本发明实施方式还提供一种机器可读存储介质，其上存储有计算机程序指令，所述计算机程序指令被处理器执行时实现上述的基因组叠阵组装方法。
本发明实施方式还提供另一种机器可读存储介质，其上存储有计算机程序指令，所述计算机程序指令被处理器执行时实现上述的基因组架构组装方法。
本发明实施方式还提供一种机器可读存储介质，其上存储有计算机程序指令，所述计算机程序指令被处理器执行时实现上述的基因组序列组装方法。
本领域技术人员可以理解实现上述实施方式的方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序存储在一个存储介质中，包括若干指令用以使得单片机、芯片或处理器(processor)执行本发明各个实施方式所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
以上结合附图详细描述了本发明的可选实施方式，但是，本发明实施方式并不限于上述实施方式中的具体细节，在本发明实施方式的技术构思范围内，可以对本发明实施方式的技术方案进行多种简单变型，这些简单变型均属于本发明实施方式的保护范围。另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明实施方式对各种可能的组合方式不再另行说明。
此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明实施方式的思想，其同样应当视为本发明实施方式所公开的内容。

Claims (21)

1.一种基因组叠阵组装方法，其特征在于，所述方法包括：
S01：抽取一定深度的第一测序序列集；
S02：通过序列比对算法对所述第一测序序列集中的测序序列进行两两比对，得到测序序列的两两比对信息；
S03：根据测序序列的两两比对信息将所有满足叠落条件的比对共同表示为第一线性回归模型；
S04：通过迭代重加权最小二乘算法对所述第一线性回归模型求解，得到全局损失函数最小的解和所述第一测序序列集中的测序序列在被测基因组上的坐标的稳健估计；
S05：根据所述全局损失函数最小的解构造的无向图将所述测序序列划分为多个连通分量，按照各个连通分量内的测序序列的坐标的稳健估计对测序序列进行排列，生成叠阵；
S06：根据测序序列的两两比对信息对所述叠阵进行合并，得到初步组装叠阵集；
S07：获取所述初步组装叠阵集中各个叠阵的一致序列，将所述测序序列比对映射到所述一致序列，得到所述测序序列向一致序列的映射结果，所述映射结果为新的组装叠阵集。
2.根据权利要求1所述的基因组叠阵组装方法，其特征在于，步骤S03包括：
S031：假定第一测序序列集中的测序序列Ri的最右端碱基对齐到测序序列Rj内部，Rj右侧悬垂部分含有yi,j个碱基，由此得到一个观测样本：yi,j＝βj-βi+εi,j；
其中，βi是Ri在其所在的DNA单链上的坐标，βj是Rj在其所在的DNA单链上的坐标，εi,j表示观测样本的观测误差；
S032：假定第二测序序列Rj的最左端碱基对齐到第一测序序列Ri内部，Ri左侧悬垂部分含有
Figure FDA0003099519490000021
个碱基，由此得到一个对偶观测样本：
Figure FDA0003099519490000022
其中，
Figure FDA0003099519490000023
是Ri的反向互补序列在互补单链上的坐标，
Figure FDA0003099519490000024
是Rj的反向互补序列在互补单链上的坐标，
Figure FDA0003099519490000025
表示对偶观测样本的观测误差；
S033：根据测序序列的两两比对信息将所有满足叠落条件的比对的观测样本和对偶观测样本整合为一个矩阵方程：Y＝Xβ+ε；
其中，每一行为一个观测样本，向量
Figure FDA0003099519490000026
为待估参数，X是一个大型稀疏矩阵，X的每一行只有两个非零元素1和-1，X储存为稀疏矩阵格式。
3.根据权利要求2所述的基因组叠阵组装方法，其特征在于，步骤S04包括：
S041：利用大型稀疏线性方程组的加速求解算法求解关于β的稀疏线性方程组：XTX·β＝XTY，得到的解作为
Figure FDA0003099519490000027
的初始值；
S042：在每一轮迭代中，根据当前的
Figure FDA0003099519490000028
计算残差向量
Figure FDA0003099519490000029
以及
Figure FDA00030995194900000210
各分量绝对值的最大值
Figure FDA00030995194900000211
S043：如果
Figure FDA00030995194900000212
则计算权重矩阵
Figure FDA00030995194900000213
其中，α为预设的收敛阈值；权重函数为
Figure FDA00030995194900000214
S044：根据计算得到的W，取互为对偶的两个样本的权重值中的较小值赋予这两个样本；
S045：利用大型稀疏线性方程组的加速求解算法求解关于β的稀疏线性方程组：XTWX·β＝XTWY；
S046：用得到的解更新
Figure FDA0003099519490000031
删除矩阵X和Y中对应权重为0的行，返回进行下一轮迭代；
S047：直到
Figure FDA0003099519490000032
停止迭代，记录最后一轮迭代的矩阵X；
S048：对最后一轮
Figure FDA0003099519490000033
的每个分量取整后输出
Figure FDA0003099519490000034
即为所述第一测序序列集中的测序序列在被测基因组上的坐标的稳健估计。
4.根据权利要求3所述的基因组叠阵组装方法，其特征在于，所述步骤S05包括：
S051：以第一测序序列集中的测序序列及所述测序序列对应的反向互补序列作为顶点，根据最后一轮迭代的矩阵X中每一行的两个非零元素确定两条测序序列，将这两条测序序列对应的顶点之间建立一条边，构造得到无向图；
S052：遍历所述无向图，将所述无向图划分为多个连通分量，得到2k个连通分量，记为
Figure FDA0003099519490000035
其中Ci与
Figure FDA0003099519490000036
所包含的测序序列一一对应，互为反向互补序列；
S053：将各个连通分量中的测序序列根据其在被测基因组上的坐标的稳健估计进行排列，得到叠阵集：
Figure FDA0003099519490000037
叠阵Li和
Figure FDA0003099519490000038
分别对应于DNA双链上互为反向互补的两条子序列。
5.根据权利要求4所述的基因组叠阵组装方法，其特征在于，所述步骤S06包括：
S061：以所述叠阵集中的叠阵作为顶点；
S062：对于任意两个叠阵，若Li3’端s个测序序列分别和Lj5’端s个测序序列具有大于t个满足叠落条件的比对，则建立一条Li对应顶点指向Lj对应顶点的有向边，构建有向图；其中，边的权重设定为满足叠落条件的比对数，s和t的取值正相关于抽样深度；
S063：从所述有向图中包含最多测序序列的叠阵Lx开始，确定从叠阵Lx出发的权重最大的边所指向的叠阵Ly；
S064：从Lx的3’端选取一条测序序列Ri，和从Ly的5’端选取一条测序序列Rj，Ri和Rj之间的比对满足叠落条件；
S065：假定Ri的最右端碱基对齐到Rj内部，Rj右侧悬垂部分含有Δ个碱基，则将Ly中所有测序序列的坐标的稳健估计都加上
Figure FDA0003099519490000041
S067：将Lx和Ly中的测序序列按照其坐标的稳健估计的大小进行排列合并成一个叠阵，同时将Lx和
Figure FDA0003099519490000042
在有向图中对应的顶点去除；
S068：从Ly开始重复步骤S063-S067，直到没有可以合并的叠阵为止，再从有向图中剩下的顶点中包含测序序列最多的叠阵开始重复步骤S063-S067，直到遍历有向图中所有顶点为止，得到的合并后的叠阵作为初步组装叠阵集，记为D1,...,Dm。
6.根据权利要求5所述的基因组叠阵组装方法，其特征在于，所述步骤S07包括：
S071：对于所述初步组装叠阵集中的任意一个叠阵Di，取每一列中频率最大的碱基生成叠阵Di对应的一致序列，记为Ei；
S072：对于叠阵Di中的任意一条测序序列Rj的最左端碱基在Di中的坐标记为a1，Rj的最右端碱基在Di中的坐标记为a2，将Ei中第a1-τ个碱基到第a2+τ个碱基之间的子序列取出与Rj进行局部序列比对；
S073：将所述局部序列比对一致的部分作为Rj向Ei的映射结果，Di中的所有测序序列向Ei的映射结果构成了一个新的叠阵
Figure FDA0003099519490000051
Figure FDA0003099519490000052
为更新后的组装叠阵集。
7.一种基因组架构组装方法，其特征在于，所述方法包括：
S11：抽取一定深度的第二测序序列集；
S12：根据所述第二测序序列集组装基因组叠阵，得到组装叠阵集
Figure FDA0003099519490000053
S13：将所述第二测序序列集与组装叠阵集
Figure FDA0003099519490000054
的一致序列集进行比对，根据比对结果筛选出单映射的测序序列及其映射上的叠阵，根据设定的条件对所述单映射的测序序列及其映射上的叠阵进行筛选，得到唯一映射的测序序列信息和筛选后的一致序列集；
S14：基于唯一映射的测序序列信息和筛选后的一致序列集构建无向有权叠阵图，将所述无向有权叠阵图划分为多个连通分量；
S15：采用组合优化模型确定各个连通分量上的每个顶点所对应的叠阵的方向；
S16：建立回归模型，估计叠阵在被测基因组上的排列坐标，根据各个顶点对应的叠阵的方向和所述排列坐标对同一连通分量中的所有叠阵进行排序，得到该连通分量的准架构，根据多个连通分量的准架构得到架构集；
S17：对所述架构集中每个架构的叠阵排列坐标进行梳理，消除坐标重叠或者冗余的叠阵的矛盾情形，得到组装基因组架构。
8.根据权利要求7所述的基因组架构组装方法，其特征在于，所述步骤S12包括：采用权利要求1-6中任一项所述的基因组叠阵组装方法根据所述第二测序序列集组装基因组叠阵。
9.根据权利要求7所述的基因组架构组装方法，其特征在于，所述根据比对结果筛选出单映射测序序列叠阵，包括：
S1301：对于双末端测序序列，从比对结果中挑选出双末端序列的两端映射到同一个叠阵一致序列上的映射结果，以估计双末端测序片段库的库长，所述映射结果包含每一端的序列比对到叠阵一致序列上的起始点坐标，分别记为s1,t1和s2,t2；库长观测值记为：in＝max(s1,t1,s2,t2)-min(s1,t1,s2,t2)；对每个片段库分别取预定数量的观测值，取观测值的中位数M作为库长的观测，并计算其标准差σ；
S1302：将分别映射到两个不同的叠阵一致序列上的双末端测序序列的映射坐标仍分别记为s1,t1和s2,t2；
S1303：根据映射坐标判断所述双末端测序序列与其映射上的叠阵一致序列的方向是否一致，根据双末端测序序列与其映射上的叠阵一致序列的方向判断结果判断两个不同的叠阵一致序列的方向是否一致；根据两个不同的叠阵一致序列的方向判断结果和所述库长计算得到两个叠阵起始点之间的距离；筛除其中一端映射到叠阵上的位置超出预设的插入长度范围的连接，得到测序序列到两两叠阵一致序列的连接信息。
10.根据权利要求9所述的基因组架构组装方法，其特征在于，所述步骤S1303包括：
定义所述双末端测序序列的左端序列映射上的叠阵方向为正，且位置在前，假定左右端映射上的叠阵起点分别为F1,F2，则：
a)若s1＜t1且s2＜t2，则两个叠阵一致序列的方向各自与双端测序序列一致，即两叠阵的方向相反，则距离为：F2-F1＝in+s1+s2，筛除len1-t1＞in+3σ或者len2-t2＞in+3σ的连接；
b)若s1＜t1且s2＞t2，则两个叠阵一致序列的方向各自与双端测序序列相反，即两叠阵的方向一致，则距离为：F2-F1＝in+s1-s2，筛除len1-t1＞in+3σ或者s2＞in+3σ的连接；
c)若s1＞t1且s2＞t2，则两个叠阵一致序列的方向各自与双端测序序列一致，即两叠阵的方向相反，则距离为：F2-F1＝in-s1-s2，筛除s1＞in+3σ或者s2＞in+3σ的连接；
d)若s1＞t1且s2＜t2，则两个叠阵一致序列的方向各自与双端测序序列相反，即两叠阵的方向一致，则距离为：F2-F1＝in-s1+s2，筛除s1＞in+3σ或者len2-t2＞in+3σ的连接。
11.根据权利要求7所述的基因组架构组装方法，其特征在于，所述根据比对结果筛选出单映射的测序序列及其映射上的叠阵，包括：
S1311：对于单分子长读长测序序列，假定比对到多个叠阵一致序列上的比对区间在测序序列上的坐标按起点排序依次为[qs1,qe1],[qs2,qe2],...,[qsk,qek]，其对应的叠阵一致序列上的坐标依次为[cs1,ce1],[cs2,ce2],...,[csk,cek]；根据上述坐标得到任意两个比对到同一测序序列上的叠阵一致序列的方向关系和起点距离作为两两叠阵一致序列的连接信息；
其中，若叠阵与测序序列方向一致，则：Fi＝qsi-csi；若叠阵与测序序列方向相反，则：Fi＝qsi-csi+leni；距离为：Fj-Fi；
S1312：筛除同一测序序列上相互交叠的连接。
12.根据权利要求9或11所述的基因组架构组装方法，其特征在于，所述步骤S14包括：
S1401：定义叠阵之间形成边的连接数阈值参数m；
S1402：以组装叠阵集的一致序列集作为顶点，单映射测序序列叠阵之间的连接数大于阈值参数m的形成一条边，构建无向有权叠阵图；
S1403：将无向有权叠阵图中顶点i和顶点j之间的同向连接数和异向连接数分别记为aij和bij，对应的距离信息分别录入同向和异向的叠阵距离列表；
S1404：将所述无向有权叠阵图划分为多个连通分量。
13.根据权利要求12所述的基因组架构组装方法，其特征在于，所述步骤S15包括：
S1501：以aij和bij中的极大值给连通分量上的每条边赋予权重并标记该边的方向，若aij≥bij，则权重wij＝aij，Dij＝1；否则权重wij＝bij，Dij＝-1；遍历所有顶点得到连通分量的每个顶点的初始方向；
S1502：迭代优化当前连通分量上的矛盾方向连接总数，每轮迭代中分别计算改变每个顶点导致目标函数的单点变化量，若存在某顶点的单点变化量为负，就选取单点变化量与总连接数比值最大的顶点改变其方向，更新方向后继续下一轮迭代，直到改变任一顶点都无法使得总矛盾数更小，停止迭代，得到每个顶点所对应的叠阵的方向。
14.根据权利要求13所述的基因组架构组装方法，其特征在于，所述步骤S16包括：
S1601：建立回归模型；
S1602：根据所述两两叠阵一致序列的连接信息以及每个顶点所对应的叠阵的方向，确定每条边的类型；
S1603：根据每条边的类型，取出每条边对应的叠阵距离列表的中位数作为回归模型的观测值；
S1604：设置加权截断最小二乘估计的残差允许阈值maxError；
S1605：以库长的方差估计的倒数与当前边所包含的符合当前顶点方向的连接数的乘积作为初始权重，即
Figure FDA0003099519490000091
计算叠阵一致序列起点坐标的加权最小二乘估计
Figure FDA0003099519490000092
其中系数矩阵X除第一行为(1,0,...,0)，其余行表示对应顶点i与j之间的连边，故第i项为-1，第j项为1；
S1606：根据上一轮的估计计算所有样本残差，当残差最大值小于残差允许阈值maxError，停止迭代，得到可行估计
Figure FDA0003099519490000093
否则选取样本中残差小于(100-iter)％分位数的样本，其中，iter表示当前迭代次数；
判断选取的样本对应的系数矩阵是否列满秩：若不满秩，根据系数矩阵将当前子图进一步划分成若干个连通分量，在子图的每个连通分量上从步骤S1501重新开始执行；若满秩，则在选取的样本子集上再次进行加权最小二乘估计，得到
Figure FDA0003099519490000094
再重复步骤S1606，直到得到可行估计
Figure FDA0003099519490000095
S1607：根据每个顶点所对应的叠阵的方向将当前连通分量中的所有叠阵按照其对应的可行估计进行排序，整理到同一个坐标系中，得到一条准架构，根据多个连通分量的准架构得到架构集。
15.根据权利要求14所述的基因组架构组装方法，其特征在于，所述步骤S17包括：
S1701：将架构中的所有叠阵一致序列按每个顶点所对应的叠阵的方向进行调整，将反向的叠阵的起点与终点坐标对换，使得架构排列中的所有叠阵起点坐标均在终点坐标之前，将叠阵按起点坐标由小到大排序；
S1702：从排序后的第一个叠阵开始，依次向后，对坐标上被当前选入的最后一个叠阵包含或交叉的叠阵进行局部序列比对，比对成功的则去除被包含的叠阵，或将被交叉的叠阵并入当前叠阵；
选择与当前叠阵终点坐标紧邻又有连接信息的叠阵加入当前架构，留置与当前叠阵终点坐标紧邻但无连接信息的叠阵；其中，
1)若两个叠阵一致序列i与j的坐标满足：
0＜startj-endi≤2kbp，且endj＞endi；则两个叠阵紧邻；
2)若两个叠阵一致序列i与j的坐标满足：
starti＜startj＜endj＜endi；则i包含j；
3)若两个叠阵一致序列i与j的坐标满足：
starti＜startj＜endi＜endj；则两个叠阵交叉；
其中，start,end分别表示叠阵一致序列在准架构中的起点和终点坐标；
S1703：更新当前叠阵为最新加入当前架构的叠阵，重复步骤S1702，直至到达该架构上的最后一个叠阵，结束当前架构，输出到最终的组装架构文件；
S1704：在被留置的叠阵中重复步骤S1702-S1703，确定下一个架构排列并输出至最终的组装架构文件；重复此步骤直至所有叠阵被安排进某一个最终的架构之中，得到组装基因组架构、叠阵在架构中的位置和方向信息、紧邻叠阵之间间隙长度的估计及估计方差、由坐标交叉且比对成功而合并得到的新的叠阵一致序列。
16.根据权利要求14所述的基因组架构组装方法，其特征在于，所述方法还包括：用所述准架构组成的准架构集替代叠阵一致序列集，重复步骤S14-S16，得到架构集。
17.一种基因组序列组装方法，其特征在于，所述方法包括：
S21：抽取一定深度的第三测序序列集，用于根据权利要求1-6中任一项所述的基因组叠阵组装方法生成组装叠阵集，所述组装叠阵集用于根据权利要求7-16中任一项所述的基因组架构组装方法生成组装基因组架构；
S22：抽取另一深度的第四测序序列集更新改进所述组装基因组架构，并对改进后的组装基因组进行质量评估，若质量达到标准，则将改进后的组装基因组储存为一条候选组装基因组；
S23：重复步骤S21-S22，直到得到预定数量的候选组装基因组，将所有候选组装基因组进行多重序列比对，将比对的每一列中频率最大的碱基连接起来得到最终的组装基因组。
18.根据权利要求17所述的基因组序列组装方法，其特征在于，所述步骤S22包括：
S2201：将第四测序序列
Figure FDA0003099519490000111
映射到组装基因组，去除其中同时映射到多个位置的测序序列，形成单映测序序列叠阵；对于叠阵中的每一列，如果其被多于d0条单映序列覆盖，则取该列中频率最大的碱基更新组装基因组的相应位置的碱基，得到改进后的组装基因组；
S2202：将第四测序序列
Figure FDA0003099519490000112
映射到改进后的组装基因组，根据映射结果对其质量进行评估，若质量达到标准，则将其储存为一条候选组装基因组。
19.一种基因组叠阵组装系统，其特征在于，所述系统应用权利要求1-6中任一项所述的基因组叠阵组装方法进行基因组叠阵组装。
20.一种基因组架构组装系统，其特征在于，所述系统应用权利要求7-16中任一项所述的基因组架构组装方法进行基因组架构组装。
21.一种基因组序列组装系统，其特征在于，所述系统应用权利要求17或18所述的基因组序列组装方法进行基因组序列组装。
CN202110620095.7A 2021-06-03 2021-06-03 基因组叠阵、基因组架构、基因组序列组装方法及系统 Pending CN115440302A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110620095.7A CN115440302A (zh) 2021-06-03 2021-06-03 基因组叠阵、基因组架构、基因组序列组装方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110620095.7A CN115440302A (zh) 2021-06-03 2021-06-03 基因组叠阵、基因组架构、基因组序列组装方法及系统

Publications (1)

Publication Number Publication Date
CN115440302A true CN115440302A (zh) 2022-12-06

Family

ID=84271563

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110620095.7A Pending CN115440302A (zh) 2021-06-03 2021-06-03 基因组叠阵、基因组架构、基因组序列组装方法及系统

Country Status (1)

Country Link
CN (1) CN115440302A (zh)

Similar Documents

Publication Publication Date Title
US10347365B2 (en) Systems and methods for visualizing a pattern in a dataset
US10600217B2 (en) Methods for the graphical representation of genomic sequence data
CN107133493B (zh) 基因组序列的组装方法、结构变异探测方法和相应的系统
CA2424031C (en) System and process for validating, aligning and reordering genetic sequence maps using ordered restriction map
US11482305B2 (en) Artificial intelligence analysis of RNA transcriptome for drug discovery
US20230222311A1 (en) Generating machine learning models using genetic data
US20120197533A1 (en) Identifying rearrangements in a sequenced genome
WO2013040583A2 (en) Determining variants in a genome of a heterogeneous sample
CN116741397B (zh) 基于多组学数据融合的癌症分型方法、系统及存储介质
CN113555062B (zh) 一种用于基因组碱基变异检测的数据分析系统及分析方法
US20120185177A1 (en) Harnessing high throughput sequencing for multiplexed specimen analysis
US20150178446A1 (en) Iterative clustering of sequence reads for error correction
CN108595912A (zh) 检测染色体非整倍性的方法、装置及系统
Genovese et al. SpeedHap: an accurate heuristic for the single individual SNP haplotyping problem with many gaps, high reading error rate and low coverage
CN115440302A (zh) 基因组叠阵、基因组架构、基因组序列组装方法及系统
JP5129809B2 (ja) ハプロタイプ推定装置、および、プログラム
WO2013191637A1 (en) Method and device for efficient calculation of allele ratio confidence intervals and uses thereof
Deshpande et al. Reconstructing and characterizing focal amplifications in cancer using AmpliconArchitect
Oehl A combinatorial approach for reconstructing rDNA repeats
Ramírez-Rafael¹ et al. Check for updates
Oluwadare Data mining and machine learning methods for chromosome conformation data analysis
Lu Big Data Analytics in Metagenomics: Integration, Representation, Management, and Visualization
CN118116462A (zh) 基于tdfps算法针对纳米孔测序中条形码的设计方法
Li Applications of Clustering Algorithms to Entity Resolution and Human Genome Variation
Linthorst et al. Computer Science Bioinformatics Track

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination