ANNOVAR 软件用法还可以更复杂
写在前面：
我多次在博客(生信菜鸟团)里面写过annovar软件的用法，而且在《直播我的基因组》里面也使用过它，https://mp.weixin.qq.com/s/GxqzKU7OPAMbjbZ5fPk7oQ 但那些都是最粗浅的使用而已，并没有深入了解ANNOVAR 软件的方方面面。
这次耗费15个小时系统性的回顾了该软件，希望可以做到教学上的最佳教程。
而且ANNOVAR软件更新频率也不容小觑，最新版是 ANNOVAR (2018Apr16) 简单注册后即可下载。
琳琅满目的人类变异注释数据库
能叫得上名字的数据库就有 dbSNP,ExAC,ESP6500,cosmic,gnomad,1000genomes,clinvar,gwas, 都被其收集整理好了，而且提供多个不同参考基因组版本的下载链接，当然， 我们默认他们做的工作都是准确无误的，毕竟自己去一个个下载数据库一个个格式化成自己需要的格式，也是不小的工作量。
官网下载中心; http://annovar.openbioinformatics.org/en/latest/user-guide/download/
# 进入软件安装目录运行
mkdir -p ~/biosoft/ANNOVAR/
nohup ./annotate_variation.pl -downdb -webfrom annovar gnomad_genome --buildver hg38 humandb/ >down.log 2>&1 &
## 或者在任意地方运行：
perl ~/biosoft/ANNOVAR/annovar/annotate_variation.pl -downdb -webfrom annovar gnomad_genome --buildver hg38 ~/biosoft/ANNOVAR/annovar/humandb/
## 可以看到需要选择参考基因组版本，这里一般是 hg38, 很少看到hg19了。
## 然后是 数据库名字，在 http://annovar.openbioinformatics.org/en/latest/user-guide/download/ 全部列出。
## 其实就是相当于执行了下面的命令而已。
wget -t 1 -T 30 -q -O hg38_gnomad_genome.txt.gz http://www.openbioinformatics.org/annovar/download/hg38_gnomad_genome.txt.gz
## 这个文件解压后有16G，注意云服务器的空间。
#
## 通常要下载的数据库非常多，包括：dbSNP,ExAC,ESP6500,cosmic,gnomad,1000genomes,clinvar,gwas等等
# 当然，注释上，并不等价于理解它们，这是很大的工作量。
如果全部下载，我挑选了一些大文件展示，大概效果如下：
1.5G Aug 26 2015 hg38_AFR.sites.2015_08.txt
3.2G Aug 26 2015 hg38_ALL.sites.2015_08.txt
1.1G Aug 26 2015 hg38_AMR.sites.2015_08.txt
5.9G Feb 8 2017 hg38_avsnp147.txt
13G Sep 30 2017 hg38_avsnp150.txt
28G Feb 22 2017 hg38_dbnsfp33a.txt
16G Mar 12 2017 hg38_gnomad_genome.txt
2.4G Sep 20 2016 hg38_hrcr1.txt
11G Mar 25 2018 hg38_intervar_20180118.txt
12G Sep 20 2016 hg38_kaviar_20150923.txt
11G Sep 20 2016 hg38_ljb26_all.txt
2.9G Nov 5 2016 hg38_mcap.txt
12G Apr 5 2018 hg38_regsnpintron.txt
2.4G Dec 6 2016 hg38_revel.txt
所以哦，能做这个分析的，计算机资源肯定是有一定保障的。
注释自己的vcf文件数据
## 可以看到 4 大参数里面，就有3个参数是注释用的。
Arguments to download databases or perform annotations
--downdb download annotation database
--geneanno annotate variants by gene-based annotation (infer functional consequence on genes)
--regionanno annotate variants by region-based annotation (find overlapped regions in database)
--filter annotate variants by filter-based annotation (find identical variants in database)
其实并不一定需要vcf文件，软件需要的只是染色体加上坐标即可，对于我们的vcf格式的变异文件， 软件通常会进行一定程度的格式化之后再进行注释 。这里的注释主要有三种方式，分别是：
-
基于基因的注释,
exonic,splicing,ncRNA,UTR5,UTR3,intronic,upstream,downstream,intergenic，使用geneanno子命令。 -
基于区域的注释,
cytoBand,TFBS,SV,bed,GWAS,ENCODE,enhancers, repressors, promoters，使用regionanno子命令。只考虑位点坐标 -
基于数据库的过滤,
dbSNP,ExAC,ESP6500,cosmic,gnomad,1000genomes,clinvar使用filter子命令。 考虑位点坐标同时关心突变碱基情况。
基于基因的注释是想搞清楚自己的vcf文件记录的突变位点，是否在基因组的哪些功能元件（主要是外显子）上面。
而且基于区域的注释
最后基于数据库的过滤，就很容易理解。
基于基因的注释
现在下载ANNOVAR最新版默认自带了hg19的数据库，所以可以很方便的注释，如果是hg38，可能得自己下载后再进行注释。
下面例子的 H3F3A.vcf 文件是 我在生信技能树直播我的基因组的时候做出来的测试数据，大家就理解为一个vcf即可。
~/biosoft/ANNOVAR/annovar/convert2annovar.pl -format vcf4old H3F3A.vcf >H3F3A.annovar 2>/dev/null
~/biosoft/ANNOVAR/annovar/annotate_variation.pl -buildver hg38 --geneanno --outfile H3F3A.anno H3F3A.annovar ~/biosoft/ANNOVAR/annovar/humandb/
~/biosoft/ANNOVAR/annovar/annotate_variation.pl -buildver hg38 --dbtype knownGene --geneanno --outfile H3F3A.anno H3F3A.annovar ~/biosoft/ANNOVAR/annovar/humandb/
一般人做到这里就结束了，因为信息量足够了，主要是对 refGene.exonic_variant_function 文件的理解，就是里面记录的变异情况：
前半部分内容如下，主要是关于每个变异位点位于基因的哪个外显子上面，造成了何种改变：
line380853 nonsynonymous SNV BRCA1:NM_007297:exon14:c.A4696G:p.S1566G
line380862 synonymous SNV BRCA1:NM_007297:exon11:c.T4167C:p.S1389S
line380863 nonsynonymous SNV BRCA1:NM_007297:exon11:c.T4070G:p.L1357R
line380866 nonsynonymous SNV BRCA1:NM_007297:exon9:c.A3407G:p.K1136R
line380867 nonsynonymous SNV BRCA1:NM_007297:exon9:c.A2972G:p.E991G
line380868 nonsynonymous SNV BRCA1:NM_007297:exon9:c.C2471T:p.P824L
line380869 nonsynonymous SNV BRCA1:NM_007297:exon9:c.T2425C:p.Y809H
line380870 synonymous SNV BRCA1:NM_007297:exon9:c.T2170C:p.L724L
line380871 synonymous SNV BRCA1:NM_007297:exon9:c.C1941T:p.S647S
line380872 nonsynonymous SNV BRCA1:NM_007297:exon9:c.G683A:p.G228D
line380873 nonsynonymous SNV BRCA1:NM_007297:exon9:c.G607A:p.E203K
line380875 nonsynonymous SNV BRCA1:NM_007297:exon8:c.G466C:p.E156Q
line380876 nonsynonymous SNV BRCA1:NM_007297:exon7:c.G430A:p.V144I
line380887 synonymous SNV BRCA1:NM_007298:exon2:c.G114A:p.K38K
通常，我们会关心 nonsynonymous ，因为它里面记录的氨基酸改变了。
但是，很明显，这样是不够的，我们只知道自己的样本是否在这个基因的这个位点突变了，而且也的确造成了氨基酸的改变，但这并不意味着这个突变是有害的，致病的。
首先需要做的就是判断这个位点是否在其它科研文献里面报道过，其他科学家使用实验验证的手段证实了这个位点是致病突变。
然后可以看这个突变位点，在一系列健康人群数据库里面是否也被发现过，一般来说人群频率很高的位点，通常就不会是致病位点。
基于区域的注释
这个功能使用频率不高，比如如果有染色体坐标对应染色体区段的数据库，就可以注释上去，每个位点属于染色体的什么区段。
cytoBand 1p36.33 chr1 69155 69155 T C hom 159.65 178 26.96 13.30
cytoBand 1p36.33 chr1 69270 69270 A G hom 59336.13 1846 27.01 32.23
cytoBand 1p36.33 chr1 69511 69511 A G hom 1279869.08 46442 46.09 27.60
cytoBand 1p36.33 chr1 69610 69610 C T hom 4538.56 6585 44.70 13.12
可以看到是一个很简单的注释，比基于基因的注释要简单太多了，仅仅是对原来的位点注释一列信息即可，就是属于哪个区域，并不需要像注释基因那些判断属于基因的什么元件，也不需要判断是造成了什么样的改变。
或者有ENCODE计划的enhancer等调控原件的坐标信息，也可以注释上去，至少在我平时数据处理经验来看，使用频率很低，当然，也欢迎生信技能树的所有其他小伙伴提出自己的意见。
基于数据库的过滤
这个使用频率非常高，而且通常是结合多个数据库信息一起过滤。
重点是检查那些后缀为 dropped 的文件，Known variants will be written to the dropped file together with allele frequencies. 如下：
238104 tmp.hg38_ALL.sites.2015_08_dropped
15575 tmp.hg38_clinvar_20170905_dropped
7527 tmp.hg38_cosmic70_dropped
172432 tmp.hg38_exac03_dropped
281449 tmp.hg38_gnomad_genome_dropped
可以看到，总共的48万位点，其中有13万是在千人基因组计划出现的，有17万是在EXAC数据库出现，但是只有区区7527个位点是在COSMIC数据库出现，，在clinvar数据库的，有15575位点。
最原始的想法，通常是找到某个基因被注释到外显子区域的nonsynonymous突变位点，然后如果其被clinvar数据库记录，那就说明找到了有证据支持的致病位点。
不过，值得一提的是 clinvar 数据库需要经常更新哦。
~/biosoft/ANNOVAR/annovar/annotate_variation.pl -downdb clinvar_20180603 humandb -buildver hg38
# 72M Jul 9 2018 hg38_clinvar_20180603.txt
~/biosoft/ANNOVAR/annovar/annotate_variation.pl -filter \
-buildver hg38 -dbtype clinvar_20180603 --outfile tmp for_annovar.input \
~/biosoft/ANNOVAR/annovar/humandb/
联合注释
参考：https://davetang.org/wiki2/index.php?title=ANNOVAR
允许在同一命令中用输出的特定顺序来对多个注释类型进行 自定义选择（custom selection）, 下面是一个最常用的联合注释情况：
dir=/home/jianmingzeng/biosoft/ANNOVAR/annovar
db=$dir/humandb/
ls $db
perl $dir/table_annovar.pl \
-buildver hg38 \
for_annovar.input $db \
-out test \
-remove -protocol \
refGene,clinvar_20180603 \
-operation g,f \
-nastring NA
## 这里只做两个数据库注释，举例说明。
‘--protocol’ 选项后跟注释来源数据库的准确名称，有多少个数据库，就对应多少种描述，即随后的operation方案。
‘--operation’ 选项后跟注释的类型: ‘g’ 表示基于基因的注释（gene-based annotation）、‘r’ 表示基于区域的注释（region-based annotation） 、‘f’ 表示基于筛选子的注释（ filter-based annotation）.
‘--outfile’ 选项是指定输出文件的前缀
关键步骤（ CR ITICAL STEP）： 1、确保注释数据库的名称正确并且是按你想要在输出文件中显示的顺序排列的； 确保 ‘--operation’指定的注释类型顺序和‘--protocol’指定的数据库顺序是一致的；确保每个protocal名称或注释类型之间只有一个逗号，并且没有空白。
通常更复杂，如下：
bin=/home/Annovar/
db=/home/Annovar/humandb
perl $bin/table_annovar.pl for_annovar.input $db -buildver hg38 -out tmp \
-protocol refGene,cytoBand,esp6500siv2_all,exac03,gnomad_genome,cosmic70,1000g2015aug_all,clinvar_20170905 \
-operation g,r,r,f,f,f,f,f -nastring NA
R包接口
值得一提的是我们生信技能树的技术骨干，也开发了 该软件的R包接口：https://github.com/JhuangLab/annovarR
annovarR: a variant annotation and visualization system based on R and Shiny framework https://life2cloud.com/tools/annovarR
致病基因分析策略
首先需要对 Pathogenic germline variants. 定义有清晰的理解。
然后需要看看权威文章通常是如何分析的，比如：12.353 Nat Commun. 2016 May
We used the ClinVar database15 to identify pathogenic germline mutations, using only those SNVs and indels recorded as being ‘probable-pathogenic’ or ‘pathogenic’, and ‘germline’, ‘inherited’, ‘paternal’, ‘maternal’, ‘biparental’ or ‘uniparental’. Variants classified as ‘germline’ by the unpaired pipeline were classified as ‘pathogenic’ using the ClinVar annotation, unless they were also present at allele frequencies >1% in the 1,000 Genomes resource.
## 挑选那些被clinVar数据库记录的致病位点，然后剔除那些人群频率较高的。
In addition, we classified SNVs absent in ClinVar but present in between one and six (1%) normal samples as ‘pathogenic’ if they were either inactivating (truncating or affecting splice sites), or identified as being ‘deleterious’ or ‘damaging’ by Provean64 Pathogenic indels present in one to six normal samples but absent from ClinVar were classified as ‘pathogenic’ if they were predicted to disrupt the reading frame or disrupt a splice junction.
## 也挑选那些并没有被clinVar数据库记录的，但是可能被预测会致病的
再比如2018年NC的一篇日本人群队列文章：Germline pathogenic variants of 11 breast cancer genes in 7,051 Japanese patients and 11,241 controls 只关心了11个基因，文章关于Pathogenic germline variants的分析方法非常值得学习。
Therefore, to maintain consistency of variant annotation across the 11 hereditary breast cancer genes, we decided to use the American College of Medical Genetics and Genomics and the Association for Molecular Pathology (ACMG/AMP) guidelines3 to assess all 11 genes analyzed in this study.
## 全套遗传咨询师流程
这个过程是非常耗费精力，需要公司团队长时间研发，通常是全套遗传咨询师流程，所以不方便分享太细致的知识点，抱歉。
而且真的是非常复杂，比如仅仅是评估一下变异位点对基因蛋白质功能影响就有：SIFT, PolyPhen2 HDIV, PolyPhen2 HVAR, LRT, MutationTaster, MutationAssessor, FATHMM, MetaSVM, MetaLR, VEST, CADD, GERP++, PhyloP and SiPhy 这些软件算法。
网友评论