需求描述
一个很简单的需求：一批水稻材料的芯片数据（位点少），想看看它们在3K Rice中处于何种亚群和位置。就需要将芯片位点与3K RG位点整合后进行分析。
已知3K Rice位点可从SNP-Seek中下载：https://snp-seek.irri.org/_download.zul;jsessionid=F2B11FD2C5BC6A9AA07D9FE198915C9E
但它是二进制Plink 格式，本身没有包含等位基因的信息：
image.png
因此当转化vcf格式时，REF和ALT是按照等位基因频率来进行分配的，部分位点的REF和ALT可能发生调换，即便是使用--keep-allele-order参数也是无用的。
这也是为什么我建议大家简化vcf信息时，不要使用plink和vcf来回转化的方式，而是直接使用bcftools，详见往期推文：如何快速简化vcf信息？
当两个vcf文件合并，位置（CHR+POS）相同时，REF和ALT可能相反，合并必然会出错。
你当然可以自己写脚本找出这些错误的位点，通过比较REF和ALT，纠正错误（可能需要用到参考基因组，找到对应位置真实的REF）。但当vcf文件很大时，处理效率还是比较低的。我的建议是用现成的工具。
$ bcftools view -R array_chr_pos.txt 3k29mio.vcf.gz -Oz -o 3k.overlap.vcf.gz
## 3K位点由于是plink1.9转化而来，ref和alt可能发生了调换，需要纠正过来方能合并。比如下面的错误是由于3k.overlap.vcf.gz中的1:715297位点REF和ALT对应的G和A 反了。
$ bcftools merge array.vcf.gz 3k.overlap.vcf.gz -Oz -o merge.vcf.gz
The REF prefixes differ: G vs A (1,1)
Failed to merge alleles at 1:715297 in 3k.overlap.vcf.gz
尝试解决
基本思路是引入参考基因组，找出REF不对应的点再纠正。
一开始尝试了bcftools的插件fixref。
$ bcftools +fixref 3k.overlap.vcf.gz -Oz -o 3k.overlap.fixref.vcf.gz -- -f msu7.fa -m top
日志显示，一些位点确实发生了调换：
# SC, guessed strand convention
SC TOP-compatible 0
SC BOT-compatible 0
# ST, substitution types
ST A>C 417 4.5%
ST A>G 1389 14.9%
ST A>T 453 4.9%
ST C>A 482 5.2%
ST C>G 310 3.3%
ST C>T 1645 17.6%
ST G>A 1538 16.5%
ST G>C 318 3.4%
ST G>T 468 5.0%
ST T>A 411 4.4%
ST T>C 1490 16.0%
ST T>G 413 4.4%
# NS, Number of sites:
NS total 9334
NS ref match 5103 54.7%
NS ref mismatch 4231 45.3%
NS flipped 334 3.6%
NS swapped 2702 28.9%
NS flip+swap 376 4.0%
NS unresolved 2817 30.2%
NS fixed pos 0 0.0%
NS skipped 0
NS non-ACGT 0
NS non-SNP 0
NS non-biallelic 0
以为解决了问题，但当我合并时，还是报错有新的位点没有纠正过来。
$ bcftools merge array.vcf.gz 3k.overlap.fixref.vcf.gz -Oz -o merge.vcf.gz
The REF prefixes differ: C vs T (1,1)
Failed to merge alleles at 1:1539076 in 3k.overlap.fixref.vcf.gz
具体原因我不知，但官方说明是不要轻易用fixref，可能会产生无意义的基因型：https://samtools.github.io/bcftools/howtos/plugin.fixref.html
image.png
正确解决
使用bcftools norm来解决，它通常用于对VCF文件进行规范化处理，包括拆分多等位位点、合并相邻位点等。
$ bcftools norm --check-ref -s -f msu7.fa 3k.overlap.vcf.gz -Oz -o 3k.overlap.fixref.vcf.gz
Lines total/split/realigned/skipped: 9334/0/0/0
REF/ALT total/modified/added: 9334/4231/0
全部纠正了，可以正确合并：
$ tabix 3k.overlap.fixref.vcf.gz
$ bcftools merge array.vcf.gz 3k.overlap.fixref.vcf.gz -Oz -o merge.vcf.gz
如上图所示，官方也提醒：不要轻易使用bcftools norm
参考： https://www.biostars.org/p/440506/；https://www.biostars.org/p/248685/；https://www.biostars.org/p/440506/ ；https://www.biostars.org/p/411202/ ；https://cloud.tencent.com/developer/ask/sof/1136765 ；https://github.com/cumc/xqtl-pipeline/issues/207 ；https://github.com/samtools/bcftools/issues/875 ；https://www.biostars.org/p/336399/
网友评论