GATK4的gvcf流程
走GVCF肯定是多个样本，比如我这里有50个病人的正常组织及肿瘤组织的WES测序数据。
得到了它们的bam文件，也是走的GATK流程，这里就不多说了。本教程首发于生信技能树VIP论坛：https://vip.biotrainee.com/d/423-gatk4-gvcf
配置GATK运行环境
参考我前面在生信菜鸟团博客分享的： https://vip.biotrainee.com/d/384-gatk4
GATK=/home/jianmingzeng/biosoft/GATK/gatk-4.0.3.0/gatk
bed=/home/jianmingzeng/annotation/CCDS/human/exon_probe.GRCh38.gene.150bp.bed
module load java/1.8.0_91
GENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
INDEX=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/bwa_index/gatk_hg38
DBSNP=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz
kgSNP=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
kgINDEL=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
其中基于hg38版本参考基因组的外显子坐标的制作方式我还是要强调一下，下载文件 CCDS.20160908.txt可以使用下面的代码：
cat CCDS.20160908.txt |perl -alne '{/\[(.*?)\]/;next unless $1;$gene=$F[2];$exons=$1;$exons=~s/\s//g;$exons=~s/-/\t/g;print "$F[0]\t$_\t$gene" foreach split/,/,$exons;}'|sort -u |bedtools sort -i >exon_probe.hg38.gene.bed
awk '{print "chr"$0"\t0\t+"}' exon_probe.hg38.gene.bed > exon_probe.GRCh38.gene.bed
awk '{print $1"\t"$2-150"\t"$3+150"\t"$4"\t"$5"\t"$6}' exon_probe.GRCh38.gene.bed > exon_probe.GRCh38.gene.150bp.bed
首先得到recal的bam文件
假设大家已经把fastq文件比对到参考基因组并且去除PCR重复了，然后：
bam=${sample}_marked_fixed.bam
time $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" BaseRecalibrator \
-I $bam -R $GENOME --output ${sample}_recal.table --known-sites $kgSNP --known-sites $kgINDEL
time $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" ApplyBQSR \
-I $bam -R $GENOME --output ${sample}_recal.bam -bqsr ${sample}_recal.table
rm ${sample}_marked_fixed.bam
然后批量对每个样本做gvcf操作
因为是外显子测序所以我才指定bed文件，其实也可以省略，我只是为了加快速度和节省空间而已! 一个WES数据大概一到两个小时运行时间吧！
for bam in /home/jianmingzeng/data/alignment/*_recal.bam
do
if((i%$1==$2))
then
sample=$(basename "$bam" _recal.bam)
echo $sample
$GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller \
-ERC GVCF -L $bed -R $GENOME -I $bam --dbsnp $DBSNP -O ${sample}_raw.vcf
fi
i=$((i+1))
done
最后把多个gvcf文件合并
因为合并只能一次给定一个区间，所以需要分染色体来做，正好相当于是并行！！！
for bed in chr{1..22} chrX chrY
do
if((i%$1==$2))
then
echo $bed
time $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" GenomicsDBImport \
-L $bed -R $GENOME \
-V 1NT-1_raw.vcf \
-V 1NT-2_raw.vcf \
--genomicsdb-workspace-path gvcfs_${bed}.db
time $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" GenotypeGVCFs \
-R $GENOME -V gendb://gvcfs_${bed}.db -O final_${bed}.vcf
fi
i=$((i+1))
done
有多少个样本，就需要在上面的命令里面添加多少个vcf文件地址。
得到的文件如下：
12M Apr 19 13:37 final_chr10.vcf
15M Apr 19 13:38 final_chr11.vcf
14M Apr 19 13:38 final_chr12.vcf
5.4M Apr 19 14:17 final_chr13.vcf
8.2M Apr 19 14:06 final_chr14.vcf
9.0M Apr 19 14:15 final_chr15.vcf
11M Apr 19 14:15 final_chr16.vcf
14M Apr 19 15:03 final_chr17.vcf
4.4M Apr 19 14:22 final_chr18.vcf
18M Apr 19 15:04 final_chr19.vcf
27M Apr 19 12:40 final_chr1.vcf
6.3M Apr 19 14:38 final_chr20.vcf
3.5M Apr 19 15:13 final_chr21.vcf
5.7M Apr 19 14:41 final_chr22.vcf
19M Apr 19 12:17 final_chr2.vcf
14M Apr 19 12:06 final_chr3.vcf
11M Apr 19 12:19 final_chr4.vcf
12M Apr 19 13:21 final_chr5.vcf
13M Apr 19 12:59 final_chr6.vcf
14M Apr 19 12:47 final_chr7.vcf
8.9M Apr 19 12:49 final_chr8.vcf
12M Apr 19 13:59 final_chr9.vcf
6.3M Apr 19 15:33 final_chrX.vcf
202K Apr 19 14:40 final_chrY.vcf
因为是区分了染色体，所以是需要把这些vcf文件合并，下次再讲咯。
生信技能树GATK4系列教程