一、简介
会得到一系列变异数据，这些变异数据只是告诉我们在基因组的某个位置发生了一段序列的改变，至于这个改变会不会影响生物学功能，我们并不清楚。而注释就是将基因组的序列变异数据转化为我们更关心的生物学功能变化的信息。
Annovar常被用在人类基因组的注释上，其实，它也可以对人类以外的基因组数据进行注释。比如，老鼠基因组的注释上。需要自己进行建立注释信息。
ANNOVAR是一个perl编写的命令行工具，能在安装了perl解释器的多种操作系统上执行。允许多种输入文件格式，包括最常被使用的VCF格式。输出文件也有多种格式，包括注释过的VCF文件、用tab或者逗号分隔的text文件。
ANNOVAR能快速注释遗传变异并预测其功能。类似的variants注释软件还有 VEP, snpEff, VAAST, AnnTools等等.
ANNOVAR 注释变异可以分成有基于基因、基于染色体区间和变异数据等三种类型. 这三种注释分别针对于每一个variant的不同方面：
-
基于基因的注释（gene-based annotation）
注释结果为突变位点位于基因的相对位置，是否改变氨基酸编码，确定受影响的氨基酸，获得变异位点的HGVS命名方式，揭示variant与已知基因直接的关系以及对其产生的功能性影响。可灵活使用RefSeq genes, UCSC genes, ENSEMBL genes, GENCODE genes或许多其他基因定义系统。 -
基于染色体区间的注释（region-based annotation）
识别特定基因组区域的变异，例如，44个物种中的保守区域，预测的转录因子结合位点, segmental duplication regions, GWAS hits, ChIP-Seq peaks, RNA-Seq peaks等等许多其他的在基因组区间的注释； -
变异数据库的注释 / 基于过滤的注释（ filter-based annotation ）
则给出这个variant的一系列信息，如： population frequency in different populations 和various types of variant-deleteriousness prediction scores, 这些可被用来过滤掉一些公共的及 probably（大概,肯定的成分较大,是most likely） nondeleterious variants. 包括Clinvar, dbSNP, Cosmic, ExAC, gnomAD等等，突变数据库整理可参考从 vcf 文件准备 ANNOVAR 数据库。鉴定特定数据库中记录的变异，例如，该变异位点是否在dbSNP中有报道，在千人基因组计划中的等位基因频率如何等等。
ANNOVAR 数据库文件实际上为特定格式的文本文件，其数据库文件命名规则为: ${path_database}/${buildver}_${database_name}.txt
ANNOVAR 所有注释结果都在 vcf 文件 INFO 列添加key-value
二、ANNOVAR的下载
填写登记表，下载ANNOVAR软件（http://annovar.openbio informatics.org/），解压文件
tar xvfz annovar.latest.tar.gz
解压后生成annovar文件夹，里面有6个perl脚本程序和2个文件夹。
### ANNOVAR 的文件结构
|-- annotate_variation.pl # 主程序，功能包括下载数据库，三种不同的注释
|-- coding_change.pl # 可用来推断蛋白质序列
|-- convert2annovar.pl # 将多种格式转为.avinput的程序
|-- retrieve_seq_from_fasta.pl # 用于自行建立其他物种的转录本
|-- table_annovar.pl # 注释程序，可一次性完成三种类型的注释
|-- variants_reduction.pl # 可用来更灵活地定制过滤注释流程
|-- example # 存放示例文件
|-- humandb # 人类注释数据库
自带了humandb是已经建立好的hg19或者GRCh37等常用的数据库，可用于人的注释。 如果要进行其他注释，需要使用 -downdb命令下载数据库到 humandb/ 目录。
三、注释
perl retrieve_seq_from_fasta.pl
--format refGene
--seqfile zunla.fasta zunla_refGene.txt
--out zunla_refGeneMrna.fa
1. 用ANNOVAR注释人类基因组variants信息
1.1 下载所有需要的注释信息库
对于基于基因的注释的数据库已经在下好的 ANNOVAR package中了。如果要进行其他注释，需要按以下命令下载数据库到 ‘humandb/’ 目录里：
perl annotate_variation.pl --downdb --buildver hg19 cytoBand humandb/
perl annotate_variation.pl --downdb --webfrom annovar --buildver hg19 1000g2014oct humandb/
perl annotate_variation.pl --downdb --webfrom annovar --buildver hg19 exac03 humandb/
perl annotate_variation.pl --downdb --webfrom annovar --buildver hg19 ljb26_all humandb/
perl annotate_variation.pl --downdb --webfrom annovar --buildver hg19 clinvar_20140929 humandb/
perl annotate_variation.pl --downdb --webfrom annovar --buildver hg19 snp138 humandb/
# -buildver: 基因组对应版本
# -webfrom annovar: 从annovar库里下载;如果annovar库中没有，则不用写该选项，会从UCSC中下载
# refGene: 数据库名称
# humandb/: 下载至该目录
这里下载的是几个通常用到的数据库：
1）cytoBand 是每个细胞间band（cytogenetic band）的染色体坐标信息 ,
2）1000g2014oct
for alternative allele frequency in the 1000 Genomes Project (version October 2014),
是2014年10版，1000基因组项目（和ExAV 外显子集合联合一样，是公开、开放的数据库）
里面供选择的等位基因频率信息
3）exac03
for the variants reported in the Exome Aggregation Consortium (version 0.3),
是0.3版外显子集合联合中报道过的variants.
4）ljb26_all
for various functional deleteriousness prediction scores from the dbNSFP database (version 2.6),
dbNSFP: A Lightweight Database of Human NonsynonymousSNPs and TheirFunctionalPredictions on ResearchGate
5）clinvar_20140929
for the variants reported in the ClinVar database (version 20140929)
ClinVar是美国国家生物技术信息中心（NCBI）于2012年11月宣布、2013年4月正式启动的公共、免费数据库。作为核心数据库，ClinVar数据库整合了十多个不同类型数据库、通过标准的命名法来描述疾病，同时支持科研人员将数据下载到本地中，开展更为个性化的研究。在遗传变异和临床表型方面，NCBI和不同的研究组已经建立了各种各样的数据库，数据信息相对比较分散，ClinVar数据库的目的在于整合这些分散的数据、将变异、临床表型、实证数据以及功能注解与分析等四个方面的信息，通过专家评审，逐步形成一个标准的、可信的、稳定的遗传变异-临床表型相关的数据库。
6）snp138
for the dbSNP database (version 138).
注意：
(i) 第一个命令中不包含 --webfrom annovar 选项, 因此是从the UCSC Genome Browser annotation database下载文件的；
(ii) --buildver hg19 选项是针对hg19这一版的基因组的；
(iii) 运行上面命令后，在 humandb/ 目录下会多几个以 hg19为前缀的文件。
1.2 table_annovar.pl注释variants
输入下列命令，用之前下载好的注释数据库来注释vcf格式文件中的variants。ANNOVAR 所有注释结果都在 vcf 文件 INFO 列添加key-value。输出的csv文件将包含输入的5列主要信息以及各个数据库里的注释。
perl table_annovar.pl <variant.vcf> humandb/ --outfile final --buildver hg19
--protocol refGene,cytoBand,1000g2014oct_eur,1000g2014oct_afr,exac03,ljb26_all,clinvar_20140929,snp138
--operation g,r,f,f,f,f,f,f
--vcfinput
<variant.vcf> 参考输入的vcf文件的名称
--protocol 选项后跟注释来源数据库的准确名称
--operation选项后跟注释的类型:
g 表示基于基因的注释（gene-based annotation）；
r 表示基于区域的注释（region-based annotation）；
f 表示基于筛选子的注释（ filter-based annotation）；
--outfile 选项是指定输出文件的前缀
关键步骤：
1、确保注释数据库的名称正确并且是按你想要在输出文件中显示的顺序排列的；
2、确保 --operation指定的注释类型顺序和--protocol指定的数据库顺序是一致的；
3、确保每个protocal名称或注释类型之间只有一个逗号，并且没有空白。
-
final.hg19_multianno.vcf.输出文件应该是以个VCF格式文件，INFO那列以key=value形式、;分割成几个小区域. eg:Func.refGene=intronic;Gene.refGene=SAMD11. 每个键值对代表一个ANNOVAR注释信息。输出文件可以用为VCF格式文件设计的基因分析软件进一步处理。 -
final.hg19_multianno.txt. 每一行代表一个variant 。用tab分隔，多余列为加上的注释信息，顺序按--protocol选项所设定的注释类型argument。
2. 用Annovar注释人类以外的基因组（Gene-based annotation）
Annovar常被用在人类基因组的注释上，其实，它也可以对人类以外的基因组数据进行注释。annovar一般只包含人类基因组注释数据库，其他的物种需要自己进行建立注释信息。
2.1 以注释小鼠基因组为例
一般如果你想看是否有某种物种,如小鼠mm9的注释库时，命令行运行
perl annotate_variation.pl -builder mm9 -downdb avdblist -webfrom annovar ./
会生成一个mm9开头的文件，里面包含小鼠mm9有多少注释数据库，然后自己可以构建一个mousedb数据库 先在annovar文件夹里面创建mousedb文件夹（名字可自取），命令mkdir mousedb 然后使用annovar文件夹下的perl程序annotate_variation.pl
perl annotate_variation.pl -downdb -buildver mm9 -webfrom annovar refGene mousedb/
这个命令能实现的是帮忙下载mm9的refGene的文件，保存在mousedb文件下，自动解压后文件名为mm9_refGene.txt。 然后程序会提示使用以下两个命令继续建库
annotate_variation.pl --buildver mm9 --downdb seq mousedb/mm9_seq
retrieve_seq_from_fasta.pl mousedb/mm9_refGene.txt -seqdir mousedb/mm9_seq -format refGene -outfile mousedb/mm9_refGeneMrna.fa
同样在annovar文件下运行这两个perl程序
perl annotate_variation.pl --buildver mm9 --downdb seq mousedb/mm9_seq
通过这个命令，会在mousedb下创建文件夹mm9_seq，并且在里面下载mm9的基因组文件chromFa.tar.gz，perl程序帮忙解压后是按染色体分开的fasta格式文件。 然后继续运行perl程序
perl retrieve_seq_from_fasta.pl mousedb/mm9_refGene.txt -seqdir mousedb/mm9_seq -format refGene -outfile mousedb/mm9_refGeneMrna.fa
该程序会会在mousedb下创建mm9_refGeneMrna.fa文件，是根据mm9_refGene.txt的信息，重新构建成的老鼠转录表达基因fasta格式文件 这样老鼠mm9 annovar gene based注释库就弄好了 以文本文件test.input为案例进行测试 生成test.input的txt格式文件，根据annovar官网介绍，只要这最基本的五列信息就可以进行注释，五列分别染色体名称，染色体上的位置，染色体上的位置，参考基因组碱基，变异碱基。
1 19215217 19215217 T C
1 33803084 33803084 A G
1 33803198 33803198 A G
1 37499237 37499237 T C
1 37499238 37499238 T C
1 37500003 37500003 T C
1 43826936 43826936 T C
1 58853960 58853960 A G
1 58854487 58854487 A G
1 60436865 60436865 T C
然后使用perl程序进行gene based的注释
perl annotate_variation.pl -out test -build mm9 test.input mousedb
注释后会生成test.variant_function，test.exonic_variant_function和test.log文件，前两个即为所需要的文件。用这个例子输出test.exonic_variant_function文件输出为空文件，因为这些位点没有在exonic区域的，所以没有结果。如果有位点在exonic中，则在test.exonic_variant_function中会更具体的描述为同义突变还是非同义突变
intronic Tfap2b 1 19215217 19215217 T C
UTR3 Bag2 1 33803084 33803084 A G
UTR3 Bag2 1 33803198 33803198 A G
UTR3 Mgat4a 1 37499237 37499237 T C
UTR3 Mgat4a 1 37499238 37499238 T C
UTR3 Mgat4a 1 37500003 37500003 T C
intronic Uxs1 1 43826936 43826936 T C
intronic Casp8 1 58853960 58853960 A G
intronic Casp8 1 58854487 58854487 A G
intronic Cyp20a1 1 60436865 60436865 T C
2.2 以注释大猩猩基因组（with the genome build identifier as panTro2.）为例。
对于gene-based annotation， ANNOVAR需要 genePred format 的 gene definition file和 FASTA format 的 transcript sequence file；
第一步：输入以下命令，下载大猩猩基因组定义文件（ gene definition file）及序列的 FASTA 文件到chimpdb/目录
perl annotate_variation.pl --downdb --buildver panTro2 gene chimpdb/
perl annotate_variation.pl --downdb --buildver panTro2 seq chimpdb/panTro2_seq
运行过程中，出现下列提示
---------------------------ADDITIONAL PROCEDURE---------------------------
--------------------------------------------------------------------------
NOTICE: the FASTA file for the genome is not available to download but can be generated by the ANNOVAR software.
PLEASE RUN THE FOLLOWING TWO COMMANDS CONSECUTIVELY TO GENERATE THE FASTA FILES (you may need to change -seqdir to -seqfile for some genomes):
annotate_variation.pl --buildver panTro2 --downdb seq chimpdb/panTro2_seq
retrieve_seq_from_fasta.pl chimpdb/panTro2_refGene.txt -seqdir chimpdb/panTro2_seq -format refGene -outfile chimpdb/panTro2_refGeneMrna.fa
第二步：注意ANNOVAR数据库中只包含人类基因组已建好的转录本，不包含其他物种的。故需要按以下命令自行建立对应物种的transcript FASTA file
perl retrieve_seq_from_fasta.pl chimpdb/panTro2_refGene.txt
--seqdir chimpdb/panTro2_seq
--format refGene
--outfile chimpdb/panTro2_refGeneMrna.fa
--seqdir说明下载的序列文件的所在目录；
--format说明 gene definition file的格式.；
--outfile 指定输出mRNA 序列文件的名称；
关键：跟在--outfile后的输出文件名应该是 <buildver>_refGeneMrna.fa这种形式，否则下一步找不到正确的 transcript FASTA sequence file.
第三步：注释variants，with the chimpanzee gene annotation:
perl table_annovar.pl <variant.vcf> chimpdb/
--vcfinput
--outfile final
--buildver panTro2
--protocol refGene
--operation g
<variant.vcf> input VCF file；
chimpdb/ directory of the downloaded data；
第四步：输出结果文件核对。 final.panTro2_multianno.txt file. The gene annotation for chimpanzee is added after the input variants.
关键：如果没有现成可用的gene definition file ，可以将基因预测工具产生的 GFF3 or GTF 文件转换成 gene definition file.
三、构建自定义注释库
ANNOVAR可以从服务器下载注释库，但是主要针对人类基因组，那么需要分析、注释其它的物种测序结果，怎么办呢，需要自建注释库。
1 以注释辣椒基因组为例
第一步：准备工作
首先，需要参考基因组的序列和注释文件，这里是名为zunla辣椒
zunla.fasta # zunla 的参考基因组序列
zunla.gff3 # zunla 的注释文件，我这里是gff3
ANNOVAR 建库需要 genePred 文件，因而需要准备几个转换 gff3 到 genePred 格式的软件
gffread # gff3 to gtf
gtfToGenePred # gtf to genePred
gffread的下载地址，需要自行编译。编译过程会自行下载https://github.com/gpertea/gclib，也可以预先下载，然后make。
git clone https://github.com/gpertea/gffread
cd gffread
make release
推荐直接conda安装
conda install gffread
gtfToGenePred的下载地址，已经编译好了，下载直接使用，他属于UCSC工具包，因而在conda环境安装gtfToGenePred的命令为
conda install ucsc-gtftogenepred
conda直接安装gtftogenepred没成功，需要加上UCSC前缀。
第二步：建立注释库
创建辣椒注释库的文件夹pepperdb，然后进去
mkdir pepperdb
cd pepperdb
转换gff3 为 gtf，推荐使用GTF格式，因为有些GFF3格式文件转换可能不正确
gffread zunla.gff -T -o zunla.gtf
转换gtf 为 GenePred
gtfToGenePred -genePredExt zunla.gtf zunla.txt
建立注释库
perl retrieve_seq_from_fasta.pl
--format refGene
--seqfile zunla.fasta zunla_refGene.txt
--out zunla_refGeneMrna.fa
这样我们的建库就完成了，获得两个重要的文件
zunla_refGeneMrna.fa
zunla_refGene.txt
强调：关于文件的命名，按照常规逻辑，这两个文件肯定不能随便命名，不然annovar无法识别！摸索了一下，前缀就是下面build参数使用的名字，这里就是zunla，下面注释就要使用“-build zunla”这个参数，对于基于基因注释的txt文件命名就是refGene，连起来就是 zunla_refGene.txt。而fa文件前缀一样，后面有稍稍差别为refGeneMrna，连起来就是zunla_refGeneMrna.fa。
第三步：转换需要注释的vcf文件
BB3.HC.vcf是HaplotypeCaller得到的vcf文件，转换成适用annovar的文件格式。得到的BB3.avinput就是我们需要注释的文件
perl convert2annovar.pl -format vcf4 BB3.HC.vcf > BB3.avinput
ANNOVAR主要使用convert2annovar.pl程序进行转换，转换后文件是精简过的，主要包含前面提到的5列内容，如果要将原格式的文件的所有内容都包含在转换后的.avinput文件中，可以使用-includeinfo参数；如果需要分开每个sample输出单一的.avinput文件，可以使用-allsample参数，等等。
ANNOVAR还主要支持以下格式转换：
- SAMtools pileup format
- Complete Genomics format
- GFF3-SOLiD calling format
- SOAPsnp calling format
- MAQ calling format
- CASAVA calling format
第四步：annotate_variation注释
perl annotate_variation.pl
-geneanno # 表示使用基于基因的注释
-dbtype refGene # 表示使用"refGene"类型的数据库
-out BB3 # 表示输出以BB3为前缀的结果文件
-build zunla
BB3.avinput
pepperdb/
-geneanno -dbtype refGene是默认值，可以省略，那么命令也可以写成
perl annotate_variation.pl
-out BB3
-build zunla
BB3.avinput
pepperdb/
得到结果文件：
BB3.exonic_variant_function # 外显子区域突变的功能、类型等
BB3.variant_function # 突变的基因及位置
BB3.log # 日志文件
第五步：table_annovar.pl注释
table_annovar.pl是ANNOVAR多个脚本的封装，可以一次性完成三种类型的注释。
我这里只有regGene类型的注释库，那么注释命令为
perl table_annovar.pl
BB3.avinput pepperdb/
-buildver zunla # 使用zunla注释库
-out BB3
-remove # 清除所有临时文件
-protocol refGene # 注释库类型为refgene
-operation g # 操作子为g
-nastring . # 缺省值用“.”代替
-csvout # 输出csv文件
image.png
table_annovar.pl也可以不经过转换，直接对vcf文件进行注释，添加-vcfinput参数就行，注释的内容将会放在vcf文件的“INFO”那一栏。但是需要注意的是不能和-csvout参数同时使用，命令如下
perl table_annovar.pl
BB3.HC.vcf
pepperdb/
-buildver zunla
-out BB3
-protocol refGene
-operation g
-nastring .
-vcfinput
-polish
所以，两种输入格式
- VCF文件：用 -vcfinput指定
- avinput
每行代表一个位点
前5列依次为：chromosome, start position, end position, the reference nucleotides, the observed nucleotides
reference nucleotides：不知道时可设置为0
observed nucleotides: insertion,deletion,block subsititution可用-表示
其余列：可有可无，如果有，在输出文件中会原样输出。
more BB3.avinput
1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
1 13211293 13211294 TC - comments: rs59770105, a 2-bp deletion
1 11403596 11403596 - AT comments: rs35561142, a 2-bp insertion
1 105492231 105492231 A ATAAA comments: rs10552169, a block substitution
1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated
cat example/ex1.avinput
1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
1 1404001 1404001 G T comments: rs149123833, a SNP in 3' UTR of ATAD3C
1 5935162 5935162 A T comments: rs1287637, a splice site variant in NPHP4
1 162736463 162736463 C T comments: rs1000050, a SNP in Illumina SNP arrays
1 84875173 84875173 C T comments: rs6576700 or SNP_A-1780419, a SNP in Affymetrix SNP arrays
1 13211293 13211294 TC - comments: rs59770105, a 2-bp deletion
1 11403596 11403596 - AT comments: rs35561142, a 2-bp insertion
1 105492231 105492231 A ATAAA comments: rs10552169, a block substitution
1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
2 234183368 234183368 A G comments: rs2241880 (T300A), a SNP in the ATG16L1 associated with Crohn's disease
该格式每列以tab分割，最重要的地方为前5列，分别是:
- 染色体(Chromosome)
- 起始位置(Start)
- 结束位置(End)
- 参考等位基因(Reference Allele)
- 替代等位基因(Alternative Allele)
- 剩下为注释部分（可选）。
ANNOVAR主要也是依靠这5处信息对数据库进行比对，进而注释变异。
报错：
chmod 777 table_annovar.pl
问题依旧！
没办法，只能跑去table_annovar.pl脚本的444行看一下具体代码
1. #generate gene anno
2. my $sc;
3. $sc = $SC_PREFIX . "annotate_variation.pl -geneanno -buildver $buildver -dbtype $protocol -outfile $tempfile.$protocol -exonsort -nofirstcodondel $queryfile $dbloc"; #20191010: add -nofirstcodondel
4. $argument and $sc .= " $argument";
5. $intronhgvs and $sc .= " -splicing_threshold $intronhgvs";
7. if ($thread) {
8. $sc .= " -thread $thread";
9. }
10. if ($maxgenethread) {
11. $sc .= " -maxgenethread $maxgenethread";
12. }
14. print STDERR "\nNOTICE: Running with system command <$sc>\n";
15. system ($sc) and die "Error running system command: <$sc>\n";
444行的内容为system ($sc) and die "Error running system command: <$sc>\n";
看得出来，这里执行了system系统命令，往上查找变量$sc的值，发现这里调用了annotate_variation.pl脚本，同样
chmod 777 annotate_variation.pl
为了避免别的脚本权限问题，干脆全部777
chmod 777 *
成功解决权限问题。
2 以构建拟南芥（Arabidopsis thaliana）的注释所需文件为例
第一步：准备
下载Arabidopsis 的 GTF file 和 genome FASTA file，到 ‘atdb’目录下。解压文件。 地址
mkdir atdb
cd atdb
wget ftp://ftp.ensemblgenomes.org/pub/release-27/plants/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.27.dna.genome.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/release-27/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.27.gtf.gz
gunzip Arabidopsis_thaliana.TAIR10.27.dna.genome.fa.gz
gunzip Arabidopsis_thaliana.TAIR10.27.gtf.gz
**第二步：
用 gtfToGenePred 工具将 GTF file 转换 GenePred file
gtfToGenePred -genePredExt Arabidopsis_thaliana.TAIR10.27.gtf AT_refGene.txt
用retrieve_seq_from_fasta.pl生成 transcript FASTA file
perl ../retrieve_seq_from_fasta.pl --format refGene --seqfile Arabidopsis_thaliana.TAIR10.27.dna.genome.fa AT_refGene.txt AT_refGeneMrna.fa
After this step, the annotation database files needed for gene-based annotation are ready. Now you can annotate a given VCF file using the procedure starting from B(iii). Please note that the ‘--buildver’ argument should be set to ‘AT’.
参考http://annovar.openbioinformatics.org/en/latest/user-guide/gene/ for more details.bases and other arguments are the same as in the human genome annotation.
总脚本如下
下载物种基因序列、注释文件
wget -c ftp://ftp.ensemblgenomes.org/pub/release-27/plants/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.27.dna.genome.fa.gz
wget -c ftp://ftp.ensemblgenomes.org/pub/release-27/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.27.gtf.gz
gzip -d Arabidopsis_thaliana.TAIR10.27.dna.genome.fa.gz
gzip -d Arabidopsis_thaliana.TAIR10.27.gtf.gz
#gtf文件格式转换
wget -c http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/gtfToGenePred
gtfToGenePred -genePredExt Arabidopsis_thaliana.TAIR10.27.gtf AT_refGene.txt
# 另一种格式转换方法，https://github.com/chengcz/pyGTF
# 使用软件包提供脚本build物种数据库，数据库buildver为AT，名称为refGene
perl retrieve_seq_from_fasta.pl --format refGene --seqfile Arabidopsis_thaliana.TAIR10.27.dna.genome.fa AT_refGene.txt --out AT_refGeneMrna.fa
五、批量注释
/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4 Sample3.varscan.snp.vcf > Sample3.annovar
/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4 Sample4.varscan.snp.vcf > Sample4.annovar
/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4 Sample5.varscan.snp.vcf > Sample5.annovar
然后用下面这个脚本批量注释
image001
Reading gene annotation from /home/jmzeng/bio-soft/annovar/humandb/hg19_refGene.txt ... Done with 50914 transcripts (including 11516 without coding sequence annotation) for 26271 unique genes
最后查看结果可知，真正在外显子上面的突变并不多
23515 Sample3.anno.exonic_variant_function
23913 Sample4.anno.exonic_variant_function
24009 Sample5.anno.exonic_variant_function
annovar软件就是把我们得到的十万多个snp分类了，看看这些snp分别是基因的哪些位置，是否引起蛋白突变
downstream
exonic
exonic;splicing
intergenic
intronic
ncRNA_exonic
ncRNA_intronic
ncRNA_splicing
ncRNA_UTR3
ncRNA_UTR5
splicing
upstream
upstream;downstream
UTR3
UTR5
UTR5;UTR3
六、一步到位：table_annovar.pl 可以同时进行3种注释
perl table_annovar.pl example/ex1.avinput humandb/ -buildver hg19 -out myanno -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation gx,r,f,f,f -nastring . -csvout -polish -xref example/gene_xref.txt
#-remove: remove all temporary files
#-operation:g,gene-based; gx,gene-based with cross-reference annotation (from -xref argument);r, region-based; f,filter-based.
#-nastring：没有对应注释，则输出`.`
#-csvout:结果用,分隔；去掉则采用默认，用Tab分隔
#-xref: whether a known genetic disease is caused by defects in this gene (this information was suffplied in the example/gene_xref.txt file in the command line) 这一项没有也OK
其中（每种数据库对应的类型参考官网）
g,gene-based,对应数据库为refGene,ensGene等
r,region-based,对应数据库为cytoBand等
f,filter-based,对应数据库为exac03,avsnp147,dbnsfp30a等
- 3种注释分开进行：annotate_variation.pl
1 gene-based
perl annotate_variation.pl -geneanno -dbtype refGene -buildver hg19 example/ex1.avinput humandb/
结果文件在example/中，ex1.avinput.variant_function和ex1.avinput.exonic_variant_function
（1）ex1.avinput.variant_function
[root@localhost example]# head ex1.avinput.variant_function
UTR5 ISG15(NM_005101:c.-33T>C) 1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
UTR3 ATAD3C(NM_001039211:c.*91G>T) 1 1404001 1404001 G T comments: rs149123833, a SNP in 3' UTR of ATAD3C
第1列：variant effects,将变异分类，如intergenic, intronic, non-synonymous SNP, frameshift deletion, large-scale duplication等
第2列：基因名，Symbol，括号中为NM_22222,为refGene编号
其余列：输入文件ex1.avinput的内容
（2）ex1.avinput.exonic_variant_function
[root@localhost example]# head ex1.avinput.exonic_variant_function
line9 nonsynonymous SNV IL23R:NM_144701:exon9:c.G1142A:p.R381Q, 1 67705958 67705958 GA comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
line10 nonsynonymous SNV ATG16L1:NM_017974:exon8:c.A841G:p.T281A,ATG16L1:NM_001190267:exon9:c.A550G:p.T184A,ATG16L1:NM_030803:exon9:c.A898G:p.T300A,ATG16L1:NM_001190266:exon9:c.A646G:p.T216A,ATG16L1:NM_198890:exon5:c.A409G:p.T137A, 2 234183368 234183368 A G comments: rs2241880 (T300A), a SNP in the ATG16L1 associated with Crohn's disease
line11 nonsynonymous SNV NOD2:NM_022162:exon4:c.C2104T:p.R702W,NOD2:NM_001293557:exon3:c.C2023T:p.R675W,16 50745926 50745926 C T comments: rs2066844 (R702W), a non-synonymous SNP in NOD2
第1列：该变异在input文件的行号
第2列：对编码基因的影响，frameshift,nonsynonymous等
第3列：被影响的基因或转录本,其中NM_22222,为refGene编号
其余列：同输入文件
用awk操作时，分隔符设定为\t；不设置时，空格也被当做分隔符，会造成错位
[root@localhost example]# head ex1.avinput.exonic_variant_function|awk -F '\t' '{print $2}'
nonsynonymous SNV
nonsynonymous SNV
nonsynonymous SNV
nonsynonymous SNV
frameshift insertion
frameshift deletion
frameshift deletion
stoploss
stopgain
frameshift substitution
[root@localhost example]# head ex1.avinput.exonic_variant_function|awk '{print $2}'
nonsynonymous
nonsynonymous
nonsynonymous
nonsynonymous
frameshift
frameshift
frameshift
stoploss
stopgain
frameshift
2 region-based
perl annotate_variation.pl -regionanno -dbtype cytoBand -buildver hg19 example/ex1.avinput humandb/
鉴定各变异的cytogenetic band,如1p36.33
结果文件在example中，ex1.avinput.hg19_cytoBand
[root@localhost example]# more ex1.avinput.hg19_cytoBand
cytoBand 1p36.33 1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
cytoBand 1p36.33 1 1404001 1404001 G T comments: rs149123833, a SNP in 3' UTR of ATAD3C
cytoBand 1p36.31 1 5935162 5935162 A T comments: rs1287637, a splice site variant in NP
HP4
cytoBand 1q23.3 1 162736463 162736463 C T comments: rs1000050, a SNP in Il
lumina SNP arrays
第1列：cytoBand
第2列：1p21.1
其余列：同输入文件
3 ilter
perl annotate_variation.pl -filter -dbtype exac03 -buildver hg19 example/ex1.avinput humandb/
结果文件在example/中,ex1.avinput.hg19_exac03_filtered(exac03中没有报道的位点）和ex1.avinput.hg19_exac03_dropped（exac03中报道的位点，包含其等位基因频率）
https://www.jianshu.com/p/84c818207240
https://www.plob.org/article/9976.html
snpEff
https://www.bioinfo-scrounger.com/archives/268/
重要
http://www.bio-info-trainee.com/2028.html
https://www.jianshu.com/p/cea3b179b8a9
ANNOVAR注释变异VCF结果说明
https://www.omicsclass.com/article/464
https://www.omicsclass.com/article/804
http://blog.sina.com.cn/s/blog_71df25810102ybtt.html
Hui Y, Kai W. Genomic variant annotation and prioritization with ANNOVAR and wANNOVAR[J]. Nature Protocols, 2015, 10(10).
网友评论