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

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

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

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人以下。