CN114242158B - ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备 - Google Patents

ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备 Download PDF

Info

Publication number
CN114242158B
CN114242158B CN202210154417.8A CN202210154417A CN114242158B CN 114242158 B CN114242158 B CN 114242158B CN 202210154417 A CN202210154417 A CN 202210154417A CN 114242158 B CN114242158 B CN 114242158B
Authority
CN
China
Prior art keywords
site
snv
reads
sites
base
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
CN202210154417.8A
Other languages
English (en)
Other versions
CN114242158A (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.)
Zhenhe Beijing Biotechnology Co ltd
Original Assignee
Zhenhe Beijing Biotechnology 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 Zhenhe Beijing Biotechnology Co ltd filed Critical Zhenhe Beijing Biotechnology Co ltd
Priority to CN202210154417.8A priority Critical patent/CN114242158B/zh
Publication of CN114242158A publication Critical patent/CN114242158A/zh
Application granted granted Critical
Publication of CN114242158B publication Critical patent/CN114242158B/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
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • 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/30—Detection of binding sites or motifs
    • 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
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/20—Polymerase chain reaction [PCR]; Primer or probe design; Probe optimisation
    • 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
    • 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

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

Abstract

本发明公开了一种ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备，属于生物医学检测技术领域。该检测方法包括接收SNV位点数据；数据预处理；提取SNV位点数据中每个位点特征，并对所提取的特征进行特征编码，将编码后的特征划分为第一特征集和第二特征集；采用Stacking策略构建SNV位点检测模型，所述Stacking策略的第一层包括两个LightGBM算法模型，第二层为逻辑回归算法学习器。所述存储装置、存储介质及设备均根据所提出的方法得以实现。本发明适用于ctDNA的单核苷酸变异检测，针对性强、特征种类多、敏感度高、结果稳定可靠。

