二代测序下机数据的数据处理

本人接触过多家从事人类基因组的二代测序公司，包括肿瘤、全外遗传病和健康人全外测序体检，很多公司的数据处理和报告解读都存在一定的问题，这个系列就作为本人在人类基因组领域的解读和分享。

目前很多公司都在使用大体的流程如下：

1）原始的下机数据fastq文件；2）读长匹配参考基因组形成BAM文件；3）GATK和/或Samtools工具call出变异；4）Annovar工具对变异做出注释；5）变异进行过滤；6）过滤后的变异进行判读，医学解读等，最后出具报告。

这是一个行业内主要使用的流程，各个环节都容易出现问题，后续将罗列各种错误出现的案例。

这种流程主导思想是先进行过滤处理，剩余不多的变异后，再进行判读。称为“先过滤后判定”。这种流程最大的问题是过滤处理的不严谨，一旦将重要变异错误过滤后，后续很容易发生错报漏报，更为严重的是，这种错误在单一流程下自身很难察觉并纠正。

因此与之对应的另一种流程就是“先判定后过滤”，也就是对所有的变异先进行致病性判定，然后在致病变异中再依据测序质量保留要报的变异。

这两种流程各有优缺点，都会在不同的情况下发生漏报错误，“先判定后过滤”的自动化程度更高，更适合批量自动处理报告。

同一个下机数据，同时采用“先过滤后判定”和“先判定后过滤”，可以最大程度避免漏报错报，2种流程相互检验校准。

“先判定后过滤”流程开发中最大的难点在于如何对一个未知的、数据库中无记载的变异做出正确的判定（而不是简单的归为意义未明的变异），并给出完整的判定证据。本人开发有VarJudg程序专门用于胚系和体系变异的判定。VarJudg变异判定有极其复杂的判定规则，说明书在百度网盘上有下载。pan.baidu.com/s/14P6Q8E95656ARHjhRLxsQw 提取码：csj9

本人独立开发有二代测序数据分析和报告出具的ETW系统（整合有VarJudg变异判定），属于“先判定后过滤”，可免费试用。

后面分享几个“先过滤后判定”的生信流程中，见到过的错误过滤情况。注意，由于每家机构所采用的软件或版本有区别，以下所列的案例可能在某些机构中并不会发生。

错误过滤案例（1）

错误发生原因：同一氨基酸的1位和3位密码子同时发生突变，Annovar软件将其注释成2个变异。

案例详细：

以下是某个下机数据的原始变异列表中的2条，检出RAD50基因的相同突变频率的2个变异：

Chr Start End Ref Alt Gene AAChange

chr5 131927638 131927638 T A RAD50 RAD50:NM_005732:exon11:c.T1705A:p.Y569N

chr5 131927640 131927640 T A RAD50 RAD50:NM_005732:exon11:c.T1707A:p.Y569X

第一条变异，RAD50 p.Y569N是数据库中无记载的一个错义突变；

第二条变异，RAD50 p.Y569X是一个显然判定为致病的截断变异；

可以看到，对于同一个氨基酸的1位和3位碱基同时发生突变，Annovar软件的注释是把它们当成2个变异来分别注释；

在常规的“先过滤后判定”流程中，第一个变异被过滤掉，第二个变异保留了下来，那么最后的人工审核只看到第二个变异，很容易就报出一个RAD50 p.Y569X致病变异。

但实际上的情形是RAD50 c.T1705C p.Y569N和RAD50 c.T1707A p.Y569X位于同一条染色体上，其氨基酸569位的密码子位TAT，最终突变为AAA即p.Y569K，并非截断突变。最终导致错报。

案例推广：

在过滤变异，包括同义突变在内的各种单碱基变异，必须排除13密码子造成的假阳变异后，才能进行过滤处理。

解决方法：

1）过滤一个变异前，需要先判定是否为1、3类密码子同时突变，否定后再进行过滤；

2）对13密码子同时发生的变异，纠正软件的注释错误；

3）ETW报告系统，进行“先判定后过滤”，这种1、3类密码子同时突变，在VarJudg变异判定系统中定义为6类假阳变异，有假阳标签后，不再发生错报。

错误过滤案例（2）

错误发生原因：数据库中有直接记载为致病的变异，忽略测序质量，跳出过滤步骤，直接达到过滤后的变异列表。

案例详细：

