SnpEff是一種變體注釋和效果預(yù)測工具，它注釋和預(yù)測遺傳變異的影響（例如氨基酸變化）。
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 構(gòu)建SnpEff數(shù)據(jù)庫
SnpEff軟件的運行，首先需要基因組fasta序列信息和GTF注釋信息粮揉，來構(gòu)建數(shù)據(jù)庫织中。
配置文件步驟如下：
1.在~/snpEff/目錄中，創(chuàng)建一個文件夾：data
2.在~/snpEFF/data目錄下病游，創(chuàng)建兩個文件夾
AT_10/ genomes/
這兩個文件夾中唇跨，分別放置了GTF文件和基因組文件
genes.gtf sequences.fa
3.編輯~/snpEff/snpEff.config文件
在文件的最后一行添加信息：
AT_10.genome: AT
構(gòu)建數(shù)據(jù)庫步驟如下：
java -jar ~/snpEff/snpEff.jar build -c ~/snpEff/snpEff.config -gtf22 -v AT_10
#參數(shù)說明
java -jar: Java環(huán)境下運行程序
-c snpEff.config配置文件路徑
-gtf22 設(shè)置輸入的基因組注釋信息是gtf2.2格式
-gff3 設(shè)置輸入基因組注釋信息是gff3格式
-v 設(shè)置在程序運行過程中輸出的日志信息
最后的AT_10參數(shù) 設(shè)置輸入的基因組版本信息，和~/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
最終會產(chǎn)生四個文件 positive.snp.eff.vcf / positive.html / positive.csv / positive.genes.txt
3. 變異注釋可視化(例如檢測某基因X的基因區(qū)變異）
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
該文件由注釋結(jié)果文件positive.snp.eff.vcf改變格式產(chǎn)生：
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的基因區(qū)變異內(nèi)容∽涛荆可以通過Tassel將VCF轉(zhuǎn)換成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")
#展示變異內(nèi)容
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
)
#展示注釋的效應(yīng)大小
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)
結(jié)果展示
geneX
