需求描述

一个很简单的需求：一批水稻材料的芯片数据（位点少），想看看它们在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/