Variant 分析阶段小结3-注释碎碎念

5300字，约120分钟

variant 注释

通过上面几步内容，我们找到了一些可信度相对高的突变位置，接下来一定会进行的一个内容就是对已有突变位点进行注释和功能预测。

注释目前常用的工具有两种，一个是snpEFF，另一个是annovar。注释的思路也可以分为两类：一类是按照基因注释，另一类是按照位置注释。

SnpEFF

关于SnpEFF这个软件的用法说难不难，但是需要注意的地方也不少。输入文件是前面生成过滤过的vcf文件（也可以是跑ChIP-seq等数据得到的peak文件），正式使用前首先需要下载好所研究物种的数据库，比对后的新vcf文件会在INFO这一列生成一个名字是ANN的tag 。软件下载安装好之后，注释过程主要有三步：

查看官方已经准备好的数据库

据说目前已经有了超过2500个物种的注释信息，植物研究常用的拟南芥和水稻自然是有的。

java -jar snpEff.jar databases | less

下载自己需要的数据库

这里我们下载水稻的数据库信息为例

java -jar /home/zf/software/snpEff/snpEff.jar Oryza_sativa

下载好之后，会在软件安装目录中的data目录下出现名字Oryza_sativa，里面会有一个后缀是bin的文件，这个文件就是后面注释时要使用的文件。

画外音：开始我使用的是官方已经准备好的数据库，但是因为这个植物相关的数据库来自Ensembl，里面的信息很多物种并不是最新的，而且会有很多我们不想要的注释内容存在。又因为水稻目前主要有来自MSU 和RAPDB的两版注释信息，所有我们选择重新构建自己的水稻注释数据库。

构建自己的数据库

# 修改 snpEff.config 文件vi snpEff.config# 添加如下两行信息，含义是构建两个版本的Rice数据库# 分别是msu和rapdb-----rapdb_rice.genome: Ricemsu_rice.genome: Rice-----mkdir ./data/msu_ricemv ./data/msu_rice# 下载相应genome,cds,protein文件和gtf文件到该文件夹# 利用gffread 把gff转换为gtf# 重命名如下data/msu_rice/├── cds.fa├── genes.gtf├── protein.fa├── sequences.fa# build databasejava -Xmx20g -jar snpEff.jar build -v msu_rice 2>&1 | tee msu_rice.build# 构建完成后结构如下data/msu_rice/├── cds.fa├── genes.gtf├── protein.fa├── sequences.fa├── snpEffectPredictor.bin # 后续注释需要的文件

对vcf结果文件进行注释

java -Xmx32g -jar /home/zf/software/snpEff/snpEff.jar msu_rice filtered_freeabyes.vcf -v > filtered_freeabyes_anno.vcf

结果文件

• snpEff_genes.txt
• snpEff_summary.html
• filteredindelsgatk_anno.vcf

其中html文件会给出主要（详细）的注释统计信息，所有关心的统计信息这个html里面都有而且还进行了简单的可视化，除了丑没啥大毛病。我这里仅用了之前freebayes产生的结果作为测试文件。例如：

整体情况

突变的影响和功能

突变分布的区域

再比如SNP的类型

snpEFF_gene.txt 这个文件里面统计了每一个基因的相关突变情况，其基因信息主要来自Ensembl数据库。

filteredindelsgatk_anno.vcf 是主要的注释vcf文件，这里面和源文件相比变化的是INFO列ANN tag , 形式如下所示，每两个竖线之间具体表示一类信息。总的来说一条注释一共有16列,有些不一定每个位点都有，如果一个位点有多条注释则使用逗号分隔。可以参考官网的详细解释。

多注释问题

通常情况下，一个位点会出现多个注释信息，比如可能是一个基因的上游，也可能是一个基因的下游，或者一个基因有过个转录本，再或者位点本身就是MNPs。既然是多个就有一个谁在前谁在后的问题，SnpEff的规则是影响大的在前，突变有害的在前，实在没得挑了，按照基因位置排序。如果一个位点有很多条注释信息并不利于我们进行后续统计，这是可以使用软件本身提供的一些脚本进行处理。

如果想保留所有的注释信息，可以让每一个注释信息独立一行

cat filtered_freeabyes_anno.vcf |/home/zf/software/snpEff/scripts/vcfEffOnePerLine.pl |less

也可以每个变异只保留第一条注释

cat filtered_freeabyes_anno.vcf |/home/zf/software/snpEff/scripts/vcfAnnFirst.py |less

annovar

这个软件也非常常用不过用起来没有SnpEff那么顺手，一是需要转换一下vcf文件变成软件支持的格式，二是生成的结果文件比较简单（是缺点但也是优点），三是生成数据库的时候需要对注释文件进行一些转换，格式但凡有些问题就会很心累。当然，如果做的是人的研究不存在这个问题，而且关于人各种注释信息 annovar 都支持的非常好。

研究植物的话，无论基于基因还是位置的注释都需要自己生成gff3 或者bed 文件。

使用过程主要有一下几个步骤。

生成注释数据库

和使用SnpEff类似的，首先准备好水稻两个版本注释的gtf文件和基因组文件。