Description

ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备
技术领域
本发明涉及生物医学检测技术领域，具体涉及一种ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备。
背景技术
循环肿瘤DNA(circulating tumor DNA，ctDNA)为肿瘤细胞破碎、凋亡或生长中的死亡后把肿瘤细胞DNA核酸物质释放到血液循环中的单链或双链DNA，其中的血浆游离DNA（cell-free DNA，cfDNA）含量非常丰富。它携带有与原发肿瘤组织相一致的分子遗传学特征，可用于无创检测肿瘤患者的基因组学状态，是一种高灵敏度、高特异性并被广泛作为多种肿瘤临床诊断的肿瘤标志物。来自肿瘤患者血液和血浆的ctDNA不仅可以用于对肿瘤进行诊断、预后以及治疗的动态检测，还可将其作为液体活检标志物进行临床应用。针对肿瘤患者血浆ctDNA的液体活检技术具有无创、便携以及高灵敏度等优点，对肿瘤早期筛查与诊断、靶向药物用药指导、预后和治疗动态监测具有重要意义。
新一代测序技术可以同时对数以万计的DNA分子进行测序，并提供突变所需的高灵敏度、易用性和准确的数据质量，常被作为分析癌症中体细胞突变的方法。单核苷酸变异(Single Nucleotide Variants，SNVs)是体细胞中最常见的单位点突变，常被作为确定蛋白质功能丧失和疾病风险的标志物，在人类各种类型的癌症的细胞增殖、肿瘤发生和肿瘤的精准治疗中起着非常重要的作用。
近年来，随着高通量测序（Next-generation sequencing, NGS）技术的迅速发展，越来越多的研究专注于SNV位点检测方法的开发。现有方法主要分为以下两类：同时使用肿瘤样本及其配对的正常样本来检测SNV位点和仅使用单个肿瘤样本来检测SNV位点。
目前体细胞变异检测方法大多数是基于肿瘤样本及其配对的正常样本进行设计开发，这类方法的核心思想是使用肿瘤样本来识别潜在的体细胞变异，并使用配对的正常样本将其中的胚系变异和杂合性丢失(Loss of heterozygosity, LoH) 变异从体细胞变异中筛选出来。这类方法的速度快，更适合高纯度的样本的SNV位点检测，但其SNV位点检测的性能有待提升。基于这种原理开发的方法以 MuTect2和VarDict软件为代表，常被用于对NGS测序进行SNV位点分析。Mutect2 先对单倍型进行局部重新组装来生成候选的变异单倍型，再结合贝叶斯体细胞似然模型来进一步从测序错误导致的假阳性变异中识别出真实的体细胞变异，然后将变异类型分为体细胞变异（配对的正常样本中没有）或胚系突变（配对的正常样本中有），最后通过硬过滤器进一步去除假阳性SNV位点。VarDict则是一个功能丰富的变异检测软件，具有多种独特的过滤方法去过滤掉假阳性位点，并支持DNA-seq与RNA-seq的数据分析。除此之外，Mutect2与VarDict虽然也可以用于分析单个肿瘤样本，但这些软件输出的变异结果中通常含有大量的假阳性位点。
由于单个的肿瘤血浆样本中会存在很多的低频位点，而这些位点常与测序误差、DNA聚合酶错误所导致的假位点相混淆，因此常使用唯一分子标记（Unique MolecularIdentifier, UMI）技术去特异性的标记样本文库中的每个分子，进而降低假阳性变异检出率，提高真实变异检测的灵敏度。这类方法的代表软件是smCounter，它通过联合考虑聚合酶链式反应（Polymerase chain reaction, PCR）和测序错误等因素，应用贝叶斯概率模型来计算每一个碱基位置发生变异的后验概率并输入到一系列的硬过滤器中进一步筛选，从而可以更加准确的检出真实的SNV位点。
由于在临床实践中常常很难获取到正常的配对样品，并且双样本检测的方法成本十分昂贵，因此基于单样本的肿瘤血浆样本的液体活检等方法被广泛的应用于肿瘤早期筛查与诊断、靶向药物用药指导、预后和治疗动态监测等领域。
尽管现有的SNV检测方法具有一定的检测率，但现有的方法是通过单个的肿瘤血浆样本去检测SNV位点，然后对检测出来的每一个SNV位点打上相应的硬过滤标签，最终还是需要人为的通过自身的经验结合不同的硬过滤标签来决定是否为一个真实的SNV变异信号，因此迫切的需要一种方法可以综合多种硬过滤指标、比对信息以及变异特有的一些特征来识别SNV位点，从而进一步提高检测的准确性。
发明内容
技术问题：本发明提供一种针对性强、特征种类多、敏感度高、结果稳定可靠的适用于ctDNA的单核苷酸变异的检测方法，以及相应的装置、存储介质和设备。
技术方案：第一方面，本发明提供一种ctDNA单核苷酸变异位点检测方法，包括：
接收SNV位点数据；
数据预处理，包括过滤测序数据、识别UMI信息并根据UMI信息得到一致性序列以及序列比对；
提取SNV位点数据中每个位点特征，并对所提取的特征进行特征编码，将编码后的特征划分为第一特征集和第二特征集；
采用Stacking策略构建SNV位点检测模型，所述Stacking策略的第一层包括两个LightGBM算法模型，第二层为逻辑回归算法学习器；
利用两个LightGBM算法模型分别对第一特征集和第二特征集进行训练，将第一层的两个不同的LightGBM模型对输入的SNV位点数据集给出的预测结果进行组合并形成新的SNV位点数据集，再作为训练集输入到第二层逻辑回归算法学习器进行训练，训练得到SNV位点检测模型，输出最终预测结果。
进一步地，所提取的特征包括：ref碱基、alt碱基、变异等位基因频率、位点总支持reads数、位点simplex支持reads数、位点duplex支持reads数、变异位点支持reads数、变异位点simplex支持reads数、变异位点duplex支持reads数、变异位点正向支持reads数、变异位点负向支持reads数、位点为A碱基的数目、位点为C碱基的数目、位点为G碱基的数目、位点为T碱基的数目、位点的碱基质量、位点的reads比对质量、位点的ref片段碱基序列、位点的比对片段碱基序列、位点的比对片段碱基质量、位点在参考基因组位置、比对片段在参考基因组的开始位置、位点UMI去重前的simplex支持reads数与错误率、位点UMI去重前的duplex支持reads数与错误率、与ref片段的编辑距离、片段的可能性、比对器输出的比对分数。
进一步地，所述的对特征进行特征编码，得到编码后的特征为：
（1）tcsV：变异位点支持reads数，包含simplex和duplex reads数；
（2）cVAF：变异等位基因频率；
（3）ref_Onehot_EIIP：ref碱基进行onehot编码成4维的特征向量，再加上该碱基对应的游离电子平均能量EIIP，一共构成5维的特征向量；其中，核苷酸A，C，G，T的EIIP值分别为0.1260，0.0806，0.1340，0.1335，其特征向量分别为 [1, 0, 0, 0, 0.1260]，[0, 1,0, 0, 0.1340]，[0, 0, 1, 0, 0.0806]，[0, 0, 0, 1, 0.1335]；
（4）alt_Onehot_EIIP：alt碱基进行onehot编码成4维的特征向量，再加上该碱基对应的游离电子平均能量EIIP，一共构成5维的特征向量；
（5）sscs/tcs：sscs与tcs的比值，其中sscs为位点simplex支持reads数，tcs为位点总支持reads数；
（6）dcsV/tcsV：dcsV与tcsV的比值，其中dcsV为变异位点duplex支持reads数；
（7）dcsV/dcs：dcsV与dcs的比值，其中dcs为位点duplex支持reads数；
（8）for_cAlt/tcsV：变异位点正向支持reads数占该变异位点支持reads数的比例，其中，for_cAlt为变异位点正向支持reads数；
（9）cA_ratio：A碱基的占比；
（10）cT_ratio：T碱基的占比；
（11）cG_ratio：G碱基的占比；
（12）Multiallelic：位点处的碱基类型数；
（13）Base_quality：位点的碱基质量值，最后取中位数、平均值、方差与标准差，一共构成4维的特征向量；
（14）Mapping_quality：位点的碱基质量值，最后取中位数、平均值、方差与标准差，一共构成4维的特征向量；
（15）Fragment_length：支持alt和ref片段长度之间的差值，取中位数、平均值、方差、标准差，一共构成4维的特征向量；
（16）Read_position：alt位点到reads首端的距离取中位数、平均值、方差、标准差，一共构成4维的特征向量；
（17）seqs：alt位点在reads中的上下游碱基，通过k-间隔核苷酸二聚体的组成编码来对reads的二聚体进行编码；
（18）family_size：通过bam文件提取该位点的reads在其同一UMI内的最大碱基深度、位点的reads在其同一UMI内的最小碱基深度、位点的reads在其同一UMI内的reads与该位点reads不一致的碱基比例，最后分别取中位数、平均值、方差与标准差，一共构成36维的特征向量；
（19）AS：比对器输出的比对分数；
（20）UM：与ref片段的编辑距离；
（21）UQ：片段的可能性。
进一步地，所述第一特征集包括编码后的特征（1）- 特征（11）；所述第二特征集包括编码后的特征（12）- 特征（21）。
进一步地，所述Stacking策略采用五折交叉验证。
进一步地，进行数据预处理时：
过滤测序数据包括：去除测序序列上的接头序列及低质量的碱基片段；
使用FGBIO软件识别序列上的UMI信息，并根据UMI信息对序列进行分组、合并得到一致性序列；
序列比对包括：将一致性序列通过软件对hg19人类参考基因组进行序列比对，并生成两个bam文件；然后通过GATK4的MergeBamAlignment模块合并上述的两个bam文件；然后对合并的bam文件进行文件排序式转换，并生成最终的bam文件。
进一步地，所述检测方法还包括：
对进行特征重要性排序，用于不同特征对SNV位点检测的贡献大小；以及，
SNV检测模型可视化解释，利用SHAP算法对特征进行排序打分，给出特征中的特征值的高低对算法预测的影响方向。
第二方面，本发明提供一种ctDNA单核苷酸变异位点检测装置，利用所提供的ctDNA单核苷酸变异位点检测方法检测循环肿瘤DNA单核苷酸变异位点，包括：
数据接收模块，其被配置为用于接收SNV位点数据；
数据预处理模块，其被配置为用于数据预处理，包括过滤测序数据、识别UMI信息并根据UMI信息得到一致性序列以及序列比对；
特征提取模块，其被配置为用于提取SNV位点数据中每个位点特征，并对所提取的特征进行特征编码，将编码后的特征划分为第一特征集和第二特征集；
模型构建模块，其被配置为采用Stacking策略构建SNV位点检测模型，所述Stacking策略的第一层包括两个LightGBM算法模型，第二层为逻辑回归算法学习器；
利用两个LightGBM算法模型分别对第一特征集和第二特征集进行训练，将第一层的两个不同的LightGBM模型对输入的SNV位点数据集给出的预测结果进行组合并形成新的SNV位点数据集，再作为训练集输入到第二层逻辑回归算法学习器进行训练，训练得到SNV位点检测模型，输出最终预测结果。
第三方面，本发明提供一种计算机可读存储介质，所述存储介质中存储有计算机指令，当所述计算机指令被处理器执行时，能够执行所述提供的ctDNA单核苷酸变异位点检测方法。
第四方面，本发明提供一种电子设备，包括：
所提出的计算机可读存储介质；
以及处理器，所述处理器能够执行存储于计算机可读存储介质中的计算机指令。
本发明与现有技术相比，具有以下优点：
（1）只需要使用UMI方式建库的单个肿瘤血浆ctDNA样本去检测低频的SNV位点，不需要正常的配对样本做参照。在临床实践中常常很难获取到正常的配对样品，并且双样本检测的方法成本十分昂贵。此外，肿瘤血浆ctDNA的含量极低，存在大量低频的SNV位点。本发明的使用UMI方式建库的单个肿瘤血浆ctDNA样本，从而可以纠正测序错误、DNA聚合酶错误所导致的低频假阳性位点，提高真实低频变异检测的灵敏度。
（2）建立了SNV位点独有的特征编码方式，包含了SNV位点的变异位点信息、序列信息以及比对信息。现有的SNV位点检测方法通常只是通过变异位点的信息去输出样本中可能的变异位点，然后再通过不同的硬过滤指标并结合自身的经验去识别真实的SNV变异位点。但这种方式需要花费大量的人力去审核位点，准确度有待提升，并且SNV位点大都不能只通过某个硬过滤的指标去区别，需要综合多个不同因素去综合的判断SNV位点的真假性。因此，本发明针对每一个SNV位点，依次提取了该位点的变异位点信息、序列信息以及比对信息进行编码，最后使用21种特征编码方式去描述该位点，综合了多个不同维度的因素去表征SNV位点。
（3）构建了基于Stacking集成学习的SNV位点检测模型，能更好地检测真实的低频SNV突变信号，具有更优的检测性能。现有的SNV位点检测方法是通过统计检验的方法或者单个的机器学习算法去构建SNV位点检测模型，虽然也具有一定的检出准确率，但泛化能力存在一定的局限性，没有利用不同机器学习算法的优势。本发明使用Stacking集成学习策略可以将三个好而不同的机器学习算法通过分层的形式进行整体模型的训练并构建了SNV位点检测的集成学习模型，然后使用三种不同的独立数据集进行验证，结果表明本发明的ctSNV检测模型相比于其他方法能更好地检测真实的低频SNV突变信号，具有更优的检测性能。
（4）使用SHAP算法对SNV位点检测集成模型进行了可视化解释，评估了每种特征维度对模型预测的影响方向。现有的SNV位点检测方法只会对SNV位点给定硬过滤标签或概率值，没有对检测方法进行可视化解释不同因素是如何影响SNV位点的阴阳性判断。本发明使用SHAP算法对21种特征的每种维度计算SHAP值，该值是基于该特征进行分裂时获得的信息增益得到的，最后对SNV位点检测集成模型进行了可视化展示每种特征维度对模型预测的影响方向，从而可以更加清晰地理解ctSNV检测模型的预测结果。
附图说明
图1为本发明的实施例中ctDNA单核苷酸变异位点检测方法的流程图；
图2为本发明的实施例中的SNV位点检测模型ctSNV的模型框架图；
图3为不同过滤方法在训练集S1上的性能比较图；
图4为不同过滤方法在独立测试集S2上的性能比较图
图5为特征的SHAP值的可视化图；
图6为本发明的实施例中ctDNA单核苷酸变异位点检测装置的框图；
图7为本发明的实施例中电子设备的框图。
具体实施方式
下面结合实施例和说明书附图对本发明作进一步的说明。说明的是，术语“第一”、“第二”等仅是出于描述的目的，并不能认为是对数量或性质等的定义。
第一方面，本发明提供一种ctDNA单核苷酸变异位点检测方法，该方法是一种基于单样本的肿瘤血浆样本的检测方法，图1示出了本发明的实施例中ctDNA单核苷酸变异位点检测方法的流程图，结合图1所示，实施例中，该方法包括：
步骤S100：接收SNV位点数据。
实施例中，为了便于说明整个方法的实现过程，结合申请人的临床数据进行说明。具体为申情人的微小残留病变（Minimal Residual Disease, MRD）流程临床报告系统提取2021年12月之前的所有样本分析结果，然后通过以下步骤来获取的阳性与阴性样本集以及阳性与阴性位点集，这样做的目的是为了能够训练SNV检测模型。
步骤S110：将MRD流程判定为阳性样本（质控合格）作为实施例的阳性样本集。
步骤S120：在阳性样本集中提取MRD流程判定为阳性SNV位点，同时这些位点通过变异支持reads数大于0且位点的p-value值 ≤ 0.05的条件进行过滤，最后得到了实施例的阳性SNV位点集。
步骤S130：提取MRD流程判定为阴性样本，同时这些样本需满足病理分期为Ⅰ、Ⅱ期，肿瘤占比 ≥ 10%以及样本质控合格等条件，最后得到了实施例的阴性样本集。
步骤S140：在阴性样本集中提取MRD流程判定为阴性SNV位点，同时这些位点需满足以下四个条件：①血浆样本中在该阴性位点的变异支持reads数≥2；②组织样本在该阴性位点的变异支持reads=0或≤3；③白细胞样本在该阴性位点的变异支持reads数=0；④这些位点必须在阳性位点集中含有。通过上述的筛选后，最后得到了实施例的阴性SNV位点集。
实施例中使用的阳性与阴性样本集以及阳性与阴性位点集的分布情况如表1所示。
表1 MRD流程基准数据集统计表
Figure 919798DEST_PATH_IMAGE001
步骤S200：数据预处理，包括过滤测序数据、识别UMI信息并根据UMI信息得到一致性序列以及序列比对。在实施例中，包括如下步骤：
步骤210：使用Trimmomatic v0.36软件去除测序序列上的接头序列及低质量的碱基片段。
步骤220：使用FGBIO v0.8软件来识别序列上的UMI信息，并根据UMI信息对序列进行分组、合并得到一致性序列simplex.consensus 和duplex consensus。
步骤230：将上述一致性序列通过BWA v0.7.12 软件对hg19（GRCh37）人类参考基因组进行序列比对，并生成simplexAligned.bam和duplexAligned.bam文件。然后通过GATK4的MergeBamAlignment模块合并上述的两个bam文件,并生成combinedConsensus.bam。再通过SAMTools v1.3软件对合并的bam文件进行文件排序式转换，并生成combinedConsensus.sorted.bam文件。
步骤S300：提取SNV位点数据中每个位点特征，并对所提取的特征进行特征编码，将编码后的特征划分为第一特征集和第二特征集。
使用pysam软件的pileup方法对combinedConsensus.sorted.bam文件提取阴阳性SNV位点集的每个位点特征，这些特征包括ref碱基、alt碱基、变异等位基因频率、位点总支持reads数（tcs）、位点simplex支持reads数（sscs）、位点duplex支持reads数（dcs）、变异位点支持reads数（tcsV）、变异位点simplex支持reads数（sscsV）、变异位点duplex支持reads数（dcsV）、变异位点正向支持reads数（for_cAlt）、变异位点负向支持reads数（rev_cAlt）、位点为A碱基的数目、位点为C碱基的数目、位点为G碱基的数目、位点为T碱基的数目、位点的碱基质量、位点的reads比对质量、位点的ref片段碱基序列、位点的比对片段碱基序列、位点的比对片段碱基质量、位点在参考基因组位置、比对片段在参考基因组的开始位置、位点UMI去重前的simplex支持reads数与错误率、位点UMI去重前的duplex支持reads数与错误率、与ref片段的编辑距离（NM）、片段的可能性（UQ）、比对器输出的比对分数（AS）。
然后将提取的特征通过以下方式进行特征编码，从而机器学习模型能够更好的针对SNV位点数据集进行训练与预测，编码后的特征为：
（1）tcsV：变异位点支持reads数，包含simplex和duplex reads数。
（2）cVAF：变异等位基因频率。
（3）ref_Onehot_EIIP：ref碱基进行onehot编码成4维的特征向量，再加上该碱基对应的游离电子平均能量EIIP，一共构成5维的特征向量；其中，核苷酸A，C，G，T的EIIP值分别为0.1260，0.0806，0.1340，0.1335，其特征向量分别为 [1, 0, 0, 0, 0.1260]，[0, 1,0, 0, 0.1340]，[0, 0, 1, 0, 0.0806]，[0, 0, 0, 1, 0.1335]。
（4）alt_Onehot_EIIP：alt碱基进行onehot编码成4维的特征向量，再加上该碱基对应的游离电子平均能量EIIP，一共构成5维的特征向量。
（5）sscs/tcs：sscs与tcs的比值，其中sscs为位点simplex支持reads数，tcs为位点总支持reads数。
（6）dcsV/tcsV：dcsV与tcsV的比值，其中dcsV为变异位点duplex支持reads数。
（7）dcsV/dcs：dcsV与dcs的比值，其中dcs为位点duplex支持reads数。
（8）for_cAlt/tcsV：变异位点正向支持reads数占该变异位点支持reads数的比例，其中，for_cAlt为变异位点正向支持reads数。
（9）cA_ratio：A碱基的占比。
（10）cT_ratio：T碱基的占比。
（11）cG_ratio：G碱基的占比。
（12）Multiallelic：位点处的碱基类型数。
（13）Base_quality：位点的碱基质量值，最后取中位数、平均值、方差与标准差，一共构成4维的特征向量。
（14）Mapping_quality：位点的碱基质量值，最后取中位数、平均值、方差与标准差，一共构成4维的特征向量。
（15）Fragment_length：支持alt和ref片段长度之间的差值，取中位数、平均值、方差、标准差，一共构成4维的特征向量。
（16）Read_position：alt位点到reads首端的距离取中位数、平均值、方差、标准差，一共构成4维的特征向量。
（17）seqs：alt位点在reads中的上下游碱基，通过k-间隔核苷酸二聚体的组成编码来对reads的二聚体进行编码。
k-间隔核苷酸二聚体的组成（Composition of k-spaced Nucleic Acid Pairs，CKSNAP）编码的定义如下：
CKSNAP特征编码计算任意k个核苷酸间隔（k = 0, 1, 2, ... , 5）的核苷酸二聚体的频率。例如当k = 0时，有16个0间隔的核苷酸二聚体（'AA', 'AC', 'AG', 'AT', 'CA', 'CC', 'CG', 'CT', 'GA', 'GC', 'GG', 'GT', 'TA', 'TC', 'TG', 'TT'）。该编码最后得到的特征向量可以定义为：
Figure 364686DEST_PATH_IMAGE002
其中，N表示对该核苷酸二聚体的计数，N all 表示核苷酸序列在不同gap下的核苷酸二聚体个数的总和。例如，如果核苷酸二聚体AA在核苷酸序列中出现N AA 次，则核苷酸二聚体AA的特征值等于N AA 除以N all 。
本发明的实施例中分别统计了gap=0,1,2,3,4,5 时的特征编码向量，此外，在实际分析中会将不确定的核苷酸标记为N，因此CKSNAP编码的最终的特征向量长度为5×5×6=150维。
（18）family_size：通过bam文件提取该位点的reads在其同一UMI内的最大碱基深度、位点的reads在其同一UMI内的最小碱基深度、位点的reads在其同一UMI内的reads与该位点reads不一致的碱基比例，最后分别取中位数、平均值、方差与标准差，一共构成36维的特征向量。
（19）AS：比对器输出的比对分数。
（20）UM：与ref片段的编辑距离。
（21）UQ：片段的可能性。
步骤S400：采用Stacking策略构建SNV位点检测模型，所述Stacking策略的第一层包括两个LightGBM算法模型，分别为LightGBM_1和LightGBM_2，第二层为逻辑回归算法学习器（Logistic Regression）。
在完成对提取的所有阳性SNV位点和阴性SNV位点进行特征编码后，SNV位点检测模型的选择变得尤为关键。针对SNV位点检测，机器学习建模方法通常比传统的统计学方法表现的更为优异。此外，由于SNV位点检测有确切的类别分类，即该点为阳性位点或阴性位点，故选择机器学习中的有监督学习算法LightGBM进行SNV位点检测的研发。梯度提升决策树(Gradient Boosting Decision Tree，GBDT) 算法是以决策树作为基函数，通过采用加法模型与不断最小化缩小训练的残差来综合多个决策树的结果进行数据分类的经典Boosting算法。LightGBM算法是对GBDT算法进行了进一步性能优化和效率提升，是一个快速高效、分布式、低内存、高性能的分类与回归算法。
除此之外，算法训练策略的选择也会十分影响算法模型的泛化能力。机器学习中常见的训练策略是集成学习，它是通过一种特定的集成学习策略将多种不同的机器学习算法（基学习器）有机的结合在一起去进行整体的训练，从而得到一个集成学习模型，相比于一种机器学习算法模型的泛化性能有不同程度的提升。根据基学习器集成策略的不同，集成学习可以分为Bagging、Boosting和Stacking三种集成学习策略。
本发明的实施例中采用Stacking集成学习策略进行SNV位点检测的研发，命名为ctSNV。由于Stacking策略通常是一种两层结构的集成学习框架，因此ctSNV在第一层使用了2个LightGBM算法模型分别针对第一特征集和第二特征集进行训练，目的是使得LightGBM算法可以从不同的特征子集的角度去检测SNV位点，从而增加SNV位点的检出率。在实施例中，第一特征集包括编码后的特征（1）-特征（11），第二特征集包括编码后的特征（12）-特征（21），利用LightGBM_1训练特征（1）-特征（11），利用LightGBM_2训练特征（12）-特征（12）。然后，将第一层的两个不同的LightGBM模型对输入的SNV位点数据集给出的预测结果进行组合并形成新的SNV位点数据集，再作为训练集输入到第二层逻辑回归（LogisticRegression）算法学习器进行训练，最终逻辑回归模型输出的预测结果即为Stacking集成学习SNV位点检测模型的最终预测结果。因此，本发明的ctSNV模型框架图2所示。
为了更准确清楚的描述ctSNV模型训练过程，使用
Figure 843072DEST_PATH_IMAGE003
表示SNV位点训练数据集，x i 表示第i个SNV位点的特征向量，y i 表示第i个SNV位点的阴阳性，y i = 1表示该SNV位点是阳性位点，y i = 0表示该SNV位点是阴性位点，n表示D数据集有n个SNV位点。ctSNV模型训练的详细步骤如下：
（1）Stacking集成学习策略采用的是五折交叉验证的方式将D划分成五个子集：
Figure 689805DEST_PATH_IMAGE004
。
（2）对于每个子集
Figure 493813DEST_PATH_IMAGE005
，分别使用剩余的4份子集去训练第一层的两个LightGBM基学习器
Figure 160418DEST_PATH_IMAGE006
，然后Ckt针对Sk子集中每个样本的特征向量x i 进行预测得到新的特征向量
Figure 708074DEST_PATH_IMAGE007
，其中
Figure 674893DEST_PATH_IMAGE008
表示基学习器Ckt对向量x i 的预测概率值，然后将所有的Ckt预测概率值组合形成二维的特征向量
Figure 649802DEST_PATH_IMAGE009
，最后将Sk转换成新的子集
Figure 69282DEST_PATH_IMAGE010
。该步骤重复迭代5次后，每个子集都被转换成新的子集，最终形成一个新的训练数据集
Figure 155050DEST_PATH_IMAGE011
，D'的特征向量是n×2的特征矩阵。
（3）将新的训练数据集D'作为第二层逻辑回归分类器的输入训练数据集，最终训练得到实施例中Stacking集成学习模型ctSNV。
利用训练好的模型ctSNV输出的最后结果即为最终预测结果。
为了验证本发明所提出的ctDNA单核苷酸变异位点检测方法相对于现有的方法具有更好的准确性，实施例中，利用不同测试数据集来评估SNV位点检测模型的性能。
实施例中使用五折交叉验证的方法来评估ctSNV模型的检测性能，为了更准确的评估不同方法检测SNV位点的性能，使用敏感度（Sensitivity, Sn）、特异度（Specificity,Sp）、精密度（Precision, PPV）、阴性预测率（NPV）、准确度（Accuracy, Acc）、F1 分数（F1-score）和马修斯相关系数（Matthews Correlation Coefficient，MCC）值等七个评价指标对不同SNV检测方法进行全面系统的评估。上述的七个指标具体定义为：
Figure 710796DEST_PATH_IMAGE012
Figure 591027DEST_PATH_IMAGE013
Figure 497804DEST_PATH_IMAGE014
Figure 387262DEST_PATH_IMAGE015
Figure 328673DEST_PATH_IMAGE016
Figure 645385DEST_PATH_IMAGE017
Figure 773878DEST_PATH_IMAGE018
其中，TP：True Positive，将正例预测为正例；FN：False Negative，将正例预测为负例；TN：True Negative，将负例预测为负例；FP：False Positive，将负例预测为正例。
测试示例1：
将步骤100中提取的SNV位点数据集按照月份进行划分，把2021年11月之前的临床数据集当作训练集并记为S1，将11月当月的临床数据集当作独立测试集并记为S2。为了更好地评估ctSNV与硬过滤方法之间的性能差异，这里选取了tcsV (variant 支持reads数≥3，则判断为SNV位点)、VAF (VAF≥0.003时，则判断为SNV位点) 、tcsV≥3或VAF≥0.003、GATK（MuTect2模块）、smCounter2、VarDict（过滤条件含variant 支持reads数≥3与VAF≥0.003）等六种过滤方法作为对比。七种不同的过滤方法在训练集S1上的性能比较结果如表2与图3所示，其中ctSNV的性能结果是其在训练集S1进行五折交叉验证的结果，从图表的结果可以很容易得出ctSNV模型在多个不同的指标中都表现最优，并且有最高的检出率，而其他过滤方法会漏检非常多阳性位点。
表2不同过滤方法在训练集S1上的性能比较表
Figure 467028DEST_PATH_IMAGE019
此外，七种不同的过滤方法在独立测试集S2上的性能比较结果如表3与图4所示。由于独立测试集中阳性位点的vaf都大于1‰，而绝大部分的阴性位点的vaf值小于1‰，因此该独立测试的结果可以看出单纯的硬过滤方法以及VarDict工具也是不错的，但是这些方法在含有大量的低频的阳性位点（如训练集）时则性能很差，而 ctSNV 在训练集和测试集上都性能稳定，检出率高。
表3不同过滤方法在独立测试集S2上的性能比较表
Figure 731787DEST_PATH_IMAGE020
测试示例2：
提取申请人内部Panel36_IVD流程八批经过实验验证的SNV位点集及其对应的样本，按照上述步骤S100-S300的流程提取对应的SNV位点的特征向量文件并将这八批数据作为独立数据集S3进行性能验证。同样的，四种不同的过滤方法在独立测试集S3上的性能比较结果如表4所示，该结果表明基于申请人MRD数据集训练的ctSNV检测模型可以较好地应用在申请人其他类似流程的SNV位点检测中。因此，当其他panel流程的测序深度、测序平台以及探针厂家与申请人的MRD流程基本一致时，基于申请人内部的MRD数据集训练的ctSNV检测模型可泛化到其他panel流程数据集进行SNV位点检测。
表4不同过滤方法在独立测试集S3上的性能比较表
Figure 484979DEST_PATH_IMAGE021
测试示例3：
采用申请人内部提取PSC1805标准品的原液用NC12878进行梯度稀释（0.05%，0.1%，0.3%，1%）后的样本和SNV位点，然后提取对应的SNV位点的特征向量文件并将这些数据作为独立数据集S4进行性能验证。同样的，四种不同的过滤方法在独立测试集S4上的性能比较结果如表5所示，该结果表明基于申请人MRD数据集训练的ctSNV检测模型在不同稀释浓度下性能都十分优异，尤其是它在超低频位点的灵敏度也是最好的。
表5不同过滤方法在独立测试集S4上的性能比较表
Figure 835189DEST_PATH_IMAGE022
测试示例4：
为了更加全面地评估ctSNV检测模型的特异性性能，采取了交叉病人特异性评估方式对ctSNV检测模型进行性能评估。同时，为了从不同视角去评估特异性，分别进行了样本水平特异性评估以及位点水平特异性评估。
在步骤S100中提取了34个阳性病人样本，每个病人样本追踪的变异位点数量不同。对于病人A的阳性样本，有B个病人样本和A病人的组织变异图谱没有交集，那就可以用A病人的组织变异图谱去B个病人样本中分别去追踪对应位点是否被判定成阳性变异位点，当有位点被判定阳性变异位点时，ctSNV检测模型将该病人样本预测为阳性样本。依次类推，重复34次直到所有阳性病人样本都进行了上述评估。
通过这种方式总共有33种病人样本组合参与样本水平特异性评估，其中有31种病人样本组合被ctSNV检测模型正确的判定为阴性样本，即样本水平特异性为93.94 %。与此同时，通过这种方式总共有37个变异位点组合参与位点水平特异性评估，其中有35个变异位点组合被ctSNV检测模型正确的判定为阴性变异位点，即位点水平特异性为94.59%。
进一步地，在本发明的实施例中，所提供的ctDNA单核苷酸变异位点检测方法还包括特征重要性分析以及对SNV位点检测模型可视化解释。
进行重要性分析的目的是探索ctSNV检测模型中不同特征对SNV位点检测的贡献大小，本发明针对ctSNV集成模型的第一层的两个基学习器分别通过调用LightGBM算法的“feature_importances_”去对各自学习器所使用的特征进行重要性排序。LightGBM_1与LightGBM_2学习器重要性排序结果如表6和表7所示。
表6 LightGBM_1学习器的特征重要性排序
Figure 66450DEST_PATH_IMAGE023
表7 LightGBM_2学习器的特征重要性排序（Top 7）
Figure 716875DEST_PATH_IMAGE024
为了进一步增强ctSNV检测模型的可解释性，本发明使用SHAP（Shapley additiveex-planations）算法对特征进行排序打分，并给出特征中的特征值的高低对算法预测的影响方向。该算法引入了特征重要性的统一度量SHAP值，该值定义如下：
Figure 109810DEST_PATH_IMAGE025
式中F表示特征集，S表示从F中移除第i个特征后得到的特征子集，
Figure 744053DEST_PATH_IMAGE026
表示该算法在包含第i个特征的特征集上训练得到的模型，f s 表示该算法模型在不包含第i个特征的特征集S上训练得到的模型，x s 表示特征子集S中的输入特征值。
因此，本发明使用SHAP算法对21种特征的每个特征维度进行了重要性评估，并将其对ctSNV检测模型的影响进行了可视化展示，如图5所示。
图5展示了最优的22个特征维度下对应的SHAP值的可视化图，图中特征的重要性从下向上依次变大，黑颜色的点代表该特征的值较大，灰颜色的点代表该特征的值较小。此外，正的SHAP值将促进ctSNV检测模型将该位点预测为阳性SNV位点，而负的SHAP值则会促进ctSNV检测模型将该位点预测为阴性SNV位点。
从图5中可以看到，排序第三重要的cG_ratio特征图表示变异位点处G碱基的占比越高，则模型越倾向于将该位点归为阴性SNV位点。反之，变异位点处G碱基的占比越低时，则模型根据单个的cG_ratio特征不能将阴阳性SNV位点区分开。同理，排序第十二重要的
dcsV/tcsV特征图表示duplex支持变异reads数占比越大，则模型越倾向于将其归为阳性SNV位点。反之，duplex支持变异reads数占比越小，则模型越倾向于将其归为阴性SNV位点。
本发明使用单个肿瘤血浆ctDNA样本，不需要正常的配对样本做参照，并且使用UMI方式建库的血浆ctDNA样本，可以纠正测序错误、DNA聚合酶错误所导致的低频假阳性位点，提高真实低频变异检测的灵敏度，降低患者的检测成本。然后对样本的比对文件进行pileup并建立SNV位点独有的21种特征编码方式，包含了SNV位点的变异位点信息、序列信息以及比对信息，综合了多个不同维度的因素去表征SNV位点，囊括了大部分SNV位点硬过滤指标相关因素。然后采用了Stacking集成学习策略将三个好而不同的机器学习算法通过分层的形式进行整体模型的训练并构建了SNV位点检测的集成学习模型，能更好地检测真实的低频SNV突变信号，具有更优的检测性能。最后使用SHAP算法对21种特征的每种维度计算SHAP值，对ctSNV检测模型进行了可视化展示每种特征维度对模型预测的影响方向，从而可以更清晰地理解ctSNV检测模型的预测结果。
第二方面，本发明的实施例中提供一种ctDNA单核苷酸变异位点检测装置，该装置可以利用本发明提供的ctDNA单核苷酸变异位点检测方法检测循环肿瘤DNA单核苷酸变异位点，如图6所示，实施例中，该装置包括：数据接收模块、数据预处理模块、特征提取模块和模型构建模块。其中，数据接收模块被配置为用于接收SNV位点数据。数据预处理模块被配置为对数据进行预处理，包括过滤测序数据、识别UMI信息并根据UMI信息得到一致性序列以及序列比对。特征提取模块被配置为提取SNV位点数据中每个位点特征，并对所提取的特征进行特征编码，将编码后的特征划分为第一特征集和第二特征集。模型构建模块被配置为采用Stacking策略构建SNV位点检测模型，所述Stacking策略的第一层包括两个LightGBM算法模型，第二层为逻辑回归算法学习器；利用两个LightGBM算法模型分别对第一特征集和第二特征集进行训练，将第一层的两个不同的LightGBM模型对输入的SNV位点数据集给出的预测结果进行组合并形成新的SNV位点数据集，再作为训练集输入到第二层逻辑回归算法学习器进行训练，训练得到SNV位点检测模型，得到预测结果。
更具体的，各个模块具体如何实现相应的功能，与上文中关于ctDNA单核苷酸变异位点检测方法的说明相对应，此处就不再赘述。
第三方面，本发明提供一种计算机可读存储介质，其中存储有计算机指令，当计算机指令由处理器执行时，能够实现本发明的实施例中所提出的ctDNA单核苷酸变异位点检测方法。本发明中所称的计算机可读介质包括各种类型的计算机存储介质，可以是通用或专用计算机能够存取的任何可用介质。举例而言，计算机可读介质可以包括RAM、ROM、EPROM、E2PROM、寄存器、硬盘、可移动盘、CD-ROM或其他光盘存储器、磁盘存储器或其他磁存储设备、或者能够用于携带或存储具有指令或数据结构形式的期望的程序代码单元并能够由通用或特定用途计算机、或者通用或特定用途处理器进行存取的任何其他临时性或者非临时性介质。如本文所使用的，盘（disk）和碟（disc）包括CD、激光碟、光碟、数字多用途光碟（DVD）、软盘和蓝光碟，其中盘通常磁性地复制数据，而碟则用激光来光学地复制数据。上述的组合也应当包括在计算机可读介质的保护范围之内。示例性存储介质耦合到处理器以使得该处理器能从/向该存储介质读写信息。在替换方案中，存储介质可以被整合到处理器。处理器和存储介质可驻留在ASIC中。ASIC可驻留在用户终端中。在替换方案中，处理器和存储介质可作为分立组件驻留在用户终端中。
本发明的第四方面，提供一种电子设备。如图7所示，电子设备包括如上文所述的任意一种计算机可读存储介质以及处理器。其中，处理器被配置成能够执行计算机可读存储介质中存储的计算机指令。需要说明的是，电子设备还可以包括其他部件，例如输入设备、显示设备等，出于清楚说明本发明的原理角度考虑，这些部件并未示出。
上述实施例仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和等同替换，这些对本发明权利要求进行改进和等同替换后的技术方案，均落入本发明的保护范围。

