这不是简单的坐标合并，多个样本的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每个步骤都是必须的吗?
网友评论