使用GATK合并多个vcf文件

这不是简单的坐标合并，多个样本的vcf文件里面都只会记录对应样本的变异位点信息，但是因为样本检测手段不一致，可能是WGS或者WES，或者不会测序，意味着基因组某些区域可能是根本就没有被测序到，也就无从得知对应位点的碱基信息了。

这个时候的合并，要从bam文件开始，同一个坐标在某个样本既可能是野生型，杂合或者纯合突变，也有可能是该位点并没有任何reads覆盖。

如果不想自己写脚本去探究，那么gatk的HaplotypeCaller的GVCF模式正好能满足要求

首先输出gvcf文件

代码如下，把所有的样本的bam文件都走一波gatk的HaplotypeCaller的GVCF模式，得到后缀为_raw.g.vcf的文件（这个文件非常占硬盘空间哦）

GENOME=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta
TMPDIR=/home/jianmingzeng/tmp/software
GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar
DBSNP=/home/jianmingzeng/annotation/variation/human/dbSNP/All_20160601.vcf.gz

ls /home/jianmingzeng/data/project/wes/bamFiles/*bam |while read id 
do 
    file=$(basename $id )
    sample=${file%%.*}
    echo $file $sample
    start=$(date +%s.%N)
    echo HaplotypeCaller `date`
    java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T HaplotypeCaller  \
    -R $GENOME -I $id --dbsnp $DBSNP -ERC GVCF -gt_mode DISCOVERY \
    -stand_emit_conf 10 -o  ${sample}_raw.g.vcf
    echo HaplotypeCaller `date`
    dur=$(echo "$(date +%s.%N) - $start" | bc)
    printf "Execution time for HaplotypeCaller : %.6f seconds" $dur
    echo
done 

理论上只要是bam文件里面表示该样本的基因组上面 某个位点被覆盖过，就会输出该位点的信息，无论其是否是突变。

这个输出的gvcf文件格式并不需要解释，也不需要理解，反正就是中间文件，当然，也欢迎有求知欲的同学继续深入了解哈。

合并多个样本的gvcf信息

因为每个样本的每个位置信息都被记录，所以这个时候的合并就很方便了，同样还是用gatk工具吧，代码如下：

GENOME=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta
TMPDIR=/home/jianmingzeng/tmp/software
GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar

java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T GenotypeGVCFs -nt 5  \
-R $GENOME   \
--variant  sample1.g.vcf   \
--variant  sample2.g.vcf   \
--variant  sample3.g.vcf   \
--variant  sample4.g.vcf \
-o output.vcf

另一种方法

以前我在直播我的基因组里面提到过，我的基因组是5条lane的独立fastq数据，期初我是先分开比对，然后把bam文件merge起来，结果发现自己在找变异的时候输出的vcf文件里面，每个lane都给出了基因型信息，也就是说根本就没有把这些lane当成是同一个样本。按照这个思路，我们可以把不同样本的bam文件合并，然后直接找变异，只有我们没有把样本信息抹掉，就仍然是独立样本独立基因型。

如果真正要合并不同的lane的数据为一个样本，需要用AddOrReplaceReadGroups来抹掉lane信息。

ls jmzeng_*.bam > files.bamlist
samtools merge -@ 5  -b files.bamlist  merged.bam
samtools index merged.bam
### AddOrReplaceReadGroups ###
java -Djava.io.tmpdir=$TMPDIR    -Xmx40g -jar $PICARD AddOrReplaceReadGroups \
INPUT=${sample}.bam OUTPUT=${sample}_tmp.bam   RGID=jmzeng  RGLB=lib_all  RGPL=illumina RGPU=x10  RGSM=jmzeng
mv ${sample}_tmp.bam ${sample}.bam
samtools index ${sample}.bam 

当然，vcf文件只是记录变异而已，变异一定要进行注释：

ls  /home/jianmingzeng/data/project/wes/variation/*_filtered_PASS.*.vcf |while read id
do
file=$(basename $id )
sample=${file%.*}
echo $sample
java -jar ~/biosoft/SnpEff/snpEff/snpEff.jar -i vcf GRCh37.75 $id >snpEFF_output/${sample}.snpEff.vcf
mv snpEff_summary.html ${sample}.snpEff_summary.html
~/biosoft/ANNOVAR/annovar/convert2annovar.pl -format vcf4old  $id  >${sample}.annovar_input
~/biosoft/ANNOVAR/annovar/annotate_variation.pl -buildver hg19 --geneanno --outfile annovar_output/${sample}_anno ${sample}.annovar_input  ~/biosoft/ANNOVAR/annovar/humandb/
done

更多gatk教程见：
一个全基因组重测序分析实战
GATK best practice每个步骤耗时如何？
GATK best practice每个步骤都是必须的吗?