肿瘤外显子数据处理系列教程（五）GATK的最佳实践

上一节我们讲到了，讲bam文件载入igv中可视化，不过是番外篇，今天我们继续上一次比对后的流程：GATK的最佳实践。

下载数据

我们在肿瘤外显子数据处理系列教程（三）就说过了GATK流程需要哪些数据，不过当时没有演示，这里就补充一下，其实很简单，就是几行命令：

nohup wget ftp://[e-mail]/bundle/hg38/dbsnp_146.hg38.vcf.gz & 
nohup wget ftp://[e-mail]/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi & 
nohup wget ftp://[e-mail]/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz & 
nohup wget ftp://[e-mail]/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi & 
nohup wget ftp://[e-mail]/bundle/hg38/Homo_sapiens_assembly38.fasta.gz & 
nohup wget ftp://[e-mail]/bundle/hg38/Homo_sapiens_assembly38.fasta.fai & 
nohup wget ftp://[e-mail]/bundle/hg38/Homo_sapiens_assembly38.dict & 
nohup wget ftp://[e-mail]/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz & 
nohup wget ftp://[e-mail]/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi & 

---

而下载这些数据的链接，是从GATK的ftp服务器上找到的，可以进入gatk官网：https://software.broadinstitute.org/gatk/，找到下面的入口：

然后拉到最后，有个ftp的链接：

ftp://[e-mail]/bundle/

找到我们要下载的数据，比如hg38，点进去

找到需要下载的数据，右键，复制链接地址，然后再服务器上用wget命令下载即可，当然这只是其中一种下载方法，比较简单，新手容易接受（墙内玩家体验感比较差）

有了数据之后，我们就可以开始GATK的最佳实践流程了，下图是官方示意图，我们已经拿到比对的bam文件了，所以接下来走Mark Duplicates和BQSR：

标记或去除重复序列

在拿到的测序结果中，有一部分序列是重复的，这是由于上机前需要进行PCR扩增，如果序列扩增次数不同，就会导致PCR重复，重复序列在igv可视化的时候就可以明显看出来。比如我们前面提取的small.bam在这个区域chr17:43,719,044-43,720,613的覆盖深度就出现了断崖式的增长。

所以我们拿到比对后的bam文件之后，要进行标记重复序列甚至去除掉重复序列。用到了GATK的MarkDuplicates工具，代码如下：

GATK=/home/jmzeng/biosoft/gatk4/gatk-4.0.6.0/gatk

cat ../case  | while read id
do
    BAM=../4.align/${id}.bam

    if [ ! -f ${id}_marked.bam ]
    then
    echo "start MarkDuplicates for ${id}" `date`
    $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates \
    -I ${BAM} \
    --REMOVE_DUPLICATES=true \
    -O ${id}_marked.bam \
    -M ${id}.metrics \
    1>${id}_log.mark 2>&1 
    echo "end MarkDuplicates for ${id}" `date`
    samtools index -@ 16 -m 4G -b ${id}_marked.bam
    fi
done

解释一下上面的代码：首先把case文件做好，也就是我们的样本名，在这里作为一个配置文件，在前面的教程中我们也一直提到它，里面记录的内容是：

$ cat ../case
case1_biorep_A_techrep
case1_biorep_B
case1_biorep_C
......
case6_biorep_C
case6_germline
case6_techrep_2

进入到while read id;do...;done循环后，每read读取一行为一个id，也就是一个样本名，对于每一个样本，再赋值bam文件的具体路径，命名为BAM。接下来就是一个if判断语句，判断是否已经存在生成的${id}_marked.bam 文件，如果存在，就不再执行，如果不存在就继续往下走。这样写控制语句的原因，是为了避免重复提交代码，如果多个样本批量跑，跑到一半的时候任务挂了，比如服务器突然停电了，那后续再次提交命令的时候，已经完成的样本就不会被再次执行，我们要做的就是重新提交这个脚本而已。紧接着就是echo start...date这是为了记录时间，和后面的echo end...date相对应，这样我们可以知道跑完一个样本需要的时间，也有的人喜欢在命令前加上一个time，也可以达到类似的效果。中间部分就是主要的命令：

-Xmx20G -Djava.io.tmpdir=./：设置调用的内存大小为20G

-I ${BAM}：指定输入文件

--REMOVE_DUPLICATES=true：删除掉重复序列，如果不加这个参数，就只是标记重复序列而不会删除。

-O ${id}_marked.bam：输出文件

1>${id}_log.mark 2>&1：这里的1和2是捕获正确和错误输出流，把内容都输出到${id}_log.mark

最后是用samtools对生成的bam文件构建索引，这是后续步骤的要求。

简单查看一下删除重复序列的结果，去除Duplicates之前的总reads数为212761899，删掉Duplicates之后剩下192699040了：

$ samtools view 4.align/case1_biorep_A_techrep.bam |wc -l
212761899
$ samtools view 5.gatk/case1_biorep_A_techrep_marked.bam |wc -l
192699040

碱基质量重校正BQSR

因为对于WGS或者WES的测序数据，我们最关心的问题就是变异检测了，而变异检测和碱基质量值密切相关。而碱基的质量值有测序仪和测序系统产生，机器带来的系统误差和噪声是不可避免的。因此，找变异之前，需要进行碱基质量重校正BQSR（Base Quality Score Recalibration）。用到了GATK的两个工具BaseRecalibrator和ApplyBQSR，代码如下：

