CN113470746A - 降低高通量测序中人工引入错误突变的方法及应用 - Google Patents

降低高通量测序中人工引入错误突变的方法及应用 Download PDF

Info

Publication number
CN113470746A
CN113470746A CN202110684865.4A CN202110684865A CN113470746A CN 113470746 A CN113470746 A CN 113470746A CN 202110684865 A CN202110684865 A CN 202110684865A CN 113470746 A CN113470746 A CN 113470746A
Authority
CN
China
Prior art keywords
mutation
score
site
throughput sequencing
mutations
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.)
Granted
Application number
CN202110684865.4A
Other languages
English (en)
Other versions
CN113470746B (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.)
Guangzhou Kingmed Diagnostics Group Co ltd
Guangzhou Jinyu Translational Medical Research Institute Co ltd
Guangzhou Kingmed Diagnostics Central Co Ltd
Original Assignee
Guangzhou Kingmed Diagnostics Group Co ltd
Guangzhou Jinyu Translational Medical Research Institute Co ltd
Guangzhou Kingmed Diagnostics Central 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 Guangzhou Kingmed Diagnostics Group Co ltd, Guangzhou Jinyu Translational Medical Research Institute Co ltd, Guangzhou Kingmed Diagnostics Central Co Ltd filed Critical Guangzhou Kingmed Diagnostics Group Co ltd
Priority to CN202110684865.4A priority Critical patent/CN113470746B/zh
Publication of CN113470746A publication Critical patent/CN113470746A/zh
Application granted granted Critical
Publication of CN113470746B publication Critical patent/CN113470746B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/50—Mutagenesis
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Bioinformatics & Cheminformatics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Genetics & Genomics (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • Chemical & Material Sciences (AREA)
  • Molecular Biology (AREA)
  • Proteomics, Peptides & Aminoacids (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)

Abstract

本发明涉及一种降低高通量测序中人工引入错误突变的方法及应用，属于生物信息学技术领域。该方法首先建立群体突变频率数据库，整合得到每个突变位点的群体突变频率，随后对各个窗口内的突变位点，按照预设的权重给定分数进行评分，得到每个突变位点的得分Svar后，考虑每个突变位点往后特定窗口大小内出现突变的次数，将这个区域内所有突变的分数相加，然后乘以这个窗口的突变数目，以达到放大的目的，从而突出存在人工引入错误突变的位点，将其识别后排除。从而可实现检出所有类型的人工引入错误突变，且具有效率高，自动化程度高，准确性高和检出全面的优势。

Description

降低高通量测序中人工引入错误突变的方法及应用
技术领域
本发明涉及生物信息学技术领域，特别是涉及一种降低高通量测序中人工引入错误突变的方法及应用。
背景技术
NGS技术为基因检测提供了便利，尤其是全外显子测序，其在价格上比全基因测序更有优势，由于性价比高，全外测序在医学上得到越来越广泛的应用。
全外测序的主要作用是寻找致病的突变，为疾病的诊疗提供重要的数据。全外筛选的突变需经过重重过滤进行质控，虽然现在对于测序的质量控制已经很严格，但是由于PCR扩增过程中的非特异性扩增等因素，部分人工引入的错误(artefacts)通过硬质量控制仍然不能排除，这对疾病的诊断产生了巨大的影响。
现阶段，虽然已有一些方法去除人工错误，但是基本上都需要以整合基因组浏览(IGV)方式，由具有丰富经验的实验人员进行人工进一步的判断，效率很低，且准确性依赖于人员的主观经验，难以标准化处理。
目前有研究推出一种利用LOFTEE排除人工引入错误突变的方法，其主要通过一组过滤器，将提前终止(stop-gained)、移码突变(frameshift)和可变剪切突变(splicing)进行评估，得到低可信和高可信的预测功能缺失突变，但是其主要在于评估和排除上述几种类型的人工错误，对于其他可能存在错误的位点并没有覆盖到，因此依然可能漏掉许多的人工错误位点。
发明内容
基于此，有必要针对上述问题，提供一种降低高通量测序中人工引入错误突变的方法，采用该方法，可通过机器自动化识别，批量排除人工错误位点，提高测序准确性。
1)建立群体突变频率数据库：获取若干样本高通量测序数据，整合得到每个突变位点的群体突变频率Fq，以及群体总数Pop；
2)位点评分：获取待测样本高通量测序数据，获得每个突变位点的信息，将各突变位点按照下述公式进行评分：
Svar＝(Fq×Pop×2)max＝M×(Score1+Score2+Score3)
其中：Svar为突变位点得分，Fq为群体突变频率，Pop为群体总数，M为预设最大值，Score1为突变类型分数，Score2为突变坐标分数，Score3为pLoF分数；
所述Score1的评分标准为：当突变为已知突变，不扣分；当突变为未知突变，扣分；
所述Score2的评分标准为：当突变为未知突变，则根据该突变所在refGene位置，按照预定规则扣分；
所述Score3的评分标准为：根据LOFTEE预测突变为功能缺失的可信度情况，按照预定规则扣分；
3)区域片段评分：从每个突变延后预定窗口大小划分，得到若干包括突变位点的序列片段；
按照下述公式对每个区域片段进行评分：
Figure BDA0003124225080000021
其中：Stotal为区域片段的总评分，N为该区域片段的总突变数；
4)位点排除：如区域片段的Stotal总评分低于预定阈值，则判定该区域片段中的突变位点为人工引入错误位点，排除该突变。
本发明人在长期实践工作的基础上，经过仔细分析、推论和求证后发现，对于人工引入的错误突变存在以下特征：插入缺失位点在捕获的测序区域两端大量聚集，导致了在某个区域内，不同的人出现不同插入缺失，导致种群中出现大量的INDEL聚集现象。以此前期研究基础出发，本发明人提出通过特定长度的窗口对上述特征进行权重评分并加权得到窗口区域片段的分数，再通过置信区间的筛选，批量排除人工错误位点的方法。
上述方法中，首先建立群体突变频率频率数据库，整合得到每个突变位点的突变频率，以及群体总数Pop，采用发明人设计的公式对突变位点进行评分，再按照预定窗口大小对序列进行划分，得到若干区域片段，对各区域片段中的突变位点得分进行整合，考虑每个突变位点往后特定窗口大小内出现突变的次数，将这个区域内所有突变的分数相加，然后乘以这个窗口的突变数目，以达到放大的目的，从而突出存在人工引入错误突变的位点，将其识别后排除。从而可实现检出所有类型的人工引入错误突变的目的，且具有效率高，自动化程度高，准确性高和检出全面的优势。
在其中一个实施例中，所述建立群体突变频率数据库步骤中，将每个突变位点按照染色体和基因组位置顺序进行排序。
在其中一个实施例中，所述位点评分步骤中，所述Score1的评分标准为：如该突变位点在dbSNP、ClinVar或gnomAD-exome数据库有注释，则判定该位点的突变为已知突变，不扣分；否则判定突变为未知突变，扣分；
所述Score2的评分标准为：当突变为未知突变，则根据该突变所在区域为refGene的exonic、splicing、UTR-3、UTR-5、upstream、downstream、intronic、intergenic或ncRNA区位置，按照预定规则扣分；
所述Score3的评分标准为：根据LOFTEE预测为突变为功能缺失(pLoF)高可信度位点或低可信位点，按照预定规则扣分。
可以理解的，上述dbSNP全称为the single nucleotide polymorphismdatabase，即单核苷酸多态性数据库，意思是“dna序列中的单一碱基对(base pair)变异”，也就是dna序列中a、t、c、g的改变，即基因组的一个特异和定位的位点出现两个或多个的核苷酸可能性，它是人类可遗传的变异中最常见的一种，因此，按照每个突变位点所在的区域及突变类型，将已报道为本领域所知悉的已知突变进行标注，有利于区分该突变是否为人工引入错误。上述dbSNP数据库优选dbSNP150。上述LOFTEE预测突变为功能缺失的pLoF突变为通过人工处理基因组聚合数据库(gnomAD)中单倍体疾病基因中的假定功能丧失(pLoF)突变，区分该突变类型，也有利于后续进行分析参考。
在其中一个实施例中，所述位点评分步骤中：
所述Score1的评分标准中，当突变为已知突变，则评分为0；当突变为未知SNV，则评分为-1；当突变为未知INDEL，则评分为-5；
所述Score2的评分标准中，当突变位于exonic区，则评分为0；当突变位于UTR-3、UTR-5、upstream、downstream、intronic、intergenic或ncRNA区，则评分为-1；当突变位于splicing区，则评分为-2；
所述Score3的评分标准中，当突变为高可信的pLoF，则评分为+3；当突变为低可信的pLoF，则评分为+2；
M＝100；
当突变位点的Svar分数大于0，则将此位点得分定义为0。
上述评分权重中，针对人工错误的主要特征(假阳性的插入缺失位点在捕获的测序区域两端大量聚集)，因此在特定窗口范围内，最大的权重应该是插入缺失(INDEL)，因此在窗口检测出未知INDEL的权重最大，其他影响因素还包括：SNP，是否为已知突变，突变的类型，是否为LOEFF注释的pLOF突变等。经过测试，我们将各因素的权重按照上述方式定义，能够较好识别出人工引入错误突变。
可以理解的，上述已知突变即为已有报道存在的突变。而对于群体频率，由于在实践中，很难得到某个位点突变的具体人数，因此一般用频率衡量这个位点的出现次数。当频率×人群×2大于100的时候，代表这个位点已经很常见，因此我们设定了一个最大的值，即M＝100。
在其中一个实施例中，所述refGene位置根据比对到NCBI refGene主转录本区域确定，所述pLoF仅保留出现在主转录本中Stop-gained、Splice site和Frameshift类型的pLoF突变位点，所述主转录本为NCBI refGene中最近更新的转录本，即距离当前时间最短的转录本。
在其中一个实施例中，所述区域片段评分步骤中，所述窗口大小通过以下方法得到：以覆盖数据库中95±5％相邻突变位点的区间长度为窗口大小。上述相邻突变位点中的相邻指两两间的区间长度，以能够覆盖绝大部分两两间突变位点区间的序列长度为窗口大小，可将人工引入错误突变正确识别。
在其中一个实施例中，所述数据库为dbSNP数据库、ECAC03全外显子数据库和/或genomeAD全外显子数据库，所述窗口大小为25±5bp。
本发明人在前期研究中，对dbSNP150，EXAC03全外显子数据库和genomeAD的全外显子数据库进行分析比较后发现，以25bp的窗口大小，可以覆盖超过95％的突变间隔。
在其中一个实施例中，所述位点排除步骤中，所述阈值为该待测样本中各区域片段按照Stotal得分由低至高排序的前5％所对应序列片段的Stotal得分。可以理解的，上述阈值也可以根据实践需要，设定某个确定的分数阈值。
本发明还公开了上述的降低高通量测序中人工引入错误突变的方法在全外显子高通量测序中的应用。
将上述的降低高通量测序中人工引入错误突变的方法应用于全外显子的测序中，能够降低人工引入突变错误，为疾病的诊疗提供准确数据，进一步提高了个性化精准医疗的可靠性。
本发明还公开了一种降低高通量测序中人工引入错误突变的装置，其特征在于，包括：
分析模块：用于按照上述的降低高通量测序中人工引入错误突变的方法，对待测样本全外显子高通量测序数据的突变数据进行分析。
可以理解的，上述装置可按照上述方法进行分析，具体呈现产品形式包括一体化设备，软件模块封装等。
与现有技术相比，本发明具有以下有益效果：
本发明的一种降低高通量测序中人工引入错误突变的方法，不仅仅针对低频及有害突变的人工错误进行识别排除，对于一般的人工引入错误突变位点也能覆盖到，具有全面、准确的优势。
并且，本发明的方法，可通过计算机自动化分析处理，克服了人工筛查存在的效率低、经验依赖性强的缺陷，适宜广泛的推广应用。
附图说明
图1为实施例中降低高通量测序中人工引入错误突变的方法流程示意图；
图2为实施例中不同窗口长度所对应Variants间隔示意图；
图3为实施例中chr9：35906583位点的IGV可视化结果；
图4为实施例中chr9：32986030位点的IGV可视化结果。
具体实施方式
为了便于理解本发明，下面将参照相关附图对本发明进行更全面的描述。附图中给出了本发明的较佳实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。本文所使用的术语“和/或”包括一个或多个相关的所列项目的任意的和所有的组合。
实施例
一种降低高通量测序中人工引入错误突变的方法，流程如图1所示，包括以下步骤：
1、建立群体突变频率数据库(S1)
获取若干样本高通量测序数据，整合得到每个突变位点的突变频率，具体如下：
1.1数据准备。
获取若干全外测序样本数据(本实施例中选择42868例样本)，得到其中所有位点的vcf合并文件，每个突变位点信息包括：位点所在染色体、具体位置、参考基因组碱基，突变的碱基类型等，并按照染色体和基因组位置进行排序。
1.2群体突变频率计算。
群体突变频率是计算的本发明突变位点评估分数的重要指标，需要通过计算得到。具体为：将上述突变位点信息数据输入为annovar的标准格式，即在vcf格式文件中加入频率数据，部分展示如下表。
表1.群体突变及频率输入文件格式列表
染色体 位置 参考碱基 突变碱基 等位基因频率 信息 纯合/杂合
chr1 14653 C T 2.33263e-05 GT 1/1
chr1 14907 A G 0.000128295 GT 1/1
chr1 14930 A G 0.000139958 GT 1/1
chr1 14932 G T 2.33263e-05 GT 1/1
chr1 14933 G A 1.16632e-05 GT 1/1
chr1 15903 G GC 6.9979e-05 GT 1/1
chr1 69331 C T 1.16632e-05 GT 1/1
chr1 69336 C T 1.16632e-05 GT 1/1
chr1 69460 C T 1.16632e-05 GT 1/1
chr1 69462 C G 2.33263e-05 GT 1/1
chr1 69474 T G 1.16632e-05 GT 1/1
chr1 69486 C A 2.33263e-05 GT 1/1
chr1 69486 C T 2.33263e-05 GT 1/1
chr1 69491 G A 1.16632e-05 GT 1/1
chr1 69495 C T 1.16632e-05 GT 1/1
chr1 69496 G A 1.16632e-05 GT 1/1
chr1 69496 G T 8.16422e-05 GT 1/1
chr1 69510 C T 1.16632e-05 GT 1/1
chr1 69511 A G 0.942186 GT 1/1
chr1 69513 A G 0.000606485 GT 1/1
chr1 69521 T C 8.16422e-05 GT 1/1
chr1 69522 T C 1.16632e-05 GT 1/1
chr1 69534 T C 0.00271752 GT 1/1
chr1 69543 C A 1.16632e-05 GT 1/1
chr1 69550 G A 9.33053e-05 GT 1/1
chr1 69550 G T 1.16632e-05 GT 1/1
chr1 69552 G A 1.16632e-05 GT 1/1
chr1 69555 T G 1.16632e-05 GT 1/1
chr1 69559 G A 0.000314906 GT 1/1
chr1 69563 A C 6.9979e-05 GT 1/1
其中：染色体和位置代表了该突变在基因组的位置定位，参考碱基为该位置上参考基因组hg19的碱基类型，突变碱基是本次测序检测出与参考基因组不同的碱基类型，等位基因频率即人群中该突变的频率，即该突变出现的次数除以(人群总人数乘以2)；信息即GT代表Genetype的缩写；1/1代表纯合、0/1代表杂合。
随后使用annovar软件自带的convert2annovar.pl-format vcf4功能模块，将vcf格式转换成annovar的标准格式，部分示例如下表。
表2.突变及频率数据注释前输入文件格式
染色体 起始位置 终止位置 参考碱基 突变碱基 纯合/杂合类型 频率
chr1 14653 14653 C T het 2.33263e-05
chr1 14907 14907 A G het 0.000128295
chr1 14930 14930 A G het 0.000139958
chr1 14932 14932 G T het 2.33263e-05
chr1 14933 14933 G A het 1.16632e-05
chr1 15903 15903 - C het 6.9979e-05
chr1 69331 69331 C T het 1.16632e-05
chr1 69336 69336 C T het 1.16632e-05
chr1 69460 69460 C T het 1.16632e-05
chr1 69462 69462 C G het 2.33263e-05
chr1 69474 69474 T G het 1.16632e-05
chr1 69486 69486 C A het 2.33263e-05
chr1 69486 69486 C T het 2.33263e-05
chr1 69491 69491 G A het 1.16632e-05
chr1 69495 69495 C T het 1.16632e-05
chr1 69496 69496 G A het 1.16632e-05
chr1 69496 69496 G T het 8.16422e-05
chr1 69510 69510 C T het 1.16632e-05
chr1 69511 69511 A G het 0.942186
chr1 69513 69513 A G het 0.000606485
chr1 69521 69521 T C het 8.16422e-05
chr1 69522 69522 T C het 1.16632e-05
chr1 69534 69534 T C het 0.00271752
chr1 69543 69543 C A het 1.16632e-05
chr1 69550 69550 G A het 9.33053e-05
chr1 69550 69550 G T het 1.16632e-05
chr1 69552 69552 G A het 1.16632e-05
chr1 69555 69555 T G het 1.16632e-05
chr1 69559 69559 G A het 0.000314906
chr1 69563 69563 A C het 6.9979e-05
2、位点评分(S2)
获取待测样本高通量测序数据，获得每个突变位点的信息，将各突变位点按照下述公式进行评分：
Svar＝(Fq×Pop×2)max＝M×(Score1+Score2+Score3)
其中：Svar为突变位点得分，Fq为群体频率，Pop为群体总数，M为预设最大值，Score1为突变类型分数，Score2为突变坐标分数，Score3为pLoF的分数。
具体步骤如下：
2.1Score1评分
对上述带有频率的人群突变数据进行annovar注释，注释依据源于dbSNP150，ClinVar数据库和gnomAD-exome数据库。
如该突变位点在dbSNP150、ClinVar或gnomAD-exome数据库有注释，则判定该位点的突变为已知突变，则评分为0；当突变为未知SNV(single nucleotide variants)，则评分为-1；当突变为未知INDEL(insertion-deletion)，则评分为-5。
2.2Score2评分
根据Score1评分结果，当突变为已知突变，Score2不评分，当突变为未知突变，则根据该突变所在区域位于refGene的不同区域评分。如：当突变位于exonic区，则评分为0；突变位于UTR-3、UTR-5、upstream、downstream、intronic、intergenic、或ncRNA区，则评分为-1；突变位于splicing区，则评分为-2。
在本实施例中，上述Score1和Score2评分中的annovar注释可在Ubuntu18.04.2LTS中采用如下具体命令行进行：able_annovar.pl all.vcf.avinput dir/humandb/-buildver hg19-out x--otherinfo-remove-protocol refGene,avsnp150,clinvar_20200316,gnomad_exome-operation g,f,f,f-nastring NA。
2.3Score3评分
标注单倍体疾病基因中的预测功能丧失pLoF突变，并对得到的pLoF突变数据进行筛选，根据该基因的主要转录本(即NCBI refGene中的基因最近更新的转录本)对pLoF位点进行annovar注释，包括pLoF位点位于基因组的位置，位于基因的具体转录本，所属变异类型，可能存在的氨基酸改变等信息。
Annovar部分注释结果如下表所示。
表3.annovar注释结果(refGene注释)
Figure BDA0003124225080000071
Figure BDA0003124225080000081
Figure BDA0003124225080000091
Figure BDA0003124225080000101
表4.annovar注释结果(dbSNP150注释结果)
染色体 起始位置 终止位置 参考碱基 突变碱基 dbSNP150注释结果
Chr Start End Ref Alt avsnp150
chr1 14653 14653 C T rs62635297
chr1 14907 14907 A G rs6682375
chr1 14930 14930 A G rs6682385
chr1 14932 14932 G T NA
chr1 14933 14933 G A rs199856693
chr1 15903 15903 - C rs557514207
chr1 69331 69331 C T NA
chr1 69336 69336 C T NA
chr1 69460 69460 C T NA
chr1 69462 69462 C G NA
chr1 69474 69474 T G rs752034042
chr1 69486 69486 C A NA
chr1 69486 69486 C T rs548369610
chr1 69491 69491 G A NA
chr1 69495 69495 C T NA
chr1 69496 69496 G A rs150690004
chr1 69496 69496 G T NA
chr1 69510 69510 C T NA
chr1 69511 69511 A G rs2691305
chr1 69513 69513 A G rs770590115
chr1 69521 69521 T C rs553724620
chr1 69522 69522 T C NA
chr1 69534 69534 T C rs190717287
chr1 69543 69543 C A NA
chr1 69550 69550 G A NA
chr1 69550 69550 G T NA
chr1 69552 69552 G A rs2531266
chr1 69555 69555 T G rs556374459
chr1 69559 69559 G A rs754025211
表5.annovar注释结果(ClinVar注释结果)
Figure BDA0003124225080000102
Figure BDA0003124225080000111
表6.annovar注释结果(gnomAD_exome注释结果)
Figure BDA0003124225080000112
Figure BDA0003124225080000121
对于pLoF位点，根据上述注释信息进行筛选，条件为：仅保留出现在主转录本中Stop-gained、Splice site和Frameshift类型的pLoF突变位点，其他不出现在主转录本而出现在非主转录本的pLoF位点忽略。
上述Stop-gained突变指使原来的氨基酸变成终止密码子，使蛋白转录提前终止的突变；Splice site指会影响转录本的正常转录，导致内含子区域也可能会转录进来，导致了整个蛋白质的结构发生变化的突变；Frameshift指移码突变，即剪接体可识别的RNA前体中内含子和外显子连接边界的序列和接头位点。假如该位置发生突变，导致mRNA的剪切出问题，蛋白质翻译就会出问题，剪切是转录的过程中将内含子剪切掉，保留外显子的过程。如果是经典剪切位点GT-AG突变则比较严重，属于致病性变异。
经过上述筛选，保留pLoF突变位点如下表所示：
表7.pLoF列表
Figure BDA0003124225080000131
Figure BDA0003124225080000141
根据最终得到pLoF位点进行Score3评分，当突变为高可信的pLoF，则评分为+3；当突变为低可信的pLoF，则评分为+2。
2.4位点评分
获得Score1-Score3评分后，将上述突变位点按照下述公式进行评分：
Svar＝(Fq×Pop×2)max＝M×(Score1+Score2+Score3)
其中：
Svar为突变位点得分；
Fq为群体频率；
Pop为群体总数，本实施例为42868；
M为预设最大值，本实施例中，M＝100。
通过这个公式，我们可以得到每一个位点的分数Svar，当Svar分数大于0，则判定为0。
3、区域片段评分(S3)
3.1窗口大小确定。
本发明人在前期研究中，对dbSNP150，EXAC03全外显子数据库和genomeAD的全外显子数据库进行分析，统计上述数据库中记载突变位点之间的距离，获得相邻两突变位点之间距离大小，发现以25bp的窗口大小，可以覆盖超过95％的突变间隔，如图2所示。因此确定切割窗口大小为25bp。
从每个突变延后25bp窗口大小进行切割，得到若干包括突变位点的区域片段。
3.2序列(区域片段)评分
按照下述公式对每个区域片段进行评分：
Figure BDA0003124225080000151
其中：Stotal为区域片段的总评分，N为该区域片段的总突变数。
4、位点排除(S4)
如区域片段的Stotal总评分低于预定阈值，则判定该区域片段中的突变位点为人工引入错误位点，排除该突变。
具体的，将所有区域片段的Stotal总评分进行比较，可通过一个具体数字(如-100或者-50)，也可以进行排序，取最低分数的前5％，这些窗口的突变则全为人口引入的错误突变。
本实施例中以-418作为筛选阈值，小于-418(即所有片段根据分数从小到大排序中的5％分数)进行筛选，得出的阳性结果(判定为人工引入错误突变)，其中部分得分负值较大及负值较小的人工引入错误突变如下表所示。
表8.人工引入错误突变及其所对应的窗口、窗口总突变数及总分数
Figure BDA0003124225080000152
Figure BDA0003124225080000161
注：每个位点均至少对应以自身位置为起点的窗口，以及可能对应其它窗口，如该位点落入任一判定为人工引入错误窗口，另一窗口未被认为属于人工引入错误，则认定为人工引入错误位点的结论不变。
上述结果看出，本实施例从8245702个突变中筛选掉155338个人工引入错误位点，筛选率为1.88％；保留了8090364个位点，占比98.12％，在正常的筛选范围(大于95％)内。
5、验证
抽取上述判定为存在人工引入错误突变的2个窗口片段进行验证。
5.1片段一
chr9：35906583位点对应的chr9:35906583-35906607窗口，其部分突变情况如下表所示。
表9.chr9：35906583位点对应的窗口(chr9:35906583-35906607)突变情况
Figure BDA0003124225080000162
Figure BDA0003124225080000171
Figure BDA0003124225080000181
Figure BDA0003124225080000191
Figure BDA0003124225080000201
Figure BDA0003124225080000211
该窗口有428个突变位点(上表示例性展示了部分位点35906583-35906589)，根据dbSNP、ClinVar和gnomAD-exome三大数据库判断，95％25bp长的窗口出现突变位点的数量不超过18个，该窗口远超过此数量，提示存在人工引入错误；且这些突变位点除了存在几个SNV外，其他均为INDEL，而INDEL的准确性不高，如此多的INDEL在同一区域内，其准确性更是大大降低。该位点所在窗口的IGV可视化结果见图3，图3为2个样本在该区域的reads比对情况，窗口长度为25bp，可以看出该区域的质量不佳，GC含量过高，可信程度不够。因此判断为人工引入的错误。与本发明的方法判断结果相吻合。
5.2片段二
再如chr9：32986030位点对应chr9:32986030-32986054窗口，其突变情况如下表和图4所示。
表10.chr9：32986030位点对应的窗口(chr9:32986030-32986054)突变情况
Figure BDA0003124225080000212
Figure BDA0003124225080000221
由上述结果可知，该窗口有45个突变位点，根据dbSNP、ClinVar和gnomAD-exome三大数据库判断，95％的25bp长窗口出现突变位点的数量不超过18个，该窗口超过此数量；且这些位点非常零散，相隔非常近，ploy A特征显著，非常容易出错；该位点iGV可视化结果见图4，图4为2个样本在该区域的reads比对情况，窗口长度为25bp，可以看出该区域的质量不佳，可信程度不够，因此判断为人工引入的错误。与本发明的方法判断结果相吻合。
验证结果证实，本发明的方法可除去绝大部分的人工引物错误突变，留下大量真实的突变。
以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种降低高通量测序中人工引入错误突变的方法，其特征在于，包括以下步骤：
1)建立群体突变频率数据库：获取若干样本高通量测序数据，整合得到每个突变位点的群体突变频率Fq，以及群体总数Pop；
2)位点评分：获取待测样本高通量测序数据，获得每个突变位点的信息，将各突变位点按照下述公式进行评分：
Svar＝(Fq×Pop×2)max＝M×(Score1+Score2+Score3)
其中：Svar为突变位点得分，Fq为群体突变频率，Pop为群体总数，M为预设最大值，Score1为突变类型分数，Score2为突变坐标分数，Score3为pLoF分数；
所述Score1的评分标准为：当突变为已知突变，不扣分；当突变为未知突变，扣分；
所述Score2的评分标准为：当突变为未知突变，则根据该突变所在refGene位置，按照预定规则扣分；
所述Score3的评分标准为：根据LOFTEE预测突变为功能缺失的可信度情况，按照预定规则扣分；
3)区域片段评分：从每个突变延后预定窗口大小划分，得到若干包括突变位点的区域片段；
按照下述公式对每个区域片段进行评分：
Figure FDA0003124225070000011
其中：Stotal为区域片段的总评分，N为该区域片段的总突变数；
4)位点排除：如区域片段的Stotal总评分低于预定阈值，则判定该区域片段中的突变位点为人工引入错误位点，排除该突变。
2.根据权利要求1所述的降低高通量测序中人工引入错误突变的方法，其特征在于，所述建立群体突变频率数据库步骤中，将每个突变位点按照染色体和基因组位置顺序进行排序。
3.根据权利要求1所述的降低高通量测序中人工引入错误突变的方法，其特征在于，所述位点评分步骤中，
所述Score1的评分标准为：如该突变位点在dbSNP、ClinVar或gnomAD-exome数据库有注释，则判定该位点的突变为已知突变，不扣分；否则判定突变为未知突变，扣分；
所述Score2的评分标准为：当突变为未知突变，则根据该突变所在区域为refGene的exonic、splicing、UTR-3、UTR-5、upstream、downstream、intronic、intergenic或ncRNA区位置，按照预定规则扣分；
所述Score3的评分标准为：根据LOFTEE预测为突变为功能缺失高可信度位点或低可信位点，按照预定规则扣分。
4.根据权利要求3所述的降低高通量测序中人工引入错误突变的方法，其特征在于，所述位点评分步骤中：
所述Score1的评分标准中，当突变为已知突变，则评分为0；当突变为未知SNV，则评分为-1；当突变为未知INDEL，则评分为-5；
所述Score2的评分标准中，当突变位于exonic区，则评分为0；当突变位于UTR-3、UTR-5、upstream、downstream、intronic、intergenic或ncRNA区，则评分为-1；当突变位于splicing区，则评分为-2；
所述Score3的评分标准中，当突变为高可信的pLoF，则评分为+3；当突变为低可信的pLoF，则评分为+2；
M＝100；
当突变位点的Svar分数大于0，则将此位点得分定义为0。
5.根据权利要求1所述的降低高通量测序中人工引入错误突变的方法，其特征在于，所述refGene位置根据比对到NCBI refGene主转录本区域确定，所述pLoF仅保留出现在主转录本中Stop-gained、Splice site和Frameshift类型的pLoF突变位点，所述主转录本为NCBI refGene中最近更新的转录本。
6.根据权利要求1所述的降低高通量测序中人工引入错误突变的方法，其特征在于，所述区域片段评分步骤中，所述窗口大小通过以下方法得到：以覆盖数据库中95±5％相邻突变位点的区间长度为窗口大小。
7.根据权利要求6所述的降低高通量测序中人工引入错误突变的方法，其特征在于，所述数据库为dbSNP数据库、ECAC03全外显子数据库和/或genomeAD全外显子数据库，所述窗口大小为25±5bp。
8.根据权利要求1所述的降低高通量测序中人工引入错误突变的方法，其特征在于，所述位点排除步骤中，所述阈值为该待测样本中各区域片段按照Stotal得分由低至高排序的前5％所对应区域片段的Stotal得分。
9.权利要求1-8任一项所述的降低高通量测序中人工引入错误突变的方法在全外显子高通量测序中的应用。
10.一种降低高通量测序中人工引入错误突变的装置，其特征在于，包括：
分析模块：用于按照权利要求1-8任一项所述的降低高通量测序中人工引入错误突变的方法，对待测样本全外显子高通量测序数据的突变数据进行分析。
CN202110684865.4A 2021-06-21 2021-06-21 降低高通量测序中人工引入错误突变的方法及应用 Active CN113470746B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110684865.4A CN113470746B (zh) 2021-06-21 2021-06-21 降低高通量测序中人工引入错误突变的方法及应用

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110684865.4A CN113470746B (zh) 2021-06-21 2021-06-21 降低高通量测序中人工引入错误突变的方法及应用

