赛福基因公开课第五节《全外显子组数据分析简介》

各位老师好，很高兴今天有机会和大家一起来探讨全外显子数据分析中的相关知识和技术。

顾名思义，全外显子测序是对基因组中所有外显子进行测序分析的方法。该方法首先利用DNA或RNA探针，将全基因组中外显子区域DNA进行捕获，然后对捕获的DNA进行PCR扩增，最后对扩增后的产物进行高通量测序。虽然外显子只占基因组的2%，却含有85%的已知致病变异， 但测序和分析价格却远比全基因组测序要低得多。所以目前很多科研和临床选用全外显子测序，而不是全基因组测序。

和WGS相比，WES有以下优势：

1.测序和数据存储的花费低。尽管WES一般比WGS测序深度高（100X vs 30X), WES的测序序列主要在占基因组2%的外显子区域，使得总的测序文件比WGS的测序文件小得多。平均测序深度为100X的全外显子测序文件的大小是6GB左右，而30X的全基因组测序文件的大小一般是90GB左右;

2.由于数据文件小，数据分析的速度就快，数据分析的花费也低； 

3.由于测序深度高，WES更容易检测出罕见变异；

4.跟非编码区相比，蛋白编码区的研究更透测，注释度更高，所以WES数据分析中的变异能够更好地被解读。 

尽管和WGS相比，WES存在很多优势，但也有不足之处：

1.WES的测序深度不均一，一些区域测序深度过高造成浪费，也有一些区域测序深度过低而无法检测到存在的变异。比如一些SNP密集的区域，富集过程中RNA或DNA探针无法和该区域杂交导致无法有效捕获，该区域的测序深度就会过低；

2.由于WES只测外显子区域（人类的外显子区域长度一般小于200bp），因而无法有效检测CNV，结构重组和其他大的结构变异；

3.由于没有捕获探针，新的编码基因无法捕获。

在高通量测序和分析中有三个常用的质量分数，了解这些质量分数将有助于对外显子数据分析的理解。 第一个是碱基质量分数，它是测序过程中碱基被测序仪错误识别的概率。 另一个是比对质量分数， 它是一个序列在参考基因组上比对正确与否的可信度。还有一个是变异质量分数，它是我们检测到的变异是否是生物变异的可信度。变异质量分数是变异位点所有序列比对质量分数的均方根。

由美国哈佛-麻省理工博德研究所 (Broad institute of MIT and Harvard) 开发的基因组分析工具包GATK（ genome analysis toolkit)是目前常用的NGS数据分析软件。今天我们以GATK的best practices 为列来介绍全外显子数据分析的流程。WES 数据分析大致分为三步：对WES数据进行前处理；利用前处理后的数据进行变异检测，最后对检测到的变异进行筛选和注释。

第一部分：WES数据的前处理流程。

WES数据分析的前处理包括对原始fastq形式的序列进行质量控制，将清理后的fastq形式的序列比对到参考基因组，对比对到参考基因组的序列进行多步质量控制，最终得到可以用来检测变异的序列。

一般情况下，我们从测序公司拿到的WES数据，是原始数据状态（raw data)。 拿到数据后的第一步一定要查看fastq文件是否完整，是否存在错误，以及数据质量的好坏。常用的软件是FASTQC. 根据FASTQC的结果报告，如果fastq文件有错误，就无法进行下面的分析。如过fastq文件没有错误，就对原始数据进行清理，包括去除低质量或者过短序列，剪切末端低质量的碱基，去除接头污染以及可能的外源污染序列。 常用的软件包括Trimmamatic, cutadapt 和fastq_mcf等。清理之后得到clean 序列。对于NGS数据质量控制的知识和技术，请参考赛福基因微课程的第四讲。

得到clean 序列后，利用比对软件，将clean 序列比对到参考基因组上去。比对简单的说就是对每一个序列， 找到它在参考基因组上对应的的位置。

