CN114649055B - 用于检测单核苷酸变异和***缺失的方法、设备和介质 - Google Patents

用于检测单核苷酸变异和***缺失的方法、设备和介质 Download PDF

Info

Publication number
CN114649055B
CN114649055B CN202210392298.XA CN202210392298A CN114649055B CN 114649055 B CN114649055 B CN 114649055B CN 202210392298 A CN202210392298 A CN 202210392298A CN 114649055 B CN114649055 B CN 114649055B
Authority
CN
China
Prior art keywords
single nucleotide
probability
activation
sequence
tandem repeat
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.)
Active
Application number
CN202210392298.XA
Other languages
English (en)
Other versions
CN114649055A (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.)
Berry Genomics Co Ltd
Original Assignee
Berry Genomics Co Ltd
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 Berry Genomics Co Ltd filed Critical Berry Genomics Co Ltd
Priority to CN202210392298.XA priority Critical patent/CN114649055B/zh
Publication of CN114649055A publication Critical patent/CN114649055A/zh
Application granted granted Critical
Publication of CN114649055B publication Critical patent/CN114649055B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • 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
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • General Health & Medical Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Biotechnology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioethics (AREA)
  • Public Health (AREA)
  • Evolutionary Computation (AREA)
  • Epidemiology (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明涉及一种用于检测单核苷酸变异和***缺失的方法、计算设备和存储介质。该方法包括：获取待测样本的比对结果数据和参考基因组序列数据；基于比对结果数据和参考基因组序列数据，生成短串联重复序列位点特征矩阵；计算待测样本的检测染色体上各个基因组位点的激活概率，以便确定多个激活区域；经由多线程技术，并行地检测每个激活区域的单核苷酸变异和***缺失，以便获取单核苷酸变异和***缺失特征值；以及基于短串联重复序列位点特征矩阵和单核苷酸变异和***缺失特征值，经由预测模型生成预测结果。本发明不仅能够显著提升检测单核苷酸变异和***缺失的速度，降低检测计算所需内存消耗，而且能够提高检测结果的准确性。

Description

用于检测单核苷酸变异和***缺失的方法、设备和介质
技术领域
本发明总体上涉及生物信息处理，并且具体地，涉及用于检测单核苷酸变异和***缺失的方法、计算设备和计算机存储介质。
背景技术
单核苷酸变异（SNP）是指基因组上单个核苷酸发生了突变，***缺失（INDEL）指的是基因组上小于50bp的序列***和缺失。SNP和INDEL变异与表型和疾病密切相关。因此，如何快速并准确地检测单核苷酸变异和***缺失变得尤为重要。
传统的检测SNP和INDEL的方法包括：基于比对序列堆叠（pileup）的方法（例如，samtools）、基于局部组装单倍体的方法（例如，platypus、GATK HaplotypeCaller）、基于概率计算单倍体的方法（例如freebayes）、结合pileup和局部单倍体组装的方法（例如strelka2）、以及基于深度学习的方法（例如deepvariant、clair）。
上述传统的检测SNP和INDEL的方法难以针对单核苷酸变异和***缺失均能快速地获得高准确率的检测结果。例如，基于比对序列堆叠的方法可以准确地检测SNP，但是检测INDEL的准确度相对较低。再例如，基于深度学习的方法虽然具有较高的准确性，但是同时也存在运行速度慢的问题。
综上，传统的检测SNP和INDEL检测单核苷酸变异和***缺失的方案存在的不足之处在于：难以同时针对单核苷酸变异和***缺失均能快速地获得高准确率的检测结果。
发明内容
本发明提供一种用于检测单核苷酸变异和***缺失的方法、计算设备和计算机存储介质，不仅能够显著提升检测单核苷酸变异和***缺失的速度，降低检测计算所需内存消耗，而且能够提高检测结果的准确性。
根据本发明的第一方面，提供了一种用于检测单核苷酸变异和***缺失的方法。该方法包括：获取待测样本的测序序列数据与参考基因组序列数据的比对结果数据和参考基因组序列数据；基于比对结果数据和参考基因组序列数据，生成短串联重复序列位点特征矩阵；计算待测样本的检测染色体上各个基因组位点的激活概率，以便确定多个激活区域；经由多线程技术，并行地检测多个激活区域中的每个激活区域的单核苷酸变异和***缺失，以便获取单核苷酸变异和***缺失特征值；以及基于短串联重复序列位点特征矩阵和单核苷酸变异和***缺失特征值，经由预测模型生成关于单核苷酸变异和***缺失的预测结果。
根据本发明的第二方面，还提供了一种计算设备，该设备包括：存储器，被配置为存储一个或多个计算机程序；以及处理器，耦合至存储器并且被配置为执行一个或多个程序使装置执行本发明的第一方面的方法。
根据本发明的第三方面，还提供了一种非瞬态计算机可读存储介质。该非瞬态计算机可读存储介质上存储有机器可执行指令，该机器可执行指令在被执行时使机器执行本发明的第一方面的方法。
在一些实施例中，并行地检测多个激活区域中的每个激活区域的单核苷酸变异和***缺失包括：将每个激活区向两侧各延伸预定数量个碱基，以便得到对应的延伸区间；以及针对每个延伸区间内所有的读长进行局部组装，以便并行地检测单核苷酸变异和***缺失，以获取单核苷酸变异和***缺失特征值。
在一些实施例中，生成短串联重复序列位点特征矩阵包括：基于参考基因组序列数据，查找参考基因组上的短串联重复序列位点数据；针对短串联重复序列位点数据进行均匀抽样，以便获得短串联重复序列位点抽样信息，短串联重复序列位点抽样信息指示短串联重复序列位点的位点位置、单元长度和单元数；基于比对结果数据和短串联重复序列位点抽样信息，经由多线程池，分别计算每个短串联重复序列位点的深度和***缺失数目，以便生成短串联重复序列位点特征数据；以及基于所生成的短串联重复序列位点特征数据，计算短串联重复序列位点特征矩阵。
在一些实施例中，基于所生成的短串联重复序列位点特征数据，计算短串联重复序列位点特征矩阵包括：确定所有短串联重复序列位点的深度和***缺失数目是否被计算完毕；以及响应于确定所有短串联重复序列位点的深度和***缺失数目被计算完毕，基于短串联重复序列位点特征数据，计算短串联重复序列位点特征矩阵，短串联重复序列位点特征矩阵包括：开空位概率矩阵、空位延伸概率矩阵和***缺失位点先验概率矩阵。
在一些实施例中，计算待测样本的检测染色体上各个基因组位点的激活概率，以便确定多个激活区域包括：配置用于存储激活概率的数组，数组的长度等于待检测染色体的长度；将待检测染色体分为预定长度的多个区间；并行地计算多个区间中的每个区间内基因组位点的激活概率，以便将所计算的激活概率存储至数组；以及使用高斯核对整条待测染色体上的各基因组位点的激活概率进行平滑，以便将激活概率大于或者等于预定阈值的连续基因组位点合并生成激活区域。
在一些实施例中，并行地计算多个区间中的每个区间内基因组位点的激活概率包括：确定当前比对序列堆叠的备选等位基因的读长数目是否大于或者等于预定读长数目阈值；响应于确定当前比对序列堆叠的备选等位基因的读长数目大于或者等于预定读长数目阈值，将当前比对序列堆叠识别为候选比对序列堆叠；响应于确定当前比对序列堆叠的备选等位基因的读长数目小于预定读长数目阈值，将当前比对序列堆叠位点的激活概率确定为零；以及使用贝叶斯算法计算候选比对序列堆叠的激活概率。
在一些实施例中，针对每个延伸区间内所有的读长进行局部组装以便并行地检测单核苷酸变异和***缺失，以获取单核苷酸变异和***缺失特征值包括：针对每个延伸区间内所有的读长进行局部组装，以便生成经组装的单倍型序列；将单倍型序列比对到参考基因组，以便检测单核苷酸变异和***缺失；保留在激活区内的、所检测的单核苷酸变异和***缺失；计算读长比对到单倍型序列的概率；以及计算单核苷酸变异位点和***缺失位点的基因型概率，以便获取单核苷酸变异和***缺失特征值。
在一些实施例中，计算读长比对到单倍型序列的概率包括：基于读长、单倍型序列、读长上每个碱基的质量值、读长上每个碱基的开空位概率、读长上每个碱基的空位延伸概率，计算读长比对到单倍型序列的概率。
提供发明内容部分是为了以简化的形式来介绍对概念的选择，它们在下文的具体实施方式中将被进一步描述。发明内容部分无意标识本发明的关键特征或主要特征，也无意限制本发明的范围。
附图说明
图1示出了根据本发明的实施例的用于实施检测单核苷酸变异和***缺失的方法的***的示意图。
图2示出了根据本发明的实施例的用于检测单核苷酸变异和***缺失的方法的流程图。
图3示出了根据本发明的实施例的用于生成短串联重复序列位点特征矩阵的方法的流程图。
图4示出了根据本发明的实施例的用于生成激活区域的方法的流程图。
图5示出了根据本发明的实施例的对序列堆叠的示意图。
图6示出了根据本发明的实施例的用于获取单核苷酸变异和***缺失特征的方法的流程图。
图7示出了根据本发明的实施例的用于生成激活区域的方法的示意图。
图8示意性示出了适于用来实现本发明实施例的电子设备的框图。
图9示出了根据本发明的实施例的用于生成短串联重复序列位点特征矩阵的方法的示意图。
在各个附图中，相同或对应的标号表示相同或对应的部分。
具体实施方式
下面将参照附图更详细地描述本发明的优选实施例。虽然附图中显示了本发明的优选实施例，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
在本文中使用的术语“包括”及其变形表示开放性包括，即“包括但不限于”。除非特别申明，术语“或”表示“和/或”。术语“基于”表示“至少部分地基于”。术语“一个示例实施例”和“一个实施例”表示“至少一个示例实施例”。术语“另一实施例”表示“至少一个另外的实施例”。术语“第一”、“第二”等等可以指代不同的或相同的对象。
如前文所描述，传统的检测单核苷酸变异和***缺失的方案存在的不足之处在于：难以同时针对单核苷酸变异和***缺失均能快速地获得高准确率的检测结果。
为了至少部分地解决上述问题以及其他潜在问题中的一个或者多个，本发明的示例实施例提出了一种用于检测单核苷酸变异和***缺失的方案。该方案包括：通过基于所生成的短串联重复序列位点特征矩阵和基于激活区域而生成的单核苷酸变异和***缺失特征值，经由预测模型生成关于单核苷酸变异和***缺失的预测结果，本发明可以基于待测样本的测序比对结果数据和参考基因组序列数据生成STR位点发生INDEL的概率模型，并使用STR位点发生INDEL的概率模型计算INDEL的后验基因型概率，提升了INDEL检测的准确性，本发明可以基于更为全面的变异特征数据来生成SNP、INDEL预测结果，因而有利于提高检测结果的准确性。另外，本发明使用多线程技术并行地计算STR位点特征值、计算参考基因组位点的激活概率和检测激活区域内的单核苷酸变异和***缺失变异，本发明可以显著提升检测SNP、INDEL的速度，减少了内存消耗。因而，本发明不仅能够显著提升检测单核苷酸变异和***缺失的速度，降低检测计算所需内存消耗，而且能够提高检测结果的准确性。
图1示出了根据本发明的实施例的用于检测单核苷酸变异和***缺失的方法的***100的示意图。如图1所示，***100包括：计算设备110、测序设备130、网络140。在一些实施例中，计算设备110、测序设备130经由网络140进行数据交互。
关于测序设备130，其例如用于针对待测对象的待测样本进行测序，以便生成关于待测样本的测序序列数据；以及将所生成的待测样本的测序序列数据发送至计算设备110。
关于计算设备110，其例如用于获取待测样本的测序序列数据与参考基因组序列数据的比对结果数据和参考基因组序列数据；以及基于比对结果数据和参考基因组序列数据，生成短串联重复序列位点特征矩阵。计算设备110还可以用于计算待测样本的检测染色体上各个基因组位点的激活概率，以便确定多个激活区域；经由多线程技术，并行地检测多个激活区域中的每个激活区域的单核苷酸变异和***缺失，以便获取单核苷酸变异和***缺失特征值；以及基于短串联重复序列位点特征矩阵和单核苷酸变异和***缺失特征值，经由预测模型生成关于单核苷酸变异和***缺失的预测结果。
在一些实施例中，计算设备110可以具有一个或多个处理单元，包括诸如GPU、FPGA和ASIC等的专用处理单元以及诸如CPU的通用处理单元。另外，在每个计算设备上也可以运行着一个或多个虚拟机。计算设备110例如包括：比对结果数据和参考基因组序列数据获取单元112、短串联重复序列位点特征矩阵生成单元114、激活区域确定单元116、单核苷酸变异和***缺失特征值获取单元118，预测结果生成单元120。上述比对结果数据和参考基因组序列数据获取单元112、短串联重复序列位点特征矩阵生成单元114、激活区域确定单元116。单核苷酸变异和***缺失特征值获取单元118，预测结果生成单元120可以配置在一个或者多个计算设备110上。
关于比对结果数据和参考基因组序列数据获取单元112，其用于获取待测样本的测序序列数据与参考基因组序列数据的比对结果数据和参考基因组序列数据。
关于短串联重复序列位点特征矩阵生成单元114，其用于基于比对结果数据和参考基因组序列数据，生成短串联重复序列位点特征矩阵。
关于激活区域确定单元116，其用于计算待测样本的检测染色体上各个基因组位点的激活概率，以便确定多个激活区域。
关于单核苷酸变异和***缺失特征值获取单元118，其用于经由多线程技术，并行地检测多个激活区域中的每个激活区域的单核苷酸变异和***缺失，以便获取单核苷酸变异和***缺失特征值。
关于预测结果生成单元120，其用于基于短串联重复序列位点特征矩阵和单核苷酸变异和***缺失特征值，经由预测模型生成关于单核苷酸变异和***缺失的预测结果。
以下将结合图2、图7、图9描述根据本发明的实施例的用于检测单核苷酸变异和***缺失的方法。图2示出了根据本发明的实施例的用于检测单核苷酸变异和***缺失的方法200的流程图。图7示出了根据本发明的实施例的用于生成激活区域的方法700的示意图。图9示出了根据本发明的实施例的用于生成短串联重复序列位点特征矩阵的方法的示意图。应当理解，方法200例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解，方法200还可以包括未示出的附加动作和/或可以省略所示出的动作，本发明的范围在此方面不受限制。
在步骤202处，计算设备110获取待测样本的测序序列数据与参考基因组序列数据的比对结果数据和参考基因组序列数据。
关于待测样本，其例如是待测对象的组织样本或者血液样本。关于测序序列数据，其例如是待测样本的DNA测序数据。
在步骤204处，计算设备110基于比对结果数据和参考基因组序列数据，生成短串联重复序列位点特征矩阵。
关于短串联重复序列位点特征矩阵，其例如包括： 开空位（Gap open）概率矩阵、空位延伸（Gap extending）概率矩阵和***缺失（INDEL）先验概率矩阵。应当理解， Gapopen概率矩阵的每个特征值为Gap open概率。Gap extending概率矩阵的每个特征值为Gapextending概率。INDEL先验概率矩阵的每个特征值为INDEL先验概率。Gap open概率矩阵、Gap extending概率矩阵和INDEL先验概率矩阵均是二维矩阵，并且形状一致，每个矩阵的行数为最大STR单元长度（例如为8），列数为最大STR单元数（例如为20）。
如图9所示，在一些实施例中，关于生成短串联重复序列位点特征矩阵的方法主要包括：计算设备110查找参考基因组上的STR位点数据；经由多线程池（例如，图9所示的线程1至线程N），分别计算每个STR位点的深度和***缺失数目，以便生成STR位点特征数据（或者STR位点特征值）；以及基于STR位点特征数据（或者STR位点特征值），计算开空位概率矩阵（即，Gap open概率矩阵，或简称为“GOP矩阵”）、 空位延伸概率矩阵（即，Gap extending概率矩阵，或简称为“GCP矩阵”）、INDEL先验概率矩阵（即，INDEL先验概率矩阵，或简称为“API矩阵”）。
应当理解，可以通过多种方式计算开空位概率矩阵、 空位延伸概率矩阵和INDEL先验概率矩阵。在一些实施例中，关于基于STR位点特征数据，计算开空位概率矩阵、 空位延伸概率矩阵和INDEL先验概率矩阵的方法，下文将结合公式（1）至（10）加以说明。例如，公式（1）示出了空位延伸概率矩阵的计算方式。
Figure 101055DEST_PATH_IMAGE002
在上述公式（1）中，M代表单元长度。N代表单元数目。
Figure 1884DEST_PATH_IMAGE003
代表单元长度为M、单元数目为N的空位延伸概率矩阵。
开空位概率矩阵和INDEL先验概率矩阵可以通过Grid-Search方法计算。Grid-Search的第一维是先验概率列表Api（
Figure 719304DEST_PATH_IMAGE004
），第二维是STR区域发生INDEL的概率列表Gp（
Figure DEST_PATH_IMAGE005
），按照i和j的步长为1进行迭代，目标是使概率函数
Figure 973568DEST_PATH_IMAGE006
最大化，L为具有相同STR单元长度和单元数目的位点数目。
Figure 974891DEST_PATH_IMAGE007
在上述公式（2）中，k代表STR位点处INDEL的个数。d代表STR位点处的总深度。以下结合公式（3）至（7）说明上文中的公式（2）中的多个变量
Figure 735037DEST_PATH_IMAGE008
、
Figure 191950DEST_PATH_IMAGE009
、
Figure 328534DEST_PATH_IMAGE010
、
Figure 856467DEST_PATH_IMAGE011
和
Figure 771202DEST_PATH_IMAGE012
的计算方式。
Figure 463215DEST_PATH_IMAGE013
求解概率函数的
Figure 183915DEST_PATH_IMAGE014
最大值为第一个局部最大值，即迭代至
Figure 644983DEST_PATH_IMAGE014
不再增大时，中止迭代并记录
Figure 999129DEST_PATH_IMAGE015
，此时，可以求得
Figure 444017DEST_PATH_IMAGE016
和
Figure 640512DEST_PATH_IMAGE017
。
开空位概率矩阵
Figure 143037DEST_PATH_IMAGE018
可以通过
Figure 743783DEST_PATH_IMAGE019
和空位延伸概率矩阵
Figure 331759DEST_PATH_IMAGE020
进行计算，在Gop（
Figure 348256DEST_PATH_IMAGE021
）范围内按照o的步长为0.25进行迭代，求解损失函数C的全局最小值，此时
Figure 758816DEST_PATH_IMAGE022
。
Figure 468146DEST_PATH_IMAGE023
在上述公式（8）至（10）中，Abs（）代表求绝对值。C代表损失函数。
Figure 605735DEST_PATH_IMAGE024
代表空位概率。
Figure 222661DEST_PATH_IMAGE025
代表非空位概率。
下文将结合图3具体说明用于生成短串联重复序列位点特征矩阵的方法300，在此，不再赘述。
在步骤206处，计算设备110计算待测样本的检测染色体上各个基因组位点的激活概率，以便确定多个激活区域。在一些实施例中，步骤206处也可以在步骤204之前进行。
关于确定激活区域的方法，其例如包括：计算设备110配置用于存储激活概率的数组，数组的长度等于待检测染色体的长度；将待检测染色体分为预定长度的多个区间；并行地计算多个区间中的每个区间内基因组位点的激活概率，以便将所计算的激活概率存储至数组；以及使用高斯核对整条待测染色体上的各基因组位点的激活概率进行平滑，以便将激活概率大于或者等于预定阈值的连续基因组位点合并生成激活区域。
例如，如图7所示，计算设备110在步骤710处配置用于存储激活概率的数组720，该数组720例如为空的数组、且长度等于待检测染色体的长度；然后，计算设备110在步骤712处将待检测染色体分为预定长度的多个区间。该多个区间例如为N个区间，N为自然数。该N个区间例如包括：标记730所指示的区间1、标记732所指示的区间2至标记734所指示的区间N。之后，计算设备110在步骤714处利用多个线程（例如线程1、线程2至线程M）并行地计算多个区间（例如N个区间）中的每个区间内基因组位点的激活概率。然后，计算设备110将所计算的激活概率存储至数组，以便形成标记722所指示的存储了整条待测染色体上的各基因组位点的激活概率的数组；之后，计算设备110在步骤716处使用高斯核对整条待测染色体上的各基因组位点的激活概率进行平滑，以便生成标记724所指示的经由平滑的激活概率数组。计算设备110在步骤718处基于经由平滑的激活概率数组确定连续基因组位点的激活概率是否大于或者等于预定阈值，以便将激活概率大于或者等于预定阈值的连续基因组位点合并生成激活区域。
下文将结合图4具体说明用于确定激活区域的方法400，在此，不再赘述。
通过确定激活区域，本发明可以确定可能存在变异的区域，以便于后续进一步针对可能存在变异的区域来检测单核苷酸变异和***缺失，由此利于提升单核苷酸变异和***缺失的检测效率和检测结果的准确率。另外，通过仅对每个区间内基因组位点的激活概率进行并行计算，而不对各基因组位点的激活概率的平滑计算也进行并行处理，本发明不仅能够提高确定激活区域的速度，而且使得计算结果与激活概率的串行计算结果相一致，因而，本发明能够在兼顾准确性的同时，提高计算速度。
在步骤208处，计算设备110经由多线程技术，并行地检测多个激活区域中的每个激活区域的单核苷酸变异和***缺失，以便获取单核苷酸变异和***缺失特征值。
例如，计算设备110将每个激活区向两侧各延伸预定数量个碱基，以便得到对应的延伸区间；以及针对每个延伸区间内所有的读长进行局部组装，以便并行地检测单核苷酸变异和***缺失，以获取单核苷酸变异和***缺失特征值。通过对每个激活区域向两侧延伸，再针对延伸之后的延伸区间并行地检测单核苷酸变异和***缺失，本发明可以解决激活区域边缘处难以检测单核苷酸变异和***缺失的问题，进而可以进一步提高检测单核苷酸变异和***缺失的准确性。
在一些实施例中，关于针对每个延伸区间内所有的读长进行局部组装，以获取单核苷酸变异和***缺失特征值的方法，其例如包括：计算设备110针对每个延伸区间内所有的读长进行局部组装，以便生成经组装的单倍型序列；将单倍型序列比对到参考基因组，以便检测单核苷酸变异和***缺失；保留在激活区内的、所检测的单核苷酸变异和***缺失；计算读长比对到单倍型序列的概率；以及计算单核苷酸变异位点和***缺失位点的基因型概率，以便获取单核苷酸变异和***缺失特征值。
下文将结合图6具体说明用于获取单核苷酸变异和***缺失特征的方法600，在此，不再赘述。
在步骤210处，计算设备110基于短串联重复序列位点特征矩阵和单核苷酸变异和***缺失特征值，经由预测模型生成关于单核苷酸变异和***缺失的预测结果。例如，计算设备110使用多个线程，经由多样本经训练的预测模型提取特征值输入矩阵M的特征，以便生成关于单核苷酸变异和***缺失的预测结果。
关于预测模型的输入特征，其例如是计算设备110基于短串联重复序列位点特征矩阵和单核苷酸变异和***缺失特征值而生成的特征值输入矩阵M。
关于预测模型，其例如而不限于是基于随机森林模型而构建的。该预测模型例如是经由关于单核苷酸变异和***缺失的多个样本进行训练的。用于训练预测模型的样本例如是已知单核苷酸变异和***缺失“真”和“假”的分类结果的特征值输入矩阵数据。
关于单核苷酸变异和***缺失的预测结果，其例如是关于单核苷酸变异和***缺失的分类和概率值。例如，单核苷酸变异“真”的概率值为0.9，单核苷酸变异“假”的概率为0.1。***缺失“真”的概率值为0.8，***缺失“假”的概率为0.1。
在一些实施例中，方法200还包括以标准的VCF格式输出SNP和INDEL变异。
经研究发现，传统的检测SNP和INDEL的方案，例如是基于DRAGEN-GATK的检测方案。在该检测方案中，用户需要将基因组拆分为多个不同的染色体，使用多个进程，每个进程处理一条染色体，然后合并每个进程输出的变异检测结果。这种多进程的方式存在一些不足之处：由于进程之间不能共享数据，每个进程都需要4-5Gb的内存，进而造成内存占用成倍的增加。例如10个进程需要40-50G的内存。另外，鉴于染色体长度不同，序列复杂度也不同，因而造成多进程之间负载不平衡，整体运行时间等于处理速度最慢的进程。此外，该检测方案数据处理较为流程复杂，需要用户处理的步骤多、不够便捷。
与DRAGEN-GATK等现有检测方法相比，本发明例如使用10Gbp数据量57Mbp捕获区间的全外显子数据进行测试。如果本发明方法使用10个线程，则生成SNP、INDEL的检测结果只需11.5分钟，所需的最大驻留内存为2.2G。以下表一示例性示出了本发明的检测SNP、INDEL的方法在运行时间和最大驻留内存方面与现有技术的对比数据。根据表一所示的对比数据可知，本发明显著提升了检测SNP、INDEL的速度，并减少了内存消耗。
表一
Figure 762096DEST_PATH_IMAGE026
在上述方案中，通过基于所生成的短串联重复序列位点特征矩阵和基于激活区域而生成的单核苷酸变异和***缺失特征值，经由预测模型生成关于单核苷酸变异和***缺失的预测结果，本发明基于待测样本的测序比对结果数据和参考基因组序列数据生成STR位点发生INDEL的概率模型，并使用STR位点发生INDEL的概率模型计算INDEL的后验基因型概率，提升了INDEL检测的准确性，本发明可以基于更为全面的变异特征数据来生成SNP、INDEL预测结果，因而有利于提高检测结果的准确性。另外，本发明使用多线程技术并行地计算STR位点特征值、计算参考基因组位点的激活概率和检测激活区域内的单核苷酸变异和***缺失变异，本发明可以显著提升检测SNP、INDEL的速度，减少了内存消耗。因而，本发明不仅能够显著提升检测单核苷酸变异和***缺失的速度，降低检测计算所需内存消耗，而且能够提高检测结果的准确性。检测
以下将结合图3描述根据本发明的实施例的用于生成短串联重复序列位点特征矩阵的方法。图3示出了根据本发明的实施例的用于生成短串联重复序列位点特征矩阵的方法300的流程图。应当理解，方法300例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解，方法300还可以包括未示出的附加动作和/或可以省略所示出的动作，本发明的范围在此方面不受限制。
在步骤302处，计算设备110基于参考基因组序列数据，查找参考基因组上的短串联重复序列位点数据。
在步骤304处，计算设备110针对短串联重复序列位点数据进行均匀抽样，以便获得短串联重复序列（或者简称为“STR”）位点抽样信息，短串联重复序列位点抽样信息指示短串联重复序列位点的位点位置、单元长度和单元数。
关于短串联重复序列位点信息，其例如是STR位点表格
Figure 642327DEST_PATH_IMAGE027
，STR位点表格
Figure 267213DEST_PATH_IMAGE027
例如至少指示位点位置，单元长度，单元数。
例如，计算设备110按照输入的抽样矩阵，针对短串联重复序列位点数据进行均匀抽样，然后输出抽样得到的STR位点表格
Figure 625513DEST_PATH_IMAGE027
。
在步骤306处，计算设备110基于比对结果数据和短串联重复序列位点抽样信息，经由多线程池，分别计算每个短串联重复序列位点的深度和***缺失数目，以便生成短串联重复序列位点特征数据。
例如，计算设备110创建一些空闲的线程（它们的集合称为线程池），然后将关于计算每个短串联重复序列位点的深度和***缺失数目的多个任务传给线程池，线程池就会对应启动多个线程来执行该多个任务，相应线程执行结束以后，该线程再次返回线程池中成为空闲状态，等待执行下一个或多个任务。本发明通过采用线程池的手段，可以避免频繁切换线程的引起的延迟，提升了运行速度。
经研究发现，生成STR位点抽样信息是预操作的过程，针对特定参考基因组序列数据，仅需执行一次计算，而生成STR位点特征数据过程需要针对每个STR位点进行计算，是方法300中最为耗时的环节，因此，本发明通过线程池的方式处理生成STR位点特征数据的步骤，可以显著地提升运行速度生成短串联重复序列位点特征数据的速度。
关于短串联重复序列位点特征数据，其例如而不限于是STR位点特征值表
Figure 287963DEST_PATH_IMAGE028
。
例如，计算设备110针对待测样本的比对结果数据和步骤304处所获得的STR位点表格
Figure 339095DEST_PATH_IMAGE027
，计算STR位点表格
Figure 451277DEST_PATH_IMAGE027
中每个STR位点的深度和***缺失数目，以便生成STR位点特征值表
Figure 613268DEST_PATH_IMAGE028
。STR位点特征值表
Figure 861715DEST_PATH_IMAGE028
相当于在STR位点表格
Figure 83749DEST_PATH_IMAGE027
的基础上，新增了深度数据和***缺失数目数据。
在步骤308处，计算设备110基于所生成的短串联重复序列位点特征数据，计算短串联重复序列位点特征矩阵。例如，计算设备110根据STR位点特征值表
Figure 683227DEST_PATH_IMAGE028
计算STR特征矩阵。
关于短串联重复序列位点特征矩阵，其例如包括：开空位概率矩阵（或简称为“GOP矩阵”）、 空位延伸概率矩阵（或简称为“GCP矩阵”）、INDEL先验概率矩阵（或简称为“API矩阵”）。短串联重复序列位点特征矩阵的行数为最大短串联重复序列单元长度，短串联重复序列位点特征矩阵的列数为最大短串联重复序列单元数。在一些实施例中，开空位概率矩阵，空位延伸概率矩阵、INDEL先验概率矩阵都是二维矩阵，并且形状一致，行数为最大STR单元长度，例如为8行，列数为最大STR单元数，例如为20列。
关于计算短串联重复序列位点特征矩阵的方法，其例如包括：确定所有短串联重复序列位点的深度和***缺失数目是否计算完毕；以及响应于确定所有短串联重复序列位点的深度和***缺失数目计算完毕，基于短串联重复序列位点特征数据，计算短串联重复序列位点特征矩阵。通过采用上述手段，可以保证所有并行计算结束后才计算短串联重复序列位点特征矩阵。
在上述方案中，本发明可以显著地提高生成短串联重复序列位点特征矩阵的速度，进而有利于进一步提升检测单核苷酸变异和***缺失的速度。
以下将结合图4描述根据本发明的实施例的用于生成激活区域的方法。图4示出了根据本发明的实施例的用于生成激活区域的方法400的流程图。应当理解，方法400例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解，方法400还可以包括未示出的附加动作和/或可以省略所示出的动作，本发明的范围在此方面不受限制。
在步骤402处，计算设备110配置用于存储激活概率的数组，数组的长度等于待检测染色体长度。例如，计算设备110针对每一个待检测染色体，分配一个长度为待检测染色体长度的数组A，以便用于存储激活概率。
在步骤404处，计算设备110将待检测染色体分为预定长度的多个区间。例如，计算设备110将待检测染色体分为长度为w（例如w=1000000）的多个区间。
在步骤406处，计算设备110并行地计算多个区间中的每个区间内基因组位点的激活概率，以便将所计算中的激活概率存储至数组。例如，计算设备110使用多线程技术针对各个区间并行地计算区间内基因组位点的激活概率。例如，计算设备110基于比对序列堆叠（pileup）的贝叶斯方法来计算基因组位点激活概率。
关于并行地计算多个区间中的每个区间内基因组位点的激活概率的方法，其例如包括：计算设备110确定当前比对序列堆叠的备选等位基因的读长数目是否大于或者等于预定读长数目阈值；响应于确定当前比对序列堆叠的备选等位基因的读长数目大于或者等于预定读长数目阈值，将当前比对序列堆叠识别为候选比对序列堆叠；响应于确定当前比对序列堆叠的备选等位基因的读长数目小于预定读长数目阈值，将当前比对序列堆叠位点的激活概率确定为零；以及使用贝叶斯算法计算候选比对序列堆叠的激活概率。
关于预定读长数目阈值，其例如而不限于是2。例如，如果计算设备110确定当前比对序列堆叠（pileup）的备选等位基因（alternative allele）reads数目大于或等于2，那么确定当前pileup为候选pileup，否则将当前pileup位点的激活概率设置为0。图5示出了根据本发明的实施例的比对序列堆叠的示意图。如图5所示，第1行为参考基因组序列，第2-9行为读长（reads）序列，每一列为一个pileup。标记510所指示的第7列的pileup，其中C碱基为参考等位基因（reference allele），T碱基为备选等位基因（alternative allele）。因此，计算设备110确定标记510所指示的第7列的pileup的备选等位基因reads数目大于或等于2，因此，可以确定标记510所指示的第7列的pileup为候选pileup。其他列的pileup，由于没有备选等位基因reads，因此，其他pileup位点的激活概率为0。关于使用贝叶斯算法计算候选比对序列堆叠的激活概率的方法，下文将结合公式（11）至（15）进行说明。
以下结合公式（11）说明原始基因型
Figure 648909DEST_PATH_IMAGE029
的概率
Figure 20372DEST_PATH_IMAGE030
的计算方法。
Figure 147728DEST_PATH_IMAGE031
在上述公式（11）中，
Figure 234501DEST_PATH_IMAGE029
代表基因型。
Figure 738295DEST_PATH_IMAGE032
和
Figure 226914DEST_PATH_IMAGE033
代表等位基因。
Figure 525171DEST_PATH_IMAGE034
代表第i条读长。
Figure 833662DEST_PATH_IMAGE035
和
Figure 406726DEST_PATH_IMAGE036
代表等位基因
Figure 221622DEST_PATH_IMAGE032
和
Figure 221939DEST_PATH_IMAGE033
的概率。
Figure 17726DEST_PATH_IMAGE035
和
Figure 128901DEST_PATH_IMAGE036
可以根据下文的公式（13）计算。
Figure 326533DEST_PATH_IMAGE037
代表基因型
Figure 966593DEST_PATH_IMAGE029
的概率。
以下结合公式（12）说明基因型
Figure 515255DEST_PATH_IMAGE029
的后验概率
Figure 164542DEST_PATH_IMAGE038
的计算方法。
Figure 231329DEST_PATH_IMAGE039
在上述公式（12）中，
Figure 307870DEST_PATH_IMAGE040
代表基因型
Figure 78249DEST_PATH_IMAGE029
的先验概率。例如，杂合SNP的先验概率为
Figure 265647DEST_PATH_IMAGE041
，杂合INDEL的先验概率为
Figure 437872DEST_PATH_IMAGE042
，纯合Reference的先验概率为1。
Figure 685313DEST_PATH_IMAGE038
代表基因型
Figure 942988DEST_PATH_IMAGE029
的后验概率。纯合SNP或INDEL的先验概率可以根据哈迪温伯格平衡或杂合纯合比（默认值为2）进行计算。
应当理解，
Figure 934078DEST_PATH_IMAGE043
通过碱基质量和预定经验值计算。以下结合公式（13）说明SNP
Figure 963738DEST_PATH_IMAGE043
的计算方式。
Figure 382081DEST_PATH_IMAGE044
在上述公式（13）中，p代表比对序列堆叠（pileup）在第i条读长
Figure 861473DEST_PATH_IMAGE045
的位置，
Figure 390675DEST_PATH_IMAGE046
为第i条读长
Figure 537491DEST_PATH_IMAGE045
的碱基质量值。
以下结合公式（14）说明于INDEL
Figure 861156DEST_PATH_IMAGE043
的计算方式。
Figure 359002DEST_PATH_IMAGE047
在上述公式（14）中，
Figure 691895DEST_PATH_IMAGE048
代表预定经验值
Figure 430568DEST_PATH_IMAGE049
。
以下结合公式（15）说明激活概率值的计算方式。
Figure 456293DEST_PATH_IMAGE050
在上述公式（15）中，
Figure 910277DEST_PATH_IMAGE051
代表纯合reference的基因型。
Figure 46860DEST_PATH_IMAGE052
代表非纯合reference的基因型概率，即激活概率。
在步骤408处，计算设备110使用高斯核对整条待测染色体上各基因组位点的激活概率进行平滑，以便将激活概率大于或者等于预定阈值的连续基因组位点合并生成激活区域。关于预定阈值，其例如是阈值k。
通过采用上述手段，本发明能够提高检测单核苷酸变异和***缺失的快速性。
以下将结合图6描述根据本发明的实施例的用于获取单核苷酸变异和***缺失特征的方法。图6示出了根据本发明的实施例的用于获取单核苷酸变异和***缺失特征的方法600的流程图。应当理解，方法600例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解，方法600还可以包括未示出的附加动作和/或可以省略所示出的动作，本发明的范围在此方面不受限制。
在步骤602处，计算设备110将每个激活区向两侧各延伸预定数量个碱基，以便得到对应的延伸区间。
例如，计算设备110将每个激活区域R0向两侧各延伸n个碱基（例如，n=100），以便得到延伸区间R1。应当理解，通过对每个激活区向两侧延伸，再针对延伸之后的延伸区间并行地检测单核苷酸变异和***缺失，本发明可以解决激活区向边缘处难以检测单核苷酸变异和***缺失的问题，进而提高检测单核苷酸变异和***缺失的准确性。
在步骤604处，计算设备110针对每个延伸区间内所有的读长进行局部组装，以便生成经组装的单倍型序列。
例如，计算设备110使用de Bruijn Graph算法对延伸区间R1内所有的reads进行局部组装，以便输出经组装的单倍型序列。
在步骤606处，计算设备110将单倍型序列比对到参考基因组，以便检测单核苷酸变异和***缺失。
例如，计算设备110使用Smith-Waterman局部比对算法将单倍型序列比对到参考基因组，以便并检测SNP、INDEL。
在步骤608处，计算设备110保留在激活区内的、所检测的单核苷酸变异和***缺失。应当理解，由于延伸区域会与其他区域存在一些干涉，因此有必要只保留在原始输入的激活区内的单核苷酸变异和***缺失。
在步骤610处，计算设备110计算读长比对到单倍型序列的概率。例如，计算设备110基于读长、单倍型序列、读长上每个碱基的质量值、读长上每个碱基的开空位概率（即，Gap open概率）、读长上每个碱基的空位延伸概率（即，Gap extending概率），计算读长比对到单倍型序列的概率。在一些实施例中，备110使用双隐马尔可夫模型（PairHMM）计算每条读长比对到单倍型序列的概率。
应当理解，每个单倍型序列与步骤606处所检测的单核苷酸变异和***缺失存在对应关系，因此，步骤610处所计算的“每条读长比对到单倍型序列的概率”可以映射到“每条读长对应到单核苷酸变异位点和***缺失的概率”。
具体而言，例如，计算设备110可以通过PairHMM使用向前算法计算reads比对到单倍型序列的概率
Figure 637110DEST_PATH_IMAGE053
。其中，PairHMM的输入为读长
Figure 568157DEST_PATH_IMAGE054
、单倍型序列
Figure 509437DEST_PATH_IMAGE055
、读长上每个碱基的质量值
Figure 184132DEST_PATH_IMAGE056
、读长上每个碱基开开空位概率（即，Gap open概率）
Figure 897398DEST_PATH_IMAGE057
、读长上每个碱基的延伸概率（即，Gap extending概率）
Figure 264925DEST_PATH_IMAGE058
。其中，开空位概率（即，Gap open概率）
Figure 427922DEST_PATH_IMAGE057
、读长上每个碱基的延伸概率（即，Gap extending概率）
Figure 171887DEST_PATH_IMAGE058
可以通过上文提及的开空位概率矩阵和空位延伸概率矩阵获取。
在步骤612处，计算设备110计算单核苷酸变异位点和***缺失位点的基因型概率，以便获取单核苷酸变异和***缺失特征值。例如，计算设备110使用贝叶斯方法计算SNP、INDEL位点的基因型概率，以用于收集SNP、INDEL的特征值。
关于单核苷酸变异和***缺失特征，其例如是特征向量矩阵，每个单核苷酸变异或***缺失变异对应于特征值矩阵中的一个特征值。
通过采用上述手段，本发明能够进一步提高检测单核苷酸变异和***缺失的速度。
图8示意性示出了适于用来实现本发明实施例的电子设备800的框图。电子设备800可以是用于实现执行图2至图4、图6所示的方法200至400、600。如图8所示，电子设备800包括中央处理单元（即，CPU 801），其可以根据存储在只读存储器（即，ROM 802）中的计算机程序指令或者从存储单元808加载到随机访问存储器（即，RAM 803）中的计算机程序指令，来执行各种适当的动作和处理。在RAM 803中，还可存储电子设备800操作所需的各种程序和数据。CPU 801、ROM 802以及RAM 803通过总线804彼此相连。输入/输出接口（即，I/O接口805）也连接至总线804。
电子设备800中的多个部件连接至I/O接口805，包括：输入单元806、输出单元807、存储单元808，CPU 801执行上文所描述的各个方法和处理，例如执行方法200至400、600。例如，在一些实施例中，方法200至400、600可被实现为计算机软件程序，其被存储于机器可读介质，例如存储单元808。在一些实施例中，计算机程序的部分或者全部可以经由ROM 802和/或通信单元809而被载入和/或安装到电子设备800上。当计算机程序加载到RAM 803并由CPU 801执行时，可以执行上文描述的方法200至400、600的一个或多个操作。备选地，在其他实施例中，CPU 801可以通过其他任何适当的方式（例如，借助于固件）而被配置为执行方法200至400、600的一个或多个动作。
需要进一步说明的是，本发明可以是方法、装置、***和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于执行本发明的各个方面的计算机可读程序指令。
计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子（非穷举的列表）包括：便携式计算机盘、硬盘、随机存取存储器（RAM）、只读存储器（ROM）、可擦式可编程只读存储器（EPROM或闪存）、静态随机存取存储器（SRAM）、便携式压缩盘只读存储器（CD-ROM）、数字多功能盘（DVD）、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波（例如，通过光纤电缆的光脉冲）、或者通过电线传输的电信号。
这里所描述的计算机可读程序指令可以从计算机可读存储介质下载到各个计算/处理设备，或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收计算机可读程序指令，并转发该计算机可读程序指令，以供存储在各个计算/处理设备中的计算机可读存储介质中。
用于执行本发明操作的计算机程序指令可以是汇编指令、指令集架构（ISA）指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码，该编程语言包括面向对象的编程语言—诸如Smalltalk、C++等，以及常规的过程式编程语言—诸如“C”语言或类似的编程语言。计算机可读程序指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络—包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机（例如利用因特网服务提供商来通过因特网连接）。在一些实施例中，通过利用计算机可读程序指令的状态信息来个性化定制电子电路，例如可编程逻辑电路、现场可编程门阵列（FPGA）或可编程逻辑阵列（PLA），该电子电路可以执行计算机可读程序指令，从而实现本发明的各个方面。
这里参照根据本发明实施例的方法、设备（***）、和计算机程序产品的流程图和/或框图描述了本发明的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。
这些计算机可读程序指令可以提供给语音交互装置中的处理器、通用计算机、专用计算机或其它可编程数据处理装置的处理单元，从而生产出一种机器，使得这些指令在通过计算机或其它可编程数据处理装置的处理单元执行时，产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中，这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作，从而，存储有指令的计算机可读介质则包括一个制造品，其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。
也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上，使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤，以产生计算机实现的过程，从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。
附图中的流程图和框图显示了根据本发明的多个实施例的设备、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分，该模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的***来实现，或者可以用专用硬件与计算机指令的组合来实现。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。
以上仅为本发明的可选实施例，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等效替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种用于检测单核苷酸变异和***缺失的方法，其特征在于，包括：
获取待测样本的测序序列数据与参考基因组序列数据的比对结果数据和参考基因组序列数据；
基于比对结果数据和参考基因组序列数据，生成短串联重复序列位点特征矩阵；
计算待测样本的检测染色体上各个基因组位点的激活概率，以便确定多个激活区域；
经由多线程技术，并行地检测多个激活区域中的每个激活区域的单核苷酸变异和***缺失，以便获取单核苷酸变异和***缺失特征值；以及
基于短串联重复序列位点特征矩阵和单核苷酸变异和***缺失特征值，经由预测模型生成关于单核苷酸变异和***缺失的预测结果；
其中，生成短串联重复序列位点特征矩阵包括：
基于参考基因组序列数据，查找参考基因组上的短串联重复序列位点数据；
针对短串联重复序列位点数据进行均匀抽样，以便获得短串联重复序列位点抽样信息，所述短串联重复序列位点抽样信息指示短串联重复序列位点的位点位置、单元长度和单元数；
基于比对结果数据和短串联重复序列位点抽样信息，经由多线程池，分别计算每个短串联重复序列位点的深度和***缺失数目，以便生成短串联重复序列位点特征数据；
确定所有短串联重复序列位点的深度和***缺失数目是否被计算完毕；以及
响应于确定所有短串联重复序列位点的深度和***缺失数目被计算完毕，基于短串联重复序列位点特征数据，计算短串联重复序列位点特征矩阵，短串联重复序列位点特征矩阵包括：开空位概率矩阵、空位延伸概率矩阵和***缺失先验概率矩阵。
2.根据权利要求1所述的方法，其特征在于，并行地检测多个激活区域中的每个激活区域的单核苷酸变异和***缺失包括：
将每个激活区向两侧各延伸预定数量个碱基，以便得到对应的延伸区间；
针对每个延伸区间内所有的读长进行局部组装，以便并行地检测单核苷酸变异和***缺失，以获取单核苷酸变异和***缺失特征值。
3.根据权利要求1所述的方法，其特征在于，计算待测样本的检测染色体上各个基因组位点的激活概率，以便确定多个激活区域包括：
配置用于存储激活概率的数组，所述数组的长度等于待检测染色体的长度；
将待检测染色体分为预定长度的多个区间；
并行地计算多个区间中的每个区间内基因组位点的激活概率，以便将所计算的激活概率存储至所述数组；以及
使用高斯核对整条待测染色体上的各基因组位点的激活概率进行平滑，以便将激活概率大于或者等于预定阈值的连续基因组位点合并生成激活区域。
4.根据权利要求3所述的方法，其特征在于，并行地计算多个区间中的每个区间内基因组位点的激活概率包括：
确定当前比对序列堆叠的备选等位基因的读长数目是否大于或者等于预定读长数目阈值；
响应于确定当前比对序列堆叠的备选等位基因的读长数目大于或者等于预定读长数目阈值，将当前比对序列堆叠识别为候选比对序列堆叠；
响应于确定当前比对序列堆叠的备选等位基因的读长数目小于预定读长数目阈值，将当前比对序列堆叠位点的激活概率确定为零；以及
使用贝叶斯算法计算候选比对序列堆叠的激活概率。
5.根据权利要求2所述的方法，其特征在于，针对每个延伸区间内所有的读长进行局部组装以便并行地检测单核苷酸变异和***缺失，以获取单核苷酸变异和***缺失特征值包括：
针对每个延伸区间内所有的读长进行局部组装，以便生成经组装的单倍型序列；
将单倍型序列比对到参考基因组，以便检测单核苷酸变异和***缺失；
保留在激活区内的、所检测的单核苷酸变异和***缺失；
计算读长比对到单倍型序列的概率；以及
计算单核苷酸变异位点和***缺失位点的基因型概率，以便获取单核苷酸变异和***缺失特征值。
6.根据权利要求5所述的方法，其特征在于，计算读长比对到单倍型序列的概率包括：
基于读长、单倍型序列、读长上每个碱基的质量值、读长上每个碱基的开空位概率、读长上每个碱基的空位延伸概率，计算读长比对到单倍型序列的概率。
7.一种计算设备，其特征在于，包括：
至少一个处理单元；
至少一个存储器，所述至少一个存储器被耦合到所述至少一个处理单元并且存储用于由所述至少一个处理单元执行的指令，所述指令当由所述至少一个处理单元执行时，使得所述设备执行根据权利要求1至6任一项所述的方法的步骤。
8.一种计算机可读存储介质，其特征在于，计算机可读存储介质上存储有计算机程序，所述计算机程序被机器执行时实现根据权利要求1至6中任一项所述的方法。
CN202210392298.XA 2022-04-15 2022-04-15 用于检测单核苷酸变异和***缺失的方法、设备和介质 Active CN114649055B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210392298.XA CN114649055B (zh) 2022-04-15 2022-04-15 用于检测单核苷酸变异和***缺失的方法、设备和介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210392298.XA CN114649055B (zh) 2022-04-15 2022-04-15 用于检测单核苷酸变异和***缺失的方法、设备和介质

Publications (2)

Publication Number Publication Date
CN114649055A CN114649055A (zh) 2022-06-21
CN114649055B true CN114649055B (zh) 2022-10-21

Family

ID=81997268

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210392298.XA Active CN114649055B (zh) 2022-04-15 2022-04-15 用于检测单核苷酸变异和***缺失的方法、设备和介质

Country Status (1)

Country Link
CN (1) CN114649055B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114974416B (zh) * 2022-07-15 2023-04-07 深圳雅济科技有限公司 一种检测相邻多核苷酸变异的方法及装置
CN115831219B (zh) * 2022-12-22 2024-05-28 郑州思昆生物工程有限公司 一种质量预测方法、装置、设备及存储介质
CN117935921B (zh) * 2024-03-21 2024-06-11 北京贝瑞和康生物技术有限公司 确定缺失/重复类型的方法、设备、介质和程序产品

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1381591A (zh) * 2001-03-30 2002-11-27 珀尔根科学公司 基因组分析方法
CN101932729A (zh) * 2007-12-05 2010-12-29 考利达基因组股份有限公司 测序反应中碱基的有效确定
CN109416928A (zh) * 2016-06-07 2019-03-01 伊路米纳有限公司 用于进行二级和/或三级处理的生物信息学***、设备和方法
CN109906276A (zh) * 2016-11-07 2019-06-18 格里尔公司 用于检测早期癌症中体细胞突变特征的识别方法
CN111292802A (zh) * 2020-02-03 2020-06-16 至本医疗科技（上海）有限公司 用于检测突变的方法、电子设备和计算机存储介质
CN111462816A (zh) * 2020-03-31 2020-07-28 至本医疗科技（上海）有限公司 用于检测胚系基因微缺失微重复的方法、电子设备和计算机存储介质
CN112687332A (zh) * 2021-03-12 2021-04-20 北京贝瑞和康生物技术有限公司 用于确定致病风险变异位点的方法、设备和存储介质
CN113168886A (zh) * 2018-08-13 2021-07-23 豪夫迈·罗氏有限公司 用于使用神经网络进行种系和体细胞变体调用的***和方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6480791B1 (en) * 1998-10-28 2002-11-12 Michael P. Strathmann Parallel methods for genomic analysis
EP1963524A2 (en) * 2005-12-01 2008-09-03 University of Bergen Diagnosis and treatment of exocrine pancreatic dysfunction and diabetes using human cel gene
CN101539967B (zh) * 2008-12-12 2010-12-01 深圳华大基因研究院 一种单核苷酸多态性检测方法
GB2506523A (en) * 2012-08-31 2014-04-02 Real Time Genomics Inc A computerised assignment of genomic sequence values based on multiple reads and probabilistic analysis
WO2017087724A1 (en) * 2015-11-17 2017-05-26 Omniome, Inc. Methods for determining sequence profiles
CN111584002B (zh) * 2020-05-22 2022-04-29 至本医疗科技（上海）有限公司 用于检测肿瘤突变负荷的方法、计算设备和计算机存储介质

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1381591A (zh) * 2001-03-30 2002-11-27 珀尔根科学公司 基因组分析方法
CN101932729A (zh) * 2007-12-05 2010-12-29 考利达基因组股份有限公司 测序反应中碱基的有效确定
CN109416928A (zh) * 2016-06-07 2019-03-01 伊路米纳有限公司 用于进行二级和/或三级处理的生物信息学***、设备和方法
CN109906276A (zh) * 2016-11-07 2019-06-18 格里尔公司 用于检测早期癌症中体细胞突变特征的识别方法
CN113168886A (zh) * 2018-08-13 2021-07-23 豪夫迈·罗氏有限公司 用于使用神经网络进行种系和体细胞变体调用的***和方法
CN111292802A (zh) * 2020-02-03 2020-06-16 至本医疗科技（上海）有限公司 用于检测突变的方法、电子设备和计算机存储介质
CN111462816A (zh) * 2020-03-31 2020-07-28 至本医疗科技（上海）有限公司 用于检测胚系基因微缺失微重复的方法、电子设备和计算机存储介质
CN112687332A (zh) * 2021-03-12 2021-04-20 北京贝瑞和康生物技术有限公司 用于确定致病风险变异位点的方法、设备和存储介质

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
Computational analysis of cancer genome sequencing data;Isidro Cortes-Ciriano 等;《Nature Reviews Genetics》;20211208(第23期);第298-314页 *
Study on genotyping and resistance evaluation for five rice blast resistance genes：web_of_science;Zhangyu 等;《Southwest China Journal of Agricultural Science》;20141231;第27卷(第5期);第1929-1936页 *
基于短序列比对的InDel检测算法研究;汪晓丹;《中国优秀硕士学位论文全文数据库 基础科学辑》;20170315;第2017年卷(第03期);第A006-432页 *
细菌VI型分泌***结构与调控的研究进展;姚丰华 等;《中国动物传染病学报》;20131231;第21卷(第6期);第80-86页 *
胰岛素受体底物-1基因5′-调控序列CAG富含区碱基缺失/***的研究;黄建军等;《湖南医科大学学报》;20010430(第02期);第103-106页 *

Also Published As

Publication number Publication date
CN114649055A (zh) 2022-06-21

Similar Documents

Publication Publication Date Title
CN114649055B (zh) 用于检测单核苷酸变异和***缺失的方法、设备和介质
Meisner et al. Inferring population structure and admixture proportions in low-depth NGS data
Pavlidis et al. SweeD: likelihood-based detection of selective sweeps in thousands of genomes
CN111292802B (zh) 用于检测突变的方法、电子设备和计算机存储介质
US20220223233A1 (en) Display of estimated parental contribution to ancestry
Browning et al. Haplotype phasing: existing methods and new developments
US20210082539A1 (en) Gene mutation identification method and apparatus, and storage medium
KR20190101966A (ko) 범-암 게놈에서 dna 접근성을 예측하기 위한 방법 및 시스템
CN114496077B (zh) 用于检测单核苷酸变异和***缺失的方法、设备和介质
Huang et al. Evaluation of variant detection software for pooled next-generation sequence data
Revell On the analysis of evolutionary change along single branches in a phylogeny
KR20220069943A (ko) 단일 세포 rna-seq 데이터 처리
Flassig et al. An effective framework for reconstructing gene regulatory networks from genetical genomics data
Pankratov et al. Prioritizing autoimmunity risk variants for functional analyses by fine-mapping mutations under natural selection
Duchemin et al. Evaluation of methods to detect shifts in directional selection at the genome scale
Zheng et al. Adaptation in structured populations and fuzzy boundaries between hard and soft sweeps
Satija et al. Combining statistical alignment and phylogenetic footprinting to detect regulatory elements
KR101770962B1 (ko) 유전자 서열 기반 개인 마커에 관한 정보를 제공하는 방법 및 이를 이용한 장치
Sitarčík et al. WarpSTR: determining tandem repeat lengths using raw nanopore signals
CN110570908B (zh) 测序序列多态识别方法及装置、存储介质、电子设备
Huang et al. Inferring sequence regions under functional divergence in duplicate genes
US11177018B2 (en) Stable genes in comparative transcriptomics
KR20200135221A (ko) Ngs 데이터를 이용하여 유전형을 예측하는 방법 및 장치
CN110021342B (zh) 用于加速变异位点的识别的方法及***
Wong et al. SpliceWiz: easy, optimized, and accurate alternative splicing analysis in R

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
GR01 Patent grant
GR01 Patent grant
REG Reference to a national code

Ref country code: HK

Ref legal event code: DE

Ref document number: 40075851

Country of ref document: HK