CN116895332A - 一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法 - Google Patents
一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法 Download PDFInfo
- Publication number
- CN116895332A CN116895332A CN202311164937.8A CN202311164937A CN116895332A CN 116895332 A CN116895332 A CN 116895332A CN 202311164937 A CN202311164937 A CN 202311164937A CN 116895332 A CN116895332 A CN 116895332A
- Authority
- CN
- China
- Prior art keywords
- mutation
- read
- reads
- site
- sequence
- 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
Links
- 230000035772 mutation Effects 0.000 title claims abstract description 181
- 238000000034 method Methods 0.000 title claims abstract description 68
- 239000012634 fragment Substances 0.000 title claims abstract description 60
- 238000001914 filtration Methods 0.000 title claims abstract description 37
- 238000001976 enzyme digestion Methods 0.000 title claims abstract description 19
- 238000010276 construction Methods 0.000 title claims abstract description 17
- 238000012163 sequencing technique Methods 0.000 claims abstract description 17
- 230000000295 complement effect Effects 0.000 claims abstract description 8
- 230000036438 mutation frequency Effects 0.000 claims abstract description 7
- 101100439666 Cupriavidus metallidurans (strain ATCC 43123 / DSM 2839 / NBRC 102507 / CH34) chrA1 gene Proteins 0.000 claims description 33
- 101100439667 Pseudomonas aeruginosa chrA gene Proteins 0.000 claims description 33
- 210000000349 chromosome Anatomy 0.000 claims description 24
- 239000011159 matrix material Substances 0.000 claims description 24
- 230000037438 passenger mutation Effects 0.000 claims description 12
- 238000011144 upstream manufacturing Methods 0.000 claims description 11
- 238000001514 detection method Methods 0.000 claims description 9
- 230000037430 deletion Effects 0.000 claims description 7
- 238000012217 deletion Methods 0.000 claims description 7
- 230000037431 insertion Effects 0.000 claims description 7
- 238000003780 insertion Methods 0.000 claims description 7
- 230000002255 enzymatic effect Effects 0.000 claims description 6
- 239000002773 nucleotide Substances 0.000 claims description 6
- 125000003729 nucleotide group Chemical group 0.000 claims description 6
- 230000000717 retained effect Effects 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000003776 cleavage reaction Methods 0.000 claims description 4
- 230000007017 scission Effects 0.000 claims description 4
- 238000007781 pre-processing Methods 0.000 claims 1
- 230000035945 sensitivity Effects 0.000 abstract description 5
- 230000007547 defect Effects 0.000 abstract description 2
- 108090000623 proteins and genes Proteins 0.000 abstract description 2
- 108020004414 DNA Proteins 0.000 description 11
- 238000005520 cutting process Methods 0.000 description 8
- 206010069754 Acquired gene mutation Diseases 0.000 description 5
- 108090000790 Enzymes Proteins 0.000 description 5
- 102000004190 Enzymes Human genes 0.000 description 5
- 230000037439 somatic mutation Effects 0.000 description 5
- 206010028980 Neoplasm Diseases 0.000 description 4
- 102000053602 DNA Human genes 0.000 description 3
- 238000007481 next generation sequencing Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011176 pooling Methods 0.000 description 2
- 108091008146 restriction endonucleases Proteins 0.000 description 2
- 238000012916 structural analysis Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 208000005623 Carcinogenesis Diseases 0.000 description 1
- 230000033616 DNA repair Effects 0.000 description 1
- 108091028043 Nucleic acid sequence Proteins 0.000 description 1
- 108020004682 Single-Stranded DNA Proteins 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 201000011510 cancer Diseases 0.000 description 1
- 230000036952 cancer formation Effects 0.000 description 1
- 231100000504 carcinogenesis Toxicity 0.000 description 1
- 230000004663 cell proliferation Effects 0.000 description 1
- 239000003153 chemical reaction reagent Substances 0.000 description 1
- 239000003795 chemical substances by application Substances 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000006862 enzymatic digestion Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 230000004853 protein function Effects 0.000 description 1
- 238000010008 shearing Methods 0.000 description 1
- 210000001082 somatic cell Anatomy 0.000 description 1
- 238000000527 sonication Methods 0.000 description 1
- 210000001519 tissue Anatomy 0.000 description 1
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
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
- G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
Landscapes
- Bioinformatics & Cheminformatics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- 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
本申请公开了一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，属于基因测序领域。该方法通过回溯比对数据寻找突变的支持reads，用支持reads的反向互补序列通过位置限定的动态规划局部比对算法与突变附近的基因组序列进行比对；分析比对结果和reads特征，将reads来源于人工片段的可能性分为6个等级，给不同等级的reads分配不同的权重计算得到支持reads中来源于人工片段的比例，通过为其设置阈值和对突变进行频率修正，过滤假阳性突变，最终精准的消除由于酶切法打断建库中产生的人工片段带来的高假阳率缺陷。该方法灵活、针对性强、敏感度高、结果稳定可靠，可以帮助酶切法打断的建库方式替代目前的超声法打断建库方式。
Description
技术领域
本申请属于基因测序技术领域，具体涉及一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法。
背景技术
新一代测序（next generation sequencing，NGS）技术可以同时对数以万计的DNA分子进行测序，并能够满足检测突变所需的高灵敏度、易用和准确的数据质量，常被作为分析癌症中体细胞突变的方法。单核苷酸变异（Single Nucleotide Variants，SNVs）是体细胞中最常见的单位点突变，常被作为确定蛋白质功能丧失和/或疾病风险的标志物，在人类各种类型的癌症的细胞增殖、肿瘤发生和肿瘤的精准治疗中起着非常重要的作用。
建库是NGS中保证数据质量的至关重要的一步，而DNA打断是建库的第一个步骤。目前，DNA打断最常用的两种方法是超声法打断和酶切法打断，其中：超声法打断的优势在于该方法将DNA暴露于均匀的超声爆发中会导致无偏剪切，从而产生高产量的均匀大小的双链DNA，文库偏好性更小，但是价格昂贵和费时；酶切法打断无需专门的仪器并且操作简单，所以更利于标准化和自动化，但是会有更高的文库偏好性，并且会产生人工片段。
人工片段产生的原理如图1所示，酶切双链DNA会在末端产生单链的DNA粘性末端，这一段粘性末端上如果存在互为反向互补的两段核苷酸序列，则这两段序列会互补结合，经过DNA修复和扩增，最终会产生和原DNA片段不同的DNA片段，这种片段就是酶切法打断过程中产生的人工片段。虽然这些人工片段比例不高，但是会产生假阳性低频突变，由于体细胞突变大多也是低频突变，所以人工片段对正确检出体细胞突变造成了显著影响，影响范围包括SNV、INDEL和complex SNV（本申请中统称为“突变”）。
人工片段出现在DNA打断阶段，无法通过唯一分子标记（Unique MolecularIdentifier，UMI）去除，只能通过生物信息方法过滤。当前的人工片段过滤方法主要是通过突变距离reads末端的距离和数据库进行匹配过滤突变，或者只分析携带softclip的reads结构，前者操作简单但是难以确定突变距离reads末端的距离的阈值且会对未记录在数据库中的新变异过于严格；后者从建库原理出发分析reads的结构是否来自于人工片段，但只针对携带softclip的reads，忽略了无softclip的reads，二者均无法确保灵敏度和精确度。
因此，迫切需要一种精准的方法识别人工片段导致的突变，从而进一步提高突变检测的准确性，使得更为简便的酶切法打断可以替代超声法打断。
发明内容
1. 要解决的问题
本申请针对酶切法打断建库中产生的人工片段会导致假阳性突变的出现，进而影响体细胞突变检测的准确性的问题，提供了一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，过滤掉人工片段产生的假阳性突变，从而克服目前酶切法打断建库中假阳性突变高的缺陷，提高体细胞突变检测的准确性。
2. 技术方案
为了解决上述问题，本申请所采用的技术方案如下：
本申请提供了一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，具体包括如下步骤：
S1，对待测样本进行测序，获得测序数据；对数据进行预处理，包括使用Trimmomatic v0.36软件去除测序序列上的接头序列及低质量的碱基片段；对序列进行比对，包括将预处理后的测序序列通过sentieon软件比对到hg19（GRCh37）人类参考基因组，经过排序和去重，生成sorted.rmdup.realign.bam文件；获取突变位点，包括用vardictv1.6软件检测突变，经过常规过滤之后得到最后的突变；
S2，支持突变位点突变的reads回溯，包括：解析S1中获得的sorted.rmdup.realign.bam文件和突变文件，获取覆盖每个突变位点的read上每个碱基与参考基因组位点一一对应关系，判断该read是否支持该突变，得到突变所有的支持reads；支持突变的判断标准如下：
（a）如果该突变是单核苷酸突变chrA:B:M→N，则当read在chrA:B的位置满足M→N，则计为chrA:B:M→N的支持read；其中A为染色体编号，B为染色体上的位点，M、N为单碱基；
（b）如果该突变是插入突变chrA:B:M→MN，则当read在chrA:B的位置对应N碱基长度的I，且read对应的碱基为N，则该read计为chrA:B:M→MN的支持read；其中A为染色体编号，B为染色体上的位点，M为单碱基，N为任意数量碱基；
（c）如果该突变是缺失突变chrA:B:MN→M，则当read在chrA:B的位置对应N长度的D，则该read计为chrA:B:MN→M的支持read；其中A为染色体编号，B为染色体上的位点，M为单碱基，N为任意长度碱基；
（d）如果该突变是复杂突变chrA:B:M→N，则当read在ref的M位点上游一个匹配碱基处正确匹配，read在ref的M位点的下游一个匹配碱基处正确匹配，且read在与M区间匹配处的read序列为N，则该read计为复杂突变chrA:B:M→N的支持read；其中A为染色体编号，B为染色体上的位点，M、N为任意数量碱基。
S3，位点限定的动态规划局部比对，对于步骤S2获取的所有支持reads，将其逐条与突变附近的参考基因组序列进行动态规划比对，包括如下步骤：
S31，取突变位点上下游各200~400 nt的序列；进一步地，取上下游各300 nt的序列；
S32，获取步骤S31中序列的反向互补序列；
S33，将S2中的reads和步骤S32的序列进行位点限定的动态规划局部比对，包括以S2中的read和S32的两个序列进行比对，构建S2中的read和S32的两个序列的矩阵，令i为序列1碱基位置的编号，xi为序列1第i个位点的碱基，j为序列2碱基位置的编号，yj为序列2的第j个碱基，F(i,j)为矩阵中i与j的分数；s(i,j)为xi与yi的匹配分数，在本申请中，每个位点的初始比对分设置为0；若xi=yi，s(i,j)=5，否则s(i,j)=-9；d为插入或缺失的得分，d=-14；利用如下公式计算比对结果得分，并填充打分矩阵；在填充后的打分矩阵中，选择最高分值进行回溯，一直到0分，获取回溯路径，从而得到最优比对结果；若匹配部分没有覆盖突变发生的位点，则将得分矩阵中的最大值替换为-1之后重复回溯，直到匹配区域覆盖突变发生的位点的碱基，得到最优比对结果；
；
S34，将支持read与步骤S32序列的比对关系定位到与步骤S31序列的比对关系，对于每条支持read，获取如下信息：
A点：read上的起始位点；
BE区域：read序列中与步骤S32序列的最佳比对区域；
C点：突变所在位点；
D点：read能够比对到参考基因组的末端位点；若read末端无softclip，则D点与E点重合；
LM：read在参考基因组序列上的比对区域；
nq：步骤S32序列上与BE区域的最佳比对区域；
r点：对应read上的D点；
N，R，Q分别为步骤S32序列上n，r，q对应的位点；
乘客突变：若BE中间存在除了C点之外的突变，且该突变能够与nq序列匹配，则为乘客突变。
S35，根据支持read来源于人工片段的可能性从低到高进行分级并统计该等级的reads数，等级及其reads数分别为：
Non-AR1：满足最佳匹配不在read两端的read数量，即E位点不在read末端的reads数量；
Non-AR2：满足最佳匹配在read两端的read数量，但匹配长度不大于9，即BE区域长度不大于9的reads数量；
Ambiguous-AR1：满足最佳匹配在read两端的read数量，但匹配长度小于18 nt，即BE区域长度小于18的reads数量；
Ambiguous-AR2：满足最佳匹配在read两端的read数量，匹配长度不小于18 nt，即BE区域长度大于或等于18的reads数量；
AR1：满足最佳匹配在read两端的read数量，且存在乘客突变的reads数量；
AR2：满足最佳匹配在read两端的read数量，且匹配区域包含read的softclip部分，即DE区域为softclip区域的reads数量；
S4，突变过滤，AR ratio为所有支持reads中来自人工片段的支持reads占比，计算公式如下：
，
其中，当AR1或者AR2不为0时，N=1，否则N=0；
为AR ratio设置阈值，不同区段的AR ratio处理方式如下：
（1）当AR ratio>=0.7时，过滤此突变；
（2）当0.7<AR ratio<=0.2时，对突变的突变频率进行修正，修正公式如下：
，
若freq>=0.02则保留突变，否则过滤掉此突变；
（3）当AR ratio<0.2时，保留突变。
3. 有益效果
本申请与现有技术相比，其有益效果在于：
（1）本申请提供的一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，基于支持reads的结构分析，现有的酶切导致的假阳突变过滤方法只针对softclipreads进行分析，在比对后分析read结构，去掉认为可能来自人工片段的reads，然而很多情况下来自人工片段的reads并不存在softclip的情况，这样就导致现有的方法难以有效地去掉所有假阳突变位点。另外现有方法在比对步骤去除reads，不能直观地看到reads去除对突变检出的影响。因此本申请针对每一个突变位点，都找到它们的支持reads，并对每条支持reads进行结构分析，判断它们是否来自于人工片段，然后配合以灵活的突变判别规则，实现对突变真假阳性的快速判断。
（2）本申请提供的一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，建立reads来源于人工片段的可能性分级方法。现有的酶切导致的假阳突变过滤方法只看reads末端与DNA负链的相似比例，对匹配长度和是否符合假阳突变来源限制较少，导致某些位于低复杂度基因组区域的突变的支持reads无法准确判断是否来源于人工序列。因此本申请对每条支持reads来源于人工片段的可能性做了6个分级，最低和最高的等级分别代表不可能来自于人工片段和非常可能来自人工片段，6个等级的reads分别分配的权重，然后计算支持reads中来自于人工片段的比例，通过给该比例设置阈值可以很好的区分真假阳突变。不同等级的可能性也可以帮助研究者对某个具体的突变了解更加充分，同时也提供了更高的灵活性，为研究者建立适当的模型提供了分类依据。
（3）本申请提供的一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，通过乘客突变对人工片段进行判断。对于酶切打断中单链互补结合所产生的人工片段，softclip是一个重要的特征，然而并非所有的人工片段都存在softclip特征，这也是现有方法的缺陷。乘客突变指的是人工序列除了携带目标突变，还同时携带的其它突变。通过判断reads上的单链互补结合区域是否包含乘客突变，给reads分配不同的人工序列可能性等级。
（4）本申请提供的一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，位点限定的动态规划比对算法。支持read的序列为A，突变上下游序列的反向互补序列为B，将序列A和序列B进行比对，若仅将匹配长度和相似比例作为标准，实际上并不符合人工序列产生假阳突变的模式。本申请中开发限定性的动态规划算法，实现序列A和序列B的匹配区域必定覆盖序列A上突变出现的位点，避免由于其它无关区域相似度过高，导致匹配结果影响对人工序列判断产生影响。最终可以得到对每一条支持reads的精准判断，进一步实现精准的假阳突变识别。
（5）本申请提供的一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，包括突变频率修正方法。本申请在处理突变同时存在真实片段和人工片段支持的情形时，不对来自人工片段的支持reads进行计数，只计算来自真实片段的支持 reads，从而对突变的突变频率进行修正，使其更符合真实情况。通过这种方式避免过滤一个真实位点，使得结果更加可靠。
附图说明
图1是酶切法建库中人工片段产生的原理图。
图2是支持read与参考序列的对照示意图。
图3是过滤前后突变统计结果。
具体实施方式
下面结合具体实施例对本申请进一步进行描述。
需要说明的是，本说明书中所引用的如“上”、“下”、“左”、“右”、“中间”等用语，亦仅为便于叙述的明了，而非用以限定可实施的范围，其相对关系的改变或调整，在无实质变更技术内容下，当亦视为本申请可实施的范畴。
除非另有定义，本文所使用的所有的技术和科学术语与属于本申请的技术领域的技术人员通常理解的含义相同；本文所使用的术语“和/或”包括一个或多个相关的所列项目的任意的和所有的组合。
实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市售购买获得的常规产品。
如本文所使用，术语“约”用于提供与给定术语、度量或值相关联的灵活性和不精确性。本领域技术人员可以容易地确定具体变量的灵活性程度。
如本文所使用，术语“......中的至少一个”旨在与“......中的一个或多个”同义。例如，“A、B和C中的至少一个”明确包括仅A、仅B、仅C以及它们各自的组合。
浓度、量和其他数值数据可以在本文中以范围格式呈现。应当理解，这样的范围格式仅是为了方便和简洁而使用，并且应当灵活地解释为不仅包括明确叙述为范围极限的数值，而且还包括涵盖在所述范围内的所有单独的数值或子范围，就如同每个数值和子范围都被明确叙述一样。例如，约1至约4.5的数值范围应当被解释为不仅包括明确叙述的1至约4.5的极限值，而且还包括单独的数字（诸如2、3、4）和子范围（诸如1至3、2至4等）。相同的原理适用于仅叙述一个数值的范围，诸如“小于约4.5”，应当将其解释为包括所有上述的值和范围。此外，无论所描述的范围或特征的广度如何，都应当适用这种解释。
实施例1
本实施例提供一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，具体包括如下步骤：
S1，对待测样本进行测序，获得测序数据；对数据进行预处理，包括使用Trimmomatic v0.36软件去除测序序列上的接头序列及低质量的碱基片段；对序列进行比对，包括将预处理后的测序序列通过sentieon软件比对到hg19（GRCh37）人类参考基因组，经过排序和去重，生成sorted.rmdup.realign.bam文件；获取突变位点，包括用vardictv1.6软件检测突变，经过常规过滤之后得到最后的突变；
S2，支持突变位点突变的reads回溯，包括：解析S1中获得的sorted.rmdup.realign.bam文件和突变文件，获取覆盖每个突变位点的read上每个碱基与参考基因组位点一一对应关系，判断该read是否支持该突变，得到突变所有的支持reads；支持突变的判断标准如下：
（a）如果该突变是单核苷酸突变chrA:B:M→N，则当read在chrA:B的位置满足M→N，则计为chrA:B:M→N的支持read；其中A为染色体编号，B为染色体上的位点，M、N为单碱基，如单核苷酸突变chr1:100:A→T，覆盖该突变位点的read上每个碱基与参考基因组位点一一对应关系如表1所示，该read在chr1:100的位置满足A→T，则该read计为chr1:100:A→T的支持 read；
（b）如果该突变是插入突变chrA:B:M→MN，则当read在chrA:B的位置对应N碱基长度的I，且read对应的碱基为N，则该read计为chrA:B:M→MN的支持read；其中A为染色体编号，B为染色体上的位点，M为单碱基，N为任意数量碱基；，如插入突变chr1:96:C→CT，覆盖该突变位点的read上每个碱基与参考基因组位点一一对应关系如表2所示，该read在chr1:96的位置对应1I，且read对应的碱基为T，则该read计为chr1:96:C→CT的支持read；
（c）如果该突变是缺失突变chrA:B:MN→M，则当read在chrA:B的位置对应N长度的D，则该read计为chrA:B:MN→M的支持read；其中A为染色体编号，B为染色体上的位点，M为单碱基，N为任意长度碱基，如缺失突变chr1:100:AC→A，覆盖该突变位点的read上每个碱基与参考基因组位点一一对应关系如表3所示，该read在chr1:100的位置对应1D，则该read计为chr1:100:AC→A的支持read；
（d）如果该突变是复杂突变chrA:B:M→N，则当read在ref的M位点上游一个匹配碱基处正确匹配，read在ref的M位点的下游一个匹配碱基处正确匹配，且read在与M区间匹配处的read序列为N，则该read计为复杂突变chrA:B:M→N的支持read；其中A为染色体编号，B为染色体上的位点，M、N为任意数量碱基，如复杂突变chr1：98 GGA→CT，覆盖该突变位点的read上每个碱基与参考基因组位点一一对应关系如表4所示，该read在ref的97 bp（突变的上游最近的一个匹配碱基）处正确匹配；该read在ref的101bp（突变的下游最近的一个匹配碱基）处正确匹配；且read在（3,6）区间的序列为CT，则该read计为复杂突变chrX:Y:MNP→OQ的支持read；
。
S3，位点限定的动态规划局部比对，位点限定的动态规划局部比对指的是在两条序列做比对时，二者能够匹配到的区域必须覆盖其中一条序列的指定位点（突变发生的位点）；对于步骤S2获取的所有支持reads，将其逐条与突变附近的参考基因组序列进行动态规划比对，包括如下步骤：
S31，取突变位点上下游各200~400 nt的序列；进一步地，取上下游各300 nt的序列；
S32，获取步骤S31序列的反向互补序列；
S33，将S2中的reads和步骤S32的序列进行位点限定的动态规划局部比对，包括以S2中的read和S32的两个序列进行比对，构建S2中的read和S32的两个序列的矩阵，令i为序列1碱基位置的编号，xi为序列1第i个位点的碱基，j为序列2碱基位置的编号，yj为序列2的第j个碱基，F(i,j)为矩阵中i与j的分数；s(i,j)为xi与yi的匹配分数，在本申请中，每个位点的初始比对分设置为0；若xi=yi，s(i,j)=5，否则s(i,j)=-9；d为插入或缺失的得分，d=-14；利用如下公式计算比对结果得分，并填充打分矩阵；在填充后的打分矩阵中，选择最高分值进行回溯，一直到0分，获取回溯路径，从而得到最优比对结果；若匹配部分没有覆盖突变发生的位点，则将得分矩阵中的最大值替换为-1之后重复回溯，直到匹配区域覆盖突变发生的位点的碱基，得到最优比对结果；
；
比如以序列1和序列2为例，序列1为AGACT，序列2为AGGCTAAAGACT，限定位点（指定位点，突变发生的位点）为序列2的第三个碱基G，比对结果的得分见表5，比对原理如下：
（a）初始打分矩阵
在位点限定的动态规划局部比对中，每个位点的初始比对分设置为0，初始打分矩阵见表6。
（b）打分公式
令i为序列1碱基位置的编号，xi为序列1第i个位点的碱基，j为序列2碱基位置的编号，yj为序列2的第j个碱基，F(i,j)为矩阵中i与j的分数；s(i,j)为xi与yi的匹配分数，在本申请中，若xi=yi，s(i,j)=5；否则s(i,j)=-9；d为插入或缺失的得分，在本发明中d=-14。
；
（c）填充打分矩阵
根据打分公式填充打分矩阵，填充后的打分矩阵见初始打分矩阵见表7。
（d）最大分值位点回溯
在填充后的打分矩阵中，选择最高分值进行回溯，一直到0分，获取回溯路径，从而得到最优比对结果，回溯路径和比对结果见表8。
（e）重复回溯
步骤d所得的比对结果中，匹配部分没有覆盖序列2的第三个碱基，并不符合预期。将得分矩阵中的最大值替换为-1之后重复回溯，直到匹配区域覆盖目标碱基。重复回溯结果见表9。
S34，将read与步骤S32序列的比对关系定位到与步骤S31序列的比对关系，示例见图2，对于每条支持read，获取如下信息：
A点：read上的起始位点；
BE区域：read序列中与步骤S32序列的最佳比对区域；
C点：突变所在位点；
D点：read能够比对到参考基因组的末端位点；若read末端无softclip，则D点与E点重合；
LM：read在参考基因组序列上的比对区域；
nq：步骤S32序列上与BE区域的最佳比对区域；
r点：对应read上的D点；
N，R，Q分别为步骤S32序列上n，r，q对应的位点；
乘客突变：若BE中间存在除了C点之外的突变，且该突变能够与nq序列匹配，则为乘客突变。
S35，根据支持read来源于人工片段的可能性从低到高进行分级并统计该等级的reads数，等级及其reads数分别为：
Non-AR1：满足最佳匹配不在read两端的read数量，即E位点不在read末端的reads数量；
Non-AR2：满足最佳匹配在read两端的read数量，但匹配长度不大于9，即BE区域长度不大于9的reads数量；
Ambiguous-AR1：满足最佳匹配在read两端的read数量，但匹配长度小于18 nt，即BE区域长度小于18的reads数量；
Ambiguous-AR2：满足最佳匹配在read两端的read数量，匹配长度不小于18 nt，即BE区域长度大于或等于18的reads数量；
AR1：满足最佳匹配在read两端的read数量，且存在乘客突变的reads数量；
AR2：满足最佳匹配在read两端的read数量，且匹配区域包含read的softclip部分，即DE区域为softclip区域的reads数量。
S4，突变过滤，AR ratio为所有支持reads中来自人工片段的支持reads占比，计算公式如下：
，
其中，当AR1或者AR2不为0时，N=1，否则N=0；
为AR ratio设置阈值，不同区段的AR ratio处理方式如下：
（1）当AR ratio>=0.7时，过滤此突变；
（2）当0.7<AR ratio<=0.2时，对突变的突变频率进行修正，修正公式如下：
，
若freq>=0.02则保留突变，否则过滤掉此突变；
（3）当AR ratio<0.2时，保留突变。
实施例2
本实施例提供一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法的验证，具体包括如下步骤：
Y1，临床数据的准备
相同的肿瘤组织样本分别使用超声打断和酶切打断建库测序，用变异检测流程分别分析获取不同打断方式的突变位点。然后通过以下步骤来获取本发明真阳位点集合与待测位点集合（超声数据检测的突变作为真阳位点；酶切数据检测的突变为待测位点，待测位点和真阳突变比较，确定待测位点的真假）。
Y11，在超声打断方法打断建库的样本集中提取变异检测流程检出的突变位点，同时按照频率>=2的条件进行过滤，最后得到了本发明的阳性突变位点集。
Y12，在酶切打断方法打断建库的样本集中提取变异检测流程检出的突变位点，同时按照频率>=2的条件进行过滤，最后得到了本发明的待测突变位点集。
本实施例使用的阳性样本共44个，共901个阳性位点；待测样本50个，共1503个待测位点。
Y2，利用实施例1中所述的一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法进行过滤，本实施例通过比较过滤前后的突变验证性能，使用敏感度（Sensitivity, Sn）、精密度（Precision, PPV）和准确度（Accuracy, Acc）3个评价指标对过滤前后的突变进行评估，计算公式如下：
其中，TP：1 Positive，将正例预测为正例；FN：0 Negative，将正例预测为负例；TN：1 Negative，将负例预测为负例FP：0 Positive，将负例预测为正例。
将过滤前后的突变位点数据集进行比较，统计结果见表10和图3，本申请在不损伤敏感度的情况下可以有效降低假阳率，提高精密度和准确度。
。/>
Claims (4)
1.一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，其特征在于，包括如下步骤：
S1，对待测样本进行测序，获得测序数据；对数据进行预处理；对序列进行比对，经过排序和去重，生成sorted.rmdup.realign.bam文件；获取突变位点；
S2，支持突变位点突变的reads回溯，包括：解析S1中获得的sorted.rmdup.realign.bam文件和突变文件，获取覆盖每个突变位点的read上每个碱基与参考基因组位点一一对应关系，判断该read是否支持该突变，得到突变所有的支持reads；
S3，位点限定的动态规划局部比对，对于步骤S2获取的所有支持reads，将其逐条与突变附近的参考基因组序列进行动态规划比对，包括如下步骤：
S31，取突变位点上下游各200~400 nt的序列；
S32，获取步骤S31中序列的反向互补序列；
S33，将S2中的reads和步骤S32的序列进行位点限定的动态规划局部比对，包括以S2中的read和S32的两个序列进行比对，构建S2中的read和S32的两个序列的矩阵，令i为序列1碱基位置的编号，xi为序列1第i个位点的碱基，j为序列2碱基位置的编号，yj为序列2的第j个碱基，F(i,j)为矩阵中i与j的分数；s(i,j)为xi与yi的匹配分数，每个位点的初始比对分设置为0；若xi=yi，s(i,j)=5，否则s(i,j)=-9；d=-14；利用如下公式计算比对结果得分，并填充打分矩阵；在填充后的打分矩阵中，选择最高分值进行回溯，一直到0分，获取回溯路径，从而得到最优比对结果；若匹配部分没有覆盖突变发生的位点，则将得分矩阵中的最大值替换为-1之后重复回溯，直到匹配区域覆盖突变发生的位点的碱基，得到最优比对结果；
；
S34，将支持read与步骤S32序列的比对关系定位到与步骤S31序列的比对关系，对于每条支持read，获取如下信息：
A点：read上的起始位点；
BE区域：read序列中与步骤S32序列的最佳比对区域；
C点：突变所在位点；
D点：read能够比对到参考基因组的末端位点；若read末端无softclip，则D点与E点重合；
LM：read在参考基因组序列上的比对区域；
nq：步骤S32序列上与BE区域的最佳比对区域；
r点：对应read上的D点；
N，R，Q分别为步骤S32序列上n，r，q对应的位点；
乘客突变：若BE中间存在除了C点之外的突变，且该突变能够与nq序列匹配，则为乘客突变；
S35，根据支持read来源于人工片段的可能性从低到高进行分级并统计该等级的reads数，等级及其reads数分别为：
Non-AR1：满足最佳匹配不在read两端的read数量，即E位点不在read末端的reads数量；
Non-AR2：满足最佳匹配在read两端的read数量，但匹配长度不大于9，即BE区域长度不大于9的reads数量；
Ambiguous-AR1：满足最佳匹配在read两端的read数量，但匹配长度小于18 nt，即BE区域长度小于18的reads数量；
Ambiguous-AR2：满足最佳匹配在read两端的read数量，匹配长度不小于18 nt，即BE区域长度大于或等于18的reads数量；
AR1：满足最佳匹配在read两端的read数量，且存在乘客突变的reads数量；
AR2：满足最佳匹配在read两端的read数量，且匹配区域包含read的softclip部分，即DE区域为softclip区域的reads数量；
S4，突变过滤，AR ratio为所有支持reads中来自人工片段的支持reads占比，计算公式如下：
，
其中，当AR1或者AR2不为0时，N=1，否则N=0；
为AR ratio设置阈值，不同区段的AR ratio处理方式如下：
（1）当AR ratio>=0.7时，过滤此突变；
（2）当0.7<AR ratio<=0.2时，对突变的突变频率进行修正，修正公式如下：
，
若freq>=0.02则保留突变，否则过滤掉此突变；
（3）当AR ratio<0.2时，保留突变。
2.根据权利要求1所述的一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，其特征在于，所述判断该read是否支持该突变的判断标准如下：
（a）如果该突变是单核苷酸突变chrA:B:M→N，则当read在chrA:B的位置满足M→N，则计为chrA:B:M→N的支持read；其中A为染色体编号，B为染色体上的位点，M、N为单碱基；
（b）如果该突变是插入突变chrA:B:M→MN，则当read在chrA:B的位置对应N碱基长度的I，且read对应的碱基为N，则该read计为chrA:B:M→MN的支持read；其中A为染色体编号，B为染色体上的位点，M为单碱基，N为任意数量碱基；
（c）如果该突变是缺失突变chrA:B:MN→M，则当read在chrA:B的位置对应N长度的D，则该read计为chrA:B:MN→M的支持read；其中A为染色体编号，B为染色体上的位点，M为单碱基，N为任意长度碱基；
（d）如果该突变是复杂突变chrA:B:M→N，则当read在ref的M位点上游一个匹配碱基处正确匹配，read在ref的M位点的下游一个匹配碱基处正确匹配，且read在与M区间匹配处的read序列为N，则该read计为复杂突变chrA:B:M→N的支持read；其中A为染色体编号，B为染色体上的位点，M、N为任意数量碱基。
3.根据权利要求2所述的一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，其特征在于，所述S31取突变位点上下游各300 nt的序列。
4.根据权利要求3所述的一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法，其特征在于，所述对数据进行预处理，包括使用Trimmomatic v0.36软件去除测序序列上的接头序列及低质量的碱基片段；对序列进行比对，包括将预处理后的测序序列通过sentieon软件比对到hg19（GRCh37）人类参考基因组，经过排序和去重，生成sorted.rmdup.realign.bam文件；获取突变位点，包括用vardict v1.6软件检测突变，经过常规过滤之后得到最后的突变。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311164937.8A CN116895332B (zh) | 2023-09-11 | 2023-09-11 | 一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311164937.8A CN116895332B (zh) | 2023-09-11 | 2023-09-11 | 一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116895332A true CN116895332A (zh) | 2023-10-17 |
| CN116895332B CN116895332B (zh) | 2023-12-05 |
Family
ID=88312488
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311164937.8A Active CN116895332B (zh) | 2023-09-11 | 2023-09-11 | 一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116895332B (zh) |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109658983A (zh) * | 2018-12-20 | 2019-04-19 | 深圳市海普洛斯生物科技有限公司 | 一种识别和消除核酸变异检测中假阳性的方法和装置 |
| CN110084314A (zh) * | 2019-05-06 | 2019-08-02 | 西安交通大学 | 一种针对靶向捕获基因测序数据的假阳性基因突变过滤方法 |
| CN110729025A (zh) * | 2019-12-17 | 2020-01-24 | 北京吉因加科技有限公司 | 基于二代测序的石蜡切片样本体细胞突变检测方法和装置 |
| CN111243664A (zh) * | 2020-03-26 | 2020-06-05 | 北京泛生子基因科技有限公司 | 一种基于高通量测序的基因变异检测方法 |
| CN112086131A (zh) * | 2020-08-18 | 2020-12-15 | 西安医学院 | 一种高通量测序中假阳性变异位点的筛选方法 |
| US20210043275A1 (en) * | 2018-02-27 | 2021-02-11 | Cornell University | Ultra-sensitive detection of circulating tumor dna through genome-wide integration |
| CN114141310A (zh) * | 2021-11-16 | 2022-03-04 | 广州燃石医学检验所有限公司 | 一种重复区域背景噪音过滤模型的构建方法及背景噪音过滤方法 |
| CN114613430A (zh) * | 2022-03-22 | 2022-06-10 | 苏州清港泉生物科技有限公司 | 一种假阳性核苷酸变异位点的过滤方法及计算设备 |
-
2023
- 2023-09-11 CN CN202311164937.8A patent/CN116895332B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20210043275A1 (en) * | 2018-02-27 | 2021-02-11 | Cornell University | Ultra-sensitive detection of circulating tumor dna through genome-wide integration |
| CN109658983A (zh) * | 2018-12-20 | 2019-04-19 | 深圳市海普洛斯生物科技有限公司 | 一种识别和消除核酸变异检测中假阳性的方法和装置 |
| CN110084314A (zh) * | 2019-05-06 | 2019-08-02 | 西安交通大学 | 一种针对靶向捕获基因测序数据的假阳性基因突变过滤方法 |
| CN110729025A (zh) * | 2019-12-17 | 2020-01-24 | 北京吉因加科技有限公司 | 基于二代测序的石蜡切片样本体细胞突变检测方法和装置 |
| CN111243664A (zh) * | 2020-03-26 | 2020-06-05 | 北京泛生子基因科技有限公司 | 一种基于高通量测序的基因变异检测方法 |
| CN112086131A (zh) * | 2020-08-18 | 2020-12-15 | 西安医学院 | 一种高通量测序中假阳性变异位点的筛选方法 |
| CN114141310A (zh) * | 2021-11-16 | 2022-03-04 | 广州燃石医学检验所有限公司 | 一种重复区域背景噪音过滤模型的构建方法及背景噪音过滤方法 |
| CN114613430A (zh) * | 2022-03-22 | 2022-06-10 | 苏州清港泉生物科技有限公司 | 一种假阳性核苷酸变异位点的过滤方法及计算设备 |
Non-Patent Citations (1)
| Title |
|---|
| 张自冠;施静艺;陈漪;孔杰;陈赛娟;黄金艳;: "基于高通量测序数据构建斑马鱼隐性遗传突变筛选系统", 现代生物医学进展, no. 09 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116895332B (zh) | 2023-12-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109767810B (zh) | 高通量测序数据分析方法及装置 | |
| CN106909806A (zh) | 定点检测变异的方法和装置 | |
| CN111755068B (zh) | 基于测序数据识别肿瘤纯度和绝对拷贝数的方法及装置 | |
| CN109801680B (zh) | 基于tcga数据库的肿瘤转移复发预测方法及系统 | |
| CN110060733B (zh) | 基于单样本的二代测序肿瘤体细胞变异检测装置 | |
| EP3729441B1 (en) | Microsatellite instability detection | |
| CN113517073B (zh) | 肺癌手术后生存率预测模型构建方法和预测模型系统 | |
| KR20180060759A (ko) | 염기서열의 변이 검출방법 및 이를 이용한 염기서열의 변이 검출 디바이스 | |
| CN116895332B (zh) | 一种酶切法打断建库中人工片段产生的假阳性突变的过滤方法 | |
| CN105528532B (zh) | 一种rna编辑位点的特征分析方法 | |
| CN105483210A (zh) | 一种rna编辑位点的检测方法 | |
| CN108460248B (zh) | 一种基于Bionano平台检测长串联重复序列的方法 | |
| CN111276189B (zh) | 基于ngs的染色体平衡易位检测分析系统及应用 | |
| CN110729025B (zh) | 基于二代测序的石蜡切片样本体细胞突变检测方法和装置 | |
| CN110164504B (zh) | 二代测序数据的处理方法、装置及电子设备 | |
| CN111763738A (zh) | 一种特征mRNA表达谱组合及肝癌早期预测方法 | |
| CN116434843A (zh) | 一种碱基测序质量评估方法 | |
| CN115954052A (zh) | 一种实体瘤微小残留病灶监控位点筛选方法及系统 | |
| WO2019132010A1 (ja) | 塩基配列における塩基種を推定する方法、装置及びプログラム | |
| CN115896256A (zh) | 基于二代测序技术的rna插入缺失突变的检测方法、装置、设备和存储介质 | |
| CN111798922B (zh) | 基于重测序数据中多态性位点密度鉴定小麦育种的基因组选择利用区间的方法 | |
| CN112382341B (zh) | 一种用于鉴定食管鳞癌预后相关的生物标志物的方法 | |
| CN111261226A (zh) | 基于ngs的微小残留病变自动化测序分析方法及装置 | |
| CN116646006B (zh) | 一种基于高通量测序和高斯混合模型的肿瘤相关基因体系突变检测方法及装置 | |
| CN116434830B (zh) | 基于ctDNA多位点甲基化的肿瘤病灶位置识别方法 |
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 |