理论上，将序列比对到参考基因组上非常简单，就是找到每个序列在参考基因组上的位置并记录下来。但事实上，很多情况下，测序数据和参考基因组间存在着各种差异，包括SNV，INDEL，CNV， SV等，这使得比对变得非常复杂。比对软件需要运用复杂的算法，将序列比对到可信度最高的一个或者少数几个位置上去。这个可信度就是我们前面提到的比对质量分数（mapping quality ）。另外，如这里的示意图，由于参考基因组中存在重复序列等原因，有些序列可以同等地比对到不同的位置，这也增加比对的难度。

将序列比对到参考基因组常用的软件包括BWA(Burrows-Wheeler Aligner), Bowtie2, Blat, SOAP 等。 比对之后的结果文件以SAM 的格式存储。 如果要对比对结果进行其他一些操作， 一般需要将SAM文件转化成压缩的二进制格式的BAM文件。

如果我们想直观地查看比对结果，可以使用可视化软件。常见的可视化软件有IGV，tablet等。用IGV或tablet可视化比对结果时，就必须要用加过索引的BAM文件。这里显示的是一个WES数据在ARX基因区的比对结果。我们可以看到，外显子区域以外也有比对的序列。 一般情况下，这些序列不要去除，更大范围内更多的比对序列有利于下游更精确的变异检测。在IGV里，我们可以对感兴趣的区域进行放大。

这个是我们对上图中中间的那个外显子区进行放大的结果。在IGV里default的设置下，比对到参考基因组上的序列有不同的颜色。灰色表示该序列在基因组上的比对结果是单一的，而且和序列对中另一个序列间的距离在正常范围内。白色说明该序列比对到2个或者2个以上的位置。如果reads是其他颜色，说明与该序列配对的序列比对到了其它的染色体， 或者两个序列间的距离大于或小于正常的插入序列长度。

另外，可以用samtools或者bamtools查看比对结果的统计值。 samtools的flagstat软件，可以查看总的比对率， 成功比对的序列里多少个是序列对，多少是单个的序列等。Samtools 的idxstat软件可以查看每个染色体的长度，该染色体上比对的序列数以及没有比对的序列数。

原始比对结果里含有一些比对质量较低或者错误比对的序列，为了提高变异识别的敏感性和精确性，要对原始比对结果进行质量控制。 这些质量控制包括：去除重复序列，对可疑的比对区域进行重新比对（Local Realignment）和重新校正碱基质量分数。 对于这些质量控制的相关知识和技术，请参考赛福基因微课堂第四讲。比对结果经过质控后，就可以进入WES数据分析的第二步了，即检测样品中的变异。

第二部分：全外显子数据变异检测流程。

比对结果通过质量控制后，我们可以开始变异识别的流程。 也就是说找出比对结果里，样品数据和参考基因组不同的位点，并计算这些位点的基因型。　

在变异识别软件尽可能多地识别真正的生物变异时，也会识别一些非生物变异，即由比对或者测序错误带来的数据和参考基因组之间的差异，这些被识别的差异称为假阳性变异。在得到原始变异文件后，需要尽可能地去除所有非生物学变异。把真正的生物学变异从来自于系统误差的非生物学变异分离出来是一个比较棘手的问题。GATK 的VQSR （Variant Quality Score Recalibration）是目前最好的去除非生物学变异的软件。

VQSR的理论基础是： 变异软件识别的每个原始变异都有一套对应的注释参数。根据这些参数做聚类分析，真正的变异趋向于聚集在一起。而这些聚类形成的组群遵循高叶斯分布。根据已知的变异建立高叶斯混合模型，然后根据建立的模型对可能的新变异进行评估，进而去除假阳性的变异。

这个是2011年发表在nature genetics上的一篇文章利用高叶斯模型评估新变异的例子。左图是利用HapMap SNP作为训练数据建立的模型，右图是根据左图的模型对新SNP评估的结果。紫色的部分和左图中橘色的部分类似，很有可能是error， 是假阳性SNP， 所以要过滤掉。 绿色部分和左图中红色部分类似，很有可能是真正的新发现SNP。

