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系列教程

我不想赚你的钱，不行吗？ （推荐阅读）