075-关于WES的简介

刘小泽写于19.1.4

今天随便写了些关于外显子的介绍，飞来珠海见了一群小伙伴，十分开心

WES（Whole Exome Sequencing）全外显子测序是对基因组中所有外显子进行测序分析的方法。首先利用DNA或RNA探针，将全基因组中外显子区域DNA进行捕获，然后对捕获的DNA进行PCR扩增，最后对扩增后的产物进行高通量测序。

什么是外显子？

外显子只占基因组的2%，却含有85%的已知致病变异，物美价廉

WES优势

• 总花费低：虽然WES一般比WGS测序深度高（100X vs 30X), 但是WES的测序序列主要在占基因组2%的外显子区域，因此总的测序文件比WGS的测序文件小得多【平均测序深度为100X的全外显子测序文件的大小是6GB左右，而30X的全基因组测序文件的大小一般是90GB左右】
• 测序深度高，WES更容易检测出罕见变异
• 数据文件小，数据分析的速度就快
• 跟非编码区相比，蛋白编码区的研究更透测，注释度更高

WES不足

• 测序深度不均一，一些区域测序深度过高造成浪费，一些区域测序深度过低而无法检测到存在的变异。例如：一些SNP密集的区域，富集过程中RNA或DNA探针无法和该区域杂交导致无法有效捕获，该区域的测序深度就会过低；
• 只测外显子区域（人类的外显子区域长度一般小于200bp），因而无法有效检测CNV、SV，一般用来检测SNV、INDEL
• 受探针限制，新的编码基因无法捕获

评价质量指标

• 碱基质量分数：测序过程中碱基被测序仪错误识别的概率
• 比对质量分数：一个序列在参考基因组上比对正确与否的可信度
• 变异质量分数：检测到的变异是否是生物变异的可信度【变异质量分数是变异位点所有序列比对质量分数的均方根】

数据分析流程

第一部分：前处理流程

原始fq质控

第一步一定要查看fastq文件是否完整，是否存在错误

过滤

过滤的话一般是：去除低质量或者过短序列，剪切末端低质量的碱基，去除接头污染以及可能的外源污染序列

比对

• 比对：序列比对到可信度最高的一个或者少数几个位置上去，可信度就是比对质量分数（mapping quality ）【问题：由于参考基因组中存在重复序列等原因，有些序列可以同等地比对到不同的位置】

• 比对后可视化：IGV【用加过索引的BAM文件】，结果显示，外显子区域以外也有比对的序列。 一般情况下，这些序列不要去除，大范围内更多的比对序列有利于下游更精确的变异检测 IGV默认颜色：灰色表示该序列在基因组上的比对结果是单一的，且和序列对中另一个序列间的距离在正常范围内；白色说明该序列比对到2个或者2个以上的位置；其他颜色，说明与该序列配对的序列比对到了其它的染色体， 或者两个序列间的距离大于或小于正常的插入序列长度

比对后质控

• 看比对情况：samtools flagstat 可以看总的比对率，成功比对的序列里多少个是pair-end序列，多少是单个的序列；idxstat 查看每个染色体的长度，该染色体上比对的序列数以及没有比对的序列数。
• 原始比对结果质控：如果有一些比对质量较低或者错误比对的序列，为了下游变异识别的敏感性和精确性，减少假阳性生物变异的识别，需要质控，比如：去除重复序列、可疑的比对区域进行重新比对（Local Realignment）、重新校正碱基质量分数
• 关于重复序列：建库过程中，需要进行多轮的PCR扩增。理想情况是每次PCR，每个区域扩增一次，得到两份。但是每次PCR中，扩增引物和不同模板结合力会存在差异，因此有的区域扩增产物大于一，就得到许多份，造成后续等位基因频率的定义以及基因型识别不准确。一般用Markduplcate功能去除重复序列
• 关于可疑的比对区域进行重新比对：最初的比对过程中，比对软件将每个序列单独和参考基因组进行比对，但有时比对软件无法确定INDEL处是INDEL还是多个SNP，导致INDEL无法检测、SNP假阳性。GATK local realignment 就是对这些区域重新比对，使所有序列错配（mismatch）的碱基数目总和最小【例如：本来参考基因组中存在7个连续的T，测序数据中在这个位置出现了一个deletion，但是比对软件偏偏认为它是多个碱基错配，导致结果中少了一个deletion，反倒多了多几个SNV】
• 关于碱基的质量分数校正：由于系统误差（可能来源测序过程化学反应试剂或者仪器报告不准），测序仪报告的碱基质量不精确，比实际质量分数偏高或者偏低。系统误差不同于随机误差，它是错误，应该被去除。 GATK recalibration 利用机器学习的方法建立误差模型，然后对碱基分数进行调整，调整后可以提高后续变异识别的准确率，减少假阳性和假阴性的变异识别。但是，这个分数校正不能纠正碱基【例如：有一个低质量的C碱基，我们不知道它是不是应该是G，但是可以知道我们相信C是正确的概率是多少】

