SnpEff软件安装以及如何从vcf文件中提取4DTV位点

前言

小编最近看到一个文章，该文章从同义突变位点中提取4dtv用来做后续分析，原文见 Genomic Consequences of Long-Term Population Decline in Brown Eared Pheasant ,所以很是好奇是怎么实现的，但文章没有说明怎么实现的，还好通过百度最终找到了破解方法，这里主要参考了大神的博客，博客原链接见从SnpEff注释得到的VCF中过滤4DTV位点

安装环境

（1）java #该软件是Java编写，小编这里用的java-1.8.0

（2）snpEff #上边提到的博主基于该软件编写的提取4DTV脚本

（3）python3以及pysam模块

pysam模块

pip3 install pysam -i https://pypi.tuna.tsinghua.edu.cn/simple

snpEff软件安装

wget  https://nchc.dl.sourceforge.net/project/snpeff/snpEff_v4_5covid19_core.zip
unzip snpEff_v4_5covid19_core.zip #解压安装
cd snpEff ; chmod 755 *jar #进入目录，可以看到两个jar执行程序，分别是SnpSift和snpEff，给执行权限

snpEff软件配置

该软件是一个很强大的软件，自带有很多基因组的数据库，小编在此不再介绍，有兴趣可以上官网查看，SnpEff & SnpSift,小编这里带你怎么根据自己有的基因组和gff3文件搭建该基因组的注释数据库，以从NCBI下载下来的苹果基因组和gff文件为例

mkdir -p data/apple
ascp -i ~/asperaweb_id_dsa.openssh  -QTr -l200m  anonftp@ftp.ncbi.nlm.nih.gov:genomes/all/GCF/002/114/115/GCF_002114115.1_ASM211411v1/GCF_002114115.1_ASM211411v1_genomic.fna.gz ./  #ascp用法见小编关于该工具的博客
ascp -i ~/asperaweb_id_dsa.openssh  -QTr -l200m  anonftp@ftp.ncbi.nlm.nih.gov:genomes/all/GCF/002/114/115/GCF_002114115.1_ASM211411v1/GCF_002114115.1_ASM211411v1_genomic.gff.gz ./

基因组和gff注释文件下载下来后，解压并命名为sequences.fa和genes.gff，这里命名错误后续会报错

cd ../../
vim snpEff.config,修改配置文件，找到Databases & Genomes，修改如下，注意命名需要严格参照小编的来，apple是上边建立在data数据库里的文件名，后缀必须有.genome

# Databases & Genomes
#Malus_domestica
apple.genome    : Malus_domestica
apple.reference : https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/114/115/GCF_002114115.1_ASM211411v1/GCF_002114115.1_ASM211411v1_genomic.fna.gz \     # Genome sequence
                 https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/114/115/GCF_002114115.1_ASM211411v1/GCF_002114115.1_ASM211411v1_genomic.gff.gz \     # gff3  

snpEff自由数据库建立

这里需要较长时间，最好是在后台nohup运行,尤其是比较大的基因组

java -jar snpEff.jar build -gff3 -v apple 

提取4DTV位点

java -jar snpEff.jar  ann apple apple.vcf.gz >apple.snpEff.vcf
python3 calc_4dTv_in_eff_vcf.py snpEff.vcf 4dTv.snpEff.vcf genome.sequences.fa

上边的apple apple.vcf.gz为最原始的输入vcf文件，4dTv.snpEff.vcf既为最终需要的4DTV输出vcf文件，这里最核心的为 calc_4dTv_in_eff_vcf.py 脚本，该脚本引自 从SnpEff注释得到的VCF中过滤4DTV位点
，脚本代码如下：

#!/usr/bin/env python3

from sys import argv
from pysam import VariantFile
from pysam import FastaFile

file_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():
    # no annotation available
    if len(rec.info.items()) == 0:
        continue
    ann = rec.info['ANN']
    info = rec.info['ANN'][0].split('|')
    # only use synonymouse variants
    if info[1] != "synonymous_variant":
        continue
    # only the 3rd position can be 4dTv
    if 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 site
    if rec.ref == info[9][-3]:
        pre = fa_in.fetch(rec.chrom, rec.pos-3, rec.pos-1)
        pre = pre.upper()
    else:
        tmp = fa_in.fetch(rec.chrom, rec.pos, rec.pos+2)
        tmp = tmp.upper()
        pre = rev_dict[tmp[1]] + rev_dict[tmp[0]]
    if pre not in codon:
        continue
    bcf_out.write(rec)