【生信】全基因组测序（WGS）
文章的文字/图片/代码部分/全部来源网络或学术论文，文章会持续修缮更新，仅供大家学习使用。
目录
1、全基因组测序(WGS) 的定义
全基因组测序（WGS, Whole Genome Sequencing）是下一代测序技术，用于快速，低成本地确定生物体的完整基因组序列。其目的是准确检测出每个样本（这里特指人）基因组中的变异集合，也就是人与人之间存在差异的那些DNA序列。
生物的主要遗传物质是DNA，细胞或生物体中一套完整的遗传物质的总和称为基因组。DNA是由外显子(Exon）和内含子(Intron)组成。外显子只占基因组的1%左右，指导体内所有蛋白的合成。内含子不编码蛋白质的合成，但绝不是无用的序列，其可以影响基因的活性和蛋白的表达。全基因组测序，即是对生物体整个基因组序列进行测序，可以获得完整的基因组信息。
?
因此有三个需要区分的概念：
- WES (whole exome sequencing, 全外显子组测序): 是指利用序列捕获技术将【全基因组的外显子区域DNA】捕获富集后进行高通量测序，能够直接发现与蛋白质功能变异相关的遗传变异SNP(单核苷酸多态性)。以人类基因组为例，虽然外显子(蛋白质编码区)只占基因组的1%，但人类基因组85%的致病突变都在外显子区域，因此具有重要意义。
WES仅是对外显子区域的测序。相比全基因组测序，外显子测序更为简便，测序成本相比也会更低，测序后数据的分析也更为简单。另外如果仅关注已知的基因，可以更依赖于靶向测序。
- WGS (whole genome sequencing, 全基因组测序): 是指对基因组整体进行高通量测序，分析不同个体间的差异，同时完成SNP及基因组结构注释。全基因组测序由于结果包含完整丰富的信息，可以得到外显子测序或靶向测序不能得到的更多信息，具有其独特的优势。且随着近年来测序技术的不断进步、测序成本的不断降低，使得全基因组测序变得触手可及。而且全基因组测序在鉴定单核苷酸变异（SNP）、插入和缺失突变（Indel）时更有优势，所以WGS逐渐成为了临床和基础研究的另一种选择。
- WGRS (whole genome re-sequencing, 全基因组重测序): 是指对已知参考基因组和注释的物种进行不同个体间的全基因组测序(WGS)，并在此基础上对个体或群体进行差异性分析，鉴定出与某类表型相关的SNP。
三个技术的覆盖程度不同：
WES：覆盖全基因组上的外显子区域
WGS：覆盖全基因组
WGRS：覆盖全基因组，是WGS在不同样本上的重复
但其实这三个测序技术其实都是在找基因组上的SNP(单核苷酸多态性)，主要分为四种：
（1）SNV(单核苷酸变异)：包括同义突变和非同义突变。
（2）InDel(插入缺失)
（3）SV(结构变异)：包括50bp以上的长片段序列插入或者删除（Big Indel）、串联重复（Tandem repeate）、染色体倒位（Inversion）、染色体内部或染色体之间的序列易位（Translocation)等等
（4）CNV(拷贝数变异)
?
?
2、GWS流程
注意，WGS本质上只是一个技术手段。重要的是，我们要明白自己所要解决的问题是什么，希望获取的结果是什么，然后再选择合适的技术。这是许多人经常忽视的一个重要问题。
下面是GWS一个简易流程图：
?下面是在网上找的一个全基因组测序的流程和应用方向：
?
2.1准备工作——分析软件
- BWA (Burrow-Wheeler Aligner)： https://github.com/lh3/bwa 这是最权威，使用最广的NGS数据比对软件，将一个序列比对到参考基因组上的工具,能够快速的进行基因组比对。
- Samtools：https://github.com/samtools/samtools 是一个专门用于处理比对数据的工具，由BWA的作者（lh3）所编写，用于操作原始数据的工具,包括堆基因数据的查看,排序和合并等等。
- Picard：http://broadinstitute.github.io/picard/ 它是目前最著名的组学研究中心-Broad研究所开发的一款强大的NGS数据处理工具，功能方面和Samtools有些重叠，但更多的是互补，它是由java编写的，我们直接下载最新的.jar包就行了。
- GATK：https://software.broadinstitute.org/gatk/download/ 同样是Broad研究所开发的，是目前业内最权威、使用最广的基因数据变异检测工具。
- Sratoolkit：用于将NCBI数据装换的工具包
2.2原始数据质控
?
数据错误率的来源：
以illumina为例目前的测序技术基本都是运用边合成边测序的技术，碱基的合成依靠的是化学反应，这使得碱基链可以不断地从5'端一直往3'端合成并延伸下去。但在这个合成的过程中随着合成链的增长，DNA聚合酶的效率会不断下降，特异性也开始变差，这会导致越到后面碱基合成的错误率就会越高，这也是为何当前NGS（下一代测序）测序读长普遍偏短的一个原因。
测序数据的质量好坏会影响我们的下游分析，但不同的测序平台其测序错误率的图谱都是有差别的，因此，非常建议在我们分析测序数据之前先搞清楚如下两个地方：
- 原始数据是通过哪种测序平台产生的，它们的错误率分布是怎么样的，是否有一定的偏向性和局限性，是否会显著受GC含量的影响等。
- 评估它们有可能影响哪些方面的分析。
下图可以看出随着测序时间的增加，测序质量是逐渐下降的。
?
?
补充，如何认识一个原始的测序数据（fastq data）：
- read各个位置的碱基质量值分布。
- 碱基的总体质量值分布。
- read各个位置上碱基分布比例，目的是为了分析碱基的分离程度。
- GC含量分布。
- read各位置的N含量。
- read是否还包含测序的接头序列。
- read重复率，这个是实验的扩增过程所引入的。
切除测序接头序列和read的低质量序列：
目前有很多工具用来切除接头序列和低质量碱基，比如SOAPnuke、cutadapt、untrimmed等不下十个，但这其中比较方便好用的是Trimmomatic（一个java程序）、sickle和seqtk。
Trimmomatic的优势是它不但可以用来切除illumina测序平台的接头序列，还可以去除由我们自己指定的特定接头序列，而且同时也能够过滤read末尾的低质量序列。
sickle和seqtk只能去除低质量碱基。具体的原理就是通过滑动一定长度的窗口，计算窗口内的碱基平均质量，如果过低，就直接往后全部切除，不是挖掉read中的这部分低质量序列，而是直接从低质量区域开始把这条read后面的所有其它碱基全部切掉否则就是在人为改变实际的基因组序列情况。
如果下机的fq数据中不含有这些测序接头，那么我们除了trimmomatic之外，也可以直接使用sickle（同时支持PE和SE数据）或者seqtk（仅支持SE），这两个处理起来会更快，消耗的计算资源也更少。
?
?2.3数据预处理
（1）序列（read）比对：
定义：测序得到的短序列（read）存储于FASTQ文件里面。
虽然它们原本都来自于有序的基因组，但在经过DNA建库和测序之后，文件中不同read之间的前后顺序关系就已经全部丢失了。
因此，FASTQ文件中紧挨着的两条read之间没有任何位置关系，它们都是随机来自于原本基因组中某个位置的短序列。
因此，需要先把这一大堆的短序列一个个去跟该物种的参考基因组比较，找到每一条read在参考基因组上的位置，然后按顺序排列好，这个过程就称为测序数据的比对。
使用BWA工具：用于流程构建的BWA就是其中最优秀的一个，它将BW(Burrows-Wheeler)压缩算法和后缀树相结合，能够让我们以较小的时间和空间代价，获得准确的序列比对结果。
?
（2）排序（sort）：
为什么要排序：BWA比对后输出的BAM文件是没顺序的。
因为FASTQ文件里面这些被测序下来的read是随机分布于基因组上面的，第一步的比对是按照FASTQ文件的顺序把read逐一定位到参考基因组上之后，随即就输出了，它不会也不可能在这一步里面自动识别比对位置的先后位置重排比对结果。
因此，比对后得到的结果文件中，每一条记录之间位置的先后顺序是乱的，我们后续去重复等步骤都需要在比对记录按照顺序从小到大排序下来才能进行，所以需要进行排序。
使用Samtools工具：用于操作原始数据的工具,包括堆基因数据的查看,排序和合并等等。
（3）去除重复序列（或者标记重复序列）：
也就是去除PCR重复序列。
重复序列来源：一般测序之前都需要先构建测序文库：通过物理（超声）打断或者化学试剂（酶切）切断原始的DNA序列，然后选择特定长度范围的序列去进行PCR扩增并上机测序。
因此，这里需要去除的重复序列来源实际上就是由PCR过程中所引入的。因为所谓的PCR扩增就是把原来的一段DNA序列复制多次。
PCR扩增的目的：很多时候构建测序文库时能用的细胞量并不会非常充足，而且在打断的步骤中也会引起部分DNA的降解，这两点会使整体或者局部的DNA浓度过低，这时如果直接从这个溶液中取样去测序就很可能漏掉原本基因组上的一些DNA片段，导致测序不全。而PCR扩增的作用就是为了把这些微弱的DNA多复制几倍乃至几十倍，以便增大它们在溶液中分布的密度，使得能够在取样时被获取到。
去除重复序列的目的：PCR扩增原本的目的是为了增大微弱DNA序列片段的密度，但由于整个反应都在一个试管中进行，因此其他一些密度并不低的DNA片段也会被同步放大，那么这时在取样去上机测序的时候，这些DNA片段就很可能会被重复取到相同的几条去进行测序。这会导致最后结果同时增大了变异检测结果的假阴和假阳率。
主要有几个原因：
- DNA在打断的那一步会发生一些损失，主要表现是会引发一些碱基发生颠换变换（嘌呤-变嘧啶或者嘧啶变嘌呤），带来假的变异。PCR过程会扩大这个信号，导致最后的检测结果中混入了假的结果；
- PCR反应过程中也会带来新的碱基错误。发生在前几轮的PCR扩增发生的错误会在后续的PCR过程中扩大，同样带来假的变异；
- 对于真实的变异，PCR反应可能会对包含某一个碱基的DNA模版扩增更加剧烈（这个现象称为PCR Bias）。因此，如果反应体系是对含有reference allele的模板扩增偏向强烈，那么变异碱基的信息会变小，从而会导致假阴。
工具：Samtools和Picard。其主要原理是对所用的序列数据进行标记（去除）或者
Samtools的rmdup是直接将这些重复序列从比对BAM文件中删除掉。使用PCR-Free的测序方案
而Picard的MarkDuplicates默认情况则只是在BAM的FLAG信息中标记出来，而不是删除，因此这些重复序列依然会被留在文件中，只是我们可以在变异检测的时候识别到它们，并进行忽略。
?
（4）局部重比对：
局部重比对通常也叫Indel局部区域重比对。
在局部重比对之前需要将同一个样本的所有比对结果合并成唯一一个大的BAM文件（merge操作）。
局部重比对目的：将BWA比对过程中所发现有潜在序列插入或者序列删除（insertion和deletion，简称Indel）的区域进行重新校正，这个过程往往还会把一些已知的Indel区域一并作为重比对的区域。
局部重比对的原因：
- 来自于参考基因组的序列特点和BWA这类比对算法本身，注意这里不是针对BWA，而是针对所有的这类比对算法，包括bowtie等。这类在全局搜索最优匹配的算法在存在Indel的区域及其附近的比对情况往往不是很准确，特别是当一些存在长Indel、重复性序列的区域或者存在长串单一碱基（比如，一长串的TTTT或者AAAAA等）的区域中。
- 在这些比对算法中，对碱基错配和开gap的容忍度是不同的。具体体现在罚分矩阵的偏向上，例如，在read比对时，如果发现碱基错配和开gap都可以的话，它们会更偏向于错配。但是这种偏向错配的方式，有时候却还会反过来引起错误的开gap。这就会导致基因组上原本应该是一个长度比较大的Indel的地方，被错误地切割成多个错配和短indel的混合集，这必然会让我们检测到很多错误的变异。这种情况还会随着所比对的read长度的增长（比如三代测序的Read，通常都有几十kbp）而变得越加严重。
Smith-Waterman算法可以极其有效地实现对全局比对结果的校正和调整，最大程度低地降低由全局比对算法的不足而带来的错误。
而且GATK的局部重比对模块，除了应用这个算法之外，还会对这个区域中的read进行一次局部组装，把它们连接成为长度更大的序列，这样能够更进一步提高局部重比对的准确性。
?
?
（5）重新校正碱基质量值（BQSR）
BQSR（Base Quality Score Recalibration）主要是通过机器学习的方法构建测序碱基的错误率模型，然后对这些碱基的质量值进行相应的调整。
在WGS分析中，变异检测是一个极度依赖测序碱基质量值的步骤。这个质量值是衡量测序出来的碱基正确率的重要（甚至是唯一）指标。它来自于测序图像数据的base calling，因此，基本上是由测序仪和测序系统来决定的。
下图中横轴（Reported quality score）是测序结果在Base calling之后报告出来的质量值，也就是我们在FASTQ文件中看到的那些；纵轴（Empirical quality score）代表的是“真实情况的质量值”（即通过统计学的技巧获得极其接近的分布结果）。
?2.4变异检测
?
变异检测的内容一般会包括：SNP、Indel，CNV和SV等。
工具：GATK HaplotypeCaller模块对样本中的变异进行检测，它也是目前最适合用于对二倍体基因组进行变异（SNP+Indel）检测的算法。HaplotypeCaller和那些直接应用贝叶斯推断的算法有所不同，它会先推断群体的单倍体组合情况，计算各个组合的几率，然后根据这些信息再反推每个样本的基因型组合。因此它不但特别适合应用到群体的变异检测中，而且还能够依据群体的信息更好地计算每个个体的变异数据和它们的基因型组合。
变异检测质控和过滤（VQSR）在获得了原始的变异检测结果之后，还需要做的就是质控和过滤。VQSR是通过构建GMM模型对好和坏的变异进行区分，从而实现对变异的质控。对于人类而言，一般来说，每个人最后检测到的变异数据大概在400万左右（包括SNP和Indel）。