第二部分：变异检测流程

目的：找出比对结果中，样品数据和参考基因组不同的位点，并计算这些位点的基因型

假阳性的存在：

变异识别软件（variant caller）不可避免会识别一些非生物变异（也就是没有生物学意义的变异，例如：比对或者测序错误带来的数据和参考基因组之间的差异）。我们要努力去除假阳性带来的影响。

但是如何区分真正的生物学变异和系统误差的非生物学变异呢？

目前最经常使用的就是：GATK VQSR （Variant Quality Score Recalibration）

关于VQSR

变异软件识别的每个原始变异都有一套对应的注释参数（如VCF的header行），根据这些参数进行聚类分析，真正的变异趋向于聚集在一起。

VQSR基于机器学习，基于群体遗传的原理，构建GMM模型，挑出和已知的变异集合Overlap的位点（如：已知且被严格验证的HapMap数据集）并分配相应的可信度权重来进行训练，训练出一个区分好变异的GMM，然后对全部数据进行打分，再把评分最低的那些拿出来，构成一个最不像正确变异的集合，用来构造一个区分坏变异的GMM，用来专门识别坏变异。最后同时用好和坏的GMM再一次同时对变异进行打分，看每个变异更像谁，就能够评判出这个变异可信的质量值。

VQSR要求比较高：好和坏变异可供训练的数目必须超过5000个，如果Overlap位点太少，无法有效训练出一个模型。因此对于WGS问题不大，但是单个的WES中加起来全部区域也就50Mb左右，外显子比较短，一般都小200bp，其中变异数目大约30K-40K，位点并不多，那么和已知的高质量变异集的overlap就更少。因此对于单个WES或者小panel（低于30 sample） 的都不适用，起码要30 sample+吧

VQSR过程一般需要SNP与Indel分别校正 ，但是也可以连续进行两次VQSR校正（第一次只对SNP校正，第二次只对Indel校正，最后再合并）

当然，如果sample不够多，数据集不够大，也可以设定阈值来手工过滤非生物学变异【注意：先分离SNP和INDEL，分别手工过滤后再重新合并】

第三部分 变异筛选、注释

筛选

因为比对会比对到参考基因组外显子区域外，因此变异识别软件识别出的变异也会有一些是分布在外显子区域之外的，可以可用bedtools和捕获区域文件，去除外显子区域外的变异，只保留外显子区域内的变异

注释

得到了外显子区域内高可信度的变异以后，可以在OMIM、HGMD、Clinvar、dbSNP等数据库进行注释

补充：

• 接头污染问题：接头序列里一般包含三个重要组成部分：区分样品的barcode序列，PCR primer序列和测序引物结合的序列。插入测序列过短时，测序反应会超过待测序列而测到3’-端的接头序列，从而造成接头序列的污染
• gVCF：gVCF 和VCF最大的区别是gVCF里含有基因组（或者感兴趣的区间）所有位点的碱基信息，不论该位点是否存在变异。

关于WES：https://ghr.nlm.nih.gov/primer/testing/sequencing

https://blueprintgenetics.com/tests/whole-exome-sequencing/

关于gVCF：https://gatkforums.broadinstitute.org/gatk/discussion/4017/what-is-a-gvcf-and-how-is-it-different-from-a-regular-vcf