SnpEff是一种变体注释和效果预测工具，它注释和预测遗传变异的影响（例如氨基酸变化）。
1. 软件安装
#方法一：软件下载
wget https://downloads.sourceforge.net/project/snpeff/snpEff_latest_core.zip
#解压缩
unzip snpEff_latest_core.zip
#方法二：conda下载
conda install -y snpeff
2. 运行
2.1 构建SnpEff数据库
SnpEff软件的运行，首先需要基因组fasta序列信息和GTF注释信息，来构建数据库。
配置文件步骤如下：
1.在~/snpEff/目录中，创建一个文件夹：data
2.在~/snpEFF/data目录下，创建两个文件夹
AT_10/ genomes/
这两个文件夹中，分别放置了GTF文件和基因组文件
genes.gtf sequences.fa
3.编辑~/snpEff/snpEff.config文件
在文件的最后一行添加信息：
AT_10.genome: AT
构建数据库步骤如下：
java -jar ~/snpEff/snpEff.jar build -c ~/snpEff/snpEff.config -gtf22 -v AT_10
#参数说明
java -jar: Java环境下运行程序
-c snpEff.config配置文件路径
-gtf22 设置输入的基因组注释信息是gtf2.2格式
-gff3 设置输入基因组注释信息是gff3格式
-v 设置在程序运行过程中输出的日志信息
最后的AT_10参数 设置输入的基因组版本信息，和~/snpEff/snpEff.config配置文件中添加的信息一致
2.2 使用SnpEff进行注释
java -Xmx10G -jar ~/snpEff/snpEff.jar eff -c ~/snpEff/snpEff.config AT_10 positive.vcf > positive.snp.eff.vcf -csvStats positive.csv -stats positive.html
最终会产生四个文件 positive.snp.eff.vcf / positive.html / positive.csv / positive.genes.txt
3. 变异注释可视化(例如检测某基因X的基因区变异）
3.1 输入文文件
.snpEff文件：
19 63539072 A C downstream_gene_variant MODIFIER
19 63539667 T C synonymous_variant LOW p.Leu63Leu
19 63539866 G C missense_variant MODERATE p.Gly129Ala
19 63539870 G A synonymous_variant LOW p.Ala130Ala
19 63539889 A C missense_variant MODERATE p.Ile137Leu
该文件由注释结果文件positive.snp.eff.vcf改变格式产生：
grep -v "#" positive.snp.eff.vcf | awk '{print $1"\t"$2"\t"$4"\t"$5}' > snpEff1
grep -v "#" positive.snp.eff.vcf | awk '{print $8}' |awk -F"ANN=" '{print $2}' | awk -F"|" '{print $2"\t"$3"\t"$11}' > snpEff2
paste -d "\t" snpEff1 snpEff2 > positive.snpEff
rm snpEff1 snpEff2
.table文件：
AFG-L1 CGGTCCCAAGGCCCCCACGNNTCNNNNNNTTGG
AFG-L3 CGGTCCCAAGGCCCCCACGAATCNCTAAATTGG
AT18487A CGGTCCNNAGGCCCCCACGGGTCCCTACNTTGG
AT18488A CGGTCCCAANGNNNNNNCGGGTCCCWAMRTTGG
AT18489A CGGTCCNNAGGCCCCCACGGGTCCNNAAATTGG
该文件包含两列，第一列是样本名，第二列是geneX的基因区变异内容。可以通过Tassel将VCF转换成table格式。
Nucleotide Codes
(Derived from IUPAC)...
A A:A
C C:C
G G:G
T T:T
R A:G
Y C:T
S C:G
W A:T
K G:T
M A:C
+ +:+ (insertion)
0 +:-
- -:- (deletion)
N Unknown
3.2 画图
#加载包
library(ggplot2)
library(ggseqlogo)
substrRight <- function(x){
num = nchar(x)-2
gsub('[0-9.]', '', substr(x, 6, nchar(x)))
}
#读取输入文件
fasta_file = read.table("geneX.table",header=F,stringsAsFactors = F)
fasta = fasta_file[,2]
snpEff <- read.table("geneX.snpEff",header=F,stringsAsFactors = F,fill=TRUE,sep="\t")
#展示变异内容
p1 <- ggseqlogo(fasta,method="prob")+theme(axis.text.x = element_blank())+labs(title = "geneX")+theme(plot.title = element_text(hjust = 0.5,size = 25))
Ref <- as.data.frame(snpEff$V3)
colnames(Ref) <- "letter"
Alt <- as.data.frame(snpEff$V4)
colnames(Alt) <- "letter"
if(dim(snpEff)[2] >6){
Old <- as.data.frame(substring(snpEff$V7,3,5))
colnames(Old) <- "letter"
Pos <- as.data.frame(gsub('[a-zA-Z.*]', '', snpEff$V7))
colnames(Pos) <- "letter"
New <- as.data.frame(substrRight(snpEff$V7))
colnames(New) <- "letter"
all <- rbind(Ref,Alt,Old,Pos,New)
num <- dim(Ref)[1]
aln <- data.frame(
letter = all,
Type=rep(c("Ref","Alt","Old","Pos","New"), each=num),
x=rep(1:num,5)
)
#展示注释信息
p2 <- ggplot(aln, aes(x, Type)) +
geom_text(aes(label=letter,)) +
scale_y_discrete(limits=c("New","Pos","Old","Alt","Ref"))+
scale_x_continuous(breaks=1:10, expand = c(0.105, 0)) + xlab('') +
theme_logo() +
theme(legend.position = 'none', axis.text.x = element_blank())
} else{
all <- rbind(Ref,Alt)
num <- dim(Ref)[1]
aln <- data.frame(
letter = all,
species=rep(c("Ref","Alt"), each=num),
x=rep(1:num,2)
)
p2 <- ggplot(aln, aes(x, Type)) +
geom_text(aes(label=letter,),check_overlap = TRUE) +
scale_y_discrete(limits=c("Alt","Ref"))+
scale_x_continuous(breaks=1:10, expand = c(0.105, 0)) + xlab('') +
theme_logo() +
theme(legend.position = 'none', axis.text.x = element_blank())
}
snpEff$impact <- NA
snpEff[which(snpEff$V6=="HIGH"),8] <- 4
snpEff[which(snpEff$V6=="MODERATE"),8] <- 3
snpEff[which(snpEff$V6=="LOW"),8] <- 2
snpEff[which(snpEff$V6=="MODIFILER"),8] <- 1
bp_data <- data.frame(
x=snpEff$V2,
Impact=snpEff$impact
)
#展示注释的效应大小
p3 <- ggplot(bp_data, aes(factor(x), Impact))+
geom_bar(stat = "identity", fill="grey")+
theme_logo()+
xlab("")+
theme(axis.text.x = element_text(angle = 90))
suppressMessages(require(cowplot))
p <- plot_grid(p1,p2,p3,ncol = 1, align = "v")
print(p)
结果展示
geneX
