GATK4 流程分析- 从fastq到vcf（样例数据处理）

（详见https://zhuanlan.zhihu.com/p/69726572?utm_id=0）

参考基因组处理

参考基因组采用在ensembl下载的鸡的参考基因组（以下命名为1.fa）
以下为对参考基因组的预处理，三个步骤一个也不能少

大于2G的参考基因组要用参数-a bwtsw

---

$ bwa index -a bwtsw 1.fa

生成1.fa.fai文件

---

$ samtools faidx 1.fa

生成参考基因组的dict⽂件

---

$ java -jar picard.jar CreateSequenceDictionary R=1.fa O=1.dict

（此处预留）

序列比对

软件：bwa （mem方法）
若一个样本有4个fastq，此处分开处理fastq，将同⼀通道的两个fq文件⼀起处理，生成⼀个sam/bam⽂件,那样子⼀个样本共生成2个sam/bam⽂件，后面还要进行merge步骤。
若一个样本只有2个fastq，此处刚好生成一个样本一个bam文件。

-R 参数涉及到后续步骤，非常重要，如果设错了需要重新跑这一步。

$ time bwa mem -t 6 -M -Y -R ‘@RG\tID:2\tSM:19P0126636\tLB:WES\tPL:Illumina’ 1.fa 2_1.fq.gz 2_2.fq.gz | samtools view -Sb - > 19P0126636WES.bam && echo “**bwa mapping done **”

（fq.gz没有的话fq也是可以的）
##查看前⾯bwamem时候设定的头⽂文件

$ samtools view -H 19P0126636WES.bam | grep ‘@RG’

