利用snpEff对基因型VCF文件进行变异注释的详细方法

利用snpEff对VCF文件进行变异注释

群体遗传研究中，在获得SNP位点后,我们需要对SNP位点进行注释，对这些SNP位点进行更深的了解。

snpEff是一个用于对基因组单核苷酸多态性(SNP)进行注释的软件，snpEff软件可以用于对VCF文件进行变异注释，使用时需要先进行安装，然后构建参考基因组数据库，即可对VCF文件进行注释，下面进行用法介绍。

---

安装方法

首先安装好java环境，通过官网下载最新版本的软件压缩包，然后解压即可，最好安装在自己熟悉的目录下。

# Download latest version
wget https://snpeff.blob.core.windows.net/versions/snpEff_latest_core.zip

# Unzip file
unzip snpEff_latest_core.zip

另外，推荐一种更加简单快捷的方法，直接使用conda安装，命令如下：

conda install snpeff -c bioconda

配置数据库

有一些物种已经有官方的注释数据库，可以直接进行下载， 比如人类基因组，用如下代码下载：

java -jar snpEff.jar download GRCh38.76

如果需要查找哪些物种有现成的数据库，可以使用如下命令：

java -jar snpEff.jar databases

如果官方没有给出数据库，就需要自行建立数据库，以下用小麦举例，通过参考基因组和注释文件建立注释数据库。

下载数据文件

需要两个主要文件，iwgsc_refseqv2.1_assembly.fa是组装好的基因组序列，iwgscRefseqv2.1HCLC.gff3是注释文件，这两个文件可以在https://wheatgenome.org/Projects/Reference-Genome-Project/RefSeq-v2.1下载，如果是研究其他物种，也可以在网络上找到这两个文件。

配置文件修改

需要修改的配置文件在snpEff目录下， 文件名为snpEff.config，打开这个文件，在最后一行添加自定义数据库信息（这里的wheat可以自定义，但是要保持一致）

echo "wheat.genome:wheat" >> snpEff.conf

数据库路径设置

按照如下设置进入snpEff的安装目录并创建文件夹：

cd snpEff  #进入 snpEff 目录下
mkdir data  #新建 data 目录
cd data  #进入 data 目录下，必须在该目录
mkdir genomes  #新建 genomes 目录，用于建立 wheat 目录中的 bin
mkdir wheat  #新建 wheat 目录,对应物种名或后续软件调用的参数名

将参考基因组.fa文件放在genomes文件夹中，并改名为wheat.fa，并将gff注释文件放在wheat文件夹中改名wheat.gff，最终形成如下文件结构。

$ tree
.
├── genomes
│   └── wheat.fa -> /NGS/Ref/IWGSC_V2.1_fa_gff/iwgsc_refseqv2.1_assembly.fa
└── wheat
    ├── genes.gff -> /NGS/Ref/IWGSC_V2.1_fa_gff/iwgscRefseqv2.1HCLC.gff3

运行脚本形成bin文件

准备好文件后，回到软件目录，并执行以下命令，自动生成参考文件。检查一下存放参考基因组注释文件的目录下是否出现一些以.bin结尾的文件（数量与参考基因组染色体数有关），有就代表构建成功。

java -jar snpEff.jar build -gff3 -v wheat -d -noCheckCds -noCheckProtein
    
    ├── sequence.1A.bin
    ├── sequence.1B.bin
    ├── sequence.1D.bin
    ├── sequence.2A.bin
    ├── sequence.2B.bin
    ├── sequence.2D.bin
    ├── sequence.3A.bin
    ├── sequence.3B.bin
    ├── sequence.3D.bin
    ├── sequence.4A.bin
    ├── sequence.4B.bin
    ├── sequence.4D.bin
    ├── sequence.5A.bin
    ├── sequence.5B.bin
    ├── sequence.5D.bin
    ├── sequence.6A.bin
    ├── sequence.6B.bin
    ├── sequence.6D.bin
    ├── sequence.7A.bin
    ├── sequence.7B.bin
    ├── sequence.7D.bin
    ├── sequence.Unknown.bin
    └── snpEffectPredictor.bin

数据库构建完成后就可以直接用了，下一次不用再重新弄。

使用方法

准备内容

• 环境：Linux 或 Ubuntu，已经安装openjdk
• 文件：基因型变异信息 VCF 格式的文件
• 参考文件：gff 或 gtf 注释文件、参考基因文件
• 软件：SnpEff

运行程序

如果是通过conda安装，直接运行以下命令即可调用，对vcf文件进行注释。

snpEff wheat ./xxx.vcf.gz > ./xxx_snpeff.vcf.gz

如果是通过本地安装，可以用java调用程序进行计算，推荐使用这种方法，更加稳定。

java -jar ~/snpeff-5.1-2/snpEff.jar 
-c ~/snpeff-5.1-2/snpEff.config wheat 
../xxx.vcf.gz > ./xxxsnp.vcf.gz

等待注释完成后会生成snpEff_genes.txt文件和snpEff_summary.html文件，记录了注释的摘要信息，另外生成一个新的vcf文件包含详细注释信息。

结果查看

运行完成后会生成一个html的网页文件，里面记录了很多重要信息，接下来进行解读（参考知乎大佬天火三玄变的帖子）

摘要信息

从上往下依次是：基因组（物种名）、注释日期、注释命令、警告信息、错误信息、输入文件行数、变异位点数（过滤之前）、非变异位点数（与参考基因组碱基一致）、变异位点数（过滤之后）、具有ID的变异位点数、非双等位基因组SNP位点数、effects个数、参考基因组总长度、参考基因组有效长度、变异率（参考基因组有效长度/变异位点数）

各染色体变异率

从左往右：染色体编号、长度、变异位点数、变异率（多少个碱基中有一个变异位点）

变异类型

包括：SNP（单核苷酸多态性）、MNP（多核苷酸多态性）、INS（插入变异）、DEL（缺失变异）、MIXED（混合变异）、INV(倒位变异）、DUP（重复变异）、BED（易位变异）、INTERVAL（间隔变异）

有效影响数量

功能分级有效数

MiSSENSE（错义突变）、NONSENSE（无义突变）、SILENT（沉默突变）

有效变异数和百分比

下图左边为按类型划分有效变异数，包括（从上往下）：3’端主要UTR变异（UTR是成熟mRNA分子5'或3'端不被翻译的部分，一般在mRNA转运、稳定性和翻译调节中起重要作用）、5’端主要UTR提前启动子获得变异、5’端主要UTR变异、下游基因变异、起始密码子编码变异、基因间隔区、内含子变异、剪接受体变异、剪接供体变异、剪接区域变异、起始缺失、起始保留变异、终止获得、终止缺失、终止保留变异、同义变异、上游基因变异。

右边为按区域划分有效变异数，包括（从上往下）：下游、外显子、间隔区、内含子、剪接位点受体、剪接位点供体、剪接位点区域、上游、3’UTR区、5’UTR区。

SNP位点碱基变异表

可以看出SNP中哪些碱基的转换比较多（A腺嘌呤、C胞嘧啶、G鸟嘌呤、T胸腺嘧啶）

总结

在使用snpEff过程中需要注意数据库的选择和构建，根据不同版本进行计算，另外尽量避免更改染色体的展示方式，防止造成识别错误。另外可以利用vcftools将vcf中的样品信息去掉，这样文件体积会大大缩小，有利用加快注释速度。

参考资料：

https://www.jianshu.com/p/77c3a2fae4ab
https://zhuanlan.zhihu.com/p/613790756
https://pcingola.github.io/SnpEff/

本文由mdnice多平台发布