CN111304308A - 一种审核高通量测序基因变异检测结果的方法 - Google Patents

一种审核高通量测序基因变异检测结果的方法 Download PDF

Info

Publication number
CN111304308A
CN111304308A CN202010135146.2A CN202010135146A CN111304308A CN 111304308 A CN111304308 A CN 111304308A CN 202010135146 A CN202010135146 A CN 202010135146A CN 111304308 A CN111304308 A CN 111304308A
Authority
CN
China
Prior art keywords
detected
mutation
reads
variation
mutation site
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
CN202010135146.2A
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.)
Genetron Health Beijing Co ltd
Original Assignee
Genetron Health Beijing 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 Genetron Health Beijing Co ltd filed Critical Genetron Health Beijing Co ltd
Priority to CN202010135146.2A priority Critical patent/CN111304308A/zh
Publication of CN111304308A publication Critical patent/CN111304308A/zh
Pending legal-status Critical Current

Links

Images

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/6869—Methods for sequencing
    • 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/6876—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes
    • C12Q1/6883—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes for diseases caused by alterations of genetic material
    • C12Q1/6886—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes for diseases caused by alterations of genetic material for cancer
    • 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
    • C12Q2600/00—Oligonucleotides characterized by their use
    • C12Q2600/106—Pharmacogenomics, i.e. genetic variability in individual responses to drugs and drug metabolism
    • 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
    • C12Q2600/00—Oligonucleotides characterized by their use
    • C12Q2600/156—Polymorphic or mutational markers