GATK=/home/jmzeng/biosoft/gatk4/gatk-4.0.6.0/gatk
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
ref=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
cat ../case  | while read id
do
    if [ ! -f ${id}_bqsr.bam ]
    then
    echo "start BQSR for ${id}" `date`
    $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./"  BaseRecalibrator \
    -R $ref  \
    -I ${id}_marked.bam  \
    --known-sites $snp \
    --known-sites $indel \
    -O ${id}_recal.table \
    1>${id}_log.recal 2>&1 

    $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./"  ApplyBQSR \
    -R $ref  \
    -I ${id}_marked.bam  \
    -bqsr ${id}_recal.table \
    -O ${id}_bqsr.bam \
    1>${id}_log.ApplyBQSR  2>&1 

    echo "end BQSR for ${id}" `date`
    fi
done

关于上面的两个工具的解释：

第一个，BaseRecalibrator，这里计算出了所有需要进行重校正的reads和特征值，然后把这些信息输出为一份校准表文件（recal.table）。--known-sites参数需要输入已知且可靠的变异位点，可以有多个。因为人与人之间基因组的差异其实是很小的，约0.1%，那么在一个群体中发现的已知变异，在某个人身上也很可能是同样存在的。所以，我们首先排除掉所有的已知变异位点，然后计算每个（报告出来的）质量值下面有多少个碱基在比对之后与参考基因组上的碱基是不同的，这些不同碱基就被我们认为是错误的碱基，它们的数目比例反映的就是真实的碱基错误率。

第二个，ApplyBQSR，这一步利用第一步得到的校准表文件（recal.table）重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新的BAM文件。

校正前后的碱基质量值：

Mutect2找变异

本来拿到了校正后的bam文件就可以走GATK的HaplotypeCaller流程，但是现在有了更好的工具，所以这里演示就跳过了HaplotypeCaller。因为要求需要有正常样本和肿瘤样本配对以便找出somatic mutations，这一步用到了GATK的Mutect2和FilterMutectCalls工具。

我们这里有normal和tumor，需要把配置文件config做成两列，如下：

$ cat config 

case1_germline  case1_biorep_A_techrep
case2_germline  case2_biorep_A
case3_germline  case3_biorep_A
......
case4_germline  case4_techrep_2
case5_germline  case5_techrep_2
case6_germline  case6_techrep_2

详细代码如下：

## mutect.sh
GATK=/home/jmzeng/biosoft/gatk4/gatk-4.0.6.0/gatk
ref=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
bed=/home/llwu/reference/index/human/CCDS/GRCh38.bed

cat config|while read id
do
    arr=(${id})
    sample=${arr[1]}
    T=../5.gatk/${arr[1]}_bqsr.bam
    N=../5.gatk/${arr[0]}_bqsr.bam
    echo "start Mutect2 for ${id}" `date`
    $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./"  Mutect2 -R $ref \
    -I ${T} -tumor $(basename "$T" _bqsr.bam) \
    -I ${N} -normal $(basename "$N" _bqsr.bam) \
    -L $bed  \
    -O ${sample}_mutect2.vcf

    $GATK  FilterMutectCalls -V ${sample}_mutect2.vcf -O ${sample}_somatic.vcf
    echo "end Mutect2 for ${id}" `date`
done

for i in *_somatic.vcf
do
j=$(basename "$i" _somatic.vcf )
echo $j
cat $i | perl -alne '{if(/^#/){print}else{next unless $F[6] eq "PASS";next if $F[0] =~/_/;print } }' > ${j}_filter.vcf
done

关于上面代码的解释：首先是Mutect2，主要是根据对正常-肿瘤样本进行位点比较寻找体细胞突变somatic mutations，如果没有正常样本，那么软件能正常使用，但假阳性会很高，而且分析流程也不一样。我们这里输入了两个配对的T-N样本    -I ${T} -tumor $(basename "$T" _bqsr.bam) 和-I ${N} -normal $(basename "$N" _bqsr.bam)。这样我们就拿到了每个tumor样本中的突变位点，${sample}_mutect2.vcf。但是这还不够，还要通过FilterMutectCalls过滤，过滤后就是可信度较高的somatic突变位点会被打上PASS标签，就是${sample}_somatic.vcf。最后我们再用一个简单的脚本把打上PASS标签的位点提取出来，得到了${j}_filter.vcf，如：

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  case1_biorep>
chr1    6146376 .       G       T       .       PASS    DP=184;ECNT=1;NLOD=22.56;N_A>
chr1    6461445 .       G       T       .       PASS    DP=21;ECNT=1;NLOD=2.70;N_ART>
chr1    31756671        .       C       A       .       PASS    DP=172;ECNT=1;NLOD=9>
chr1    32672798        .       A       T       .       PASS    DP=448;ECNT=1;NLOD=2>
chr1    39441098        .       G       T       .       PASS    DP=411;ECNT=1;NLOD=3>
......

到这里拿到了最后的vcf文件，我们就可以进一步分析，而GATK的流程也到这里差不多结束，希望你也能学会这个分析流程。下一节我们将会把vcf文件和bam文件关联在一起以进一步了解突变位点，同时会用ANNOVAR软件对突变位点进行注释，希望你不要错过。最后留一个彩蛋，公众号后台回复GATK，谢谢观看！我们下周再见~

■   ■   ■ 

全国巡讲约你

第一站-重庆  （已结束）

粤港澳大湾区专场 （已结束）

第二站-济南 （已结束）

千呼万唤进北京（已结束）

七月份我们不外出，只专注单细胞！

系统学习单细胞分析，报名生信技能树的线下培训，手慢无。

一年一度的生信技能树单细胞线下培训班（已经结束）

下一场由你指定！！