参数option：
-R 设定头文件 必填，ID：通道名或样本名,通过这个信息分组，必须唯一（但是可以自定义，自定义必须唯一）；SM：样本名；LB：⽂库名；PL：测序平台信息[COMPLETE,ILLUMINA,SANGER]。以上这些信息后续GATK和markduplicate会用到，不可出错
-t 核数
-M :-M 将 shorter split hits 标记为次优，以兼容Picard’s markDuplicates 软件 关于alignment，
由于比对算法的区别，DNA数据一般用bwa和bowtie2，RNA一般用STAR，如果用bwa，R的设定一定（一定！）要对好，ID和Sample我一般直接用${prefix}
(prefix是一直设定的样本前缀，一般为样本名）来运行，避免出错，如果这里没有设定好，一路运行下去也不会报错，就是出来的vcf没有任何位点，让人无法debug…

sam/bam文件预处理

软件：samtools，picard（基于java）

#如果前面一个样本有4个fastq，此处需要先将bam文件合并,也可以指定合并区域

---

$ time samtools merge -@ 6 -h sample_1.bam output.bam sample_1.bam sample_2.bam

---

#-h FILE 指定FILE内的’@’头复制到输出bam文件中并替换输出文件的文件头。否则，输出文件的文件头将从第一个输入文件复制过来；

#Sortsam,非常占用内存和IO读写，耗时,-@ 核数

---

$ time samtools sort -@ 30 -o 19P0126636WES.sorted.bam 19P0126636WES.bam

#markduplicate，占用内存，耗时
#error： open too much file 打开太多文件； 解决方法 $ ulimit -n 10240
#error: GC-overhead limit exceeded 内存不足；解决方法，运行是加参数 -Xmx64g，64g是指64g内存，请按照自己电脑情况修改数字,笔记本可以删去
#因为这个数据集之前的duplicate并不过关，所以加了一个设定REMOVE_DUPLICATES=true，默认只是标记而不删除
#SE可以用samtools或picard去重，PE一般用picard，效果好。

---

$ time java -jar picard.jar.1 MarkDuplicates I=19P0126636WES.sorted.bam O=19P0126636WES.sorted.markdup.bam M=19P0126636WES.sorted.markdup.txt REMOVE_DUPLICATES=true

---

#index
$ time java -jar picard.jar.1 BuildBamIndex -Xmx64g I=19P0126636WES.sorted.markdup.bam

---

##检查一下生成的bam文件质量是否合格（optional）
$ time samtools flagstat -@ 30 19P0126636WES.sorted.markdup.bam > 19P0126636WES.sorted.markdup.stat

---

（非必须）
$ mkdir flagstat
$ multiqc flagstat/
$ cat 19P0126636WES.sorted.markdup.stat

生成vcf文件

单样本

单样本到这一步就可以生成VCF了

使用HaplotypeCaller对上述的bam文件进行call variant，就是寻找变异的过程，生成一个VCF文件用于后续位点的质控和注释,从一个bam到一个vcf文件

$ time gatk HaplotypeCaller -R 1.fa -I 19P0126636WES.sorted.markdup.BQSR.bam -O 19P0126636WES.HC.vcf

之后就可以用不同的方法过滤掉置信度不高的位点，然后进行注释，注释过程见下文。

多样本

单独为每个样本生成后续分析所需的中间文件——gVCF文件。这一步中包含了对原始fastq数据的质控、比对、排序、标记重复序列、BQSR和HaplotypeCaller gVCF等过程。这些过程全部都适合在单样本维度下独立完成。值得注意的是，与单样本模式不同，该模式中每个样本的gVCF应该成为这类流程的标配，在后续的步骤中我们可以通过gVCF很方便地完成群体的Joint Calling；
第二步，依据第一步完成的gVCF对这个群体进行Joint Calling，从而得到这个群体的变异结果和每个人准确的基因型（Genotype），最后使用VQSR完成变异的质控。这两个步骤其实还包含了许多细节，具体可见我在流程中的注释。

要生成gvcf前期步骤与单样品分析步骤一致，只是在最后HaplotypeCaller有差别，多了-ERC GVCF 参数

$ time gatk HaplotypeCaller -R /home/zhn/5/1.fa -I 19P0126636WES.sorted.markdup.bam -ERC GVCF -O 19P0126636WES2.g.vcf

多个gvcf整合(optional)

#CombineGVCFs：旧方法，速度慢，但是可以一次全部合并（合并不同样本的文件）
$ gatk CombineGVCFs
-R /path/to/hg38/hg38.fa
-V sample1.g.vcf
-V sample2.g.vcf
-O combined.g.vcf

---

#样本较多时，按如下操作
$ find gvcf/ -name “*.g.vcf” >input.list
$ nohup gatk -Xmx20g CombineGVCFs -R /data/all_data/ref/hg38/hg38.fa --variant
$ input.list -o combined.g.vcf & #Xmx20g 使⽤用20G内存

---

#joint genotyping
$ gatk GenotypeGVCFs
-R /path/to/hg38/hg38.fa
-V combined.g.vcf
-G StandardAnnotation
-O raw_variants.vcf (这个就是后续命令行中的19P0126636WES.HC.vcf，VQSR的输入文件）

变异质控
VQSR中参考的指标阈值有6个，分别是：
QualByDepth（QD）
FisherStrand (FS)
StrandOddsRatio (SOR)
RMSMappingQuality (MQ)
MappingQualityRankSumTest (MQRankSum)
ReadPosRankSumTest (ReadPosRankSum)
（详见https://blog.csdn.net/weixin_33378570/article/details/112818035
以及https://blog.csdn.net/Gossie/article/details/109320960）

参考程序：

gatk VariantFiltration -R 1.fa\
    -V all.vcf \
    -filter "QD > 2.0" --filter-name "QD2" \
    -filter "FS > 60.0" --filter-name "FS60" \
    -filter "MQ > 40.0" --filter-name "MQ40" \
    -filter "MQRankSum > -12.5" --filter-name "MQRankSum-12.5" \
    -filter "ReadPosRankSum > -8.0" --filter-name "ReadPosRankSum-8" \
    -O all.gatk.vcf