以下是某个下机数据的原始变异列表中的2条，检出BRCA2基因的相同突变频率的2个变异：

Chr Start End Ref Alt Gene AAChange

chr13 32907421 32907421 A - BRCA2 BRCA2:NM_000059:exon10:c.1806delA:p.G602fs

chr13 32907421 32907421 - A BRCA2 BRCA2:NM_000059:exon10:c.1806dupA:p.G602fs

这个突变发生的位点是BRCA2基因的8个A碱基连续位置，测序错误或mapping错误常常发生，属于不可靠的变异。2个变异同时出现，很容易发现是假阳性变异。

但是第一个变异rs80359307，属于Clinvar数据库中有直接记载的变异，记载为致病。

而第二个变异尽管也判定为致病变异，但是数据库中并无直接记载。

在某些流程当中，数据库中有直接记载为致病的变异跳出过滤步骤，直接达到过滤后的变异列表。

这样就导致过滤后的变异列表，只看到第一个变异，而没有看到第二个变异，误报第一个致病变异。

案例推广：

数据库中有直接记载为致病的高质量可靠变异，并不能直接绕开过滤步骤，还是要做一定的假阳性判定后才能到达最后的变异列表。

解决方法：

1）流程中，尽量不要有绕开过滤步骤的变异，有绕开过滤步骤的变异在报之前需谨慎复核。

2）ETW报告系统，进行“先判定后过滤”（合计有10类假阳性变异），所有的变异在过滤前均需要做假阳判定打上假阳标签，不再发生假阳性变异错报。

错误注释过滤案例（3）

错误发生原因：Annovar软件对邻近多个碱基同时发生突变且跨越多个氨基酸的变异，不能给出正确的注释。

案例详细1：

以下是某个下机数据的原始变异列表中的1条，检出KRAS基因的2个碱基同时发生变异：

Chr Start End Ref Alt Gene AAChange

chr12 25398285 25398286 CA AG KRAS KRAS:NM_004985:exon2:c.33_34CT

该变异实际上是KRAS基因的11位氨基酸密码子的最后1个碱基和12位氨基酸密码子的第一个碱基同时发生突变，最终导致11位氨基酸不发生改变，而12位氨基酸p.G12C。Annovar软件没有给出正确的注释，导致这个致病变异被漏报。

案例详细2：

以下是某个下机数据的原始变异列表中的1条，检出ERBB2基因的1个碱基发生复杂插入突变：

Chr Start End Ref Alt Gene AAChange

chr17 37880997 37880997 G TTGT ERBB2 ERBB2:NM_004448:exon20:c.2326_2326delinsTTGT

同样的，这个Annovar软件没有给出这个复杂变异的正确注释，实质上手工拼接后发现这个变异是p.G776delinsLC（属于20号外显子插入突变），导致这个致病变异被漏报。

案例推广：

对邻近多个碱基同时发生突变且跨越多个氨基酸的变异，Annovar软件没有给出正确的氨基酸注释的变异，这类变异不能直接过滤。

解决方法：

1）修复Annovar软件bug，对于不能注释的变异进行特殊处理不能直接过滤。

2）ETW报告系统，对于Annovar不能正确注释的变异判定10类假阳性变异，常见的10类假阳性变异已经能够正确注释，一旦检出新发未注释的10类假阳变异，系统报警提醒人工拼接序列。

错误注释过滤案例（4）

错误发生原因：GATK或Samtools软件对多个碱基同时发生插入缺失的复杂突变，不能正确的mapping。

案例详细1：

以下是某个下机数据的原始变异列表中的2条，检出EGFR基因的多碱基发生复杂缺失突变：

Chr Start End Ref Alt Gene AAChange

chr7 55242467 55242483 AATTAAGAGAAGCAACA - EGFR EGFR:NM_005228:exon19:c.2237_2253del:p.E746fs

chr7 55242487 55242487 C - EGFR EGFR:NM_005228:exon19:c.2257delC:p.P753fs

这2个移码变异经手工重新组装以后，实际为EGFR p.E746_S751delinsVS，软件并不能给出变异的真实结果。

案例详细2：

以下是某个下机数据的原始变异列表中的6条，检出EGFR基因的多碱基发生复杂缺失突变：

Chr Start End Ref Alt Gene AAChange

chr7 55242467 55242468 AA - EGFR EGFR:NM_005228:exon19:c.2237_2238del:p.E746fs