Claims (8)

1.一种ctDNA单核苷酸变异位点检测方法，其特征在于，包括：
接收SNV位点数据；
数据预处理，包括过滤测序数据、识别UMI信息并根据UMI信息得到一致性序列以及序列比对；
提取SNV位点数据中每个位点特征，并对所提取的特征进行特征编码，将编码后的特征划分为第一特征集和第二特征集；所提取的特征包括：ref碱基、alt碱基、变异等位基因频率、位点总支持reads数、位点simplex支持reads数、位点duplex支持reads数、变异位点支持reads数、变异位点simplex支持reads数、变异位点duplex支持reads数、变异位点正向支持reads数、变异位点负向支持reads数、位点为A碱基的数目、位点为C碱基的数目、位点为G碱基的数目、位点为T碱基的数目、位点的碱基质量、位点的reads比对质量、位点的ref片段碱基序列、位点的比对片段碱基序列、位点的比对片段碱基质量、位点在参考基因组位置、比对片段在参考基因组的开始位置、位点UMI去重前的simplex支持reads数与错误率、位点UMI去重前的duplex支持reads数与错误率、与ref片段的编辑距离、片段的可能性、比对器输出的比对分数；
所述的对特征进行特征编码，得到编码后的特征为：
（1）tcsV：变异位点支持reads数，包含simplex和duplex reads数；
（2）cVAF：变异等位基因频率；
（3）ref_Onehot_EIIP：ref碱基进行onehot编码成4维的特征向量，再加上该碱基对应的游离电子平均能量EIIP，一共构成5维的特征向量；其中，核苷酸A，C，G，T的EIIP值分别为0.1260，0.0806，0.1340，0.1335，其特征向量分别为 [1, 0, 0, 0, 0.1260]，[0, 1, 0,0, 0.1340]，[0, 0, 1, 0, 0.0806]，[0, 0, 0, 1, 0.1335]；
（4）alt_Onehot_EIIP：alt碱基进行onehot编码成4维的特征向量，再加上该碱基对应的游离电子平均能量EIIP，一共构成5维的特征向量；
（5）sscs/tcs：sscs与tcs的比值，其中sscs为位点simplex支持reads数，tcs为位点总支持reads数；
（6）dcsV/tcsV：dcsV与tcsV的比值，其中dcsV为变异位点duplex支持reads数；
（7）dcsV/dcs：dcsV与dcs的比值，其中dcs为位点duplex支持reads数；
（8）for_cAlt/tcsV：变异位点正向支持reads数占该变异位点支持reads数的比例，其中，for_cAlt为变异位点正向支持reads数；
（9）cA_ratio：A碱基的占比；
（10）cT_ratio：T碱基的占比；
（11）cG_ratio：G碱基的占比；
（12）Multiallelic：位点处的碱基类型数；
（13）Base_quality：位点的碱基质量值，最后取中位数、平均值、方差与标准差，一共构成4维的特征向量；
（14）Mapping_quality：位点的碱基质量值，最后取中位数、平均值、方差与标准差，一共构成4维的特征向量；
（15）Fragment_length：支持alt和ref片段长度之间的差值，取中位数、平均值、方差、标准差，一共构成4维的特征向量；
（16）Read_position：alt位点到reads首端的距离取中位数、平均值、方差、标准差，一共构成4维的特征向量；
（17）seqs：alt位点在reads中的上下游碱基，通过k-间隔核苷酸二聚体的组成编码来对reads的二聚体进行编码；
（18）family_size：通过bam文件提取该位点的reads在其同一UMI内的最大碱基深度、位点的reads在其同一UMI内的最小碱基深度、位点的reads在其同一UMI内的reads与该位点reads不一致的碱基比例，最后分别取中位数、平均值、方差与标准差，一共构成36维的特征向量；
（19）AS：比对器输出的比对分数；
（20）UM：与ref片段的编辑距离；
（21）UQ：片段的可能性；
采用Stacking策略构建SNV位点检测模型，所述Stacking策略的第一层包括两个LightGBM算法模型，第二层为逻辑回归算法学习器；
利用两个LightGBM算法模型分别对第一特征集和第二特征集进行训练，将第一层的两个不同的LightGBM模型对输入的SNV位点数据集给出的预测结果进行组合并形成新的SNV位点数据集，再作为训练集输入到第二层逻辑回归算法学习器进行训练，训练得到SNV位点检测模型，输出最终预测结果。
2.根据权利要求1所述的ctDNA单核苷酸变异位点检测方法，其特征在于，所述第一特征集包括编码后的特征（1）- 特征（11）；所述第二特征集包括编码后的特征（12）-特征（21）。
3.根据权利要求2所述的ctDNA单核苷酸变异位点检测方法，其特征在于，所述Stacking策略采用五折交叉验证。
4.根据权利要求1所述的ctDNA单核苷酸变异位点检测方法，其特征在于，数据预处理时：
过滤测序数据包括：去除测序序列上的接头序列及低质量的碱基片段；
使用FGBIO软件识别序列上的UMI信息，并根据UMI信息对序列进行分组、合并得到一致性序列；
序列比对包括：将一致性序列通过软件对hg19人类参考基因组进行序列比对，并生成两个bam文件；然后通过GATK4的MergeBamAlignment模块合并上述的两个bam文件；然后对合并的bam文件进行文件排序式转换，并生成最终的bam文件。
5.根据权利要求1-4任一项所述的ctDNA单核苷酸变异位点检测方法，其特征在于，所述检测方法还包括：
对特征进行重要性排序，用于不同特征对SNV位点检测的贡献大小；以及，
SNV检测模型可视化解释，利用SHAP算法对特征进行排序打分，给出特征中的特征值的高低对算法预测的影响方向。
6.一种ctDNA单核苷酸变异位点检测装置，利用权利要求1-5任一项所述的ctDNA单核苷酸变异位点检测方法检测循环肿瘤DNA单核苷酸变异位点，其特征在于，包括：
数据接收模块，其被配置为用于接收SNV位点数据；
数据预处理模块，其被配置为用于数据预处理，包括过滤测序数据、识别UMI信息并根据UMI信息得到一致性序列以及序列比对；
特征提取模块，其被配置为用于提取SNV位点数据中每个位点特征，并对所提取的特征进行特征编码，将编码后的特征划分为第一特征集和第二特征集；
模型构建模块，其被配置为采用Stacking策略构建SNV位点检测模型，所述Stacking策略的第一层包括两个LightGBM算法模型，第二层为逻辑回归算法学习器；
利用两个LightGBM算法模型分别对第一特征集和第二特征集进行训练，将第一层的两个不同的LightGBM模型对输入的SNV位点数据集给出的预测结果进行组合并形成新的SNV位点数据集，再作为训练集输入到第二层逻辑回归算法学习器进行训练，训练得到SNV位点检测模型，输出最终预测结果。
7.一种计算机可读存储介质，其特征在于，所述存储介质中存储有计算机指令，当所述计算机指令被处理器执行时，能够执行权利要求1-5任一项所述的ctDNA单核苷酸变异位点检测方法。
8.一种电子设备，其特征在于，包括：
权利要求7所述的计算机可读存储介质；
以及处理器，所述处理器能够执行存储于计算机可读存储介质中的计算机指令。
CN202210154417.8A 2022-02-21 2022-02-21 ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备 Active CN114242158B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210154417.8A CN114242158B (zh) 2022-02-21 2022-02-21 ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210154417.8A CN114242158B (zh) 2022-02-21 2022-02-21 ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备