首先要把gtf文件（如果是gff文件需要先用gffread转换成gtf文件）转换为genepred文件，转换的工具是gtfToGenePred，要想用这个软件，得先安装kentUtils。多少有点麻烦。

# 转换 msu文件gffread genes.gff3 -T -o msu.gtfgtfToGenePred -genePredExt msu.gtf msu_refGene.txt

接下来是生成mRNA fasta文件

retrieve_seq_from_fasta.pl --format refGene --seqfile ../msu_rice/sequences.fa msu_refGene.txt --out msu_refGeneMrna.fa

然后把这两个文件分别放在annovar下的msu文件夹即可。

生成输入文件

convert2annovar.pl -format vcf4 filtered_samtools.vcf > filtered_samtools.annovar.txt

根据基因信息进行注释

annotate_variation.pl -buildver msu filtered_samtools.annovar.txt /home/zf/software/annovar/msu/ --geneanno --outfil filtered_samtools.annovar

结果文件

filtered samtools.annovar.variantfunction

intergenic LOC_Os01g01100(dist=1986),LOC_Os01g01110(dist=1684) Chr1 55398 55398 A G hom 38.415 2exonic LOC_Os01g01380 Chr1 194123 194123 C T hom 39.4149 2

filtered samtools.annovar.exonicvariant_function

line2 nonsynonymous SNV LOC_Os01g01380:LOC_Os01g01380.1:exon6:c.C260T:p.A87V, Chr1 194123 194123 C T hom 39.4149 2

关于结果的解读可以参考官网，和SnpEff相比，variantfunction文件第一列是突变的位置，第二列是相关基因，后面就是输入文件。而exonicvariant_function则是专门针对外显子区域的突变进行注释。第一列的信息是在原始注释文件的行数，第二列是突变的种类，和SnpEff不同的是，如果有多个突变注释信息，annovar只会根据突变的权重列出最重要的一个。

关于突变区域和类型的解释参考官网即可，这里不在着墨。

基于位置注释

基于位置注释有多种具体的思路，比如可以是组蛋白修饰区域，转录因子集合区域，或者是miRAN区域，也可以是基因组上的重复片段，还可以是基因组上一些特征区域比如启动子。

其中重复区域的注释比较重要，可以参考下面这段解释

Genetic variants that are mapped to segmental duplications are most likely sequence alignment errors and should be treated with extreme caution. Sometimes they manifest as SNPs with high fold coverage and probably high confidence score, but they may actually represent two non-polymorphic sites in the genomes that happen to have the same flanking sequence.

针对植物研究而言没有现成的数据下载，所以通常使用我们自己分析生成的gff3文件或者bed文件来进行注释。

比如我们利用repeatmasker 软件找到了水稻中的简单重复区域，可以利用其生成的注释文件对我们的vcf结果进行基于位置信息的进一步注释。

说到RepeatMasker这个软件就再说两句，软件安装完成后利用以下命令，可以完成repeat的寻找

RepeatMasker -parallel 10 -species rice -gff -dir repeat rice.fa

会得到四个结果文件，其中gff和.out文件是需要的，但是这个gff文件是gff2格式，并不是annovar 要求的gff3 格式（每一列必须有ID=），虽然这个软件有一个脚本可以out文件改为gff3，但还是不符合要求，需要讲结果文件的Target=改为ID=

cat /home/zf/software/annovar/msu/rice_repeat.gff3|sed 's/Target=/ID=/' |cut -d' ' -f1 >rice_repeat.gffannotate_variation.pl -regionanno -buildver msu -dbtype gff3 -gff3dbfile rice_repeat.gff filtered_samtools.annovar.txt /home/zf/software/annovar/msu/ --outfile filtered_samtools.annovar

注释结果: 在测试文件的9万个变异位点中，有1088个注释到了重复区域。

gff3 Score=14;Name=(AAT)n Chr1 444478 444478 - GTA hom 121 3gff3 Score=14;Name=(AGGAA)n Chr1 528020 528020 - AAAGGGAAGAG hom 30.4183 1gff3 Score=60;Name=(AT)n Chr1 563936 563936 A T hom 36.4154 2gff3 Score=26;Name=(CAATAC)n Chr1 574725 574725 A G hom 36.4154 2

可以看出富集最多的重复区域是A-rich

cut -f2 filtered_samtools.annovar.msu_gff3|cut -d';' -f2|sort \|uniq -c |sed 's/ *//'|tr ' ' '\t' |sort -k1,1nr|head118 Name=A-rich57 Name=(AT)n40 Name=(TA)n28 Name=GA-rich16 Name=(TC)n14 Name=(ATAT)n14 Name=(ATT)n12 Name=(AG)n

• Variant 分析阶段小结2- 变异寻找碎碎念
• Variant 分析阶段小结1-基础碎碎念
• 谁来拯救你 我的屁屁踢
• RNA-seq 从原理到应用
• 生物统计学与R极简手册
• 用 Excel 怎么了，你咬我啊？
• 高效使用云笔记只需“一五一十”

作者博客地址：http://kaopubear.top