如何快速纠正VCF文件中REF和ALT的位置错误?

需求描述

一个很简单的需求:一批水稻材料的芯片数据(位点少),想看看它们在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/1136765https://github.com/cumc/xqtl-pipeline/issues/207https://github.com/samtools/bcftools/issues/875https://www.biostars.org/p/336399/

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,904评论 6 497
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,581评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 160,527评论 0 350
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,463评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,546评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,572评论 1 293
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,582评论 3 414
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,330评论 0 270
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,776评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,087评论 2 330
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,257评论 1 344
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,923评论 5 338
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,571评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,192评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,436评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,145评论 2 366
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,127评论 2 352

推荐阅读更多精彩内容