「生信练习题」从SnpEff注释得到的VCF中过滤4DTV位点

在2019年的一篇NG文章"Resequencing of 414 cultivated and wild watermelon accessions identifies selection for fruit quality traits"的方法部分，作者提到他们使用位于四倍兼并位点( fourfold degenerated sites)上的SNP构建系统发育树。原本的19,725,853个SNP经过这一过滤，最终只剩下89,914 个SNP。

Phylogenetic and population analyses. Bi-allelic SNPs with a missing data rate less than 15% and a minor allele count greater than three were kept for population genomic analyses. Additionally, only SNPs at fourfold degenerated sites (89,914 SNPs) were used to construct a neighbor-joining phylogenetic tree using MEGA7 with 500 bootstraps61.

所谓四倍兼并位点(4DTV), 指无论怎么突变，最终氨基酸都一样的位点，例如"TCA/TCT/TCC/TCG"中的第三位就是一个4DTV，只要在这个位点的碱基，无论怎么突变最终都是丝氨酸(Ser)。

我稍微检索了一下，发现并没有相应的脚本，于是就有了这道练习题。

给定一个用SnpEff注释的VCF文件，请过滤出其中的位于四倍兼并位点上的SNP。

我花了一个晚上时间写了代码，使用方法为python calc_4dTv_in_eff_vcf.py input.vcf output.vcf ref.fa，三个参数分别是输入的vcf，输出的vcf文件名，以及参考基因组序列。

需要注意的是，通常SnpEff会对一个位点注释多个信息，而我只选择其中第一条注释信息进行过滤。

#!/usr/bin/env python3from sys import argv
from pysam import VariantFile
from pysam import FastaFilefile_in = argv[1]
file_out = argv[2]
fafile = argv[3]codon = set(["TC", "CT", "CC", "CG", "AC", "GT", "GC", "GG"])
rev_dict = dict(A='T',T='A', C='G', G='C')bcf_in = VariantFile(file_in)
bcf_out = VariantFile(file_out, "w", header = bcf_in.header)
fa_in = FastaFile(fafile)for rec in bcf_in.fetch():ann = rec.info['ANN']info = rec.info['ANN'][0].split('|')# only use synonymouse variantsif info[1] != "synonymous_variant":continue# only the 3rd position can be 4dTvif int(info[9][2:-3]) % 3 != 0:continue# determine the strand by the REF column and mutation# if the ref is not same as the mutation siteif rec.ref == info[9][-3]:pre = fa_in.fetch(rec.chrom, rec.pos-3, rec.pos-1)else:tmp = fa_in.fetch(rec.chrom, rec.pos, rec.pos+2)tmp.upper()pre = rev_dict[tmp[1]] + rev_dict[tmp[0]]if pre not in codon:continuebcf_out.write(rec)

「生信练习题」从SnpEff注释得到的VCF中过滤4DTV位点相关推荐

「一道面试题」输入URL到渲染全面梳理中-页面渲染篇
前置知识此文是一道面试题,又不仅仅是一道面试题,不过这道题共分了三篇来说,嗯..可想而知接上文,上文我们讲了网络通信的部分,详细请看「一道面试题」输入URL到渲染全面梳理上-网络通信篇, 那么该说 ...
「2022 最新版」未认证微信公众号图文中插入外部链接教程
如何在微信公众号图文中插入外部链接呢?作为一名公众号小编,领导经常要求在图文中直接访问外部链接,但是由于微信平台的限制,公众号图文中不允许直接访问外链,只能插入其他公众号文章的链接. 现在,可以通过小 ...
中国二维码--汉信码(中国主导的首个二维码码制国际标准「汉信码」ISO/IEC 20830:2021《信息技术自动识别与数据采集技术汉信码条码符号规范》)
国际标准化组织(ISO)和国际电工协会(IEC)正式发布汉信码 ISO/IEC 国际标准 --ISO/IEC 20830:2021<信息技术自动识别与数据采集技术汉信码条码符号规范>. ...
「点击领取」数据社免费红包封面发放中
祝大家新的一年兔步青云前兔似锦兔飞猛进大展宏兔兔年大吉今年的微信红包封面来了,老读者应该都知道这个封面的来源<我是如何设计这款数据人专属红包封面的>
生信宝典教程大放送，一站式学习生信技术
生物信息学包含生物数据分析.数据可视化.重复工作程序化,是生物.医学科研必备的技能之一.生信宝典精心组织生信学习系列教程.生信工具精品教程,通过大量的生信例子.关键的注释.浓缩的语句和录制的视频帮助快 ...
生信宝典文章集锦，一站式学习生信！众多干货，有趣有料
生信的作用越来越大,想学的人越来越多,不管是为了以后发展,还是为了解决眼下的问题.但生信学习不是一朝一夕就可以完成的事情,也许你可以很短时间学会一个交互式软件的操作,却不能看完程序教学视频后就直接写程 ...
告别「灭霸式审稿」，IJCAI-21 的投稿者爽到家！
点击上方"视学算法",选择加"星标"或"置顶" 重磅干货,第一时间送达转自 | 新智元编辑 | Q 没有「灭霸式审稿」的IJCAI-21 ...
你想要的生信知识全在这——生信宝典目录 (181202)
生信的作用越来越大,想学的人越来越多,不管是为了以后发展,还是为了解决眼下的问题.但生信学习不是一朝一夕就可以完成的事情,也许你可以很短时间学会一个交互式软件的操作,却不能看完程序教学视频后就直接写程 ...
生信宝典文章集锦，你想看的都在
本文转载自"生信宝典",己获授权. 生信的作用越来越大,想学的人越来越多,不管是为了以后发展,还是为了解决眼下的问题.但生信学习不是一朝一夕就可以完成的事情,也许你可以很短时间学会 ...

「生信练习题」从SnpEff注释得到的VCF中过滤4DTV位点

「生信练习题」从SnpEff注释得到的VCF中过滤4DTV位点相关推荐

最新文章

热门文章