CN114758720B - 用于检测拷贝数变异的方法、设备和介质 - Google Patents
用于检测拷贝数变异的方法、设备和介质 Download PDFInfo
- Publication number
- CN114758720B CN114758720B CN202210664065.0A CN202210664065A CN114758720B CN 114758720 B CN114758720 B CN 114758720B CN 202210664065 A CN202210664065 A CN 202210664065A CN 114758720 B CN114758720 B CN 114758720B
- Authority
- CN
- China
- Prior art keywords
- window
- determining
- missing
- type
- window position
- 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
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/10—Ploidy or copy number 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
- 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
Landscapes
- Bioinformatics & Cheminformatics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Biotechnology (AREA)
- Biophysics (AREA)
- Genetics & Genomics (AREA)
- Molecular Biology (AREA)
- Chemical & Material Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Analytical Chemistry (AREA)
- Evolutionary Biology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及一种用于检测拷贝数变异的方法、设备和介质。该方法包括:基于比对结果数据,计算每一条读长在参考基因组中的位置;将基因组按照预定尺寸的窗口进行划分,以便统计每一个窗口内的唯一比对数;针对每一个窗口内的唯一比对数进行预处理;基于经由预处理的唯一比对数,分别针对每条染色体上的可能添加的断点,进行第一类型片段拟合,以便确定所划分的第一类型片段和关于第一类型片段的断点集合;以及针对所划分的第一类型片段,计算每个窗口的重复表征数据和缺失表征数据,以便确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异。本发明针对低深度全基因组测序,显著提高检测拷贝数变异的准确性。
Description
技术领域
本发明总体上涉及生物信息处理,并且具体地,涉及用于检测拷贝数变异的方法、计算设备和计算机存储介质。
背景技术
随着测序成本快速下降,测序技术逐渐从实验室走向了临床。例如,低深度全基因组测序技术可以作为一线产前诊断技术,对可能存在胎儿染色体异常的孕妇进行产前诊断。在实践中,低深度全基因组测序技术往往需要搭配优质地拷贝数变异(CNV)鉴定算法,才能更好地服务于临床。拷贝数变异(copy number variation,CNV)是指染色体上大于1kb的DNA片段的增加或者减少,主要表现为亚显微水平的缺失和重复。 应当理解,由于CNV的种类、片段大小等千差万别,且观测测序数据可能受到多方因素(例如,基因组GC含量、重复区域等因素)影响,进而导致算法处理某些特殊的CNV片段不能达到理想的效果。因此,一套能够高准确地同时鉴定各自不同类型的CNV(例如,能够准确鉴定长短片段重复/缺失/嵌合、染色体重复/缺失/嵌合、高GC区域CNV等)的算法是至关重要的。
在传统的用于检测拷贝数变异的方法中,由于CNV可能出现的情况复杂,因此会在一个或几个方面存在缺陷,例如,对于低比例的嵌合鉴定不够准确,或者对于纯合缺失或多拷贝重复出错等等。
综上,传统的用于检测拷贝数变异的方案存在的不足之处在于:无法区分正常和低比例嵌合,单独针对低深度全基因组测序的准确性不高。
发明内容
本发明提供一种用于检测拷贝数变异的方法、计算设备和计算机存储介质,针对低深度全基因组测序,显著提高检测拷贝数变异的准确性。
根据本发明的第一方面,提供了一种用于检测拷贝数变异的方法。该方法包括:基于待测样本的测序数据与参考基因组序列的比对结果数据,计算每一条读长在参考基因组中的位置;将基因组按照预定尺寸的窗口进行划分,以便统计每一个窗口内的唯一比对数;针对每一个窗口内的唯一比对数进行预处理;基于经由预处理的唯一比对数,分别针对每条染色体上的可能添加的断点,进行第一类型片段拟合,以便确定所划分的第一类型片段和关于第一类型片段的断点集合;以及针对所划分的第一类型片段,计算每个窗口的重复表征数据和缺失表征数据,以便确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异,第一类型片段的尺寸大于预定阈值。
根据本发明的第二方面,还提供了一种计算设备,该设备包括:存储器,被配置为存储一个或多个计算机程序;以及处理器,耦合至存储器并且被配置为执行一个或多个程序使装置执行本发明的第一方面的方法。
根据本发明的第三方面,还提供了一种非瞬态计算机可读存储介质。该非瞬态计算机可读存储介质上存储有机器可执行指令,该机器可执行指令在被执行时使机器执行本发明的第一方面的方法。
在一些实施例中,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:针对所划分的第一类型片段,将每个窗口的重复表征数据进行累加,以生成连续窗口的重复表征数据累计值;确定连续窗口的重复表征数据累计值是否大于0;响应于确定连续窗口的重复表征数据累计值大于0,在连续窗口的重复表征数据累计值上继续累加下一窗口的重复表征数据,直至累加之后的连续窗口的重复表征数据累计值大于重复阈值;以及将重复表征数据累计值大于重复阈值时的连续窗口所在片段确定为重复第二类型片段。
在一些实施例中,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:针对所划分的第一类型片段,将每个窗口的缺失数据表征数据分别进行累加,以生成连续窗口的缺失表征数据累计值;确定连续窗口的缺失表征数据累计值是否大于0;响应于确定连续窗口的缺失表征数据累计值大于0,在连续窗口的缺失表征数据累计值上继续累加下一窗口的缺失表征数据,直至累加之后的连续窗口的缺失表征数据累计值大于缺失阈值;以及将缺失表征数据累计值大于缺失阈值时的连续窗口所在片段确定为缺失第二类型片段。
在一些实施例中,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:基于所确定的重复第二类型片段,确定关于重复第二类型片段的断点集合;基于所确定的缺失第二类型片段,确定关于缺失第二类型片段的断点集合;针对关于重复第二类型片段的断点集合和关于缺失第二类型片段的断点集合进行去重,以便基于去重后的断点集合确定拷贝数变异。
在一些实施例中,确定所划分的第一类型片段和关于第一类型片段的断点集合包括:将染色体上所有窗口作为一个整体片段,设置断点集合为空;在断点集合的基础上,遍历所有其他可能添加断点的位置,以便依次构建临时断点集合;分别基于断点集合和临时断点集合将染色体分为多个片段,以便计算每个片段的平均拷贝数;针对断点集合和临时断点集合,分别计算每个窗口的拷贝数、每个窗口的拷贝数到所在片段的平均拷贝数的距离,以便将所计算的距离的平均值作为误差;确定针对断点集合的初始误差和临时断点集合的最小误差,以便确定初始误差与最小误差之间的差值是否小于第一类型片段拟合距离阈值;响应于确定初始误差与最小误差之间的差值小于第一类型片段拟合距离阈值,以断点集合进行片段的划分;以及响应于确定初始误差与最小误差之间的差值大于或者等于第一类型片段拟合距离阈值,以与最小误差对应的临时断点集合替换断点集合,以便在替换后的断点集合上依次构建临时断点集合。
在一些实施例中,针对每一个窗口内的唯一比对数进行预处理包括:针对每一个窗口内的唯一比对数进行归一化处理;基于归一化后的每一个窗口内的唯一比对数,计算关于待测样本的Y染色体唯一比对数的占比,以便确定关于待测样本的所属性别;分别构建常染色体的阴性参考集和与所确定的性别相对应的性染色体的阴性参考集;分别对每个样本上所有的窗口进行GC矫正;以及基于矫正后的每个窗口的唯一比对数和所构建的阴性参考集,确定拷贝数的观测值。
在一些实施例中,针对每一个窗口内的唯一比对数进行归一化处理包括:基于单条染色体最大窗口数目、样本数量、染色体数目,确定归一化比例;基于归一化比例和每一个窗口内的唯一比对数,计算每一个窗口内的归一化后的唯一比对数;确定当前窗口内的归一化后的唯一比对数是否小于唯一比对数平均值的预定比例;以及响应于确定当前窗口内的归一化后的唯一比对数小于唯一比对数平均值的预定比例,确定当前窗口为检测盲区;以及将当前窗口内的归一化后的唯一比对数替换为空缺罚分值。
在一些实施例中,确定关于待测样本的所属性别包括:计算当前待测样本的Y染色体唯一比对数的占比;确定所计算的Y染色体唯一比对数的占比是否小于或者等于预定占比阈值;响应于确定所计算的Y染色体唯一比对数的占比小于或者等于预定占比阈值,确定当前样本所属的性别为女性;以及响应于确定所计算的Y染色体唯一比对数的占比大于预定占比阈值,确定当前样本所属的性别为男性。
在一些实施例中,分别构建常染色体的阴性参考集和与所确定的性别相对应的性染色体的阴性参考集包括:分别计算每个窗口内的唯一比对数的均值和标准差;基于所计算的均值和标准差,针对每个窗口,计算保留区间;以及基于每个窗口的归一化后的唯一比对数与保留区间的比较,确定针对每个样本的、每窗口的保留系数,以便分别构建常染色体的阴性参考集和与所确定的性别相对应的性染色体的阴性参考集。
在一些实施例中,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:针对经由第一类型片段拟合后的染色体,确定起始窗口位置、断点集合,以便计算起始窗口位置的重复表征数据累计值和重复阈值;响应于确定起始窗口位置的重复表征数据累计值大于0,使得起始窗口位置叠加至下一窗口位置,以便获得自起始窗口位置至下一窗口位置的重复表征数据累计值;响应于确定自起始窗口位置至下一窗口位置的重复表征数据累计值大于0,确定自起始窗口位置至下一窗口位置的重复表征数据累计值是否大于重复阈值;响应于确定自起始窗口位置至下一窗口位置的重复表征数据累计值大于重复阈值,将自起始窗口位置至下一窗口位置的重复表征数据累计值作为重复阈值;确定下一窗口位置是否为第一类型片段的最终窗口位置;响应于确定下一窗口位置为第一类型片段的最终窗口位置,确定重复阈值是否大于第二类型片段阈值;以及响应于确定重复阈值大于第二类型片段阈值,将起始窗口位置和重复阈值处对应窗口添加入断点集合,以便形成关于重复第二类型片段的断点集合和重复第二类型片段。
在一些实施例中,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:针对经由第一类型片段拟合后的染色体,确定起始窗口位置、断点集合,以便计算起始窗口位置的缺失表征数据累计值和缺失阈值;响应于确定起始窗口位置的缺失表征数据累计值大于0,使得起始窗口位置叠加至下一窗口位置,以便获得自起始窗口位置至下一窗口位置的缺失表征数据累计值;响应于确定自起始窗口位置至下一窗口位置的缺失表征数据累计值大于0,确定自起始窗口位置至下一窗口位置的缺失表征数据累计值是否大于缺失阈值;响应于确定自起始窗口位置至下一窗口位置的缺失表征数据累计值大于缺失阈值,将自起始窗口位置至下一窗口位置的缺失表征数据累计值作为缺失阈值;确定下一窗口位置是否为第一类型片段的最终窗口位置;响应于确定下一窗口位置为第一类型片段的最终窗口位置,确定缺失阈值是否大于第二类型片段阈值;以及响应于确定缺失阈值大于第二类型片段阈值,将起始窗口位置和缺失阈值处对应窗口添加入断点集合,以便形成关于缺失第二类型片段的断点集合和缺失第二类型片段。
计算每个窗口的重复表征数据和缺失表征数据包括:基于每个窗口所在片段的平均拷贝数、重复状态平均拷贝数、缺失状态平均拷贝数,计算重复比率与缺失比率;基于所计算的重复比率与缺失比率,计算每个窗口的观测数据正常概率密度、重复概率密度和缺失概率密度;以及基于所计算的每个窗口的观测数据正常概率密度、重复概率密度和缺失概率密度,计算每个窗口的重复表征数据和缺失表征数据。
提供发明内容部分是为了以简化的形式来介绍对概念的选择,它们在下文的具体实施方式中将被进一步描述。发明内容部分无意标识本发明的关键特征或主要特征,也无意限制本发明的范围。
附图说明
图1示出了根据本发明的实施例的用于实施检测拷贝数变异的方法的系统的示意图。
图2示出了根据本发明的实施例的用于检测拷贝数变异的方法的流程图。
图3示出了根据本发明的实施例的用于针对每一个窗口内的唯一比对数进行预处理的方法的流程图。
图4示出了根据本发明的实施例的用于确定所划分的第一类型片段和关于第一类型片段的断点集合的方法的流程图。
图5示出了根据本发明的用于计算每个窗口的重复表征数据和缺失表征数据的方法的流程图。
图6示出了根据本发明的实施例的用于识别重复小片段的断点的方法的流程图。
图7示出了根据本发明的实施例的用于识别缺失小片段的断点的方法的流程图。
图8示意性示出了适于用来实现本发明实施例的电子设备的框图。
在各个附图中,相同或对应的标号表示相同或对应的部分。
具体实施方式
下面将参照附图更详细地描述本发明的优选实施例。虽然附图中显示了本发明的优选实施例,然而应该理解,可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反,提供这些实施例是为了使本发明更加透彻和完整,并且能够将本发明的范围完整地传达给本领域的技术人员。
在本文中使用的术语“包括”及其变形表示开放性包括,即“包括但不限于”。除非特别申明,术语“或”表示“和/或”。术语“基于”表示“至少部分地基于”。术语“一个示例实施例”和“一个实施例”表示“至少一个示例实施例”。术语“另一实施例”表示“至少一个另外的实施例”。术语“第一”、“第二”等等可以指代不同的或相同的对象。
如前文所描述,传统的用于检测拷贝数变异的方案存在的不足之处在于:无法区分正常和低比例嵌合,单独针对低深度全基因组测序的准确性不高。
为了至少部分地解决上述问题以及其他潜在问题中的一个或者多个,本发明的示例实施例提出了一种用于检测拷贝数变异的方案。通过基于经由预处理的每个窗口的唯一比对数,分别针对每条染色体上的可能添加的断点,进行第一类型片段拟合,以便确定所划分的第一类型片段和关于第一类型片段的断点集合; 以及针对所划分的第一类型片段,计算每个窗口的重复表征数据和缺失表征数据,以便确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异;本发明可以通过基于关于第一类型片段(即,大片段)的断点集合和关于重复第二类型片段(即,重复小片段)和缺失第二类型片段(即,缺失小片段)的断点集合,确定拷贝数变异;使得对于大片段CNV和小片段CNV采用了不同的鉴别流程,进而既能保证大片段CNV整体不过于碎片化,且对于低比例的嵌合足够敏感,又能够保证针对小片段CNV有足够的灵敏度和准确性,因此,本发明针对低深度全基因组测序,显著提高检测拷贝数变异的准确性。图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描述根据本发明的实施例的用于检测拷贝数变异的方法。图2示出了根据本发明的实施例的用于检测拷贝数变异的方法200的流程图。应当理解,方法200例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解,方法200还可以包括未示出的附加动作和/或可以省略所示出的动作,本发明的范围在此方面不受限制。
在步骤202处,计算设备110基于待测样本的测序数据与参考基因组序列的比对结果数据,计算每一条读长在参考基因组中的位置。例如,计算设备110基于比对结果数据,针对单个样本,计算出每一条read在参考基因组中的位置。
关于待测样本的测序数据,其例如是待测样本的符合预定测序深度范围(例如,预定测序深度范围例如为低深度范围,如0.06X~0.2X)的全基因组测序数据。例如是单张芯片上所有低深度全基因组测序样本数据。
关于比对结果数据,其例如是采用BWA或者其他任意比对算法,将单张芯片上所有低深度全基因组测序样本数据分别与hg19参考基因组进行比对而生成的比对结果数据。
在步骤204处,计算设备110将基因组按照预定尺寸的窗口进行划分,以便统计每一个窗口内的唯一比对数。在一些实施例中,计算设备110还统计唯一比对的读长上出现碱基G或C的总数目,以用于GC矫正。
关于预定尺寸,其例如而不限于是20Kb。例如,计算设备110将将基因组按照20Kb的窗口(bin)大小进行划分。统计每一个窗口内唯一比对数,以及唯一比对的read上出现碱基G或C的总数目。
在步骤206处,计算设备110针对每一个窗口内的唯一比对数进行预处理。
关于针对每一个窗口内的唯一比对数进行预处理的方法,其例如包括:计算设备110针对每一个窗口内的唯一比对数进行归一化处理;基于归一化后的每一个窗口内的唯一比对数,计算关于待测样本的Y染色体唯一比对数的占比,以便确定关于待测样本的所属性别;分别构建常染色体的阴性参考集和与所确定的性别相对应的性染色体的阴性参考集;分别对每个样本上所有的窗口进行GC矫正;以及基于矫正后的每个窗口的唯一比对数和所构建的阴性参考集,确定拷贝数的观测值。下文将结合图3具体说明用于针对每一个窗口内的唯一比对数进行预处理的方法300,在此,不再赘述。
在步骤208处,计算设备110基于经由预处理的唯一比对数,分别针对每条染色体上的可能添加的断点,进行第一类型片段拟合,以便确定所划分的第一类型片段和关于第一类型片段的断点集合。
应当理解,每条染色体由两条染色单体组成,中间狭窄处称为着丝粒(centromere),它将染色体分为短臂(即,p臂)与长臂(即,q臂)。如果计算设备110染色体存在p臂与q臂,则分别对p臂与q臂进行第一类型片段拟合;以及将针对p臂的第一类型片段拟合结果和针对q臂的第一类型片段拟合结果进行拼接,以便确定所划分的第一类型片段和关于第一类型片段的断点集合。
例如,计算设备110通过迭代的方式,在大片段的每两个窗口之间设置断点,以便确定所设置的断点是否成立,如果所设置的断点成立,则所设置的断点为“真”,如果所设置的断点不成立,则所设置的断点为“假”。如果则所设置的断点为“真”,则将片段划分为两个片段,进而分别针对划分后的两个片段确定每两个窗口之间设置的断点是否成立,以此类推,不断迭代,直至满足第一类型片段拟合距离阈值,将所划分的片段确定为第一类型片段。
具体而言,关于确定所划分的第一类型片段和关于第一类型片段的断点集合的方法,其例如包括:将染色体上所有窗口作为一个整体片段,设置断点集合为空;在断点集合的基础上,遍历所有其他可能添加断点的位置,以便依次构建临时断点集合;分别基于断点集合和临时断点集合将染色体分为多个片段,以便计算每个片段的平均拷贝数;针对断点集合和临时断点集合,分别计算每个窗口的拷贝数、每个窗口的拷贝数到所在片段的平均拷贝数的距离,以便将所计算的距离的平均值作为误差;确定针对断点集合的初始误差和临时断点集合的最小误差,以便确定初始误差与最小误差之间的差值是否小于第一类型片段拟合距离阈值;响应于确定初始误差与最小误差之间的差值小于第一类型片段拟合距离阈值,以断点集合进行片段的划分;以及响应于确定初始误差与最小误差之间的差值大于或者等于第一类型片段拟合距离阈值,以与最小误差对应的临时断点集合替换断点集合,以便在替换后的断点集合上依次构建临时断点集合。下文将结合图4具体说明用于确定所划分的第一类型片段和关于第一类型片段的断点集合的方法400,在此,不再赘述。
在步骤210处,计算设备110针对所划分的第一类型片段,计算每个窗口的重复表征数据和缺失表征数据,以便确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异,第一类型片段的尺寸大于预定阈值。
重复表征数据例如是重复得分。缺失表征数据例如是缺失得分。例如,计算设备110基于所划分的大片段,针对每个窗口,通过高斯混合模型确定重复得分和缺失得分。
关于确定关于重复第二类型片段的方法,其例如包括:针对所划分的第一类型片段(大片段),将每个窗口的重复表征数据(例如,重复得分)进行累加,以生成连续窗口的重复表征数据累计值;确定连续窗口的重复表征数据累计值是否大于0;响应于确定连续窗口的重复表征数据累计值大于0,在连续窗口的重复表征数据累计值上继续累加下一窗口的重复表征数据,直至累加之后的连续窗口的重复表征数据累计值大于重复阈值;以及将重复表征数据累计值大于重复阈值时的连续窗口所在片段确定为重复第二类型片段。如果连续窗口的重复表征数据累计值或缺失表征数据累计值为非正值(小于或者等于0),则连续窗口所在片段不存在重复或者缺失,则将连续窗口的重复表征数据累计值或缺失表征数据累计值替换为0。
类似的,确定关于缺失第二类型片段的断点集合的方法,其例如包括:针对所划分的第一类型片段,将每个窗口的缺失数据表征数据(例如,缺失得分)分别进行累加,以生成连续窗口的缺失表征数据累计值;确定连续窗口的缺失表征数据累计值是否大于0;响应于确定连续窗口的缺失表征数据累计值大于0,在连续窗口的缺失表征数据累计值上继续累加下一窗口的缺失表征数据,直至累加之后的连续窗口的缺失表征数据累计值大于缺失阈值;以及将缺失表征数据累计值大于缺失阈值时的连续窗口所在片段确定为缺失第二类型片段。
然后,计算设备110基于所确定的重复第二类型片段,确定关于重复第二类型片段的断点集合;基于所确定的缺失第二类型片段,确定关于缺失第二类型片段的断点集合;针对关于重复第二类型片段的断点集合和关于缺失第二类型片段的断点集合进行去重,以便基于去重后的断点集合确定拷贝数变异。
在上述方案中,通过基于关于第一类型片段(即,大片段)的断点集合和关于重复第二类型片段(即,重复小片段)和缺失第二类型片段(即,缺失小片段)的断点集合,确定拷贝数变异;使得对于大片段CNV和小片段CNV采用了不同的鉴别流程,进而既能保证大片段CNV整体不过于碎片化,且对于低比例的嵌合足够敏感,又能够保证针对小片段CNV有足够的灵敏度和准确性,因此,本发明针对低深度全基因组测序,显著提高检测拷贝数变异的准确性。
以下将结合图3描述根据本发明的实施例的用于针对每一个窗口内的唯一比对数进行预处理的方法。图3示出了根据本发明的实施例的用于针对每一个窗口内的唯一比对数进行预处理的方法300的流程图。应当理解,方法300例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解,方法300还可以包括未示出的附加动作和/或可以省略所示出的动作,本发明的范围在此方面不受限制。
在步骤302处,计算设备110针对每一个窗口内的唯一比对数进行归一化处理。例如,计算设备110将长染色体上的比对上的读长数目归一化到例如5兆。进行归一化处理的原因在于,有时测序深度较高,有时测序深度较低,通过归一化处理可以将测序深度统一到相同测序水平上,以利于更准确地进行分析。
关于针对每一个窗口内的唯一比对数进行归一化处理的方法,其例如包括:计算设备110 基于单条染色体最大窗口数目、样本数量、染色体数目,确定归一化比例;基于归一化比例和每一个窗口内的唯一比对数,计算每一个窗口内的归一化后的唯一比对数;确定当前窗口内的归一化后的唯一比对数是否小于唯一比对数平均值的预定比例;以及响应于确定当前窗口内的归一化后的唯一比对数小于唯一比对数平均值的预定比例,确定当前窗口为检测盲区;以及将当前窗口内的归一化后的唯一比对数替换为空缺罚分值。
关于预定比例,其例如而不限于为20%。
关于空缺罚分值,其例如而不限于为NA_penalty。
应当理解,如果当前窗口为检测盲区,则在该区域比对上的读长数目会非常有限。因此,如果某个窗口,其归一化后的唯一比对数低于唯一比对数平均值的20%,则认为当前窗口为检测盲区,将唯一比对数由空缺罚分值替代。通过采用上述手段,即,对盲区进行空缺罚分,本发明可以灵活地控制CNV区间是否应该跨越盲区,对于CNV区间的边界更为精确。
以下结合公式(1)说明确定归一化比例的算法。
在上述公式(1)中,代表第k个待测样本归一化比例。i代表第i条染色体,其中i取值从1取到22。j代表单条染色体第j个窗口。k代表第k个待测样本。代表第k个待测样本、第i条染色体、第j个窗口的唯一比对数。12500为单条染色体最大窗口数目。22为常染色体条数。如果单张芯片的待测样本数目为k。唯一比对数矩阵X为24*12500*n的矩阵。其中24表示24条染色体(包括22条常染色体+X染色体+Y染色体)。
以下结合公式(2)说明计算每一个窗口内的归一化后的唯一比对数。
在步骤304处,计算设备110基于归一化后的每一个窗口内的唯一比对数,计算关于待测样本的Y染色体唯一比对数的占比,以便确定关于待测样本的所属性别。
由于男性与女性的差异主要在于Y染色体,因此,可以通过判定Y染色体上匹配上的reads数目是否较多来确定关于待测样本的所属性别。如果Y染色体上匹配上的reads数目较多,则确定待测样本属于男性;如果Y染色体上匹配上的reads数目较少,则确定当前样本所属的性别为女性。
具体而言,关于确定关于待测样本的性别的方法,其例如包括:计算设备110计算当前样本的Y染色体唯一比对数的占比;确定所计算的Y染色体唯一比对数的占比是否小于或者等于预定占比阈值;响应于确定所计算的Y染色体唯一比对数的占比小于或者等于预定占比阈值,确定当前样本所属的性别为女性;响应于确定所计算的Y染色体唯一比对数的占比大于预定占比阈值,确定当前样本所属的性别为男性。
以下结合公式(3)说明当前样本的Y染色体唯一比对数的占比的计算方式。
以下结合公式(4)说明确定关于待测样本的性别的计算方式。
在步骤306处,计算设备110分别构建常染色体的阴性参考集和与所确定的性别相对应的性染色体的阴性参考集。例如,计算设备110针对常染色体构建阴性参考集;以及针对性染色体,按照步骤304处所确定的关于待测样本的性别,构建对应的男性或者女性的阴性参考集。 应当理解,女性没有Y染色体,X染色体有两条。男性X染色体有一条,Y染色体有一条。因此,女性与男性的待测样本所测得比对上的reads的数目在染色体的分布是不一致的。因而,需要针对常染色体和性染色体按照不同方式构建阴性参考集。
关于确定阴性参考集的方法,其例如包括:分别计算每个窗口内的唯一比对数的均值和标准差;基于所计算的均值和标准差,针对每个窗口,计算保留区间;以及基于每个窗口的归一化后的唯一比对数与保留区间的比较,确定针对每个样本的、每窗口的保留系数,以便分别构建常染色体的阴性参考集和与所确定的性别相对应的阴性参考集
以下结合公式(5)和(6)说明计算每个窗口上唯一比对数的均值和标准差的算法。
以下结合公式(7)说明计算保留区间的算法。
以下结合公式(8)说明计算保留系数的算法。
以下结合公式(9)说明阴性参考集的确定算法。
针对常染色体,直接通过以上公式构建阴性参考集;针对性染色体,按照性别分别构建男性和女性的阴性参考集。
在步骤308处,计算设备110分别对每个样本上所有的窗口进行GC矫正。应当理解,在二代测序仪上测出的数据,通常都会表现出测序深度与GC 含量的相关性,称为GC bias。在染色体GC含量高或者染色体GC含量的位置,对比上的reads数目比较少。因此,通过GC矫正,可以去除由于染色体GC含量差异而造成的测序结果GC偏差。
以下结合公式(10)说明对每个样本上所有的窗口进行GC矫正的方法。
在上述公式(10)中,代表预测出每个窗口上的预测比对数。表示常染色体上所有存在比对的窗口的归一化后的唯一比对数的中位数。代表第k个待测样本、第i条染色体、第j个窗口的归一化后的唯一比对数。代表经由GC矫正后的每个窗口的唯一比对数。
例如,以代表k个待测样本、第i条染色体、第j个窗口的碱基G或C的总数目。以Y(即,GC含量)为自变量,X(即,唯一比对数,或称“比对上的reads数目”)为因变量,基于每个窗口的GC含量和比对上的reads数目,画出散点图。采用smooth.spline算法,拟合曲线,以便训练GC矫正模型。再通过经由训练的GC矫正模型,输入Y,预测出每个窗口上的预测比对数,即,预测出每个窗口上的经由GC矫正后的唯一比对数。
在步骤310处,计算设备110基于矫正后的每个窗口的唯一比对数和所构建的阴性参考集,确定拷贝数的观测值。
关于确定拷贝数的观测值的方法,其例如包括:计算设备110基于所预测出的每个窗口上的预测比对数和所构建的阴性参考集,计算fold change值;以及基于所确定的性别和所计算的fold change值,计算拷贝数的观测值。
以下结合公式(11)说明计算fold change值的方法。
以下结合公式(12)说明计算每窗口的观测拷贝数的方法。
以下将结合图4描述根据本发明的实施例的用于确定所划分的第一类型片段和关于第一类型片段的断点集合的方法。图4示出了根据本发明的实施例的用于确定所划分的第一类型片段和关于第一类型片段的断点集合的方法400的流程图。应当理解,方法400例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解,方法400还可以包括未示出的附加动作和/或可以省略所示出的动作,本发明的范围在此方面不受限制。
在步骤402处,计算设备110将染色体上所有窗口作为一个整体片段,设置断点集合为空。
在步骤404处,计算设备110在断点集合的基础上,遍历所有其他可能添加断点的位置,以便依次构建临时断点集合。
在步骤406处,计算设备110分别基于断点集合和临时断点集合将染色体分为多个片段,以便计算每个片段的平均拷贝数。例如,所计算的每个片段的平均拷贝数为SCN。
在步骤408处,计算设备110针对断点集合和临时断点集合,分别计算每个窗口的拷贝数、每个窗口的拷贝数到所在片段的平均拷贝数的距离,以便将所计算的距离的平均值作为误差。
以下结合公式(13)说明每个窗口的拷贝数到所在片段的平均拷贝数的距离的方法。
在上述公式(13)中,CN代表每个窗口的拷贝数。SCN代表所在片段的平均拷贝数。error代表误差。
在步骤410处,计算设备110确定针对断点集合的初始误差和临时断点集合的最小误差。
在步骤412处,计算设备110确定初始误差与最小误差之间的差值是否小于第一类型片段拟合距离阈值。例如,针对断点集合,令误差error 为初始误差 error_begin。针对临时断点集合,选取误差最小的情形,该最小的误差为error_end。与最小误差所对应的临时断点集合为最终临时断点集合。
在上述公式(14)中,CD_cutoff代表第一类型片段拟合距离阈值。
在步骤414处,如果计算设备110确定初始误差与最小误差之间的差值小于第一类型片段拟合距离阈值,以断点集合进行片段的划分。
在步骤416处,如果计算设备110确定初始误差与最小误差之间的差值大于或者等于第一类型片段拟合距离阈值,以与最小误差对应的临时断点集合替换断点集合,以便在替换后的断点集合上依次构建临时断点集合。
通过采用上述手段,本发明能够区分整个大片段是因为数据波动而存在关于拷贝数变异的噪音,还是真实地存在大片段的拷贝数变异,进而能够准确鉴别大片段的拷贝数变异。
以下将结合图5描述根据本发明的实施例的用于计算每个窗口的重复表征数据和缺失表征数据包括的方法。图5示出了根据本发明的实施例的用于计算每个窗口的重复表征数据和缺失表征数据包括的方法500的流程图。应当理解,方法500例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解,方法500还可以包括未示出的附加动作和/或可以省略所示出的动作,本发明的范围在此方面不受限制。
在步骤502处,计算设备110基于每个窗口所在片段的平均拷贝数、重复状态平均拷贝数、缺失状态平均拷贝数,计算重复比率与缺失比率。
关于重复比率,其例如指示重复状态下的拷贝数相对于标准状态下拷贝数的变化情况。
关于缺失比率,其例如指示缺失状态下的拷贝数相对于标准状态下拷贝数的变化情况。
以下结合公式(15)和(16)分别说明计算重复比率与缺失比率的算法。
在上述公式(15)和(16)中,代表重复比率。代表缺失比率。代表第k个待测样本、第i条染色体、第j个窗口所在片段的平均拷贝数。代表第k个待测样本、第i条染色体、第j个窗口的重复状态平均拷贝数。代表第k个待测样本、第i条染色体、第j个窗口的缺失状态的平均拷贝数。
在步骤504处,计算设备110基于所计算的重复比率与缺失比率,计算每个窗口的观测数据正常概率密度、重复概率密度和缺失概率密度。
应当理解,由于二代测序服从泊松分布,泊松分布与高斯分布存在近似关系,可近似认为每个窗口上的观测数据服从正常状态、重复状态和缺失状态下的高斯分布,如以下公式(17)至(19)所示。
在上述公式(17)至(19)中,代表第k个待测样本、第i条染色体、第j个窗口的观测数据。代表重复状态下的高斯分布。代表缺失状态下的高斯分布。代表正常状态下的高斯分布。其中,和满足以下公式(20)和(21)
在上述公式(20)和(21)中,代表阴性参考集的观测数据的重复比率,其指示如果第k个待测样本、第i条染色体、第j个窗口存在重复状态下的reads的数目。代表阴性参考集的观测数据的缺失比率,其指示如果第k个待测样本、第i条染色体、第j个窗口存在缺失状态下的reads的数目。代表重复比率。代表缺失比率。代表阴性参考集的第k个待测样本、第i条染色体、第j个窗口比对上的reads的数目。
以下结合公式(22)至(24)分别说明用于计算正常概率密度、重复概率密度和缺失概率密度的算法。
在步骤506处,计算设备110基于所计算的每个窗口的观测数据正常概率密度、重复概率密度和缺失概率密度,计算每个窗口的重复表征数据和缺失表征数据。
以下结合公式(25)至(26)分别说明用于计算每个窗口的重复表征数据和缺失表征数据的算法。
应当理解,对于盲区的窗口,该窗口的观测数据的重复表征数据score_dup和窗口的观测数据的缺失表征数据score_del均为NA_penalty。NA_penalty代表空缺(测序盲区bin)罚分值。通过空缺(测序盲区bin)罚分值可以灵活调整不同测序平台导致的差异。
应当理解,方法500用于计算每个窗口的重复表征数据和缺失表征数据(即,重复得分和缺失得分)。例如,每个窗口的重复得分和缺失得分的取值区间为“-1”至“1”。如果重复得分越接近“1”,则对应的窗口越可能存在重复。如果重复得分越接近“-1”,则对应的窗口越趋近于不存在重复。
通过采用上述手段,本发明能够准确计算每个窗口的重复表征数据和缺失表征数据。以下将结合图6描述根据本发明的实施例的识别重复小片段的和缺失小片段的断点的方法。图6示出了根据本发明的实施例的用于识别重复小片段的和缺失小片段的断点的方法600的流程图。应当理解,方法600例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解,方法600还可以包括未示出的附加动作和/或可以省略所示出的动作,本发明的范围在此方面不受限制。
应当理解,方法600用于将每个窗口的观测数据的重复得分和缺失得分通过累加的方式扩展至片段上连续窗口的重复表征数据累计值和缺失表征数据累计值。如果片段上的窗口的重复得分连续为正值(大于0),则继续将窗口的重复得分进行累加,直至重复表征数据累计值大于重复阈值,则确定连续窗口所在片段为重复片段。如果连续窗口的重复表征数据累计值为非正值(小于或者等于0),则连续窗口所在片段不存在重复,则将窗口的重复表征数据累计值替换为0。同理,如果连续窗口的缺失表征数据累计值为正值(大于0),则将窗口的缺失得分继续进行累加,直至缺失表征数据累计值大于缺失阈值,则确定连续窗口所在片段为缺失片段。方法600例如具体包括以下步骤。
在步骤602处,计算设备110针对所划分的第一类型片段,确定起始窗口位置、断点集合,以便计算起始窗口位置的重复表征数据累计值和重复阈值。
例如,计算设备110针对所划分的第一类型片段,确定起始窗口位置i为1,再令下一窗口位置j,k=i。断点集合初始值为空。通过以下公式(27)分别计算起始窗口位置的重复表征数据累计值score_dup。
在步骤604处,计算设备110确定起始窗口位置的重复表征数据累计值是否大于0。
在步骤606处,如果计算设备110确定起始窗口位置的重复表征数据累计值大于0,使得起始窗口位置叠加至下一窗口位置,以便获得自起始窗口位置至下一窗口位置的重复表征数据累计值。
获得自起始窗口位置至下一窗口位置的重复表征数据累计值的方法例如包括:在起始窗口位置的重复表征数据累计值上叠加下一窗口位置的重复表征数据,以更新重复表征数据累计值。例如如以下公式(28)所示。
在步骤608处,计算设备110确定自起始窗口位置至下一窗口位置的重复表征数据累计值是否大于0。
在步骤610处,如果计算设备110确定自起始窗口位置至下一窗口位置的重复表征数据累计值大于0,确定自起始窗口位置至下一窗口位置的重复表征数据累计值是否大于重复阈值。
在步骤612处,如果计算设备110确定自起始窗口位置至下一窗口位置的重复表征数据累计值大于重复阈值,将自起始窗口位置至下一窗口位置的重复表征数据累计值作为重复阈值。如果计算设备110确定自起始窗口位置至下一窗口位置的重复表征数据累计值小于或者等于重复阈值,跳转至步骤620,使得下一窗口位置迭代至再下一窗口位置,以便继续确认自起始窗口位置至再下一窗口位置的重复表征数据累计值是否大于0。
在步骤614处,计算设备110确定下一窗口位置是否为第一类型片段的最终窗口位置。
在步骤616处,如果计算设备110确定下一窗口位置为第一类型片段的最终窗口位置,确定重复阈值是否大于第二类型片段阈值。如果计算设备110确定下一窗口位置并非为第一类型片段的最终窗口位置,跳转至步骤620,使得下一窗口位置迭代至再下一窗口位置,以便继续确认自起始窗口位置至再下一窗口位置的重复表征数据累计值是否大于0。
在步骤618处,如果计算设备110确定重复阈值大于第二类型片段阈值,将起始窗口位置和重复阈值处对应窗口位置添加入断点集合,以便形成关于重复第二类型片段的断点集合和重复第二类型片段。如果计算设备110确定重复阈值小于或者等于第二类型片段阈值,确定起始窗口位置是否为第一类型片段的最终窗口位置。
通过采用上述手段,本发明可以准确地识别重复小片段的断点集合和重复小片段。
以下将结合图7描述根据本发明的实施例的识别重复小片段的和缺失小片段的断点的方法。图7示出了根据本发明的实施例的用于识别重复小片段的和缺失小片段的断点的方法700的流程图。应当理解,方法700例如可以在图8所描述的电子设备800处执行。也可以在图1所描述的计算设备110处执行。应当理解,方法700还可以包括未示出的附加动作和/或可以省略所示出的动作,本发明的范围在此方面不受限制。
在步骤702处,计算设备110针对经由第一类型片段拟合后的染色体,确定起始窗口位置、断点集合,以便计算起始窗口位置的缺失表征数据累计值和缺失阈值。
例如,计算设备110针对所划分的第一类型片段(即,经由大片段拟合后的大片段),确定起始窗口位置i为1,再令下一窗口位置j,k=i。初始断裂点集合为空。通过以下公式(29)分别计算缺失表征数据累计值score_del。
在步骤704处,计算设备110确定始窗口位置的缺失表征数据累计值是否大于0。
在步骤706处,如果计算设备110确定起始窗口位置的缺失表征数据累计值大于0,使得起始窗口位置叠加至下一窗口位置,以便获得自起始窗口位置至下一窗口位置的缺失表征数据累计值。
获得下一窗口位置的缺失表征数据累计值的方法例如包括:在起始窗口位置的缺失表征数据累计值上叠加下一窗口位置的缺失表征数据,以更新缺失表征数据累计值。例如以下公式(30)所示。
在步骤708处,计算设备110确定自起始窗口位置至下一窗口位置的缺失表征数据累计值是否大于0。
在步骤710处,如果计算设备110确定自起始窗口位置至下一窗口位置的缺失表征数据累计值大于0,确定自起始窗口位置至下一窗口位置的缺失表征数据累计值是否大于缺失阈值。
在步骤712处,如果计算设备110确定自起始窗口位置至下一窗口位置的缺失表征数据累计值大于缺失阈值,将自起始窗口位置至下一窗口位置的缺失表征数据累计值作为缺失阈值。如果计算设备110确定自起始窗口位置至下一窗口位置的缺失表征数据累计值小于或者等于缺失概率阈值,跳转至步骤720,使得下一窗口位置迭代至再下一窗口位置,以便继续确认自起始窗口位置至再下一窗口位置的缺失表征数据累计值是否大于0。
在步骤714处,计算设备110确定下一窗口位置是否为第一类型片段的最终窗口位置。
在步骤716处,如果计算设备110确定下一窗口位置为第一类型片段的最终窗口位置,确定缺失阈值是否大于第二类型片段阈值。如果计算设备110确定下一窗口位置并非为第一类型片段的最终窗口位置,跳转至步骤720,使得下一窗口位置迭代至再下一窗口位置,以便继续确认自起始窗口位置至再下一窗口位置的缺失表征数据累计值是否大于0。
在步骤718处,如果计算设备110确定缺失阈值大于第二类型片段阈值,将起始窗口位置和缺失阈值处对应窗口位置添加入断点集合,以便形成关于缺失第二类型片段的断点集合和缺失第二类型片段。如果计算设备110确定缺失概率阈值小于或者等于第二类型片段阈值,确定起始窗口位置是否为第一类型片段的最终窗口位置。
通过采用上述手段,本发明可以准确地识别缺失小片段的断点集合和缺失小片段。
图8示意性示出了适于用来实现本发明实施例的电子设备800的框图。电子设备800可以是用于实现执行图2至图7所示的方法200至700。如图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至700。例如,在一些实施例中,方法200至700可被实现为计算机软件程序,其被存储于机器可读介质,例如存储单元808。在一些实施例中,计算机程序的部分或者全部可以经由ROM 802和/或通信单元809而被载入和/或安装到电子设备800上。当计算机程序加载到RAM 803并由CPU 801执行时,可以执行上文描述的方法200至700的一个或多个操作。备选地,在其他实施例中,CPU 801可以通过其他任何适当的方式(例如,借助于固件)而被配置为执行方法200至700的一个或多个动作。
需要进一步说明的是,本发明可以是方法、装置、系统和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质,其上载有用于执行本发明的各个方面的计算机可读程序指令。
计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括:便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身,诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波(例如,通过光纤电缆的光脉冲)、或者通过电线传输的电信号。
这里所描述的计算机可读程序指令可以从计算机可读存储介质下载到各个计算/处理设备,或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收计算机可读程序指令,并转发该计算机可读程序指令,以供存储在各个计算/处理设备中的计算机可读存储介质中。
用于执行本发明操作的计算机程序指令可以是汇编指令、指令集架构(ISA)指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码,该编程语言包括面向对象的编程语言—诸如Smalltalk、C++等,以及常规的过程式编程语言—诸如“C”语言或类似的编程语言。计算机可读程序指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中,远程计算机可以通过任意种类的网络—包括局域网(LAN)或广域网(WAN)—连接到用户计算机,或者,可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。在一些实施例中,通过利用计算机可读程序指令的状态信息来个性化定制电子电路,例如可编程逻辑电路、现场可编程门阵列(FPGA)或可编程逻辑阵列(PLA),该电子电路可以执行计算机可读程序指令,从而实现本发明的各个方面。
这里参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或框图描述了本发明的各个方面。应当理解,流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合,都可以由计算机可读程序指令实现。
这些计算机可读程序指令可以提供给语音交互装置中的处理器、通用计算机、专用计算机或其它可编程数据处理装置的处理单元,从而生产出一种机器,使得这些指令在通过计算机或其它可编程数据处理装置的处理单元执行时,产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中,这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作,从而,存储有指令的计算机可读介质则包括一个制造品,其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。
也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上,使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤,以产生计算机实现的过程,从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。
附图中的流程图和框图显示了根据本发明的多个实施例的设备、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上,流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分,该模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中,方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如,两个连续的方框实际上可以基本并行地执行,它们有时也可以按相反的顺序执行,这依所涉及的功能而定。也要注意的是,框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合,可以用执行规定的功能或动作的专用的基于硬件的系统来实现,或者可以用专用硬件与计算机指令的组合来实现。
以上已经描述了本发明的各实施例,上述说明是示例性的,并非穷尽性的,并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下,对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择,旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进,或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。
以上仅为本发明的可选实施例,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等效替换、改进等,均应包含在本发明的保护范围之内。
Claims (12)
1.一种用于检测拷贝数变异的方法,其特征在于,包括:
基于待测样本的测序数据与参考基因组序列的比对结果数据,计算每一条读长在参考基因组中的位置;
将基因组按照预定尺寸的窗口进行划分,以便统计每一个窗口内的唯一比对数;
针对每一个窗口内的唯一比对数进行预处理;
基于经由预处理的唯一比对数,分别针对每条染色体上的可能添加的断点,进行第一类型片段拟合,以便确定所划分的第一类型片段和关于第一类型片段的断点集合; 以及
针对所划分的第一类型片段,计算每个窗口的重复表征数据和缺失表征数据,以便确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异,第一类型片段的尺寸大于预定阈值;
确定所划分的第一类型片段和关于第一类型片段的断点集合包括:
将染色体上所有窗口作为一个整体片段,设置断点集合为空;
在断点集合的基础上,遍历所有其他可能添加断点的位置,以便依次构建临时断点集合;
分别基于断点集合和临时断点集合将染色体分为多个片段,以便计算每个片段的平均拷贝数;
针对断点集合和临时断点集合,分别计算每个窗口的拷贝数、每个窗口的拷贝数到所在片段的平均拷贝数的距离,以便将所计算的距离的平均值作为误差;
确定针对断点集合的初始误差和临时断点集合的最小误差,以便确定初始误差与最小误差之间的差值是否小于第一类型片段拟合距离阈值;
响应于确定初始误差与最小误差之间的差值小于第一类型片段拟合距离阈值,以断点集合进行片段的划分;以及
响应于确定初始误差与最小误差之间的差值大于或者等于第一类型片段拟合距离阈值,以与最小误差对应的临时断点集合替换断点集合,以便在替换后的断点集合上依次构建临时断点集合;
计算每个窗口的重复表征数据和缺失表征数据包括:
基于每个窗口所在片段的平均拷贝数、重复状态平均拷贝数、缺失状态平均拷贝数,计算重复比率与缺失比率;
基于所计算的重复比率与缺失比率,计算每个窗口的观测数据正常概率密度、重复概率密度和缺失概率密度;以及
基于所计算的每个窗口的观测数据正常概率密度、重复概率密度和缺失概率密度,计算每个窗口的重复表征数据和缺失表征数据。
2.根据权利要求1所述的方法,其特征在于,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:
针对所划分的第一类型片段,将每个窗口的重复表征数据进行累加,以生成连续窗口的重复表征数据累计值;
确定连续窗口的重复表征数据累计值是否大于0;
响应于确定连续窗口的重复表征数据累计值大于0,在连续窗口的重复表征数据累计值上继续累加下一窗口的重复表征数据,直至累加之后的连续窗口的重复表征数据累计值大于重复阈值;以及
将重复表征数据累计值大于重复阈值时的连续窗口所在片段确定为重复第二类型片段。
3.根据权利要求2所述的方法,其特征在于,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:
针对所划分的第一类型片段,将每个窗口的缺失表征数据分别进行累加,以生成连续窗口的缺失表征数据累计值;
确定连续窗口的缺失表征数据累计值是否大于0;
响应于确定连续窗口的缺失表征数据累计值大于0,在连续窗口的缺失表征数据累计值上继续累加下一窗口的缺失表征数据,直至累加之后的连续窗口的缺失表征数据累计值大于缺失阈值;以及
将缺失表征数据累计值大于缺失阈值时的连续窗口所在片段确定为缺失第二类型片段。
4.根据权利要求3所述的方法,其特征在于,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:
基于所确定的重复第二类型片段,确定关于重复第二类型片段的断点集合;
基于所确定的缺失第二类型片段,确定关于缺失第二类型片段的断点集合;
针对关于重复第二类型片段的断点集合和关于缺失第二类型片段的断点集合进行去重,以便基于去重后的断点集合确定拷贝数变异。
5.根据权利要求1所述的方法,其特征在于,针对每一个窗口内的唯一比对数进行预处理包括:
针对每一个窗口内的唯一比对数进行归一化处理;
基于归一化后的每一个窗口内的唯一比对数,计算关于待测样本的Y染色体唯一比对数的占比,以便确定关于待测样本的所属性别;
分别构建常染色体的阴性参考集和与所确定的性别相对应的阴性参考集;
分别对每个样本上所有的窗口进行GC矫正;以及
基于矫正后的每个窗口的唯一比对数和所构建的阴性参考集,确定拷贝数的观测值。
6.根据权利要求5所述的方法,其特征在于,针对每一个窗口内的唯一比对数进行归一化处理包括:
基于单条染色体最大窗口数目、样本数量、染色体数目,确定归一化比例;
基于归一化比例和每一个窗口内的唯一比对数,计算每一个窗口内的归一化后的唯一比对数;
确定当前窗口内的归一化后的唯一比对数是否小于唯一比对数平均值的预定比例;以及
响应于确定当前窗口内的归一化后的唯一比对数小于唯一比对数平均值的预定比例,确定当前窗口为检测盲区;以及
将当前窗口内的归一化后的唯一比对数替换为空缺罚分值。
7.根据权利要求5所述的方法,其特征在于,确定关于待测样本的所属性别包括:
计算当前待测样本的Y染色体唯一比对数的占比;
确定所计算的Y染色体唯一比对数的占比是否小于或者等于预定占比阈值;
响应于确定所计算的Y染色体唯一比对数的占比小于或者等于预定占比阈值,确定当前样本所属的性别为女性;以及
响应于确定所计算的Y染色体唯一比对数的占比大于预定占比阈值,确定当前样本所属的性别为男性。
8.根据权利要求5所述的方法,其特征在于,分别构建常染色体的阴性参考集和与所确定的性别相对应的阴性参考集包括:
分别计算每个窗口内的唯一比对数的均值和标准差;
基于所计算的均值和标准差,针对每个窗口,计算保留区间;以及
基于每个窗口的归一化后的唯一比对数与保留区间的比较,确定针对每个样本的、每窗口的保留系数,以便分别构建常染色体的阴性参考集和与所确定的性别相对应的阴性参考集。
9.根据权利要求1所述的方法,其特征在于,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:
针对所划分的第一类型片段,确定起始窗口位置、断点集合,以便计算起始窗口位置的重复表征数据累计值和重复阈值;
响应于确定起始窗口位置的重复表征数据累计值大于0,使得起始窗口位置叠加至下一窗口位置,以便获得自起始窗口位置至下一窗口位置的重复表征数据累计值;
响应于确定自起始窗口位置至下一窗口位置的重复表征数据累计值大于0,确定自起始窗口位置至下一窗口位置的重复表征数据累计值是否大于重复阈值;
响应于确定自起始窗口位置至下一窗口位置的重复表征数据累计值大于重复阈值,将自起始窗口位置至下一窗口位置的重复表征数据累计值作为重复阈值;
确定下一窗口位置是否为第一类型片段的最终窗口位置;
响应于确定下一窗口位置为第一类型片段的最终窗口位置,确定重复阈值是否大于第二类型片段阈值;以及
响应于确定重复阈值大于第二类型片段阈值,将起始窗口位置和再下一窗口位置添加入断点集合,以便形成关于重复第二类型片段的断点集合和重复第二类型片段。
10.根据权利要求1所述的方法,其特征在于,确定关于重复第二类型片段和缺失第二类型片段的断点集合,以用于确定拷贝数变异包括:
针对所划分的第一类型片段,确定起始窗口位置、断点集合,以便计算起始窗口位置的缺失表征数据累计值和缺失阈值;
响应于确定起始窗口位置的缺失表征数据累计值大于0,使得起始窗口位置叠加至下一窗口位置,以便获得自起始窗口位置至下一窗口位置的缺失表征数据累计值;
响应于确定自起始窗口位置至下一窗口位置的缺失表征数据累计值大于0,确定自起始窗口位置至下一窗口位置的缺失表征数据累计值是否大于缺失阈值;
响应于确定自起始窗口位置至下一窗口位置的缺失表征数据累计值大于缺失阈值,将自起始窗口位置至下一窗口位置的缺失表征数据累计值作为缺失阈值;
确定下一窗口位置是否为第一类型片段的最终窗口位置;
响应于确定下一窗口位置为第一类型片段的最终窗口位置,确定缺失阈值是否大于第二类型片段阈值;以及
响应于确定缺失阈值大于第二类型片段阈值,将起始窗口位置和再下一窗口位置添加入断点集合,以便形成关于缺失第二类型片段的断点集合和缺失第二类型片段。
11.一种计算设备,其特征在于,包括:
至少一个处理单元;
至少一个存储器,所述至少一个存储器被耦合到所述至少一个处理单元并且存储用于由所述至少一个处理单元执行的指令,所述指令当由所述至少一个处理单元执行时,使得所述设备执行根据权利要求1至10任一项所述的方法的步骤。
12.一种计算机可读存储介质,其特征在于,计算机可读存储介质上存储有计算机程序,所述计算机程序被机器执行时实现根据权利要求1至10中任一项所述的方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210664065.0A CN114758720B (zh) | 2022-06-14 | 2022-06-14 | 用于检测拷贝数变异的方法、设备和介质 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210664065.0A CN114758720B (zh) | 2022-06-14 | 2022-06-14 | 用于检测拷贝数变异的方法、设备和介质 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114758720A CN114758720A (zh) | 2022-07-15 |
CN114758720B true CN114758720B (zh) | 2022-09-02 |
Family
ID=82336629
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210664065.0A Active CN114758720B (zh) | 2022-06-14 | 2022-06-14 | 用于检测拷贝数变异的方法、设备和介质 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114758720B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115579054B (zh) * | 2022-11-17 | 2023-06-02 | 北京大学 | 单细胞拷贝数变异探测方法、装置、设备及介质 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104221022A (zh) * | 2012-04-05 | 2014-12-17 | 深圳华大基因医学有限公司 | 一种拷贝数变异检测方法和系统 |
CN111916150A (zh) * | 2019-05-10 | 2020-11-10 | 北京贝瑞和康生物技术有限公司 | 一种基因组拷贝数变异的检测方法和装置 |
CN113270141A (zh) * | 2021-06-10 | 2021-08-17 | 哈尔滨因极科技有限公司 | 一种基因组拷贝数变异检测整合算法 |
Family Cites Families (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
SG11201404079SA (en) * | 2012-01-20 | 2014-10-30 | Bgi Diagnosis Co Ltd | Method and system for determining whether copy number variation exists in sample genome, and computer readable medium |
CN105574361B (zh) * | 2015-11-05 | 2018-11-02 | 上海序康医疗科技有限公司 | 一种检测基因组拷贝数变异的方法 |
CN105760712B (zh) * | 2016-03-01 | 2019-03-26 | 西安电子科技大学 | 一种基于新一代测序的拷贝数变异检测方法 |
CN108427864B (zh) * | 2018-02-14 | 2019-01-29 | 南京世和基因生物技术有限公司 | 一种拷贝数变异的检测方法、装置以及计算机可读介质 |
CN108573125B (zh) * | 2018-04-19 | 2022-05-13 | 上海亿康医学检验所有限公司 | 一种基因组拷贝数变异的检测方法及包含该方法的装置 |
-
2022
- 2022-06-14 CN CN202210664065.0A patent/CN114758720B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104221022A (zh) * | 2012-04-05 | 2014-12-17 | 深圳华大基因医学有限公司 | 一种拷贝数变异检测方法和系统 |
CN111916150A (zh) * | 2019-05-10 | 2020-11-10 | 北京贝瑞和康生物技术有限公司 | 一种基因组拷贝数变异的检测方法和装置 |
CN113270141A (zh) * | 2021-06-10 | 2021-08-17 | 哈尔滨因极科技有限公司 | 一种基因组拷贝数变异检测整合算法 |
Also Published As
Publication number | Publication date |
---|---|
CN114758720A (zh) | 2022-07-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US8412463B2 (en) | Methods of associating an unknown biological specimen with a family | |
CN114496077B (zh) | 用于检测单核苷酸变异和插入缺失的方法、设备和介质 | |
CN111210326A (zh) | 一种用于构建用户画像的方法及系统 | |
CN107480470B (zh) | 基于贝叶斯与泊松分布检验的已知变异检出方法和装置 | |
CN111462816B (zh) | 用于检测胚系基因微缺失微重复的方法、电子设备和计算机存储介质 | |
US20150286702A1 (en) | Adaptive variable selection for data clustering | |
CN111916150A (zh) | 一种基因组拷贝数变异的检测方法和装置 | |
CN114758720B (zh) | 用于检测拷贝数变异的方法、设备和介质 | |
CN111584002B (zh) | 用于检测肿瘤突变负荷的方法、计算设备和计算机存储介质 | |
KR20140006846A (ko) | Dna 서열의 데이터 분석 | |
CN114649055A (zh) | 用于检测单核苷酸变异和插入缺失的方法、设备和介质 | |
CN107451422A (zh) | 一种基因序列数据分析与在线交互可视化的方法 | |
CN111861328B (zh) | 建立物流识别库的方法、物流轨迹查询更新方法及设备 | |
US20220207302A1 (en) | Machine learning method and machine learning apparatus | |
CN110570908B (zh) | 测序序列多态识别方法及装置、存储介质、电子设备 | |
CN113535458B (zh) | 异常误报的处理方法及装置、存储介质、终端 | |
CN110797081B (zh) | 激活区域识别方法及装置、存储介质及电子设备 | |
CN113299342B (zh) | 基于芯片数据的拷贝数变异检测方法及检测装置 | |
WO2019213810A1 (zh) | 检测染色体非整倍性的方法、装置及系统 | |
CN113327646B (zh) | 测序序列的处理方法及装置、存储介质、电子设备 | |
CN114785616A (zh) | 数据风险检测方法、装置、计算机设备及存储介质 | |
KR102349023B1 (ko) | 뉴클레오티드 서열 변이의 빈도 분포 결정 | |
CN114792548B (zh) | 校正测序数据、检测拷贝数变异的方法、设备和介质 | |
CN107844478B (zh) | 一种专利文件的处理方法及装置 | |
US8359329B2 (en) | Method, computer apparatus and computer program for identifying unusual combinations of values in data |
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: 40077961 Country of ref document: HK |