CN112111565A - 一种细胞游离dna测序数据的突变分析方法和装置 - Google Patents

一种细胞游离dna测序数据的突变分析方法和装置 Download PDF

Info

Publication number
CN112111565A
CN112111565A CN201910538698.5A CN201910538698A CN112111565A CN 112111565 A CN112111565 A CN 112111565A CN 201910538698 A CN201910538698 A CN 201910538698A CN 112111565 A CN112111565 A CN 112111565A
Authority
CN
China
Prior art keywords
mutation
site
reads
cfdna
total
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
CN201910538698.5A
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.)
SHANGHAI GENMINIX INFORMATICS CO Ltd
Original Assignee
SHANGHAI GENMINIX INFORMATICS 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 SHANGHAI GENMINIX INFORMATICS CO Ltd filed Critical SHANGHAI GENMINIX INFORMATICS CO Ltd
Priority to CN201910538698.5A priority Critical patent/CN112111565A/zh
Publication of CN112111565A publication Critical patent/CN112111565A/zh
Pending legal-status Critical Current

Links

Classifications

    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
    • C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
    • C12Q1/6844—Nucleic acid amplification reactions
    • C12Q1/6858—Allele-specific amplification
    • 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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/50—Mutagenesis
    • 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/10—Sequence alignment; Homology search

Landscapes

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

Abstract

本发明涉及一种细胞游离DNA测序数据突变分析方法及装置。所述方法包括以下步骤：1.预处理原始的测序数据；2.比对到参考标准基因组，对PCR重复序列进行鉴定与排除，并校正序列的比对误差，得到准确的序列比对文件；3.进行肿瘤相关体细胞变异位点的鉴定，采用Sentieon软件TNScope工具，设置参数确保能检测到高于万分之五变异丰度的cfDNA变异位点，得到潜在的肿瘤相关体细胞变异位点表；4.分析和鉴定高度可信的低频变异，最终得到期望的阳性变异分析结果。本发明的方法能够达到超高的突变灵敏度检测以及高验证准确率。

Description

