如何合并多个体细胞突变检测工具的输出结果

简介

在肿瘤基因组数据分析方法中，很多文章会使用多个软件 call 突变，然后进行合并，常用见的方法就是多个软件得到的结果中，一个突变在两个或两个以上的软件检测到即保留，认为其可信度较高。合并策略通常有两种：取交集（只保留两个工具都检测到的变异）或取并集（保留任一工具检测到的变异）。但是不同的突变检测工具，得到的 vcf 文件并不完全一致，因为不同工具定义的 format tag 不太一致。所以直接进行合并可能导致格式有问题。这里简单简介一种比较粗暴的vcf 取交集的方法。主要是针对体细胞突变检测工具 Mutect2 和 Strelka2 的结果。

vcf 文件格式

既然是要对 vcf 文件进行合并，就要深入了解 vcf 文件格式。vcf(Variant Call Format)：保存基因组变异数据的文本文件，分为头部区和主体区。

头部区：以 ## 开头的注释信息，注明了vcf的版本，生成的命令，以及对主体区各部分信息的解释。

主体区：除了列名之外，每一行即为一个突变信息。

一个实际的 vcf 文件（体细胞突变）主体区的截图如下：

对以上信息的每一列解释为：

如果想对vcf文件有一个更加深入的了解，可以查看官方帮助文档：https://samtools.github.io/hts-specs/VCFv4.2.pdf

Mutect2 的 vcf 文件

Mutect2 的体细胞突变检测方法在GATK 的Somatic Mutation流程--肿瘤基因组测序数据分析专栏 已经有详细介绍了。如果不看 vcf 的 header 的话，主要突变信息的内容是：

1. less -S 6.mutect/case1_biorep_B_filter.vcf |grep -v '##'|less -S

Strelka2 的 vcf 文件

Strelka2 的体细胞突变检测方法在Strelka2 call Somatic 流程--肿瘤基因组测序数据分析专栏 已经有详细介绍了。其体细胞突变检测的结果分为 SNV 和 INDELs 两个 vcf 文件，需要先进行合并：

1. zcat 6.strelka/case1_biorep_B/results/variants/somatic.snvs.vcf.gz | \
2. grep -v '^#'>6.strelka/case1_biorep_B_filter.vcf
4. zcat 6.strelka/case1_biorep_B/results/variants/*gz | grep -v '^#' | \
5. grep 'PASS' |sort -V >>6.strelka/case1_biorep_B_filter.vcf
7. less -S 6.strelka/case1_biorep_B_filter.vcf |grep -v '##' |less -S

合并两个 vcf 文件

如果突变检测只用了 Mutect2 和 Strelka2 ，可以用下面代码进行合并。首先使用 bcftools index 对文件进行构建索引：

1. bgzip -c 6.mutect/case1_biorep_B_filter.vcf \
2. >6.mutect/case1_biorep_B_filter.vcf.gz
4. bcftools index 6.mutect/case1_biorep_B_filter.vcf.gz
6. bgzip -c 6.strelka/case1_biorep_B_filter.vcf \
7. >6.strelka/case1_biorep_B_filter.vcf.gz
9. bcftools index 6.strelka/case1_biorep_B_filter.vcf.gz

然后使用 bcftools isec 进行取交集：

1. bcftools isec -p 6.merge-n=2-w1 \
2. 6.mutect/case1_biorep_B_filter.vcf.gz \
3. 6.strelka/case1_biorep_B_filter.vcf.gz
5. mv 6.merge/0000.vcf6.merge/case1_biorep_B_filter.vcf.gz

这样最后得到的 6.merge/case1_biorep_B_filter.vcf就是 Mutect2 和 Strelka2 的 overlap 了。可以看到合并后的 vcf 文件 header 部分有多出来 bcftools 的处理信息。

这种方法比较简单粗暴，实际上也有其他工具可以实现，欢迎大家在评论区提出自己的方法。