首页 > 其他分享 >如何快速纠正VCF文件中REF和ALT的位置错误?

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

时间:2023-11-11 16:47:17浏览次数:40  
标签:VCF vcf gz overlap ST ALT NS REF

目录

需求描述

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

更多信息请关注微信公众号:生物信息与育种

标签:VCF,vcf,gz,overlap,ST,ALT,NS,REF
From: https://www.cnblogs.com/miyuanbiotech/p/17826042.html

相关文章

  • Altium Designer自学笔记
    本次使用AD20为基础进行练习。1.1新建工程包括:原理图、PCB、原理图库、PCB库。 1.2新建元器件 点击右下角的“Panels"面板,调出新建元器件界面1.3视图---->栅格------->切换捕捉栅格右边DesignerltemID可修改器件名称,绘制状态下Tap键可暂停修改。 11.4复制元器件按......
  • CF1485F Copy or Prefix Sum 题解
    思路考虑\(a_i\)要么是\(b_i\)要么是\(b_i-s\)。考虑\(s\)代表着什么。它是\(a\)的前缀和。那么必然是往前一段\(b\)的和。因为每个\(b\)代表着要么是这一位的\(a\)或者前面所有的\(a\)。考虑设\(f_i\)为这个位置填\(b_i\)的方案数。\(g_i\)为这个......
  • [20231109]bash shell快捷键alt+number的问题.txt
    [20231109]bashshell快捷键alt+number的问题.txt--//前一阵子,我想实现12行合并1行的输出,理论讲要使用paste命令加入12个-.输入命令时候要数输入了多少-.我知道bashshell有一--//个快捷键alt+number可以产生连续输入某个字符,但是我一直不知道如何关掉这个功能.有时候误触发这......
  • [BalticOI 2019 Day2] 汤姆的餐厅
    [BalticOI2019Day2]汤姆的餐厅题目背景译自BalticOI2019Day2T1.Tom'sKitchen题目描述Tom'sKitchen是一家非常受欢迎的餐厅,其受欢迎的原因之一是每份菜都由至少$K$名厨师进行准备。今天有$N$份菜需要准备,第$i$份菜的准备时间是$A_i$小时。Tom可以......
  • 以下代码执行后,输出结果为 抛出ReferenceError
    letx=10;letfoo=()=>{console.log(x);letx=20;x++;}foo();使用let声明的变量,既不会发生变量提升,同时又存在“暂时性死区”,所以在块级作用域内,如果使用let声明一个变量,那么该变量在声明之前是不可用的,否则会抛出ReferenceError异常一楼的回答说:”l......
  • 【Qt初入江湖】Qt QSqlRelationalTableModel 底层架构、原理详细描述
    鱼弦:内容合伙人、新星导师、全栈领域创作新星创作者、51CTO(Top红人+专家博主)、github开源爱好者(go-zero源码二次开发、游戏后端架构https://github.com/Peakchen) QtQSqlRelationalTableModel是Qt中用于实现具有关联表格的模型类,它继承自QSqlTableModel。QSqlRelationalTable......
  • RLHF · PBRL | PEBBLE:通过 human preference 学习 reward model
    论文题目:PEBBLE:Feedback-EfficientInteractiveReinforcementLearningviaRelabelingExperienceandUnsupervisedPre-training,貌似是ICML2021的文章。本博客为论文阅读笔记,【不能代替】阅读原文的工作量。原文写的也很好,是AI顶会的风格,相对容易读懂。阅读材料:p......
  • salt自定义模块内使用日志例子
    如果你想要在你的SaltMinion中使用自定义的Salt模块并且记录日志,你可以创建一个自定义Salt模块,并在模块中使用Python的标准`logging`库来记录日志。以下是一个示例:首先,在SaltMaster上创建一个自定义模块的目录,例如`/srv/salt/_modules/`。然后在该目录中创建一个Python文件,例......
  • DataGrip连接MySql数据库失败:dataGrip java.net.ConnectException: Connection refuse
    1.问题报错:dataGripjava.net.ConnectException:Connectionrefused:connect.详细错误:[08S01]CommunicationslinkfailureThelastpacketsentsuccessfullytotheserverwas0millisecondsago.Thedriverhasnotreceivedanypacketsfromtheserver.Communica......
  • ALLEGRO导网表报错This reference has already been assigned to a different package
     (1)QUESTION(ORCAP-1589):Nethastwoormorealiases-possibleshort?原因:器件默认管脚命名(NET名称)与所连接网络的NET名称不一致导致的措施:可忽略。或关闭Tools->DesignRulesCheck->PhysicalRules->Checkpowergroundshort(2)ReportforInvalidReferencesERROR(ORCAP-......