一种细胞游离DNA测序数据的突变分析方法和装置
技术领域
本申请属于生物信息学领域，具体地，涉及一种细胞游离DNA测序数据的突变分析方法和装置。
背景技术
随着液体活检技术在精准肿瘤学领域的快速发展，形成了一系列以细胞游离DNA(cell-free DNA，以下简称cfDNA)为核心分析物(analyte)的相关检测技术与分析方法来辅助确认肿瘤进展状态以及相关用药指导。由于具有无创伤性(直接抽提外周血)且能够按需多次采样等优势，使得cfDNA检测分析，成为肿瘤学研究与应用中不可代替的技术方案。但是，应用高通量并行短序列测序(即二代测序)技术，来鉴定细胞游离DNA中循环肿瘤(ctDNA)来源的相关低频变异仍然存在诸多难点。由于循环肿瘤DNA占整个细胞游离DNA的比例较低，所以在测序文库的制备过程中引入的较小变化(如DNA片段化、末端修复以及PCR扩增误差等)都会导致准确检测变异的难度增大；另外，对于测序产出数据的分析，由于cfDNA的长度较短、血液中自身存在的背景亚克隆非肿瘤突变以及人基因组区域的复杂度(特别是部分同源序列以及短重复序列区域)等因素，加之前面提到的测序技术因素都增加了通过生物信息分析流程来灵敏检测微量突变以及排除非肿瘤突变的难度，更多相关细节可参考文献[1,2]。随着分子标签(Unique molecule identifiers,UMI)技术在二代基因组测序上的应用，一定程度上增加了cfDNA中的低频(<1％)变异检测的可靠性以及准确性，但是目前仍然没有有效的方案来准确检测这部分低频突变。当前可以实施参考的方案可见于参考文献[2,3,4,5,6]，其中由于无法有效的鉴定测序过程中引入的背景噪音(如参考文献[2,3])或血浆cfDNA中本身存在的许多与肿瘤发展无关的背景噪音(如参考文献[4])，所以上述方案要么无法检测低频突变(<0.5％)，要么分析得到的突变位点中存在较多的假阳性结果；而参考文献[5,6]的方法由于无法较好的排除测序技术以及分析过程中引入的错误，故导致超低频突变结果不仅会遗漏部分突变，而且会引入部分假阳性突变。
发明内容
技术目的
本发明的一个目的是提供一种cfDNA测序数据的突变分析方法，其能够充分提高cfDNA突变检测的极限与准确率，并具有较强的临床可操作性。
本发明的另一个目的是提供一种cfDNA测序数据突变分析装置。
技术方案
根据一个方面，本发明提供一种cfDNA测序数据突变分析方法，其包括以下步骤：
1)利用分子标签(UMI)技术构建待测样本文库，然后用二代测序方法得到原始数据，对原始数据进行预处理，所述预处理包括：去除测序数据中包含的接头(adapter)序列、排除包括低复杂度序列以及长度较短的序列的部分低可信序列、将双端分子标签信息(Dual UMI)提取到Reads的名称中、比较双端Reads中的存在交叠的重复测序序列并矫正低质量的不一致碱基；
2)将步骤1得到的测序数据比对到人类标准参考基因组，而后借助之前提取的分子标签信息鉴定出同一类分子标签的PCR重复序列，对局部区域进行重新比对来校正插入缺失频繁区域容易出现的比对误差，最后针对测序数据进行单个碱基层面的质量校正，得到准确的序列比对文件；
3)上述步骤1和2在实际处理中需要重复两次，一次是待测cfDNA样本，另一次是同一个个体的血液对照样本，得到两个序列比对信息文件(bam)，利用这两个文件进行肿瘤相关体细胞变异(Somatic Variant)位点的鉴定，采用Sentieon软件TNScope工具，设置参数确保能检测到高于万分之五变异丰度的cfDNA变异位点以达到预期的漏检率，得到潜在的肿瘤相关体细胞变异(Somatic Variant)位点表；
4)分析和鉴定高度可信的低频变异，利用健康对照个体的背景信息以及分子标签(UMI)标记的序列信息来确定阴性变异与阳性变异，然后以此拟合机器学习模型的参数来判定全部的变异结果，最终得到期望的阳性变异分析结果。
具体地，上述步骤4)可包括以下三个步骤：
步骤4.1：利用健康的个体分别采用上述步骤1、2与3提到的方法进行变异位点分析，用TNER软件综合多个健康个体样本来构建正常样本的背景信息；
步骤4.2：对待测个体采用上述步骤1、2与3提到的方法进行肿瘤相关体细胞变异(Somatic Variant)位点突变分析，得到原始的突变位点结果；针对待测cfDNA样本经上述步骤2处理之后的序列比对信息文件，提取原始的突变位点表中每个突变位点相关的分子标签信息，记录Dual分子标签分别出现在双端Reads的次数，记为DMC1和DMC2，对于DMC1与DMC2都大于等于3的突变位点利用氧化DNA鉴定软件排除氧化导致的变异后，记上述得到位点集合为TP；另外，如果原始突变位点中存在与步骤4.1中的背景噪音位点列表相同的突变，则标记原始突变位点表中这部分突变位点为FP；将TP与FP中的交集位点以及不在TP或FP中的其他突变位点，记为UD，这部分位点是需要进行后续判断的位点列表，另外，分别记TP与FP中不在TP与FP交集中的位点为TMP和FMP；
步骤4.3：提取步骤4.2中的TMP、FMP、UD突变中如下所述的16个突变特征，该16个突变特征分为基因组信息特征与测序Reads信息特征，其中测序Reads信息特征又分类为二维特征与一维特征，
其中所述二维特征针对待测cfDNA样本与血液对照样本分别计算，其具体计算方式如下：提取每个突变的测序特征以及基因组信息特征，定义单个突变位点相关的测序片段集合为所有在基因组比对区域上与突变位置有重叠的Reads，称为覆盖总Reads，其总计数为TD；覆盖总Reads中该突变位点发生变异的Reads的集合记为突变总Reads，其总数为MD；然后计算得到突变等位基因的突变比例，记为MF；另外，分别计算覆盖总Reads与突变总Reads的平均比对质量，分别记为TMQ与MMQ；然后统计突变总Reads中每条read的碱基编辑距离，取平均值记为MED；统计突变总Reads中与评估突变一致的碱基的质量，取平均值记为MBQ，
所述一维特征仅针对cfDNA样本计算，其计算方法如下：利用cfDNA中突变Reads与对照胚系突变Reads的信息进行Fisher精确检验，并计算其p值，记为MFP；cfDNA样本中突变总Reads中每条read与评估突变一致的碱基距离read末端的长度，并求得其平均值与标准差，分别记为TML与TSD；cfDNA样本中突变总Reads的链偏差，采用TNscope计算得出的相关记录，记为SOR；另外记每个突变位点的插入缺失的序列长度为IDL，其中单核苷酸变异长度为0，
剩余的基因组信息特征计算方式如下：突变位点上下游50bp范围内重复碱基对的数量，记为NHP；突变位点上下游50bp范围内GC碱基对占总碱基对的比例，记为PGC；突变位点上下游50bp范围内碱基序列的复杂度得分，记为LCS；突变位点上下游10bp碱基对范围内鉴定为重复元素的比例，记为FRE，
每个位点的16个特征构建完成后，通过设置步骤4.2中TMP中位点为阳性位点，FMP中位点为阴性位点，利用随机森林模型为上述位点进行机器学习训练，得到拟合的模型，利用该拟合的模型判定UD中位点为阳性或阴性位点，从而得到待测个体的cfDNA变异位点突变分析结果。
本发明方法步骤1)中，优选地，采用fastp软件[7]进行处理，具体的参数如下--low_complexity_filter--length_required 60--correction--umi_loc＝per_read--umi_len＝3--umi_skip＝2。
本发明方法步骤2)中，优选地，采用参考基因组hd37d5，并另外设置BWA MEM的相关参数-M-Y-K 10000000，并添加分子标签(UMI)序列到RX字段中；另外在重复Reads标记过程中，通过PICARD中的子工具UmiAwareMarkDuplicatesWithMateCigar来实现带分子标签Reads的重复标记。
本发明方法步骤3)中，优选地，采用Sentieon软件的体细胞鉴定工具TNscope进行原始变异位点的鉴定，并设置相关参数为--sv_mask_ext 10--max_fisher_pv_active0.05--min_tumor_allele_frac 0.0005--no_mapq_cap 1--clip_by_minbq 1--max_error_per_read 3--min_init_tumor_lod 1.0--assemble_mode 4。
根据另一个方面，本发明提供了一种cfDNA测序数据的突变分析装置，其包括处理器，所述处理器被配置为执行上述方法。
技术效果
本发明充分利用对照样本、正常群体的cfDNA变异信息以及带分子标签测序技术进行超低频突变分析，依据样本的对照信息以及分子标签提取初步的阳性突变位点以及依据正常群体cfDNA变异信息提取阴性突变位点，并采用16个人工筛选的突变相关的位点特征来构造训练机器学习模型，由此，本发明的方法能够达到超高的突变灵敏度(检测低频突变达到0.1％)检测以及高验证准确率。
具体实施方式
下文中，通过具体实施方式来更好地说明本发明以使其被更好地理解，但这些实施方式不用于限制本发明的范围。
定义
在本发明中，“细胞游离DNA(circulating cell-free DNA,cfDNA)”是对游离在血液循环中的DNA片段的总称，包含了正常组织细胞以及肿瘤组织细胞凋亡或分泌产生的DNA片段。其中，肿瘤组织细胞凋亡或分泌产生的DNA片段特指循环肿瘤DNA(cell-free tumorDNA或ctDNA)，在血液中的半衰期短，可以实时反映肿瘤的动态变化。
术语“人类标准参考基因组”指通用的人类单倍体细胞核以及线粒体的DNA分子序列参照库，以FASTA格式存在，是生物信息学研究的基础，本文使用hs37d5人类参考基因组。
术语“同一个个体的血液对照样本”是指采用全血制备文库作为同一个检测个体的对照遗传DNA样本，使血细胞裂解从而使细胞核内的DNA被片段化后制备成文库进行后续测序处理，该步骤是常规的胚系遗传物质检测方法。而cfDNA样本制备特殊地从人类血清或血浆中纯化提取出游离的DNA片段，这些片段大多以短片段存在，浓度极低，长度在170bp左右，来源于组织细胞的凋亡或分泌。
具体地，本发明的经优化的检测cfDNA突变的方法可通过如下的步骤来实现。
步骤1.cfDNA样本用分子标签(UMI)技术构建文库后测序，对其原始双端数据fastq进行分析预处理：其中包括去除测序文库构建过程中引入的接头序列，过滤低复杂度Reads以及长度较短的Reads，并校正双端Reads中的存在交叠的重复测序序列的碱基内容以及将Dual分子标签信息提取到fastq文件中每条Reads以‘@’开头的第一行序列标识相关的描述信息中，此处可以采用fastp软件[7]进行处理，具体的参数如下--low_complexity_filter--length_required 60--correction--umi_loc＝per_read--umi_len＝3--umi_skip＝2。
步骤2.将步骤1得到的Reads比对到人类标准参考基因组，然后基于分子标签信息进行PCR重复序列的标记与去除，最后对剩余Reads进行局部区域的插入缺失重比对与碱基质量的校正，最后得到无冗余信息的序列比对结果文件(bam)。上述处理过程与BoardInstitute的GATK最佳实践中相关的处理过程基本一致，具体操作包含比对到参考基因组、重复Reads标记、局部INDEL重矫正和碱基质量矫正，具体可参考文献[8,9,10]，但是在比对以及重复序列标记上采用的具体参数上有所差异。比对过程的具体差异内容如下，采用参考基因组hd37d5，并另外设置BWA MEM的相关参数-M-Y-K10000000并添加分子标签(UMI)序列到RX字段中；另外在重复Reads标记过程中，通过PICARD中的子工具UmiAwareMarkDuplicatesWithMateCigar来实现带分子标签Reads的重复标记。其余方法与GATK最佳实践一致，包括：局部区域的插入缺失重比对，采用根据千人基因组数据库及Mills数据库的已知插入缺失信息来生成重比对的区域，然后对目标区域进行重比对。碱基质量的校正，是对Reads的碱基质量值进行重新校正，通过机器学习建模的方法，使输出结果的Reads碱基质量值能够更加接近真实情况与参考基因组之间错配的概率。比对步骤采用BWA软件处理，重复Reads标记采用PICARD软件，插入缺失重比对与碱基质量的校正采用GATK软件或者Sentieon软件的相关工具，最终得到带有分子标签信息的序列比对结果文件(bam)。
步骤3.上述步骤1和2在实际处理中需要重复两次，一次是待测个体的cfDNA样本，另一次是同一个体的胚系细胞对照样本(一般为血液样本)，得到两个序列比对信息文件(bam)。默认cfDNA样本需要用分子标签技术处理，按上述步骤1和2进行处理。胚系细胞对照样本如果采用带分子标签文库构建方式，则按照本文档提供的步骤1与步骤2进行处理，如果没有采用分子标签技术，那么可以采用参考文献[8,9,10]中GATK最佳实践的处理包含比对到人类标准参考基因组、重复Reads标记、局部INDEL重矫正和碱基质量矫正、生成序列比对信息文件(bam)。考虑到整个细胞游离DNA中循环肿瘤DNA(ctDNA)本身丰度可能非常低，采用Sentieon软件的体细胞鉴定工具TNscope[11,12]进行原始变异位点的鉴定，并设置相关参数为--sv_mask_ext10--max_fisher_pv_active 0.05--min_tumor_allele_frac0.0005--no_mapq_cap1--clip_by_minbq 1--max_error_per_read 3--min_init_tumor_lod 1.0--assemble_mode 4，上述设置能够保证高于万分之五的变异能够鉴定出来。初步变异结果相关的bam参数头文件对上述得到的初步变异结果进行保守的假阳性过滤，具体实施方法为过滤掉VCF中满足如下任意一个条件的位点：
PV>0.3(Fisher检验ref与alt在肿瘤及其配对样本中是否显著差异的P值大于0.3的位点)、PV2>0.3(高质量Reads做Fisher检验ref与alt在肿瘤及其配对样本中是否显著差异的P值大于0.3的位点)、ECNT>10(一个单倍型发现大于10个潜在的体细胞突变事件的位点)、FILTER字段包含germline_risk的位点(疑似遗传变异的位点)或FILTER字段包含triallelic_site的位点(同一位点两个以上变异方向的位点)。
步骤4.分析和鉴定高度可信的低频变异，具体的分析过程可包括下述步骤4.1、4.2以及4.3：
步骤4.1.利用健康的对照个体(特别是样本量较少时，一般少于6个大于2个)的突变分析结果来构建正常样本的背景信息，测序数据的分析采用本文档提到的步骤1、2与3来得到原始的突变结果；然后利用参考文献[13]的TNER方法(相关源代码可参考链接https://github.com/ctDNA/TNER)对所有的对照样本进行背景噪音分析，设置背景检测显著性水平阈值0.01，得到一系列背景噪音位点，这些位点将视为具有低频突变的噪音位点，需要在后续突变位点结果中予以排除。
步骤4.2.对待测个体的cfDNA样本及胚系细胞对照样本的测序结果采用步骤1、2与3提到的方法进行突变分析，得到原始的突变位点结果；针对cfDNA样本在步骤2的序列比对文件提取原始的突变位点表中每个突变位点相关的分子标签信息(UMI)，记录Dual分子标签(UMI)中分别出现在双端Reads的次数，记为(DMC1，DMC2)，对于DMC1与DMC2都大于等于3的突变位点利用参考文献[14]氧化DNA鉴定软件排除DNA损伤导致的变异后，记上述得到位点集合为TP；另外，如果原始突变位点位于步骤4.1中的背景噪音位点列表中并且与背景位点记录的突变一致(突变发生位置一致，突变变异内容一致)，则标记这部分突变位点为FP；将TP与FP中的交集位点以及不在TP或FP中的其他突变位点，记为UD(待定位点)，这部分位点即是需要进行后续判断的位点列表，另外分别记TP与FP中不在TP与FP交集中的位点为TMP与FMP。
步骤4.3.提取步骤4.2中TMP、FMP及UD的突变的如下16个突变特征，其中16个特征分为基因组信息特征与测序Reads信息特征，其中测序Reads信息特征又分类为二维特征(针对cfDNA样本与配对血液样本分别计算)与一维特征(仅针对cfDNA样本计算)。二维特征的具体计算方式如下：提取每个突变的测序特征以及基因组信息特征，定义单个突变位点相关的测序片段集合为所有在基因组比对区域上与突变位点有重叠的Reads，称为覆盖总Reads，其总计数为TD；覆盖总Reads中该突变位点发生变异的Reads的集合记为突变总Reads，其总数为MD；然后计算得到突变等位基因的突变比例，记为MF；另外，分别计算覆盖总Reads与突变总Reads的平均比对质量，分别记为TMQ与MMQ；然后统计突变总Reads中每条read的碱基编辑距离，取平均值记为MED；统计突变总Reads中与评估突变一致的碱基的质量，取平均值记为MBQ。所述一维特征仅针对cfDNA样本计算，其计算方法如下：利用cfDNA中突变Reads与对照白细胞突变Reads的信息进行Fisher精确检验，并计算其p值，记为MFP；cfDNA样本中突变总Reads中每条read与评估突变一致的碱基距离read末端的长度，并求得其平均值与标准差，分别记为TML与TSD；cfDNA样本中突变总Reads的链偏差，采用TNscope计算得出的相关记录，记为SOR；另外记每个突变位点的插入缺失的序列长度为IDL，其中单核苷酸变异长度为0。剩余的基因组信息特征计算方式如下：突变位点上下游50bp范围内重复碱基对的数量，记为NHP；突变位点上下游50bp范围内GC碱基对占总碱基对的比例，记为PGC；突变位点上下游50bp范围内碱基序列的复杂度得分(DUST)，记为LCS；突变位点上下游10bp碱基对范围内鉴定为重复元素的比例(RepeatMasker)，记为FRE。至此每个位点的16个特征构建完成，然后通过设置步骤4.2中TMP中位点为阳性位点，FMP中位点为阴性位点，利用随机森林模型(设定随机种子为1234)为上述位点进行机器学习训练，得到拟合的模型，这样可以利用此模型来判定UD中位点为阳性还是阴性位点，从而得到最终的cfDNA突变分析结果(包括TMP中位点以及上述模型判定为阳性的位点)。
关于一些边缘情况的说明，FP的位点通常会比较多，一般上百个，而TP位点较少，如果导致TMP位点少于10个的时候，可以考虑利用同批次的所有cfDNA样本的TMP位点以及FMP位点构造训练集进行训练，得到统一的拟合模型，然后对单个样本的cfDNA预测模型进行微调，结合人工审核随机森林结果模型中每个决策分支的生物学意义；如果同一批次上机的所有样本TMP位点小于20的情况，可以采用降低(DMC1，DMC2)阈值的策略来重新构建训练集，并进行模型训练；最后，如果总体TMP位点数量较少，可以考虑报告UD中的所有变异位点。
上述步骤所用变异位点相关特征信息示例如下：
CHR POS REF ALT TD MD MF MFP SOR TMQ MMQ MBQ IDL MED TML TSD NHP PGCLCS FRE Dual_UMI DMC1 DMC2
chr2 141283568 A G(1483,703) (72,0)(0.057,0.0) 0 0.956 (60.0,60.0)(60.0,0) (36.975,0)0 (1.99375,0) (56.78908742820676,64.92761394101876)(35.52531112192436,39.75588081822495)0 0.336633663 0 0 GCT-ATA 6 3
chr2 141299473G A(1898,894) (4,0) (0.002,0.0) 0.2783 1.153(60.003358119453104,60.0)(60.0,0) (40.2,0) 0 (1.2666666666666666,0)(59.75509714559846,70.2799740764744)(37.34407883559418,40.61024337862157)00.475247525 0 0 TGA-CTG 4 7
chr6 32821444G A(1689,834) (4,0) (0.003,0.0) 0.2015 1.044(59.999665831244776,60.0) (60.0,0) (37.75,0) 0 (1.0833333333333333,0)(53.11345029239766,74.51835443037974) (41.74388627466131,44.24873929207103)00.742574257 0 0 AGA-CTG 3 3
chr7 140420288 T A(912,864) (2,0) (0.002,0.0) 0.2821 2.811 (60.0,59.99552715654952) (60.0,60.0) (35.0,12.0) 0 (1.6363636363636365,16.0)(59.93936207500876,75.20575079872205)(39.453658075267825,42.422006979791625) 00.306930693 0 0 GAC-TTG5 5
chr9 98231172 C T (3314,868) (5,0) (0.002,0.0) 0.2996 0.996 (60.0,59.96615905245347) (60.0,0) (39.142857142857146,0)0 (1.0,0)(74.30174096973641,75.45910885504794) (45.16028856299894,41.19806399352321)00.653465347 0 0 CAG-TCT 5 4
chr12 25398284 C T (1781,726) (122,0) (0.076,0.0) 0 1.282(59.972667542706965,59.995077932731746) (60.0,0)(37.943820224719104,0)0(1.1647940074906367,0)(68.0712220762155,73.01558654634947)(38.177969370922305,40.35560454986908) 0 0.425742574 0 0 GAG-TAG4 6
chr12 111884735 C T (2878,617) (37,0) (0.013,0.0) 0.0017 0.749(59.99561073488839,60.0) (60.0,0) (32.76136363636363,0) 0(1.2613636363636365,0)(90.66666666666667,77.97672162948594)(42.817576657608704,40.48986463350886) 0 0.534653465 0 0 AAA-CCG 5 4
1.Chan,L.L.,Jiang,P.(2015).Bioinformatics analysis of circulatingcell-free DNA sequencing data.Clin Biochem,48(15),962-975.
2.Chen,S.,Liu,M.,Zhou,Y.(2018).Bioinformatics Analysis for Cell-FreeTumor DNA Sequencing Data.Methods Mol Biol,1754,67-95.
3.Mansukhani,S.et al.(2018).Ultra-Sensitive Mutation Detection andGenome-Wide DNA Copy Number Reconstruction by Error-Corrected CirculatingTumor DNA Sequencing.Clin Chem,64(11),1626-1635.
4.Xu,C.et al.(2019).smCounter2:an accurate low-frequency variantcaller for targeted sequencing data with unique molecularidentifiers..Bioinformatics.35(8),1299–1309
5.Newman,A.M.et al.(2014).An ultrasensitive method for quantitatingcirculating tumor DNA with broad patient coverage.Nat Med,20(5),548-554.
6.Newman,A.M.et al.(2016).Integrated digital error suppression forimproved detection of circulating tumor DNA.Nat Biotechnol,34(5),547-555.
7.Chen,S.Zhou,Y.Chen,Y.Gu,J.(2018).fastp:an ultra-fast all-in-oneFASTQ preprocessor.Bioinformatics,34(17),884-890.
8.McKenna,A.et al.(2010).The Genome Analysis Toolkit:a MapReduceframework for analyzing next-generation DNA sequencing data.Genome Res,20(9),1297-1303.
9.DePristo,M.A.et al.(2011).Aframework for variation discovery andgenotyping using next-generation DNAsequencing data.Nat Genet.43(5),491-498.
10.Van der Auwera,G.A.et al.(2013).From FastQ data to high confidencevariant calls:the Genome Analysis Toolkit best practices pipeline.Curr.Protoc.Bioinformatics,43(1),11.10.1-11.10.33.
11.Freed,Donald et al.(2018).The Sentieon Genomics Tools-Afast andaccurate solution to variant calling from next-generation sequence data.doi:10.1101/115717.
12.Freed,Donald et al.(2018).TNscope:Accurate Detection of SomaticMutations with Haplotype-based Variant Candidate Detection and MachineLearning Filtering.doi:10.1101/250647.
13.Deng,S.et al.(2018).TNER:ANovel Bayesian Background ErrorSuppression Method for Mutation Detection in Circulating Tumor DNA.BMCBioinformatics.19:387.
14.Costello,M.et al.(2013).Discovery and characterization ofartifactual mutations in deep coverage targeted capture sequencing data dueto oxidative DNAdamage during sample preparation.Nucleic Acids Research,41(6).e67.

Claims (6)

1.一种cfDNA测序数据突变分析方法，包括以下步骤：
1)利用分子标签技术构建待测样本文库，然后用二代测序方法得到原始数据，对原始数据进行预处理，所述预处理包括：去除测序数据中包含的接头序列、排除包括低复杂度序列以及长度较短的序列的部分低可信序列、将双端分子标签信息提取到Reads的名称中、比较双端Reads中的存在交叠的重复测序序列并矫正低质量的不一致碱基；
2)将步骤1)得到的测序数据比对到人类标准参考基因组，而后借助之前提取的分子标签信息鉴定出同一类分子标签的PCR重复序列，对局部区域进行重新比对来校正插入缺失频繁区域容易出现的比对误差，最后针对测序数据进行单个碱基层面的质量校正，得到准确的序列比对文件；
3)将上述步骤1)和2)重复两次，一次是针对待测cfDNA样本，另一次是针对同一个个体的血液对照样本，得到两个序列比对信息文件，利用这两个文件进行肿瘤相关体细胞变异位点的鉴定，采用Sentieon软件TNScope工具，设置参数确保能检测到高于万分之五变异丰度的cfDNA变异位点从而达到预期的漏检率，得到潜在的肿瘤相关体细胞变异位点表；
4)分析和鉴定高度可信的低频变异，利用健康对照个体的背景信息以及分子标签标记的序列信息来确定阴性变异与阳性变异，然后以此拟合机器学习模型的参数来判定全部的变异结果，最终得到期望的阳性变异分析结果。
2.根据权利要求1所述的分析方法，其中，步骤4)包括以下三个步骤：
步骤4.1：利用健康的个体分别采用步骤1)、2)与3)提到的方法进行变异位点分析，用TNER软件综合多个健康个体样本来构建正常样本的背景信息；
步骤4.2：对待测个体采用步骤1)、2)与3)提到的方法进行肿瘤相关体细胞变异位点突变分析，得到原始的突变位点结果；针对待测cfDNA样本经上述步骤2)处理之后的序列比对信息文件，提取原始的突变位点表中每个突变位点相关的分子标签信息，记录Dual分子标签分别出现在双端Reads的次数，记为DMC1和DMC2，对于DMC1与DMC2都大于等于3的突变位点利用氧化DNA鉴定软件排除氧化导致的变异后，记上述得到位点集合为TP；另外，如果原始突变位点中存在与步骤4.1中的背景噪音位点列表相同的突变，即，突变发生位置一致，突变Alt内容一致，则标记原始突变位点表中这部分突变位点为FP；将TP与FP中的交集位点以及不在TP或FP中的其他突变位点，记为UD，这部分位点是需要进行后续判断的位点列表，另外，分别记TP与FP中不在TP与FP交集中的位点为TMP和FMP；
步骤4.3：提取步骤4.2中的TMP、FMP、UD突变中如下所述的16个突变特征，该16个突变特征分为基因组信息特征与测序Reads信息特征，其中测序Reads信息特征又分类为二维特征与一维特征，
其中所述二维特征针对待测cfDNA样本与血液对照样本分别计算，其具体计算方式如下：提取每个突变的测序特征以及基因组信息特征，定义单个突变位点相关的测序片段集合为所有在基因组比对区域上与突变位置有重叠的Reads，称为覆盖总Reads，其总计数为TD；覆盖总Reads中该突变位点发生变异的Reads的集合记为突变总Reads，其总数为MD；然后计算得到突变等位基因的突变比例，记为MF；另外，分别计算覆盖总Reads与突变总Reads的平均比对质量，分别记为TMQ与MMQ；然后统计突变总Reads中每条read的碱基编辑距离，取平均值记为MED；统计突变总Reads中与评估突变一致的碱基的质量，取平均值记为MBQ，
所述一维特征仅针对cfDNA样本计算，其计算方法如下：利用cfDNA中突变Reads与对照胚系突变Reads的信息进行Fisher精确检验，并计算其p值，记为MFP；cfDNA样本中突变总Reads中每条read与评估突变一致的碱基距离read末端的长度，并求得其平均值与标准差，分别记为TML与TSD；cfDNA样本中突变总Reads的链偏差，采用TNscope计算得出的相关记录，记为SOR；另外记每个突变位点的插入缺失的序列长度为IDL，其中单核苷酸变异长度为0，
剩余的基因组信息特征计算方式如下：突变位点上下游50bp范围内重复碱基对的数量，记为NHP；突变位点上下游50bp范围内GC碱基对占总碱基对的比例，记为PGC；突变位点上下游50bp范围内碱基序列的复杂度得分，记为LCS；突变位点上下游10bp碱基对范围内鉴定为重复元素的比例，记为FRE，
每个位点的16个特征构建完成后，通过设置步骤4.2中TMP中位点为阳性位点，FMP中位点为阴性位点，利用随机森林模型为上述位点进行机器学习训练，得到拟合的模型，利用该拟合的模型判定UD中位点为阳性或阴性位点，从而得到待测个体的cfDNA变异位点突变分析结果。
3.根据权利要求1所述的分析方法，其中，在步骤1)中，采用fastp软件进行处理，具体的参数如下--low_complexity_filter--length_required 60--correction--umi_loc＝per_read--umi_len＝3--umi_skip＝2。
4.根据权利要求1所述的分析方法，其中，在步骤2)中，采用参考基因组hd37d5，并另外设置BWA MEM的相关参数-M-Y-K 10000000，并添加分子标签序列到RX字段中；在重复Reads标记过程中，通过PICARD中的子工具UmiAwareMarkDuplicatesWithMateCigar来实现带分子标签Reads的重复标记。
5.根据权利要求1所述的分析方法，其中，在步骤3)中，采用Sentieon软件的体细胞鉴定工具TNscope进行原始变异位点的鉴定，并设置相关参数为--sv_mask_ext 10--max_fisher_pv_active 0.05--min_tumor_allele_frac 0.0005--no_mapq_cap 1--clip_by_minbq 1--max_error_per_read 3--min_init_tumor_lod 1.0--assemble_mode 4。
6.一种cfDNA测序数据的突变分析装置，其包括处理器，其中，所述处理器被配置为执行权利要求1至5中任一项所述的方法。
CN201910538698.5A 2019-06-20 2019-06-20 一种细胞游离dna测序数据的突变分析方法和装置 Pending CN112111565A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910538698.5A CN112111565A (zh) 2019-06-20 2019-06-20 一种细胞游离dna测序数据的突变分析方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910538698.5A CN112111565A (zh) 2019-06-20 2019-06-20 一种细胞游离dna测序数据的突变分析方法和装置

