基本概念
SNP:单核苷酸多态性主要是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态性。
SNP位点注释:顾名思义也就是把SNP位点信息和基因组信息相关联,比如有多少SNP位点落在编码区,这些位点有多大比例是非同义突变。
数据准备
做SNP位点注释需要三样原件,分别是:参考基因组序列文件,参考基因组序GFF注释文件和之前Call好的SNP数据VCF文件。
使用SnpEff注释VCF文件
推荐使用snpEFF来做突变注释[1]。
软件安装:
wget https://snpeff.blob.core.windows.net/versions/snpEff_latest_core.zip
unzip snpEff_latest_core.zip
构建参考数据库
首先,我们需要编辑snpEff文件夹下的snpEff.config, 在# Databases & Genomes后增加一个新的物种信息。
#Databases & Genomes
物种名.genome : 物种名
这里的物种名是物种注释信息的名称,必须要和输入文件名对应。这里在snpEff/data/路径下建一个以你的物种名命名的子文件夹,用来放构建参考数据库所需要的文件。
mkdir -p snpEff/data/物种名
在data/物种名/文件夹下面存放两个文件
参考基因组:sequences.fa.gz:
注释文件,GFF3格式(也可以是GFF2格式):genes.gff.gz
之后就可以用build子命令进行构建参考数据库。
java -jar snpEff.jar build -gff3 -v 物种名
注释:
进入工作路径,比如在自己的文件夹下的snpEff_results文件夹,执行以下命令,即可获得注释文件。
java -jar snpEff.jar ann 物种名 input.vcf.gz > snpeff.vcf
默认情况下snpEff的注释信息会很多,都存放在snpeff.vcf文件中,其实它还是一个vcf格式文件,只是在原来的输出文件基础上增加了一个tag:ANN
我可以将这个vcf格式文件稍微处理下,保留原来的vcf文件的前5列,再加上ANN列形成一个新文件来查看:
perl -alne 'next if $_ =~ /^#/;$F[7] =~ /(ANN=\S+)/;print "$F[1]\t$F[2]\t$F[3]\t$F[4]\t$F[5]\t$1"' snpeff.vcf >snpeff.anntag.vcf
ANN列里的信息包括很多,大致列举几个,具体的内容大家看看官网官方解释
Allele :字母表示该突变在目的基因组上的碱基
-
Annotation :由sequence ontology定义的突变类型,第一个突变位点的downstream_gene_variant在SO系统中的定位如下:
如果变异位点属于多个类型时,多个类型之间用&符号连接。
-
Annotation_Impact :表示snpEFF对这个突变的影响的预测,有4个程度(HIGH, MODERATE, LOW, MODIFIER)
Gene Name :表示该突变所在基因的基因名,如果是这个突变位于intergenic,则使用该突变离的最近的一个基因的名称
Gene ID :表示gene id
Feature type :表示突变所在区域的类型,比如transcript, motif, miRNA等
Feature ID :表示Feature type对应的id
Transcript biotype :转录本类型, 通常采用Ensembl数据库的转录本类型
Rank / total :只有当变异位点位于基因区域时才有值,会给出变异位点所处的exon/intron的编号和该基因的exon/intron的总数,比如一个突变位点位于基因的第3个exon上,该基因一共有12个exon, 对应的Rank的值为3/12
当变异位点位于基因区域以外时,该字段的值为空HGVS.c :采用HGVS标准命名的基因水平的变异情况
HGVS.p:采用HGVS标准命名的蛋白质水平的变异情况,只有当突变位点位于编码区是才会有值
cDNA.pos/cDNA.length:突变位点在cDNA上的位置/cDNA的总长度
CDS.pos/CDS.length:突变位点在CDS上的位置/CDS的总长度
AA.pos/AA.length:突变位点在氨基酸序列上的位置/氨基酸序列的总长度
Distance:变异位点与最近的特征的距离,当变异位点位于基因间区时,会给出与最近的基因之间的距离;当变异位点位于exon区域时,会给出与最近的内含子边界的距离,不同的情况,距离的定义不同。
-
ERRORS/WARNINGS/INFO:对注释结果的可靠程度进行评估,各种取值代表的含义如下图:
我们可以用如下几个参数来简化输出
-no-downstream
-no-upstream
-no-utr
-no-intergenic
-no-intron
比如说我们只关注CDS中的注释信息,不考虑上游、下游、UTR、基因间区等信息
java -jar snpEff.jar ann -no-utr -no-downstream -no-upstream -no-intergenic 物种名 input.vcf.gz > snpeff.vcf
最终除了输出的vcf文件外,我们还会有额外两个文件,记录总结性信息:
snpEff_genes.txt: 总结每个基因的突变位点数
snpEff_summary.html: 总结突变的类型数,可视化数据