Landscapes

  • Chemical & Material Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Organic Chemistry (AREA)
  • Health & Medical Sciences (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Engineering & Computer Science (AREA)
  • Zoology (AREA)
  • Wood Science & Technology (AREA)
  • Immunology (AREA)
  • Analytical Chemistry (AREA)
  • Genetics & Genomics (AREA)
  • Physics & Mathematics (AREA)
  • Biophysics (AREA)
  • Biotechnology (AREA)
  • Microbiology (AREA)
  • Molecular Biology (AREA)
  • Pathology (AREA)
  • Biochemistry (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Hospice & Palliative Care (AREA)
  • Oncology (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明公开了一种审核高通量测序基因变异检测结果的方法。该方法包括如下步骤：构建训练集，包括若干阳性变异位点和阴性非变异位点的测序数据；提取并向量化变异位点特征；采用随机森林法构建模型，然后利用模型判断待测位点变异是否为变异位点。本发明方法能快速、准确完成原需人工检视突变的环节，大大节约人力成本、提高整体变异检测效率；达到甚至在一定程度上超越了人工审核的准确度，并能大大降低人工检验的主观性；使用机器学习算法，在该领域具有一定创新；基于真实的临床数据开发模型，更贴近于真实案例，并在金标准数据集上表现与人工结果无出其二，并能适配各种不同的细胞变异检测软件，具有良好的实用性。

Description

一种审核高通量测序基因变异检测结果的方法
技术领域
本发明涉及生物信息学领域，具体涉及一种审核高通量测序基因变异检测结果的方法。
背景技术
目前，在体细胞变异检测中，二代测序技术已得到了极为广泛的应用，并为临床中的用药指导、预后监测等提供了必要的基础。因而变异检测的准确性成为精准医疗的一个必要前提。然而，精确地检测体细胞变异(尤其是低频变异)面临着样本制备损坏、测序错误、参考基因组序列比对错误等一系列问题的严峻挑战。用于体细胞变异检测的生物信息学软件则是下游、同时也是十分重要的一环。目前众多优秀的变异检测软件(如Mutect,Strelka,Varscan等)在体细胞变异检测市场上得到了广泛的应用。然而，面对复杂的变异情况，一般常用的变异检测软件检出来的变异结果都存在一定比例的错误(从15％到30％不等)，故都需要通过进一步的可视化的变异呈现软件(例如最为常用软件为IGV(Integrative Genome Viewer)[http://www.igv.org/])来进行人工的审核、排除假阳并给出最终的真实变异结果。然而，该人工审核过程存在三方面的不足：一、使用可视化软件审核变异并给出较为准确的结果，需要比较深厚的生物信息学专业知识背景和长期的训练，同时也存在较强的主观判断因素；二、在测序数据越来越庞大的情况下，人工审核过程效率和人工成本上在一定程度上制约着变异检测的时效性和经济性；三、人工审核的主观性、审核标准的不统一性也对变异检测的准确性构成了一定的威胁。由于IGV的审核规则存在较大的主观性，目前Erica K.等对其判定准则的规范化进行了系统性的研究整理和测试(Standard operating procedure for somatic variant refinement of sequencingdata with paired tumor and normal samples DOI:10.1038/s41436-018-0278-z)。该研究通过一批金标准数据集，确立并验证了一系列关于人工审核变异的标准。实验中专业人员通过学习该套标准，其对体细胞变异判读的准确性从77.4％提高到了94.1％。尽管IGV审核规则得到了专业性的确立和验证，然而，人工审核时效低的问题还未得到有效的解决。
发明内容
针对上述问题，本发明的目的是利用机器学习算法及人工标注的金标准数据集，提供一种用来替代或者自动化人工审核体细胞变异的方法，即一种审核基因变异检测结果的方法。
第一方面，本发明要求保护一种审核高通量测序基因变异检测结果的方法。
本发明所要求保护的审核高通量测序基因变异检测结果的方法，可包括如下步骤：
(A)构建训练集，包括若干阳性变异位点和阴性非变异位点的测序数据。
(B)从训练集中提取并向量化变异位点的特征；所述变异位点的特征包括以下任意6种或6种以上，例如任意7种、8种、9种、10种或11种：
突变支持数：肿瘤组织中支持待检测突变位点的reads数；
突变频率：肿瘤组织中待检测突变位点的频率；
碱基质量：肿瘤组织中支持待检测突变位点的reads中待检测突变位点的碱基质量的均值；
错误比对率：肿瘤组织中支持待检测突变位点的reads中待检测突变位点的上下游50bp序列错误比对率的均值；
HDR值：肿瘤组织中支持待检测突变位点的reads中疑似同源性比对错误得分；
SideB值：肿瘤组织中支持待检测突变位点的reads中边缘偏好性得分；
StrandB值：肿瘤组织中支持待检测突变位点的reads的链偏好性得分；
碱基比率：正常对照组织和肿瘤组织中支持待检测突变位点reads数的比值；
其他变异类型比率：肿瘤组织中如有某种变异类型外的其他变异类型，所述其他变异类型与所述某种变异类型的支持reads数之比；
比对长度：肿瘤组织中支持待检测突变位点reads的平均长度；
InDel长度：Indel的缺失或插入长度；
(C)利用步骤(B)得到的特征结果，采用随机森林法构建模型，然后利用模型判断待测位点是否为变异位点。
上述方法中，所述变异位点的特征还可包括以下任意一种或一种以上，例如任意两种、三种、四种、五种或六种：
比对质量：肿瘤组织中支持待检测突变位点的reads的比对质量的均值；
环境质量：肿瘤组织中支持待检测突变位点的reads中待检测突变位点的上下游各50bp序列的碱基质量的均值；
插入片段大小：肿瘤组织中支持待检测突变位点的reads的插入片段大小的均值；
基因组复杂度得分：待检测突变位点上下游各20bp的参考基因组序列的复杂度得分；
正常覆盖深度：正常对照组织中待检测突变位点的覆盖深度；
正常InDel存在率：肿瘤组织和/或正常组织中待检测突变位点的上下游各50bp是否有InDel；如有，则计算该InDel的变异频率和长度的乘积，或进一步将乘积相加；如无则为0。例如，如果肿瘤组织或正常组织中待检测突变位点的上下游各50bp有InDel，则计算该InDel的变异频率和长度的乘积，如果肿瘤组织和正常组织中待检测突变位点的上下游各50bp都有InDel，则计算该InDel的变异频率和长度的乘积后进行相加。
上述方法中，所述基因变异的类型为碱基置换或InDel；所述模型根据变异类型构建，利用同一变异类型的变异位点的测序数据构建训练集，提取并向量化特征并构建模型。例如，利用变异类型为碱基置换的变异位点的测序数据构建训练集，提取并向量化特征后，构建模型。
在本发明的具体实施方式中，当所述基因变异的类型为碱基置换时，所述变异位点的特征，优选的，包括：突变支持数、突变频率、碱基质量、错误比对率、HDR值、SideB值、StrandB值和碱基比率；更优选的，包括突变支持数、突变频率、碱基质量、错误比对率、HDR值、SideB值、StrandB值、碱基比率和比对质量。
在本发明的具体实施方式中，当所述基因变异的类型为InDel时，所述变异位点的特征，优选的，包括突变支持数、突变频率、错误比对率、HDR值、SideB值、StrandB值、碱基比率、其他变异类型比率、InDel长度和比对质量。
其中所述的碱基置换可以是单碱基(SNV)、双碱基(DNV)和/或三碱基(TNV)的置换。
上述方法中，所述HDR值的计算方法为：
当支持待检测突变位点的reads数小于5时，或者当支持待检测突变位点的reads数大于等于5、支持待检测突变位点的reads上下游某个变异在正常对照组织中与肿瘤组织中的变异频率之差的绝对值小于0.7时，HDR值为0；
当支持待检测突变位点的reads数大于等于5、支持待检测突变位点的reads上下游某个变异在正常对照组织中与肿瘤组织中的变异频率之差的绝对值大于等于0.7时，HDR值计算方式为：
Figure BDA0002397010850000031
其中，Ii等于1；n为符合HDR条件的位点个数。
HDR值越小，表明疑似同源性的比对错误概率越高。
上述方法中，所述SideB值的计算方式为：
当支持待检测突变位点的reads数大于等于10时，所述SideB值按以下公式计算：
Figure BDA0002397010850000032
其中，n＝50，Dti是表示突变位点上游第i个bp的支持reads的深度；Dni是突变位点下游第i个bp的支持reads的深度。abs表示取绝对值。
当支持待检测突变位点的reads数小于10时，所述SideB值为0。
SideB值得分越低，则说明该位点越具有边缘偏好的特点。
上述方法中，所述StrandB值的计算方式为：
当支持待检测突变位点的reads数大于等于10时，所述StrandB值按以下公式计算：
Figure BDA0002397010850000041
其中，R为待检测突变位点支持reads中正链的比例。abs表示取绝对值。
当支持待检测突变位点的reads数小于10时，所述StrandB值为0。
StrandB值越高，则表明该突变越不具备链偏好的特点。
进一步地，所述复杂度得分的计算方法为：如果参考基因组序列上待检测突变位点的上下游20bp范围内出现6个或6个以上连续的单碱基重复序列或串联重复序列(例如‘AAAAAA’或‘ATCATCATCATCATCATCATC’)，则该特征值为1，反之则为0,；所述连续的单碱基重复序列即一个相同碱基连续重复出现，例如‘AAAAAA’；所述连续的串联重复序列中，所述串联序列指两个或两个以上碱基串联的序列，例如‘AT’或‘ATC’，所述连续的串联重复序列，例如‘ATCATCATCATCATCATCATC’。
在所述方法在步骤(B)从训练集中提取并向量化变异位点的特征前，还可包括按照如下对训练集中的支持待检测突变位点reads进行过滤的步骤：过滤掉具有如下三种情况中至少一种的reads：reads长度低于55、reads比对质量小于10、软截断(Soft Clips)长度大于20bp。
在本发明的具体实施方式中，是采用python2.7.8的pysam(v0.11)对候选变异位点的成对样本中的reads进行的过滤排除。
通过上述步骤，可过滤掉不可靠的reads，防止对输入的训练集引入不必要的噪音。
进一步地，模型的构建可利用R语言的randomForest程序包(如R3.5.1中的randomForest_4.16-14包)完成。
更进一步地，使用十折交叉验证法对于模型的主要参数(如树的个数(ntree)及随机抽取特征数(mtry)构建梯度)进行优化。
更加具体地，对于碱基置换变异而言，优化后的主要参数如下：ntree设定为900；mtry设定为5；对于InDel变异而言，优化后参数如下：ntree设定为2000；mtry设定为4。
在所述方法的步骤(A)中，所述阳性变异位点和阴性非变异位点是指按照现有软件进行检测后，通过注释，人工审核后的阳性变异位点或阴性非变异位点。
其中，所述检测包括，使用现有的变异检测软件对测序数据进行检测，例如，对于碱基置换变异使用GATK3.1中的Mutect1(版本3.1)进行检测，对于InDel变异(插入和缺失)使用Strelka(v1.0.14)进行检测(此处使用上述两个软件进行变异检测，但该算法在适用方面不限于上述两个软件的检测结果)。
对变异进行检测后，例如，使用VEP(v83)对变异结果进行注释。
对变异结果进行注释后，可选的，使用下列过滤条件对候选变异进行过滤：过滤WES(全外显子)区域外的变异；过滤掉同义突变；对照COSMIC数据库(https://cancer.sanger.ac.uk/cosmic/)的非热点变异，过滤掉变异频率低于5％的变异；对于COSMIC数据库中收录的热点变异，过滤掉变异频率低于1％的变异；过滤掉常用千人基因数据库中位点变异频率大于1％的变异(疑似common variants)；过滤掉低覆盖位点(具体为tumor中位点覆盖小于40X，或normal中覆盖小于10X)；过滤掉绝对支持数低于或等于6条reads的变异。
其中，所述常用千人基因数据库中的位点变异频率大于1％的变异位点包括：
NHLBI-ESP数据库中人群频率大于1％的变异位点(https://evs.gs.washington.edu/EVS/)；
千人基因组计划数据库中人群频率大于1％的变异位点(http://phase3browser.1000genomes.org/index.html)；
千人基因组计划数据的东亚人群中频率大于1％的变异位点(http://phase3browser.1000genomes.org/index.html)。
经过上述过滤环节，大量难以确定的低频变异及临床意义不明显的变异将被过滤掉。
最后，通过人工审核得出真、假变异，作为阳性变异位点和阴性非变异位点；例如，通过IGV对候选变异进行真、假标签的评估的。
所述基因变异指的是人类肿瘤组织和正常对照组织相比所发生的基因变异，可为碱基置换或者InDel。
在本发明中，所述碱基置换包括SNV变异、DNP变异和TNP变异的情况，即单碱基变异、双碱基变异和/或三碱基变异。
第二方面，本发明要求保护一种审核高通量测序基因变异检测结果的系统。
本发明所要求保护的审核基因变异检测结果的系统，包括装置A、装置B和装置C；
所述装置A能够按照前文所述方法中的步骤(A)构建包括若干阳性变异位点和阴性非变异位点的测序数据的训练集；
所述装置B能够按照前文所述方法中的步骤(B)从所述训练集中提取并向量化所述变异位点的特征；
所述装置C能够按照前文所述方法中的步骤(C)用所述装置B提取并向量化的所述变异位点的特征，采用随机森林法构建模型，然后利用模型判断待测位点是否为变异位点。
第三方面，本发明要求保护如下任一应用：
(I)照前文第一方面所述的方法或照前文第二方面所述的系统在制备用于肿瘤早筛、肿瘤预后、肿瘤分类和/或肿瘤用药指导的产品中的应用；
所述产品能够按照前文第一方面所述方法的步骤实现对高通量测序基因变异检测结果进行审核。
(II)前文第一方面所述的方法或照前文第二方面所述的系统在肿瘤早筛、肿瘤预后、肿瘤分类和/或肿瘤用药指导中的应用；
(III)前文第一方面所述的方法或照前文第二方面所述的系统在基因变异检测中的应用。
本发明由于采取以上技术方案，其具有以下优点：
(1)该算法能够快速、准确的完成原需人工检视突变的环节，可大大节约人力成本、提高整体变异检测流程的效率；
(2)该算法在一定程度上超越了人工审核的准确度，并能够大大降低人工检验的主观性；
(3)该算法使用机器学习算法，在该领域上具有一定的创新；
(4)该算法基于真实的临床数据开发模型，更贴近于真实案例，并在金标准数据集上表现与人工结果无出其二，并能适配各种不同的细胞变异检测软件，具有良好的实用性。
附图说明
图1为本发明审核基因变异检测结果方法的算法流程图。
图2为筛选特征时在训练集上进行交叉验证计算的错误率曲线图。
图3为碱基置换和InDel模型的AUC指标。
图4为碱基置换和InDel模型在独立验证集上的预测。
具体实施方式
下述实施例中所使用的实验方法如无特殊说明，均为常规方法。
下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。
实施例1、本发明审核基因变异检测结果方法的建立及应用
根据如图1所示的流程图，模型构建方法包括以下步骤：
一、训练数据准备
为了采用监督学习算法(随机森林)，将人工审核体细胞变异(肿瘤组织、正常对照组织样本)这一过程自动化，选取了94例癌症的全外显子测序样本，进行Illumina测序，平均深度约为Tumor 200X，Normal 100X。
使用比对软件bwa-0.7.10对测序后的fastq文件进行序列比对，然后对输出的bam数据进行变异检测，例如，其中，对于SNP、DNP、TNP(单碱基、双碱基和三碱基变异)使用GATK3.1中的Mutect1(版本3.1)进行检测，对InDel(插入和缺失)使用Strelka(v1.0.14)进行检测(此处使用上述两个软件进行变异检测，但该算法在适用方面不限于上述两个软件的检测结果)。之后使用VEP(v83)对变异结果进行注释。在对注释结果进行最终人工审核之前，根据突变的临床意义使用了下列过滤条件对候选变异进行了过滤：
1、过滤WES(全外显子)区域外的变异；
2、过滤掉同义突变；
3、对照COSMIC数据库(https://cancer.sanger.ac.uk/cosmic/)的非热点变异，过滤掉变异频率低于5％的变异；
4、对于COSMIC数据库中收录的热点变异，过滤掉变异频率低于1％的变异；
5、过滤掉常用千人基因数据库中位点变异频率大于1％的变异(疑似commonvariants)。
其中，所述常用千人基因数据库中的位点变异频率大于1％的变异位点包括：
NHLBI-ESP数据库中人群频率大于1％的变异位点(https://evs.gs.washington.edu/EVS/)；
千人基因组计划数据库中人群频率大于1％的变异位点(http://phase3browser.1000genomes.org/index.html)；
千人基因组计划数据的东亚人群中频率大于1％的变异位点(http://phase3browser.1000genomes.org/index.html)。
6、过滤掉低覆盖位点，具体为tumor中位点覆盖小于40X，或normal中覆盖小于10X的；
7、过滤掉绝对支持数低于或等于6条reads的变异。
经过上述过滤环节，大量难以确定的低频变异及临床意义不明显的变异将被过滤掉。
上述候选变异真实标签的确定是通过两位研究人员确定的，其具有两年或两年以上相关生物信息学工作经验的一批研究者(本领域普通技术人员)，通过IGV SOP的学习[Standard operating procedure for somatic variant refinement of sequencingdata with paired tumor and normal samples DOI:10.1038/s41436-018-0278-z]，规范并统一了IGV评判规则。之后对上述候选变异通过IGV进行了真、假标签的评估。评估结束后将两位人员不一致的位点去除，剩余的作为用来训练监督学习模型的金标准位点。该步骤结束后，共有11223个变异，其中SNP为9181个，Indel为2042个。
真实的体细胞变异应具有足够的分子数(reads)支持(指扩增、比对、去重后具有唯一标识的分子)，一般要求要≥7条reads；而常见的错误类型变异则会有以下几种情形：
(1)疑似胚系变异。当正常对照细胞(通常为白细胞)和肿瘤组织中同时检测出相同的变异类型，且能够较大概率下排除肿瘤细胞污染正常组织细胞的可能(该变异在正常组织细胞中的变异频率也偏高)，则该变异是胚细胞变异的可能性大。
(2)正常组织细胞中无reads覆盖或覆盖深度低。当某候选变异位点对应的正常细胞中没有测序覆盖或覆盖深度小于10层时，因无法判定该变异在正常组织细胞中是否有支持，故该种情况的位点一般会被归类为假阳。
(3)癌组织细胞中无reads覆盖或覆盖深度低时(小于40层)，位点一般被归类为假阳。
(4)支持reads有高度的链偏好性(比如大于90％的reads都是正链或者负链)；或(且)有一定比例的reads上存在比对错误的碱基。此种情况下的位点标签较难以判定。
(5)支持的reads上同时有众多变异碱基，但对应的正常组织细胞中并不存在该变异。此种情况多为比对错误或参考基因组序列有误引起。
(6)多种变异类型同时出现且都占有一定的比例。
(7)变异频率太低。一般变异频率在5％以下的都存疑，但这个阈值与测序深度有较大的相关性，因实验而异。本发明肿瘤平均测序深度约200X，统一按照5％，突变热点1％的标准过滤。
(8)支持的每条reads上都有不统一的、零星的碱基变异。
(9)比对质量过低。例如超过60％的reads的比对质量都小于30。
(10)碱基质量过低。例如支持reads的平均碱基质量不到20。
(11)支持的reads上邻近部位(上下游50bp)有Indel。此种情况增大了比对错误的概率。
(12)边缘偏好性：即所有的支持reads中，变异位点都处于reads的同侧末端。
(13)邻近序列复杂性低。
二、过滤reads
该过滤步骤使用python2.7.8的pysam(v0.11)对候选变异位点的成对样本中的reads进行了过滤排除。其过滤条件包括如下情形：
(1)Reads长度低于55的；
(2)Reads比对质量小于10的；
(3)软截断长度大于20bp的。
通过上述步骤，可过滤掉不可靠的reads，防止对输入的训练集引入不必要的噪音。
三、提取并向量化变异位点的特征
对每个变异位点的下列相关特征进行了提取并向量化。即：基于上述过滤步骤后的reads，结合人工审核时常用的位点判读考量因素，计算单个位点及其上下游各50bp位点的相关特征。相关特征及计算方式如下：
突变支持数(mut_Count)：肿瘤组织中支持突变位点的reads数。
突变频率(mut_Freq)：肿瘤组织中待检测突变位点的频率(突变reads支持数除以该位点总深度)。
比对质量(mapping_Quality)：肿瘤组织中支持突变reads的比对质量的均值。
碱基质量(base_Quality)：肿瘤组织中支持突变reads中突变位点碱基质量的均值。
环境质量(BQ_average)：肿瘤组织中支持待检测突变位点的reads中待检测突变位点的上下游各50bp序列的碱基质量的均值。
错误比对率(misMatch_Average)：肿瘤组织中支持突变reads中突变位点上下游50bp序列错误比对率的均值。
HDR值(HDR_score)：肿瘤组织中支持突变reads中疑似同源性比对错误得分。
同源性比对错误常导致在支持突变的reads上下游其他碱基位置上有数个突变频率相同或相近的变异，进而易导致该位点假阳性的判断。只有支持突变的reads数大于等于5时，才计算该值，否则HDR为0。
根据这种错误特点，可先统计出现上下游位置变异的个数，其后计算相应的HDR值。当上下游某个变异在正常对照组织中与肿瘤组织中的变异频率之差的绝对值大于等于0.7时，认为这个变异位置符合HDR的情况，HDR计算方式为：
Figure BDA0002397010850000091
其中，Ii等于1，n为符合HDR条件的位点个数。
当上下游某个变异在正常对照组织中与肿瘤组织中的变异频率之差的绝对值小于0.7时，认为这个变异位置不符合HDR的情况，HDR为0。
HDR值越小，表面疑似同源性的比对错误概率越高。
SideB值(sideBias)：肿瘤组织中支持突变reads中边缘偏好性得分，计算方式为：
当支持突变的reads数大于等于10时：
Figure BDA0002397010850000101
其中，n＝50，Dti是表示突变位点上游第i个bp的支持reads的深度；Dni是突变位点下游第i个bp的支持reads的深度。abs表示取绝对值。
当支持突变的reads数小于10时，SideB为0。
SideB得分越低说明该位点越具有边缘偏好的特点。
插入片段大小(insertSize)：肿瘤组织中支持待检测突变位点的reads的插入片段大小的均值；
StrandB值(strandBias)：肿瘤组织中支持突变reads的链偏好性得分，统计这些reads中正链的比例R，计算链偏好性得分值。
对于支持突变的reads数大于等于10条的变异位点：
Figure BDA0002397010850000102
其中，R为支持突变reads中正链的比例。abs表示取绝对值。
当支持突变的reads数小于10时，StrandB为0。
StrandB得分越高表明该突变越不具备链偏好的特点。
基因组复杂度得分(repeative_Flag)：突变位置上下游各20bp的参考基因组序列的复杂度得分。
如果参考基因组序列上待检测突变位点的上下游20bp范围内出现6个或6个以上连续的单碱基重复序列或串联重复序列(例如‘AAAAAA’或‘ATCATCATCATCATCATCATC’)，则该特征值为1，反之则为0；所述连续的单碱基重复序列即一个相同碱基连续重复出现，例如‘AAAAAA’；所述连续的串联重复序列中，所述串联序列指两个或两个以上碱基串联的序列，例如‘AT’或‘ATC’，所述连续的串联重复序列，例如‘ATCATCATCATCATCATCATC’。
碱基比率(ntRatio)：正常对照组织和肿瘤组织中支持突变reads的比值。
其他变异类型比率(var_TypeRatio)：肿瘤组织中如有某种变异类型外的其他变异类型，所述其他变异类型与所述某种变异类型的支持reads数之比。
正常覆盖深度(normal_Coverage)：正常对照组织突变位点的覆盖深度。
比对长度(query_Length)：肿瘤组织中支持突变reads的平均长度。
正常InDel存在率(normal_Indels)：肿瘤组织和/或正常组织中待检测突变位点的上下游各50bp是否有InDel；如有，则计算该InDel的变异频率和长度的乘积，或进一步将乘积相加；如无则为0。如果肿瘤组织或正常组织中待检测突变位点的上下游各50bp有InDel，则计算该InDel的变异频率和长度的乘积，如果肿瘤组织和正常组织中待检测突变位点的上下游各50bp都有InDel，则计算该InDel的变异频率和长度的乘积后进行相加。
InDel长度(indel_Lengths)：InDel的缺失或插入长度。
四、模型训练、评估及预测
由于碱基置换和InDel变异特征具有较明显的差异性，例如InDel的长短和本身序列的复杂度对位点本身的真假影响较大、串联重复区域的InDel错误率更高等，在衡量位点真假中碱基置换和InDel具有不同的参考标准和阈值，所以我们将碱基置换和InDel训练集划分为两套独立数据分别训练了模型。
1、模型构建过程
将步骤三中提取的相关特征的数据进行训练集和验证集的划分。具体为：在碱基置换(包括SNV、DNV和TNV)和InDel中分别随机抽选了6181例和1430个变异位点作为训练集，其余位点作为测试集，如表1所示。
表1模型构建训练集和测试集的划分
模型 训练集 测试集
碱基置换 6181 3000
Indel 1430 612
接下来碱基置换和Indel的模型构建和优化分别在6181和1430个训练集上进行，模型性能的评估在测试集上进行。模型的构建使用R3.5.1中的randomForest_4.16-14包进行。
(1)碱基置换模型构建中，使用的特征如表2所示。将表2中的特征从上往下进行数量上的递进(例如，模型取3个特征时，所选取的特征为“mismatch_Average”、“ntRatio”、“mutCount”)，在训练集上使用5折交叉验证法对各梯度下模型的错误率进行评估(使用函数为“randomForest::rfcv”)，结果如图2中左图所示，横坐标为特征个数，纵坐标为交叉验证错误率。可见，对于碱基置换模型来说，当特征个数增长到8个时，错误率趋势曲线趋于平缓。
使用上述不同个数特征下产生的模型对测试集进行预测，评估曲线下面积AUC并计算模型的准确度，梯度1、2特征个数太少未纳入分析，结果如表3所示。
表2碱基置换模型构建的特征
序号 特征名称
1 misMatch_Average
2 ntRatio
3 mut_Count
4 HDR_score
5 Base_Quality
6 strandBias
7 sideBias
8 Mut_Freq
9 Mapping_Quality
10 BQ_average
11 insertSize
12 normal_Coverage
13 normal_Indels
14 repeative_Flag
表3不同特征个数的模型效果
特征个数 测试集AUC 测试集准确度
3 0.981±0.005 0.939
4 0.982±0.004 0.939
5 0.989±0.003 0.955
6 0.992±0.003 0.96
7 0.992±0.003 0.96
8 0.992±0.003 0.961
<u>9</u> <u>0.993±0.003</u> <u>0.963</u>
10 0.994±0.003 0.964
11 0.992±0.003 0.965
12 0.992±0.003 0.967
13 0.995±0.002 0.968
14 0.995±0.002 0.967
(2)InDel模型构建中，使用的特征如表4所示。将表4中的特征从上往下进行数量上的递进，在训练集上使用5折交叉验证法对各梯度下模型的错误率进行评估(使用函数为“randomForest::rfcv”)，结果如图2中右图所示，横坐标为特征个数，纵坐标为交叉验证错误率。可见，对于InDel模型来说，当特征个数增长到10个时，错误率趋势曲线趋于平缓。
使用上述不同个数特征下产生的模型对测试集进行预测，评估曲线下面积AUC并计算模型的准确度，梯度1、2特征个数太少未纳入分析，结果如表5所示。
表4 InDel模型构建的特征
Figure BDA0002397010850000121
Figure BDA0002397010850000131
表5不同特征个数的模型效果
特征个数 测试集AUC 测试集准确度
3 0.986±0.009 0.943
4 0.991±0.007 0.948
5 0.991±0.007 0.951
6 0.992±0.007 0.948
7 0.992±0.007 0.953
8 0.992±0.007 0.956
9 0.992±0.007 0.951
<u>10</u> <u>0.994±0.006</u> <u>0.962</u>
11 0.993±0.007 0.964
12 0.993±0.006 0.966
13 0.994±0.006 0.971
14 0.994±0.006 0.969
15 0.994±0.006 0.969
16 0.994±0.006 0.972
对于上述模型构建的主要参数如树的个数(ntree)及随机抽取特征数(mtry)构建梯度，在训练集上使用十折交叉验证法对其进行了优化。优化后最佳参数详细见表6。
表6 SNV和Indel模型最优参数
模型 最佳ntree 最佳mtry
碱基置换 900 5
Indel 2000 4
2、测试集上模型的表现
用上述步骤中训练好的模型对独立测试集进行了预测和模型评估。
将待审核位点按照步骤三提取并向量化位点特征，然后用预测模型进行预测，根据阈值判断阴性或阳性，例如，取阈值为0.5时，低于阈值(0.5)为阴性，高于或等于阈值(0.5)为阳性。
依据图2中的交叉验证错误率曲线、测试集上的AUC指标以及参考特征本身的意义三方面因素，选择碱基置换的特征个数为9、InDel特征个数梯度为10的模型，如表3和5所示，其AUC分别为0.993±0.003、0.994±0.006。各梯度下独立测试集上模型的准确度如表4和6所示。其中，准确率＝(预测真阳数+预测真阴数)/所有样本数。
在取阈值为0.5时，灵敏度和特异度详见表7。碱基置换和InDel模型的AUC值分别为0.993±0.002和0.994±0.006，详见图3。测试集预测得分值的分别线箱图见图4。图4所示，两个模型在独立测试集上对位点阴、阳性的鉴别效果良好。
表7碱基置换和Indel模型在测试集上的表现指标
Figure BDA0002397010850000141
注：表中标注的真实阳性和真实阴性是以前文所述IGV法进行评估的。
上述各实施例仅用于说明本发明，其中关于模型特征计算、算法调优都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (10)

1.一种审核高通量测序基因变异检测结果的方法，包括如下步骤：
(A)构建训练集，包括若干阳性变异位点和阴性非变异位点的测序数据；
(B)从训练集中提取并向量化变异位点的特征；所述变异位点的特征包括以下任意6种或6种以上：
突变支持数：肿瘤组织中支持待检测突变位点的reads数；
突变频率：肿瘤组织中待检测突变位点的频率；
碱基质量：肿瘤组织中支持待检测突变位点的reads中待检测突变位点的碱基质量的均值；
错误比对率：肿瘤组织中支持待检测突变位点的reads中待检测突变位点的上下游50bp序列错误比对率的均值；
HDR值：肿瘤组织中支持待检测突变位点的reads中疑似同源性比对错误得分；
SideB值：肿瘤组织中支持待检测突变位点的reads中边缘偏好性得分；
StrandB值：肿瘤组织中支持待检测突变位点的reads的链偏好性得分；
碱基比率：正常对照组织和肿瘤组织中支持待检测突变位点reads数的比值；
其他变异类型比率：肿瘤组织中如有某种变异类型外的其他变异类型，所述其他变异类型与所述某种变异类型的支持reads数之比；
比对长度：肿瘤组织中支持待检测突变位点reads的平均长度；
InDel长度：Indel的缺失或插入长度；
(C)利用步骤(B)得到的特征结果，采用随机森林法构建模型，然后利用模型判断待测位点变异是否为变异位点。
2.根据权利要求1所述的方法，其特征在于，所述变异位点的特征还包括以下任意一种或一种以上：
比对质量：肿瘤组织中支持待检测突变位点的reads的比对质量的均值；
环境质量：肿瘤组织中支持待检测突变位点的reads中待检测突变位点的上下游各50bp序列的碱基质量的均值；
插入片段大小：肿瘤组织中支持待检测突变位点的reads的插入片段大小的均值；
基因组复杂度得分：待检测突变位点上下游各20bp的参考基因组序列的复杂度得分；
正常覆盖深度：正常对照组织中待检测突变位点的覆盖深度；
正常InDel存在率：肿瘤组织和/或正常组织中待检测突变位点的上下游各50bp是否有InDel；如有，则计算该InDel的变异频率和长度的乘积，或进一步将乘积相加；如无则为0。
3.根据权利要求1或2所述的方法，其特征在于：所述基因变异的类型为碱基置换或InDel；所述模型根据变异类型构建，利用同一变异类型的变异位点的测序数据构建训练集，提取并向量化特征并构建模型。
4.根据权利要求1-3中任一所述的方法，其特征在于：所述HDR值的计算方法为：
当支持待检测突变位点的reads数小于5时，或者当支持待检测突变位点的reads数大于等于5、支持待检测突变位点的reads上下游某个变异在正常对照组织中与肿瘤组织中的变异频率之差的绝对值小于0.7时，HDR值为0；
当支持待检测突变位点的reads数大于等于5、支持待检测突变位点的reads上下游某个变异在正常对照组织中与肿瘤组织中的变异频率之差的绝对值大于等于0.7时，HDR值计算方式为：
Figure FDA0002397010840000021
其中，Ii等于1；n为符合HDR条件的位点个数。
5.根据权利要求1-4中任一所述的方法，其特征在于：所述SideB值的计算方式为：
当支持待检测突变位点的reads数大于等于10时，所述SideB值按以下公式计算：
Figure FDA0002397010840000022
其中，n＝50，Dti是表示突变位点上游第i个bp的支持reads的深度；Dni是突变位点下游第i个bp的支持reads的深度；
当支持待检测突变位点的reads数小于10时，所述SideB值为0。
6.根据权利要求1-5中任一所述的方法，其特征在于：所述StrandB值的计算方式为：
当支持待检测突变位点的reads数大于等于10时，所述StrandB值按以下公式计算：
Figure FDA0002397010840000023
其中，R为待检测突变位点支持reads中正链的比例；
当支持待检测突变位点的reads数小于10时，所述StrandB值为0。
7.根据权利要求2-6中任一所述的方法，其特征在于：所述基因组复杂度得分的计算方法为：如果参考基因组序列上待检测突变位点的上下游20bp范围内出现6个或6个以上连续的单碱基重复序列或串联重复序列，则该特征值为1，反之则为0。
8.根据权利要求1-7中任一所述的方法，其特征在于：所述方法在步骤(B)从训练集中提取并向量化变异位点的特征前，还包括按照如下对训练集中的支持待检测突变位点reads进行过滤的步骤：过滤掉具有如下三种情况中至少一种的reads：reads长度低于55、reads比对质量小于10、软截断长度大于20bp。
9.一种审核基因变异检测结果的系统，包括装置A、装置B和装置C；
所述装置A能够按照权利要求1-8任一中的步骤(A)构建包括若干阳性变异位点和阴性非变异位点的测序数据的训练集；
所述装置B能够按照权利要求1-8任一中的步骤(B)从所述训练集中提取并向量化所述变异位点的特征；
所述装置C能够按照权利要求1-8任一中的步骤(C)用所述装置B提取并向量化的所述变异位点特征，采用随机森林法构建模型，然后利用模型判断待测位点变异是否为变异位点。
10.如下任一应用：
(I)权利要求1-8中任一所述的方法或权利要求9所述的系统在制备用于肿瘤早筛、肿瘤预后、肿瘤分类和/或肿瘤用药指导的产品中的应用；
(II)权利要求1-8中任一所述的方法或权利要求9所述的系统在肿瘤早筛、肿瘤预后、肿瘤分类和/或肿瘤用药指导中的应用；
(III)权利要求1-8中任一所述的方法或权利要求9所述的系统在基因变异检测中的应用。
CN202010135146.2A 2020-03-02 2020-03-02 一种审核高通量测序基因变异检测结果的方法 Pending CN111304308A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010135146.2A CN111304308A (zh) 2020-03-02 2020-03-02 一种审核高通量测序基因变异检测结果的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010135146.2A CN111304308A (zh) 2020-03-02 2020-03-02 一种审核高通量测序基因变异检测结果的方法

Publications (1)

Publication Number Publication Date
CN111304308A true CN111304308A (zh) 2020-06-19

Family

ID=71149431

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010135146.2A Pending CN111304308A (zh) 2020-03-02 2020-03-02 一种审核高通量测序基因变异检测结果的方法

Country Status (1)

Country Link
CN (1) CN111304308A (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111899790A (zh) * 2020-08-17 2020-11-06 天津诺禾医学检验所有限公司 测序数据的处理方法及装置
CN112634988A (zh) * 2021-01-07 2021-04-09 内江师范学院 基于Python语言的基因变异检测方法及系统
CN113470746A (zh) * 2021-06-21 2021-10-01 广州市金域转化医学研究院有限公司 降低高通量测序中人工引入错误突变的方法及应用
CN115171781A (zh) * 2022-07-13 2022-10-11 广州市金圻睿生物科技有限责任公司 肿瘤变异位点是否为噪音的识别方法、系统、装置和介质
WO2023207396A1 (zh) * 2022-04-25 2023-11-02 天津华大基因科技有限公司 用于分析变异检测结果的模型的构建方法
CN117577182A (zh) * 2024-01-15 2024-02-20 迈杰转化医学研究(苏州)有限公司 一种快速识别药物标识位点的系统及其应用

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108690871A (zh) * 2018-03-29 2018-10-23 深圳裕策生物科技有限公司 基于二代测序的插入缺失突变检测方法、装置和存储介质
CN109994155A (zh) * 2019-03-29 2019-07-09 北京市商汤科技开发有限公司 一种基因变异识别方法、装置和存储介质
CN110010195A (zh) * 2018-12-04 2019-07-12 志诺维思（北京）基因科技有限公司 一种探测单核苷酸突变的方法及装置
CN110846411A (zh) * 2019-11-21 2020-02-28 上海仁东医学检验所有限公司 一种基于二代测序的单独肿瘤样本区分基因突变类型的方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108690871A (zh) * 2018-03-29 2018-10-23 深圳裕策生物科技有限公司 基于二代测序的插入缺失突变检测方法、装置和存储介质
CN110010195A (zh) * 2018-12-04 2019-07-12 志诺维思（北京）基因科技有限公司 一种探测单核苷酸突变的方法及装置
CN109994155A (zh) * 2019-03-29 2019-07-09 北京市商汤科技开发有限公司 一种基因变异识别方法、装置和存储介质
CN110846411A (zh) * 2019-11-21 2020-02-28 上海仁东医学检验所有限公司 一种基于二代测序的单独肿瘤样本区分基因突变类型的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
BARNELL E.K.等: "Standard operating procedure for somatic variant refinement of sequencing data with paired tumor and normal samples", 《GENETICS IN MEDICINE》 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111899790A (zh) * 2020-08-17 2020-11-06 天津诺禾医学检验所有限公司 测序数据的处理方法及装置
CN112634988A (zh) * 2021-01-07 2021-04-09 内江师范学院 基于Python语言的基因变异检测方法及系统
CN112634988B (zh) * 2021-01-07 2021-10-08 内江师范学院 基于Python语言的基因变异检测方法及系统
CN113470746A (zh) * 2021-06-21 2021-10-01 广州市金域转化医学研究院有限公司 降低高通量测序中人工引入错误突变的方法及应用
CN113470746B (zh) * 2021-06-21 2023-11-21 广州市金域转化医学研究院有限公司 降低高通量测序中人工引入错误突变的方法及应用
WO2023207396A1 (zh) * 2022-04-25 2023-11-02 天津华大基因科技有限公司 用于分析变异检测结果的模型的构建方法
CN115171781A (zh) * 2022-07-13 2022-10-11 广州市金圻睿生物科技有限责任公司 肿瘤变异位点是否为噪音的识别方法、系统、装置和介质
CN115171781B (zh) * 2022-07-13 2023-04-07 广州市金圻睿生物科技有限责任公司 肿瘤变异位点是否为噪音的识别方法、系统、装置和介质
CN117577182A (zh) * 2024-01-15 2024-02-20 迈杰转化医学研究(苏州)有限公司 一种快速识别药物标识位点的系统及其应用
CN117577182B (zh) * 2024-01-15 2024-04-02 迈杰转化医学研究(苏州)有限公司 一种快速识别药物标识位点的系统及其应用

Similar Documents

Publication Publication Date Title
CN111304308A (zh) 一种审核高通量测序基因变异检测结果的方法
CN106909806A (zh) 定点检测变异的方法和装置
CN109880910A (zh) 一种肿瘤突变负荷的检测位点组合、检测方法、检测试剂盒及系统
WO2022141775A1 (zh) 基于dna甲基化谱的肿瘤免疫检查点抑制剂治疗有效性评估模型的构建方法
CN109767810B (zh) 高通量测序数据分析方法及装置
CN112802548A (zh) 单样本全基因组预测等位基因特异性拷贝数变异的方法
CN111968701B (zh) 检测指定基因组区域体细胞拷贝数变异的方法和装置
CN104794371B (zh) 检测逆转座子插入多态性的方法和装置
CN105986008A (zh) Cnv检测方法和装置
WO2022170909A1 (zh) 药物敏感预测方法、电子设备及计算机可读存储介质
WO2020237184A1 (en) Systems and methods for determining whether a subject has a cancer condition using transfer learning
CN105404793A (zh) 基于概率框架和重测序技术快速发现表型相关基因的方法
CN108304694B (zh) 基于二代测序数据分析基因突变的方法
CN108647495A (zh) 身份关系鉴定方法、装置、设备及存储介质
CN109461473B (zh) 胎儿游离dna浓度获取方法和装置
CN112837748A (zh) 一种区分不同解剖学起源肿瘤的系统及其方法
CN116356001B (zh) 一种基于血液循环肿瘤dna的双重背景噪声突变去除方法
KR102217272B1 (ko) 유전체 변이 정보를 이용한 질병 진단 바이오마커 추출 방법
CN109390034B (zh) 一种检测肿瘤组织中正常组织含量和肿瘤拷贝数的方法
CN116312779A (zh) 检测样本污染和识别样本错配的方法和装置
CN116287204A (zh) 检测特征基因的突变情况在制备静脉血栓栓塞症风险检测产品中的应用
CN114067908B (zh) 一种评估单样本同源重组缺陷的方法、装置和存储介质
US11535896B2 (en) Method for analysing cell-free nucleic acids
Dueker et al. Analysis of genetic linkage data for Mendelian traits
CN107208152B (zh) 检测突变簇的方法和装置

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