2018-09-24 wes流程文件解读2
上次我们整理到bwa比对后得到bam文件，下一步我们要通过GATK流程从bam文件中call variant。
一、使用GATK前须知事项：
- 对GATK的测试主要使用的是人类全基因组和外显子组的测序数据，而且全部是基于illumina数据格式，目前还没有提供其他格式文件（如Ion Torrent）或者实验设计（RNA-Seq）的分析方法。
- GATK是一个应用于前沿科学研究的软件，不断在更新和修正，因此，在使用GATK进行变异检测时，最好是下载最新的版本。
- 在GATK使用过程中，有些步骤需要用到已知变异信息，对于这些已知变异，GATK只提供了人类的已知变异信息 ，可以在GATK的FTP站点下载（GATK resource bundle）。如果要研究的不是人类基因组，需要自行构建已知变异，GATK提供了详细的构建方法。
- GATK在进行BQSR和VQSR的过程中会使用到R软件绘制一些图，因此，在运行GATK之前最好先检查一下是否正确安装了R和所需要的包，所需要的包大概包括ggplot2、gplots、bitops、caTools、colorspace、gdata、gsalib、reshape、RColorBrewer等。如果画图时出现错误，会提示需要安装的包的名称。
第一步：原始数据的处理：对原始下机fastq文件进行过滤和比对（mapping）对于Illumina下机数据推荐使用bwa进行mapping。前边已讲。
bwa中 -r str：定义头文件。‘@RG\tID:foo\tSM:bar’，如果在此步骤不进行头文件定义，在后续GATK分析中还是需要重新增加头文件。GATK2.0以上版本将不再支持无头文件的变异检测。 - ID str：输入reads集ID号；LB：read集文库名；PL：测序平台（illunima或solid）；PU：测序平台下级单位名称（run的名称）；SM：样本名称
- 对于最后得到的sam文件，将比对上的结果提取出来（awk即可处理），即可直接用于GATK的分析。
注意：由于GATK在下游的snpcalling时，是按染色体进行callsnp的。因此，在准备原始sam文件时，可以先按染色体将文件分开，这样会提高运行速度。但是当数据量不足时，可能会影响后续的VQSR分析，这是需要注意的。 - 对sam文件进行进行重新排序（reorder）
由BWA生成的sam文件时按字典式排序法进行的排序（lexicographically）进行排序的（chr10，chr11…chr19，chr1，chr20…chr22，chr2，chr3…chrM，chrX，chrY），但是GATK在进行callsnp的时候是按照染色体组型（karyotypic）进行的（chrM，chr1，chr2…chr22，chrX，chrY），因此要对原始sam文件进行reorder。可以使用picard-tools中的ReorderSam完成。 - 对bam文件进行sort排序处理 一步是将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序。可以使用picard-tools中SortSam完成。
第一步： Duplicates Marking - 在制备文库的过程中，由于PCR扩增过程中会存在一些偏差，也就是说有的序列会被过量扩增。这样，在比对的时候，这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。而这些过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由PCR扩增所形成的duplicates，这一步可以使用picard-tools来完成。
- 去重复的过程是给这些序列设置一个flag以标志它们，方便GATK的识别。还可以设置 REMOVE_DUPLICATES=true 来丢弃duplicated序列。对于是否选择标记或者删除，对结果应该没有什么影响，GATK官方流程里面给出的例子是仅做标记不删除。这里定义的重复序列是这样的：如果两条reads具有相同的长度而且比对到了基因组的同一位置，那么就认为这样的reads是由PCR扩增而来，就会被GATK标记。
最主要目的就是尽量减小文库构建时引入文库的PCR bias
标记后对结果文件生成索引文件。
GATK=/home/jmzeng/biosoft/gatk4/gatk-4.0.6.0/gatk
ref=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
snp=/public/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz
indel=/public/biosoft/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
for sample in {7E5239.L1,7E5240,7E5241.L1}
do
echo $sample
$GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates \
-I $sample.bam \
-O ${sample}_marked.bam \
-M $sample.metrics \
1>${sample}_log.mark 2>&1
samtools index ${sample}_marked_fixed.bam
第二步：Local realignment around indels
这一步的目的就是将比对到indel附近的reads进行局部重新比对，将比对的错误率降到最低。一般来说，绝大部分需要进行重新比对的基因组区域，都是因为插入/缺失的存在，因为在indel附近的比对会出现大量的碱基错配，这些碱基的错配很容易被误认为SNP。还有，在比对过程中，比对算法对于每一条read的处理都是独立的，不可能同时把多条reads与参考基因组比对来排错。因此，即使有一些reads能够正确的比对到indel，但那些恰恰比对到indel开始或者结束位置的read也会有很高的比对错误率，这都是需要重新比对的。
主要分为两步：
- 通过运行RealignerTargetCreator来确定要进行重新比对的区域。
java -jar GenomeAnalysisTK.jar
-R hg19.fa
-T RealignerTargetCreator
-I hg19.reorder.sort.addhead.dedup_04.bam
-o hg19.dedup.realn_06.intervals
-known Mills_and_1000G_gold_standard.indels.hg38.vcf
-known 1000G_phase1.indels.hg38.vcf
参数说明：
-R： 参考基因组；
-T： 选择的GATK工具；
-I： 输入上一步所得bam文件；
-o： 输出的需要重新比对的基因组区域结果；
-maxInterval:允许进行重新比对的基因组区域的最大值，不能太大，太大耗费会很长时间，默认值500；
-known:已知的可靠的indel位点，指定已知的可靠的indel位点，重比对将主要围绕这些位点进行，对于人类基因组数据而言，可以直接指定GATK resource bundle里面的indel文件（必须是vcf文件）。
- 通过运行IndelRealigner在这些区域内进行重新比对
-R hg19.fa
-T IndelRealigner
-targetIntervals hg19.dedup.realn_06.intervals
-I hg19.reorder.sort.addhead.dedup_04.bam
-o hg19.dedup.realn_07.bam
-known Mills_and_1000G_gold_standard.indels.hg38.vcf
-known 1000G_phase1.indels.hg38.vcf
注意：1. 第一步和第二步中使用的输入文件（bam文件）、参考基因组和已知indel文件必须是相同的文件。
2. 当在相同的基因组区域发现多个indel存在时，这个工具会从其中选择一个最有可能存在比对错误的indel进行重新比对，剩余的其他indel不予考虑。
第三步：.Base quality score recalibration
这一步是对bam文件里reads的碱基质量值进行重新校正，使最后输出的bam文件中reads中碱基的质量值能够更加接近真实的与参考基因组之间错配的概率。这一步适用于多种数据类型，包括illunima、solid、454、CG等数据格式。在GATK2.0以上版本中还可以对indel的质量值进行校正，这一步对indel calling非常有帮助。
举例说明，在reads碱基质量值被校正之前，我们要保留质量值在Q25以上的碱基，但是实际上质量值在Q25的这些碱基的错误率在1%，也就是说质量值只有Q20，这样就会对后续的变异检测的可信度造成影响。还有，在边合成边测序的测序过程中，在reads末端碱基的错误率往往要比起始部位更高。另外，AC的质量值往往要低于TG。BQSR的就是要对这些质量值进行校正。
- 利用工具BaseRecalibrator，根据一些known sites，生成一个校正质量值所需要的数据文件，GATK网站以“.grp”为后缀命名。
-R $ref \
-I ${sample}_marked_fixed.bam \
--known-sites $snp \
--known-sites $indel \
-O ${sample}_recal.table \
1>${sample}_log.recal 2>&1
第四步：ApplyBQSR
利用已知snp数据库，通过机器学习方法调整碱基质量分数
-R $ref \
-I ${sample}_marked_fixed.bam \
-bqsr ${sample}_recal.table \
-O ${sample}_bqsr.bam \
1>${sample}_log.ApplyBQSR 2>&1
这时候，GATK流程文件准备工作结束，完成了variants calling所需要的所有准备，生成了用于下一步变异检测的bam文件。
第五步：Variant Calling
GATK在这一步里面提供了两个工具进行变异检测——UnifiedGenotyper和HaplotypeCaller
- HaplotypeCaller不支持Reduce之后的bam文件，因此，当选择使用HaplotypeCaller进行变异检测时，不需要进行Reduce reads
-R $ref \
-I ${sample}_bqsr.bam \
--dbsnp $snp \
-O ${sample}_raw.vcf \
1>${sample}_log.HC 2>&1
对输入的bam文件中的所有样本进行变异检测，最后生成一个vcf文件，vcf文件中会包含所有样本的变异位点和基因型信息。但是现在所得到的结果是最原始的、没有经过任何过滤和校正的Variants集合。这一步产生的变异位点会有很高的假阳性，尤其是indel，因此，必须要进行进一步的筛选过滤。这一步还可以指定对基因组的某一区域进行变异检测，只需要增加一个参数 -L：target_interval.list，target_interval.list 格式是bed格式文件。
- bed文件格式，三列必须
TCGAGA
#对应bed文件的坐标应为
#chrome start end
chr1 0 5
- 注意：GATK进行变异检测的时候，是按照染色体排序顺序进行的（先call chr1，然后chr2，然后chr3…最后chrY），并非多条染色体并行检测的，因此，如果数据量比较大的话，建议分染色体分别进行，对性染色体的变异检测可以同常染色体方法。
大多数参数的默认值可以满足大多数研究的需求，因此，在做变异检测过程中，如果对参数意义不是很明确，不建议修改。
第六步： 对原始变异检测结果进行过滤（hard filter and VQSR）
这一步的目的就是对上一步call出来的变异位点进行过滤，去掉不可信的位点。这一步可以有两种方法，一种是通过GATK的VariantFiltration，另一种是通过GATK的VQSR（变异位点质量值重新校正）进行过滤。
GATK是推荐使用VQSR的，但使用VQSR数据量一定要达到要求，数据量太小无法使用高斯模型。在使用VAQR时，indel和snp要分别进行。 - VQSR原理介绍：
这个模型是根据已有的真实变异位点（人类基因组一般使用HapMap3中的位点，以及这些位点在Omni 2.5M SNP芯片中出现的多态位点）来训练，最后得到一个训练好的能够很好的评估真伪的错误评估模型，可以叫他适应性错误评估模型。这个适应性的错误评估模型可以应用到call出来的原始变异集合中已知的变异位点和新发现的变异位点，进而去评估每一个变异位点发生错误的概率，最终会给出一个得分。这个得分最后会被写入vcf文件的INFO信息里，叫做VQSLOD，就是在训练好的混合高斯模型下，一个位点是真实的概率比上这个位点可能是假阳性的概率的log odds ratio（对数差异比），因此，可以定性的认为，这个值越大就越好。
VQSR主要分两个步骤，这两个步骤会使用两个不同的工具：VariantRecalibrator和ApplyRecalibration。
VariantRecalibrator：通过大量的高质量的已知变异集合的各个注释（包括很多种，后面介绍）的值来创建一个高斯混合模型，然后用于评估所有的变异位点。这个文件最后将生成一个recalibration文件。
原理简单介绍： 这个模型首先要拿到真实变异数据集和上一步骤中得到的原始变异数据集的交集，然后对这些SNP值相对于具体注释信息的分布情况进行模拟，将这些变异位点进行聚类，最后根据聚类结果赋予所有变异位点相应的VQSLOD值。越接近聚类核心的变异位点得到的VQSLOD值越高。
ApplyRecalibration：这一步将模型的各个参数应用于原始vcf文件中的每一个变异位点，这时，每一个变异位点的注释信息列中都会出现一个VQSLOD值，然后模型会根据这个值对变异位点进行过滤，过滤后的信息会写在vcf文件的filter一列中。
原理简单介绍：在VariantRecalibrator这一步中，每个变异位点已经得到了一个VQSLOD值了，同时，这些LOD值在训练集里也进行了排序。当你在这一步中设置一个tranche sensitivity 的阈值（这个阈值一般是一个百分数，如设置成99%），那么，如果LOD值从大到小排序的话，这个程序就会认为在这个训练集中，LOD值在前99%的是可信的，当这个值低于这个阈值，就认为是错误的。最后，程序就会用这个标准来过滤上一步call出来的原始变异集合。如果LOD值超过这个阈值，在filter那一列就会显示PASS，如果低于这个值就会被过滤掉，但是这些位点仍然会显示在结果里面，只不过会在filter那一列标示出他所属于的tranche sensitivity 的名称。在设置tranche sensitivity 的阈值时，要兼顾敏感度和质量值。
tranche值的设定
前面提到了，这个值得设定是用来在后续的ApplyRecalibration中如何根据这个阈值来过滤变异位点的，也就是说，如果这个值设定的比较高的话，那么最后留下来的变异位点就会多，但同时假阳性的位点也会相应增加；如果设定的低的话，虽然假阳性会减少，但是会丢失很多真实的位点。因此，跟选择注释时一样，可以run两遍VariantRecalibrator，第一遍的时候多写几个阈值，第一遍跑完之后看结果，看那个阈值好，选择一个最好的阈值，再run一遍VariantRecalibrator。至于说怎么区分好坏，有几个标准：
- 看结果中已知变异位点与新发现变异位点之间的比例，这个比例不要太大，因为大多数新发现的变异都是假阳性，如果太多的话，可能假阳性的比例就比较大；
- 看保留的变异数目，这个就要根据具体的需求进行选择了。
- 看TI/TV值，对于人类全基因组，这个值应该在2.15左右，对于外显子组，这个值应该在3.2左右，不要太小或太大，越接近这个数值越好，这个值如果太小，说明可能存在比较多的假阳性。
千人中所选择的tranche值是99
注意：Indel不支持tranche值的选择，另外，一部分注释类型在做indel的校正时也不支持，具体信息可以详查GATK网站。
当数据量太小时，可能高斯模型不会运行，因为变异位点数满足不了模型的统计需求。这时候可以通过降低--maxGaussian的值，让程序运行。这个值表示的是程序将变异位点分成的最大的组数，降低这个值让程序把变异位点聚类到更少的组里面，使每个组中的变异位点数增加来满足统计需求，但是这样做降低程序分辨真伪的能力。因此，在运行程序的时候，要对各方面进行权衡。
过滤后生成的vcf文件的格式说明，即每一列所代表的的内容，可参考下面的网站，有详细的说明
http://www.broadinstitute.org/gatk/guide/article?id=1268
第七步：初步注释分析利用GATK中的VariantEval来完成，或者用annovar注释
ANNOVAR是由王凯编写的一个注释软件，可以对SNP和indel进行注释，也可以进行变异的过滤筛选。
ANNOVAR能够利用最新的数据来分析各种基因组中 的遗传变异。主要包含三种不同的注释方法，Gene-based Annotation（基于基因的注释）、Region-based Annotation（基于区域的注释）、Filter-based Annotation（基于筛选的注释）。ANNOVAR由Perl编写。
优点：提供多个数据可直接下载、支持多种格式、注释直观；缺点：没有数据库的物种无法注释。
│ annotate_variation.pl #主程序，功能包括下载数据库，三种不同的注释
│ coding_change.pl #可用来推断蛋白质序列
│ convert2annovar.pl #将多种格式转为.avinput的程序
│ retrieve_seq_from_fasta.pl #用于自行建立其他物种的转录本
│ table_annovar.pl #注释程序，可一次性完成三种类型的注释
│ variants_reduction.pl #可用来更灵活地定制过滤注释流程
│
├─example #存放示例文件
│
└─humandb #人类注释数据库
ANNOVAR下载数据库
命令示例
# -buildver 表示version
# -downdb 下载数据库的指令
# -webfrom annovar 从annovar提供的镜像下载，不加此参数将寻找数据库本身的源
# humandb/ 存放于humandb/目录下
输出的csv文件将包含输入的5列主要信息以及各个数据库里的注释，此外，table_annoval.pl可以直接对vcf文件进行注释（不需要转换格式），注释的内容将会放在vcf文件的“INFO”那一栏。
说明：本文主要说明wes分析流程的理论，代码仅作示例。
参考文献
1.https://www.jianshu.com/p/49d035b121b8
2.https://blog.csdn.net/zhu_si_tao/article/details/53321374
3.https://www.plob.org/article/9976.html
