目录
需求描述
一个很简单的需求:一批水稻材料的芯片数据(位点少),想看看它们在3K Rice中处于何种亚群和位置。就需要将芯片位点与3K RG位点整合后进行分析。
已知3K Rice位点可从SNP-Seek中下载:https://snp-seek.irri.org/_download.zul;jsessionid=F2B11FD2C5BC6A9AA07D9FE198915C9E
但它是二进制Plink 格式,本身没有包含等位基因的信息:
因此当转化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
正确解决
使用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/
更多信息请关注微信公众号:生物信息与育种