chr7 55242470 55242485 TAAGAGAAGCAACATC - EGFR EGFR:NM_005228:exon19:c.2240_2255del:p.L747fs

chr7 55242471 55242488 AAGAGAAGCAACATCTCC CA EGFR EGFR:NM_005228:exon19:c.2241_2258CA

chr7 55242471 55242478 AAGAGAAG - EGFR EGFR:NM_005228:exon19:c.2241_2248del:p.L747fs

chr7 55242472 55242487 AGAGAAGCAACATCTC - EGFR EGFR:NM_005228:exon19:c.2242_2257del:p.R748fs

chr7 55242481 55242488 ACATCTCC - EGFR EGFR:NM_005228:exon19:c.2251_2258del:p.T751fs

这6个移码变异经手工重新组装以后，实际仅为1条变异EGFR p.E746_S752delinsV（chr7_55242467_55242485_AATTAAGAGAAGCAACATC_T； c.2237_2255>T），需要重新手工从原始读长去拼接序列。

案例推广：

多个碱基复杂的插入缺失变异，GATK或Samtools软件都不能很好的正确组装，经常call出一堆错误的移码变异。不能把软件call出的变异列表当成真正的变异列表。

解决方法：

1）修复GATK或Samtools软件软件bug，对于错误,mapping的变异进行修复。同时引入其它的软件，例如mutscan等。

2）ETW报告系统，对于指定基因和指定区域的已知的复杂插入缺失突变，进行了组装修复，并设定有报警机制，一旦出现新出现的未知的复杂移码突变，人工干预手工组装。

3）对于肿瘤体细胞检测，最明显的是EGFR和ERBB2这2个基因的19和20号外显子的插入缺失突变，此区域进行扫描报警，即该区域出现未知的移码或插入缺失突变，需要人工复核变异；同样的对于遗传病检测，根据遗传病的种类，也需要对检测范围内的基因特定区域设置报警。

4）多碱基的复杂变异，几乎在每份下机数据中都有出现的概率，生信需要定期统计各种碱基变异出现的频谱并与参考值做比对，从统计学上防止软件漏报的情形。某些机构的变异call出频率明显不正常，下机数据中几乎没有复杂变异，过度依赖相信国外软件。

案例（5）

错误发生原因：多个氨基酸同时发生突变，可以有多种注释方式，需要匹配文献报道或数据库记载的注释方式。

案例详细：

以下是某个下机数据的原始变异列表中的2条，检出EGFR基因的2个相邻的错义突变：

Chr Start End Ref Alt Gene AAChange

chr7 55249020 55249020 A T EGFR EGFR:NM_005228:exon20:c.A2318T:p.H773L

chr7 55249022 55249022 G A EGFR EGFR:NM_005228:exon20:c.G2320A:p.V774M

在同一条染色体上的2个相邻的错义突变p.H773L和p.V774M，还有另一种写法即插入缺失突变，即p.H773_V774delinsLM，核查文献和数据库中的记载及解析后，发现文献中对该变异的描述均为p.H773_V774delinsLM，因此该变异用药解析调用的是p.H773_V774delinsLM。

案例推广：

同一条染色体上多个相邻氨基酸同时发生突变，调用变异的解析数据库需要2种不同的方式：单个错义突变逐个匹配和把这些变异当成一个插入缺失去匹配。

案例（6）

错误发生原因：很多机构都把数据库中无记载同义突变和人群频率高的突变过滤掉。

同义突变可能致病的原因

1）同义突变正好发生在GT-AG剪切位点上面，如果预测为影响剪切，那么这个未知的同义突变则视为剪切失活突变，大概率会判为致病。

2）该基因有多个功能转录本，在其中一个转录本中为同义突变，但是在另一个转录本中为错义或stopgain突变，但某些非常特殊的情况下Annovar软件注释只写了一种同义突变。

在Clinvar数据库中记载有很多的已知的致病的同义突变，以上2个致病因素可以解释大约70%的同义突变致病原因。

案例推广：

“先过滤后判定”的流程中，对每一步过滤都需要谨慎是否有错误过滤，试跑1份10份100份报告均没有发生错误过滤并不代表流程的正确性。如何在流程中发现新的过滤错误原因并纠正出来是最大的难点。

“先判定后过滤”流程中，对每一个未知的变异（包含同义突变和基因间区域变异）都进行完整的判定，确认是非致病后再进行过滤，从另一个角度来进行数据处理。

