使用GATK合并比较多个vcf文件

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

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

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

首先输出gvcf文件

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

  1. GENOME=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta

  2. TMPDIR=/home/jianmingzeng/tmp/software

  3. GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar

  4. DBSNP=/home/jianmingzeng/annotation/variation/human/dbSNP/All_20160601.vcf.gz

  5. ls /home/jianmingzeng/data/project/wes/bamFiles/*bam |while read id

  6. do

  7.    file=$(basename $id )

  8.    sample=${file%%.*}

  9.    echo $file $sample

  10.    start=$(date +%s.%N)

  11.    echo HaplotypeCaller `date`

  12.    java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T HaplotypeCaller  \

  13.    -R $GENOME -I $id --dbsnp $DBSNP -ERC GVCF -gt_mode DISCOVERY \

  14.    -stand_emit_conf 10 -o  ${sample}_raw.g.vcf

  15.    echo HaplotypeCaller `date`

  16.    dur=$(echo "$(date +%s.%N) - $start" | bc)

  17.    printf "Execution time for HaplotypeCaller : %.6f seconds" $dur

  18.    echo

  19. done

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

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

合并多个样本的gvcf信息

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

  1. GENOME=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta

  2. TMPDIR=/home/jianmingzeng/tmp/software

  3. GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar

  4. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T GenotypeGVCFs -nt 5  \

  5. -R $GENOME   \

  6. --variant  sample1.g.vcf   \

  7. --variant  sample2.g.vcf   \

  8. --variant  sample3.g.vcf   \

  9. --variant  sample4.g.vcf \

  10. -o output.vcf

另一种方法

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

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

  1. ls jmzeng_*.bam > files.bamlist

  2. samtools merge -@ 5  -b files.bamlist  merged.bam

  3. samtools index merged.bam

  4. ### AddOrReplaceReadGroups ###

  5. java -Djava.io.tmpdir=$TMPDIR    -Xmx40g -jar $PICARD AddOrReplaceReadGroups \

  6. INPUT=${sample}.bam OUTPUT=${sample}_tmp.bam   RGID=jmzeng  RGLB=lib_all  RGPL=illumina RGPU=x10  RGSM=jmzeng

  7. mv ${sample}_tmp.bam ${sample}.bam

  8. samtools index ${sample}.bam

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

  1. ls  /home/jianmingzeng/data/project/wes/variation/*_filtered_PASS.*.vcf |while read id

  2. do

  3. file=$(basename $id )

  4. sample=${file%.*}

  5. echo $sample

  6. java -jar ~/biosoft/SnpEff/snpEff/snpEff.jar -i vcf GRCh37.75 $id >snpEFF_output/${sample}.snpEff.vcf

  7. mv snpEff_summary.html ${sample}.snpEff_summary.html

  8. ~/biosoft/ANNOVAR/annovar/convert2annovar.pl -format vcf4old  $id  >${sample}.annovar_input

  9. ~/biosoft/ANNOVAR/annovar/annotate_variation.pl -buildver hg19 --geneanno --outfile annovar_output/${sample}_anno ${sample}.annovar_input  ~/biosoft/ANNOVAR/annovar/humandb/

  10. done

更多gatk教程见：

一个全基因组重测序分析实战

GATK best practice每个步骤耗时如何？

GATK best practice每个步骤都是必须的吗？

本文只是合并了，至于如何比较，请听下回分解！

赞 (0)

相关推荐