参考文章:
http://annovar.openbioinformatics.org/en/latest/
https://cloud.tencent.com/developer/article/1054541
1.下载ANNOVAR软件包
需要使用.edu邮箱进行注册,然后作者会发最新的软件包到你邮箱
解压之后是有6个perl脚本,和两个文件夹humandb/, example/
drwxr-xr-x 3 1001 fqq1024 4.0K Jun 19 18:32 humandb
drwxr-xr-x 2 1001 fqq1024 4.0K Apr 17 2018 example
-rwxr-xr-x 1 1001 fqq1024 217K Apr 17 2018 annotate_variation.pl
-rwxr-xr-x 1 1001 fqq1024 27K Apr 17 2018 coding_change.pl
-rwxr-xr-x 1 1001 fqq1024 167K Apr 17 2018 convert2annovar.pl
-rwxr-xr-x 1 1001 fqq1024 19K Apr 17 2018 retrieve_seq_from_fasta.pl
-rwxr-xr-x 1 1001 fqq1024 39K Apr 17 2018 table_annovar.pl
-rwxr-xr-x 1 1001 fqq1024 22K Apr 17 2018 variants_reduction.pl
2.获取wegene 测序得到的SNP
cat wegene/20190615wegene_core_results.txt |less -SN
1 # This data file generated by WeGene at: Sat, 15 Jun 2019 12:04:02
2 #
3 # This file contains the genotype called by WeGene with our
4 # internal quality control pipeline. The low quality sites
5 # were discarded. If we could not determine the genotypes,
6 # it will be —-. As such, the call rate and accuracy may not
7 # be one hundred percent. So this data is suitable only for
8 # research, educational, and informational use and not for
9 # medical or other use.
10
11 # This text file is a list of your data which are TAB-
12 # separated. Each line corresponds to a single SNP or short
13 # InDel (insertion or deletion).
14 # For each SNP or short InDel, we provide its identifier (an
15 # rsid or an internal id), its location on the reference human
16 # genome (human assembly build 37, GRCh37) and the genotype
17 # call oriented with respect to the plus strand on the human
18 # reference sequence. For consistency, the genotypes are always
19 # two base pairs, including hemizygous calls.
20
21 # Please note, as our ability to call genotypes improves, it is
22 # possible that your data may be slightly different at different
23 # times.
24
25 # rsid chromosome position genotype
26 rs9701055 1 565433 CC
27 rs9651229 1 567667 CC
提取SNP信息,也就是把空格和#开头的删除
cat wegene/20190615wegene_core_results.txt |grep -v "#"|grep -v "^$"|less -SN
1 rs9701055 1 565433 CC
2 rs9651229 1 567667 CC
3 rs9701872 1 568208 TT
4 rs11497407 1 568527 GG
cat wegene/20190615wegene_core_results.txt |grep -v "#"|grep -v "^$" > mysnpID.txt
3. 将SNP文件转换为ANNOVAR的输入文件.avinput
perl convert2annovar.pl -format rsid example/snplist.txt -dbsnpfile humandb/hg19_snp138.txt > snplist.avinput
这个时候你可能需要下载hg19_snp138数据至本地的humandb/
参考http://annovar.openbioinformatics.org/en/latest/user-guide/download/
#约12G大小
perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar snp138 humandb/
4. 基于GWAS对SNP进行注释(Identify variants reported in previously published GWAS)
#下载GWAS数据60M
perl annotate_variation.pl -build hg19 -downdb gwasCatalog humandb/
#注释
perl annotate_variation.pl -regionanno -build hg19 -out ex1 -dbtype gwasCatalog snplist.avinput humandb/
##注释很快3min之内完成
ls -lht annovar
-rw-r--r-- 1 root root 3.0M Jun 20 07:45 ex1.hg19_gwasCatalog
-rw-r--r-- 1 root root 45M Jun 20 07:34 snp.avinput
drwxr-xr-x 3 1001 fqq1024 4.0K Jun 19 18:32 humandb
drwxr-xr-x 2 1001 fqq1024 4.0K Apr 17 2018 example
-rwxr-xr-x 1 1001 fqq1024 217K Apr 17 2018 annotate_variation.pl
-rwxr-xr-x 1 1001 fqq1024 27K Apr 17 2018 coding_change.pl
-rwxr-xr-x 1 1001 fqq1024 167K Apr 17 2018 convert2annovar.pl
-rwxr-xr-x 1 1001 fqq1024 19K Apr 17 2018 retrieve_seq_from_fasta.pl
-rwxr-xr-x 1 1001 fqq1024 39K Apr 17 2018 table_annovar.pl
-rwxr-xr-x 1 1001 fqq1024 22K Apr 17 2018 variants_reduction.pl
5. 查看GWAS注释结果
cat annovar/ex1.hg19_gwasCatalog |less -SN
1 gwasCatalog Name=Pancreatic cancer chr1 894573 894573 G A rs13303010
2 gwasCatalog Name=Body mass index chr1 1005806 1005806 C T rs3934834
3 gwasCatalog Name=IgG glycosylation chr1 1079198 1079198 T C rs11260603
4 gwasCatalog Name=Blood protein levels chr1 1162435 1162435 C A rs3766186
5 gwasCatalog Name=Inflammatory bowel disease,Ulcerative colitis chr1 1247494 1247494 T C rs12103
或者
6. 基于基因对SNP进行注释
#ANNOVAR压缩包里已经自带,太好了
perl annotate_variation.pl -downdb -buildver hg19 -webfrom annovar refGene humandb/
#基于hg19基因进行注释,注释很快2min之内完成
perl annotate_variation.pl -out ex1 -build hg19 snp.avinput humandb/
ls -lht annovar
-rw-r--r-- 1 root root 1.1K Jun 20 13:18 ex1.log
-rw-r--r-- 1 root root 11M Jun 20 13:18 ex1.exonic_variant_function
-rw-r--r-- 1 root root 16K Jun 20 13:18 ex1.invalid_input
-rw-r--r-- 1 root root 85M Jun 20 13:18 ex1.variant_function
7.查看基因的注释结果
有两种注释结果,具体可参考http://annovar.openbioinformatics.org/en/latest/user-guide/gene/
#全部的gene variant
cat annovar/ex1.variant_function |less -SN
#外显子上的gene variant
cat annovar/ex1.exonic_variant_function |less -SN
#查看多少种gene variant
cat ex1.exonic_variant_function |cut -f2|cut -d" " -f1|sort |uniq -c |sort -nr
49585 nonsynonymous
12001 synonymous
2743 stopgain
541 unknown
385 frameshift
46 stoploss
26 nonframeshift