Publications (1)

Publication Number Publication Date
CN112111565A true CN112111565A (zh) 2020-12-22

Family

ID=73796034

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910538698.5A Pending CN112111565A (zh) 2019-06-20 2019-06-20 一种细胞游离dna测序数据的突变分析方法和装置

Country Status (1)

Country Link
CN (1) CN112111565A (zh)

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113257360A (zh) * 2021-06-24 2021-08-13 北京橡鑫生物科技有限公司 癌症筛查模型、癌症筛查模型的构建方法及构建装置
CN113257350A (zh) * 2021-06-10 2021-08-13 臻和(北京)生物科技有限公司 基于液体活检的ctDNA突变程度分析方法和装置、ctDNA性能分析装置
CN114187964A (zh) * 2021-12-13 2022-03-15 深圳市海普洛斯生物科技有限公司 一种肺癌围手术期分子残留病灶基因检测panel及检测模型的构建方法
CN114596918A (zh) * 2022-03-11 2022-06-07 苏州吉因加生物医学工程有限公司 一种检测突变的方法及装置
CN115083521A (zh) * 2022-07-22 2022-09-20 角井(北京)生物技术有限公司 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及系统
CN115458052A (zh) * 2022-08-16 2022-12-09 珠海横琴铂华医学检验有限公司 基于一代测序的基因突变分析方法、设备和存储介质
CN115458051A (zh) * 2022-09-28 2022-12-09 北京泛生子基因科技有限公司 一种可保留分子标签信息的在测序数据中模拟小变异的方法、装置及计算机可读存储介质
CN115831233A (zh) * 2023-02-07 2023-03-21 杭州联川基因诊断技术有限公司 一种基于mTag的靶向测序数据预处理的方法、设备和介质
WO2023115662A1 (zh) * 2021-12-24 2023-06-29 广州燃石医学检验所有限公司 一种变体核酸的检测方法
CN116364178A (zh) * 2023-04-18 2023-06-30 哈尔滨星云生物信息技术开发有限公司 一种体细胞序列数据分类方法及相关设备
CN116504318A (zh) * 2023-06-25 2023-07-28 西安交通大学医学院第一附属医院 一种基于机器学习的肿瘤ctDNA信息统计处理方法
WO2023207396A1 (zh) * 2022-04-25 2023-11-02 天津华大基因科技有限公司 用于分析变异检测结果的模型的构建方法
CN117079720A (zh) * 2023-10-16 2023-11-17 北京诺禾致源科技股份有限公司 高通量测序数据的处理方法和装置

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113257350A (zh) * 2021-06-10 2021-08-13 臻和(北京)生物科技有限公司 基于液体活检的ctDNA突变程度分析方法和装置、ctDNA性能分析装置
CN113257350B (zh) * 2021-06-10 2021-10-08 臻和(北京)生物科技有限公司 基于液体活检的ctDNA突变程度分析方法和装置、ctDNA性能分析装置
CN113257360A (zh) * 2021-06-24 2021-08-13 北京橡鑫生物科技有限公司 癌症筛查模型、癌症筛查模型的构建方法及构建装置
CN114187964A (zh) * 2021-12-13 2022-03-15 深圳市海普洛斯生物科技有限公司 一种肺癌围手术期分子残留病灶基因检测panel及检测模型的构建方法
WO2023115662A1 (zh) * 2021-12-24 2023-06-29 广州燃石医学检验所有限公司 一种变体核酸的检测方法
CN114596918A (zh) * 2022-03-11 2022-06-07 苏州吉因加生物医学工程有限公司 一种检测突变的方法及装置
WO2023207396A1 (zh) * 2022-04-25 2023-11-02 天津华大基因科技有限公司 用于分析变异检测结果的模型的构建方法
CN115083521A (zh) * 2022-07-22 2022-09-20 角井(北京)生物技术有限公司 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及系统
CN115458052A (zh) * 2022-08-16 2022-12-09 珠海横琴铂华医学检验有限公司 基于一代测序的基因突变分析方法、设备和存储介质
CN115458051A (zh) * 2022-09-28 2022-12-09 北京泛生子基因科技有限公司 一种可保留分子标签信息的在测序数据中模拟小变异的方法、装置及计算机可读存储介质
CN115831233A (zh) * 2023-02-07 2023-03-21 杭州联川基因诊断技术有限公司 一种基于mTag的靶向测序数据预处理的方法、设备和介质
CN115831233B (zh) * 2023-02-07 2023-05-16 杭州联川基因诊断技术有限公司 一种基于mTag的靶向测序数据预处理的方法、设备和介质
CN116364178A (zh) * 2023-04-18 2023-06-30 哈尔滨星云生物信息技术开发有限公司 一种体细胞序列数据分类方法及相关设备
CN116364178B (zh) * 2023-04-18 2024-01-30 哈尔滨星云生物信息技术开发有限公司 一种体细胞序列数据分类方法及相关设备
CN116504318A (zh) * 2023-06-25 2023-07-28 西安交通大学医学院第一附属医院 一种基于机器学习的肿瘤ctDNA信息统计处理方法
CN116504318B (zh) * 2023-06-25 2023-08-25 西安交通大学医学院第一附属医院 一种基于机器学习的肿瘤ctDNA信息统计处理方法
CN117079720A (zh) * 2023-10-16 2023-11-17 北京诺禾致源科技股份有限公司 高通量测序数据的处理方法和装置
CN117079720B (zh) * 2023-10-16 2024-01-30 北京诺禾致源科技股份有限公司 高通量测序数据的处理方法和装置

