annovar的avinput转为vcf

最近遇到问题,要将annovar的结果转为vcf,然后用这个vcf对比另一个vcf,找二者共有的变异。其实就是一个bedtools intersect的事,但是annovar的start,end,Ref,Alt和vcf的pos,Ref,Alt没法直接对比,vcf中的Ref,Alt是不能出现横杠(-)的。
在biostars上有人回答了这个问题Convert annovar file to vcf file (biostars.org),应该说除了没给代码,该有的都有了。
本着能不写代码就不写代码的原则,github上刚好有这这个项目,虽然只有一个星。😂(用annovar注释结果做vcf取交集这种操作,用annovar直接就输出vcf就完事了,也就我这种拿别家分析结果搞事情,才有这种需求)
KijinKims/Annovar2VCF: table annovar format into VCF format (github.com)
下载这个项目

git clone https://github.com/KijinKims/Annovar2VCF.git

由于这个脚本实际调用了annovar附带的retrieve_seq_from_fasta.pl,根据输入的参考基因组文件,去获取indel位置的碱基。因此,如果你存放annovar相关脚本的目录,没有加入环境变量中。那么就需要指定这个脚本的绝对路径。

cd Annovar2VCF/
vi Annovar2VCF.py  #编辑此脚本
#找到
retrieve_args=["retrieve_seq_from_fasta.pl", "-format",...]
#将retrieve_seq_from_fasta.pl,修改为绝对路径
#例如我
retrieve_args=["/Data2/soft/annovar/retrieve_seq_from_fasta.pl", "-format",...]

修改完成后,需要安装这个脚本的依赖包numpy,pandas注意,这个脚本是python2写的,由于我个人常年用conda的python3,这几个包都是家中常备。几乎没用过或linux自带的python2,因此我这里直接将此脚本转为了python3脚本。

which python #查看我当前python的路径
#例如我,返回 /home/tanqiang/miniconda3/bin/python
#有python的地方就会有2to3,因此2to3也在/home/tanqiang/miniconda3/bin/里
2to3 -w Annovar2VCF.py #2to3会直接在原始文件上修改,但会生成一个xxx.py.bak的原始文件
#修改后的文件有个空格和tab的缩进小问题

完成这些后,可以运行示例试试,如果成功运行,说明没有问题啦

python Annovar2VCF.py --input sample.txt  --output sample.vcf -r /Data2/ref_tq/ucsc.hg19.fasta

查看vcf可以看到,表头(姑且这么叫)没有FORMAT列,不能作为常规的vcf文件,如果有其它数据构造需求,就得自己动手了。

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    Func.refGene    Gene.refGene    GeneDetail.refGene      ExonicFunc.refGene      AAChange.refGene

如果用这样生成的vcf文件,作为bedtools的输入

bedtools intersect -a sample.vcf  -b demo_tar.vcf 

会报错,因为bedtools必须识别到vcf的header,即首行需要添加

##fileformat=VCF

Trouble auto-detecting certaain VCF files · Issue #131 · arq5x/bedtools2 (github.com)

©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容