“先过滤后判定”和“先判定后过滤”相互验证可以最大程度保障没有漏报错报。

最后说几点心得：

1，没有完美的算法和软件，过度依赖国外软件的结果就是一大堆错报漏报，必须有持续纠错的思想准备并具有发现问题能力。

2，人类基因组的解读极其复杂，问题层出不穷。由医学人员制定规则，软件人员根据规则来开发算法的思路是有问题的，仅变异判定这一块的逻辑节点就超过百万，这么多规则很难完整传递。

3，对错报漏报的零容忍，ETW系统需要每天更新补丁，一个版本的系统在使用中很难撑过3个月，就必须推倒重新开发。有些机构开发出一套系统可以用好几年，纠错能力非常有限，或者是面对各种漏洞采用人工方法去填补。

4，基因检测行业中两类不同的机构：有些只用一套报告系统的，很难发现自身错误，经比对发现有不少漏报错报，但这类机构的人员普遍都非常的自信基本不接受外界的纠错。另一类机构采用两套报告系统的，相互印证，经常能发现系统的各种错误及漏洞，对报告系统反复进行大版本的升级，这类机构的人员反而非常谨慎。

5，评判一套数据分析和报告系统的好坏的标准：1）能直接生成最终的报告无需任何人工干预的报告占比多少；2）系统生成一份报警的报告后，后续人工的工作量究竟多大；3）好的系统可以让单人每天出具50-100份报告；4）从下机数据到报告出具，包括数据库的搭建，整条生产线人员可以在5人以下。

二代测序下机数据的数据处理相关推荐

关于文献中二代测序数据下载（NCBI）的问题
关于文献中二代测序数据下载(NCBI)的问题现在二代测序用于生物学研究非常广泛,大部分文章的序列会上传到Sequence Read Archive(SRA)上,这东西也属于NCBI数据库中的吧,我理 ...
PacBio HiFi测序介绍及百迈客最新下机数据公布
PacBio HiFi测序介绍及百迈客最新下机数据公布百迈客生物已认证账号已关注 3 人赞同了该文章众所周知,要获得基因组的完整图片,就必须组装reads,以目前主要的测序技术来看,短读长 ...
宏基因组数据处理 - Nanopore下机数据fast5格式
过年期间,我的三代Nanopore测序数据回来了.本来期待的是几十G的数据吧,结果人家寄来的硬盘上来就是两三T,人直接傻了.经过整理,发现测序公司送来的数据分有两种类型,一种就是我们熟悉的FASTQ格 ...
木桶排序算法_【生信常识】二代测序的比对算法浅析
前言本来我只打算将孟大哥的视频内容做一个文字版的概述,然后孟大哥说,不如再加一个算法推导吧,然后我就开始看多一些东西,然后就想着把孟大哥视频里面大概提及然后没有仔细讲的部分做一些补充,完善整个体系的 ...
二代测序技术之illumina测序技术原理简介
现今的生信领域几乎就是和无数的序列打交道,而这些序列的来源就是如今风靡的高通量测序技术,现今的测序不论是测RNA.DNA.miRNA还是ChIP-Seq等等,都是基于NGS(二代测序,next-gen ...
MPB：微生物所蔡磊组-基于二代测序的真菌基因组组装和注释
为进一步提高<微生物组实验手册>稿件质量,本项目新增大众评审环节.文章在通过同行评审后,采用公众号推送方式分享全文,任何人均可在线提交修改意见.公众号格式显示略有问题,建议电脑端点击文末阅 ...
二代测序linux软件,二代测序数据分析软件包大全
二代测序数据分析软件包大全 Integrated solutions*CLCbio Genomics Workbench-de novoand reference assembly of Sanger ...
二代测序技术中生物信息学的应用
随着科学技术的巨大进步,产生了大量的"组学"数据.理解生物系统各个层次产生的大量序列和结构数据是关键,由此产生了"生物信息学". "生物信息学&quo ...
和dump文件_SRA数据库及下载二代测序原始数据转换为fastq文件
一.SRA数据库: NCBI网站储存二代测序原始数据的数据库. (一)SRA数据类型: 1.Studies:研究课题 2.Experiments:实验设计 3.Samples:样品信息 4.Runs: ...

二代测序下机数据的数据处理

二代测序下机数据的数据处理相关推荐

最新文章

热门文章