Similar Documents

Publication Publication Date Title
CN112111565A (zh) 一种细胞游离dna测序数据的突变分析方法和装置
CN111341383B (zh) 一种检测拷贝数变异的方法、装置和存储介质
CN107423578B (zh) 检测体细胞突变的装置
CN109767810B (zh) 高通量测序数据分析方法及装置
CN110739027B (zh) 一种基于染色质区域覆盖深度的癌症组织定位方法及系统
CA3044231A1 (en) Validation methods and systems for sequence variant calls
US20190228131A1 (en) Novel method capable of differentiating fetal sex and fetal sex chromosome abnormality on various platforms
CN108319813A (zh) 循环肿瘤dna拷贝数变异的检测方法和装置
CN113724791B (zh) Cyp21a2基因ngs数据分析的方法、装置及应用
KR20190085667A (ko) 무세포 dna를 포함하는 샘플에서 순환 종양 dna를 검출하는 방법 및 그 용도
CN111326212A (zh) 一种结构变异的检测方法
EP3513345A1 (en) Methods and systems for ultra-sensitive detection of genomic alterations
CN114023381B (zh) 一种肺癌mrd融合基因判定方法、装置、存储介质及设备
WO2019213811A1 (zh) 检测染色体非整倍性的方法、装置及系统
Li et al. Comparative sequence alignment reveals River Buffalo genomic structural differences compared with cattle
CN116013419A (zh) 检测染色体拷贝数变异的方法
Zamperin et al. Sequencing of animal viruses: quality data assurance for NGS bioinformatics
Smith et al. Benchmarking splice variant prediction algorithms using massively parallel splicing assays
CN113373234A (zh) 一种基于突变特征的小细胞肺癌分子分型确定方法及应用
CN112837748A (zh) 一种区分不同解剖学起源肿瘤的系统及其方法
KR20210040714A (ko) 핵산 서열 분석에서 위양성 변이를 검출하는 방법 및 장치
CN116564406A (zh) 一种遗传变异自动化解读方法及设备
CN115132276A (zh) 一种实体瘤突变基因检测分析方法及系统
CN114093428B (zh) 一种ctDNA超高测序深度下低丰度突变的检测系统和方法
WO2018186687A1 (ko) 생물학적 시료의 핵산 품질을 결정하는 방법

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
WD01 Invention patent application deemed withdrawn after publication
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20201222