Publications (2)

Publication Number Publication Date
CN114242158A CN114242158A (zh) 2022-03-25
CN114242158B true CN114242158B (zh) 2022-05-13

Family

ID=80747558

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210154417.8A Active CN114242158B (zh) 2022-02-21 2022-02-21 ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备

Country Status (1)

Country Link
CN (1) CN114242158B (zh)

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11443832B2 (en) * 2019-03-07 2022-09-13 Nvidia Corporation Genetic mutation detection using deep learning
CN110010197B (zh) * 2019-03-29 2021-07-20 深圳裕策生物科技有限公司 基于血液循环肿瘤dna的单核苷酸变异检测方法、装置和存储介质
US20200342956A1 (en) * 2019-04-25 2020-10-29 Carnegie Mellon University Methods and Systems for Use in Cancer Prediction
CN113223613A (zh) * 2021-05-14 2021-08-06 西安电子科技大学 基于多维度单核苷酸变异特征的癌症检测方法

Also Published As

Publication number Publication date
CN114242158A (zh) 2022-03-25

Similar Documents

Publication Publication Date Title
CN112888459B (zh) 卷积神经网络系统及数据分类方法
US20210257050A1 (en) Systems and methods for using neural networks for germline and somatic variant calling
CN109767810B (zh) 高通量测序数据分析方法及装置
US20230222311A1 (en) Generating machine learning models using genetic data
US20210090686A1 (en) Single cell rna-seq data processing
US20220277811A1 (en) Detecting False Positive Variant Calls In Next-Generation Sequencing
CN113674803A (zh) 一种拷贝数变异的检测方法及其应用
CN112951327A (zh) 药物敏感预测方法、电子设备及计算机可读存储介质
CN111180013B (zh) 检测血液病融合基因的装置
CN115715416A (zh) 基于机器学习的医学数据检查器
US20220101135A1 (en) Systems and methods for using a convolutional neural network to detect contamination
CN114242158B (zh) ctDNA单核苷酸变异位点检测方法、装置、存储介质及设备
US20200105374A1 (en) Mixture model for targeted sequencing
El Qumsan Assessment of supervised classification methods for the analysis of RNA-seq data
CN117935933B (zh) Cdkn2a/b纯合性缺失的分析方法及系统
CN117153392A (zh) 一种胃癌预后预测的标志物、评估模型及其构建方法
EP4138003A1 (en) Neural network for variant calling
Lam Characterizing Sequencing Artifacts
Wang et al. Shenjie Wang1, 2, Yuqian Liu1, 2, Juan Wang1, 2, 3*, Xiaoyan Zhu 1, 2, Yuzhi Shi3, Xuwen Wang1, 2, Tao Liu3, Xiao Xiao2, 4 and Jiayin Wang1, 2
Niehus Multi-Sample Approaches and Applications for Structural Variant Detection
Chieruzzi Identification of RAS co-occurrent mutations in colorectal cancer patients: workflow assessment and enhancement
Persson Comparing Two Algorithms for the Detection of Cross-Contamination in Simulated Tumor Next-Generation Sequencing Data
Wang et al. TMBstable: a variant caller controls performance variation across heterogeneous sequencing samples
Kariotis Unsupervised machine learning of high dimensional data for patient stratification
Miller A Method for Identification of Pancreatic Cancer Through Methylation Signatures in Cell-Free DNA

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