Publications (2)

Publication Number Publication Date
CN113470746A true CN113470746A (zh) 2021-10-01
CN113470746B CN113470746B (zh) 2023-11-21

Family

ID=77868784

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110684865.4A Active CN113470746B (zh) 2021-06-21 2021-06-21 降低高通量测序中人工引入错误突变的方法及应用

Country Status (1)

Country Link
CN (1) CN113470746B (zh)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014152990A1 (en) * 2013-03-14 2014-09-25 University Of Rochester System and method for detecting population variation from nucleic acid sequencing data
CN105861710A (zh) * 2016-05-20 2016-08-17 北京科迅生物技术有限公司 测序接头、其制备方法及其在超低频变异检测中的应用
CN107229841A (zh) * 2017-05-24 2017-10-03 重庆金域医学检验所有限公司 一种基因变异评估方法及系统
US20190017119A1 (en) * 2017-07-12 2019-01-17 The General Hospital Corporation Genetic Risk Predictor
WO2019169042A1 (en) * 2018-02-27 2019-09-06 Cornell University Ultra-sensitive detection of circulating tumor dna through genome-wide integration
CN110997936A (zh) * 2017-09-08 2020-04-10 深圳华大生命科学研究院 基于低深度基因组测序进行基因分型的方法、装置及其用途
CN111304308A (zh) * 2020-03-02 2020-06-19 北京泛生子基因科技有限公司 一种审核高通量测序基因变异检测结果的方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014152990A1 (en) * 2013-03-14 2014-09-25 University Of Rochester System and method for detecting population variation from nucleic acid sequencing data
CN105861710A (zh) * 2016-05-20 2016-08-17 北京科迅生物技术有限公司 测序接头、其制备方法及其在超低频变异检测中的应用
CN107229841A (zh) * 2017-05-24 2017-10-03 重庆金域医学检验所有限公司 一种基因变异评估方法及系统
US20190017119A1 (en) * 2017-07-12 2019-01-17 The General Hospital Corporation Genetic Risk Predictor
CN110997936A (zh) * 2017-09-08 2020-04-10 深圳华大生命科学研究院 基于低深度基因组测序进行基因分型的方法、装置及其用途
WO2019169042A1 (en) * 2018-02-27 2019-09-06 Cornell University Ultra-sensitive detection of circulating tumor dna through genome-wide integration
CN111304308A (zh) * 2020-03-02 2020-06-19 北京泛生子基因科技有限公司 一种审核高通量测序基因变异检测结果的方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
M.A. KAUNISTO, ET AL: "Novel splice site CACNA1A mutation causing episodic ataxia type 2", 《NEUROGENETICS 5》, pages 69 - 73 *
刘思捷，等: "CITED2基因在内脏反位患者中的突变分析", 《上海交通大学学报（医学版）》, pages 500 - 504 *