在做VQSR时，SNP和INDEL要分别进行质量校正。但可以将VQSR连续run两次，而不必将SNP和INDEL分离、分别run VQSR之后再合并。以原始的包括SNP和INDEL的文件作为输入，在run第一步时，只对SNP进行质量分时校正，而INDEL不受影响。第二步对INDEL做质量分数校正，而SNP不受影响。

VQSR是目前最好的识别并过滤假阳性变异的方法，但只有原始变异集足够大的情况下，VQSR才能正常运行并过滤假阳性变异。在原始变异集不够大，不能够运行VQSR的情况下，可以通过对注释参数设定阈值来手工过滤非生物学变异。这种情况下，就必须先分离SNP和INDEL，分别手工过滤后再重新合并。

从前面IGV可视化比对结果的图里可以看出，比对到参考基因组的序列不只是在外显子区，很多在外显子以外的区域。那么变异软件识别的变异里有些也会分步在外显子区域以外。利用bedtools和捕获区域文件，可以去除外显子区以外的变异，只保留外显子区内的变异。

在得到外显子区域内高可信度的变异后，我们可以选取可靠的遗传致病数据库OMIM、HGMD、Clinvar、dbSNP等作为主要检索数据库对变异进行注释，然后结合样本来源的遗传信息或科研需要，筛选变异。这些内容我们会在下一节课里进行讲述。

感谢各位老师在百忙之中抽出时间收看及收听我们这次课程。欢迎各位老师加入我们的基因大课堂共同探讨、共同交流。

现场问答环节

1.什么是gVCF？gVCF和普通的VCF有什么区别？

答：VCF是Variant Call Format的缩写。是标准化的存储SNP，INDEL，和结构变异的text文件。VCF里只含有变异的位点信息。  gVCF是genomic VCF。是VCF的一种。gVCF 和VCF最大的区别是gVCF里含有基因组（或者感兴趣的区间）所有位点的碱基信息，不论该位点是否存在变异。

2.在数据量多大的情况下不可以用VQSR，而只能用手工过滤？

答：一般情况下，全外显子测序样品少于30个，或者全基因组测序样品少于2个，VQSR都不能很好的运行。即使运行通过不出错，结果也不可信。

3.怎样知道bam文件是否已经排序了？

答：可以用samtools查看bam文件的header。 如果在header里有SO:coordiate的字样，就是已经sort过的。

4.我先想问下WES为啥不能有效检测CNV，这个不是特别理解。这个人是纯合CNV或杂合CNV都不能有效检测吗？

答：是的，利用WES检测纯合和杂合的CNV都比较困难。这是因为人类的外显子都比较短，一般都小200bp。另外，WES测序的测序深度不均一， 用reads count检测CNV的话也会比较困难。

5.老师，因为外显子短于200，所以不好看CNV这个不明白。是因为缺的太多MAPPING不上了吗。是缺的比如100bp  剩下80bp是不能比对到参考基因吗？

答：不是mapping不上。而是说一般的CNV或者其他的SV都比较大，一般大于外显子。

6.有reads  您觉得counts数不准。 老师，这个测序深度不均一性是测序建库技术操作的问题 ，还是每个人的个体化差异导致不均一的现象。比如特定位置A基因是每个人的样本都会那不均一吗，还是有的人均一有的人不均一？

答：不均一的原因个体DNA和测序应该都有，看不同的情况。比如我们曾经分析过一个WES, 在一个本应该检测到变异的基因，我们无论如何检测不到。可视化发现该区域内没有reads。而查看这段基因序列，发现是100% 的G。而对测序仪来说，如果G含量超过80%，就很难成功测序。而如果个体某段DNA含有比较多的SNP，捕获探针就无法很好和DNA杂交，而不能有效的捕获。