Also Published As

Publication number Publication date
CN113470746B (zh) 2023-11-21

Similar Documents

Publication Publication Date Title
CN109767810B (zh) 高通量测序数据分析方法及装置
CN108624650B (zh) 判断实体瘤是否适合免疫治疗的方法和检测试剂盒
US10496679B2 (en) Computer algorithm for automatic allele determination from fluorometer genotyping device
KR20210113237A (ko) 무 세포 dna 말단 특성
CN109545281B (zh) 一种基于二代高通量测序的trio家系遗传突变模式的分析方法
CN111139291A (zh) 一种单基因遗传性疾病高通量测序分析方法
CN110838341B (zh) 一种ATAC-seq测序数据的生物信息分析方法
CN105483210A (zh) 一种rna编辑位点的检测方法
CN111276189B (zh) 基于ngs的染色体平衡易位检测分析系统及应用
CN113470746B (zh) 降低高通量测序中人工引入错误突变的方法及应用
CN111662973A (zh) 与慢性阻塞性肺疾病易感性辅助诊断相关的snp位点及其应用
CN116741272A (zh) 基于基因组突变特征及基因集表达特征的卵巢癌hrd分型系统及方法
CN110373458A (zh) 一种地中海贫血检测的试剂盒及分析系统
WO2022266790A1 (zh) 降低高通量测序中人工引入错误突变的方法及应用
CN113462783B (zh) 基于MassArray核酸质谱的脑胶质瘤染色体lp/19q检测方法及应用
CN102154452B (zh) 一种鉴定顺式和反式调控作用的方法和系统
CN109033752A (zh) 一种基于长读长测序的多基因融合检测方法
TW202300656A (zh) 基因組序列上之拷貝數變異之候選斷點之機械性檢測
Walne et al. Genome-wide whole-blood transcriptome profiling across inherited bone marrow failure subtypes
CN115820654A (zh) Loxhd1基因突变体及其应用
CN113764044B (zh) 一种构建骨髓增生异常综合征进展基因预测模型的方法
US20220267836A1 (en) Method of haplotyping
CN117577177A (zh) 一种用于慢阻肺病多基因风险评分的snp标志物、装置及其应用
Quinones-Valdez et al. Long-read RNA-seq demarcates cis-and trans-directed alternative RNA splicing
CN115109847A (zh) 一组近视、高度近视相关的snp标志物及其应用

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