Annovar各程序的功能,自建avdb,构建索引

ANNOVAR的程序模块(本人目录)
├── annotate_variation.pl   //annovar主程序,功能包括下载数据库,三种不同的注释
├── annovar_index.pl   //index构建的程序,对于染色体是第一列的文件有效,修改网上程序
├── coding_change.pl   //用来推断蛋白质序列的程序
├── convert2annovar.pl   //可以转换各种不同格式的文件到.avinput,适合后期输出txt文件
├── example   //例子和一些有用的文件
├── humandb   //人类注释所用到的数据库的文件夹目录
├── idx_annovar.pl   //index构建的程序,对于染色体是第一列的文件有效,自写
├── index_annovar.pl   //来自annovar作者王凯老师的index构建程序
├── prepare_annovar_user.pl   //重新格式化并准备annovar注释数据库,比如COSMIC、CLINVAR
├── retrieve_seq_from_fasta.pl   //从全基因组fasta文件重新格式化特定基因组位置的序列
├── summarize_annovar.pl   //自动运行的流程,注释突变文件并获得以逗号分隔的文件中
├── table_annovar.pl   //注释程序,可一次性完成三种类型的注释,输出可以是csv/txt
└── variants_reduction.pl   //可用来更灵活地定制过滤注释流程

PS:

  1. 上述有些程序来自于很早以前的收藏,还有两个来自于网上和自写;
  2. convert2annovar.pl 支持的格式:pileup, cg, cgmastervar, gff3-solid, soap, maq, casava, vcf4, vcf4old, rsid, maq, annovar, annovar2vcf, bed, region, transcript.
自建avdb

因为一些特定的注释数据库在annovar中是不一定提供的,那么如果想和annovar一起注释,要如何操作呢?举个例子,想注释CCRs(constrained coding regions)到突变结果中,这个数据库是基于区域的有害性预测打分工具。首先下载数据文件,还有X染色体

CCRs的数据库构建

先看一下下载的文件的格式:



这个数据库中,其实有用的自由第四列,但是前三列代表region,也是需要的,所以,只需要把钱四列截取出来:

gzip -dc ccrs.autosomes.v2.20180420.bed.gz | grep -v '^#' cut -f 1-4 > hg19_CCRs.txt

这一步的目的有三个,一是留取有用的信息,二是减小数据库的大小,三是从bed文件转为txt的文件后缀,annovar的数据库的命名格式,一般是hg19_dbname.txt/hg18_dbname.txt/hg38_dbname.txt等等,便于程序识别读取数据库。
单单修改成这样是不足以用于注释,在注释的主程序中,对于dbtype也有一个设定,如果想使用CCRs,这里有两种方法,第一种修改annotate_variation.pl这个主程序,第二种就是修改这个文件的格式。
先看一下主程序关于dbtype的选择,也就是文件格式的匹配,在annotate_variation.pl的第2998行,有一个预设,如下图:



上图中有三个红框,可以发现关于dbtype的预设有很多匹配,第二个红框里面是我修改了,加入了ccrs这个匹配,因为这个文件格式其实和和CCRs的格式是一致,position是1、2、3列,第四列是需要的输出列,这就可以开始注释了;
第二种方法就是我们看第三个红框,这个红框是对于其他所有的未知的格式来进行匹配,我们先看一下源码中的格式类型,第2、3、4列是position,第五列是需要的输出列,按照这个格式,我们可以修改获得的hg19_CCRs.txt这个文件,只需要给这个文件加第一列,这里我选择原文件的最后一列:

gzip -dc ccrs.autosomes.v2.20180420.bed.gz | cut -f 13 grep -v "unique" > first.txt
paste first.txt hg19_CCRs.txt > hg19_ccrs.txt

修改完的文件就可以用于注释:

perl convert2annovar.pl -format vcf4old samp.raw.vcf > samp.avinput
perl table_annovar.pl samp.avinput humandb/ -outfile samp.anno -buildver hg19 -protocol ccrs -operation r -nastring . -remove --otherinfo

这里的数据库不大,可以构建hg19_ccrs.txt的索引,也可以不用构建;构建索引的程序可以给王凯老师发邮件获得,我以前不知道,自己写了一个能用的idx_annovar.pl:

#!/usr/bin/perl -w
use strict;
use Getopt::Long;
use File::Basename;
use POSIX qw(strftime);
my ($in,$bins,$help);
GetOptions
(
        "i|in=s"=>\$in, "b|bins=i"=>\$bins,     "help|?"=>\$help,
);
my $usage=<<INFO;
Usage:
        perl $0 [options]
Options:
                -i <file> <in.input the database file for indexing>
                -b <INT> <bin.size the size of the index's bin>
INFO
die $usage if ($help || !$in || !$bins);
open IN,"$in" || die $!;
my $out = ${in} . ".idx";
my $size_in = -s $in;
my $len = 0;
my %idxLast=();
my %idxFirst=();
my $key_bin = 0;
$idxLast{$key_bin} = 0;
while(<IN>){
    my ($chr, $start, $end) = (split /\t/, $_)[0, 1, 2];
    $chr =~ s/^chr//i;
    my $bin_start =$bins * int($start/$bins);
    $len += length $_;
    $key_bin = "$chr-$bin_start";
    if(exists $idxLast{$key_bin}){
        $idxLast{$key_bin} = $len;
    }
    else{
        $idxFirst{$key_bin} = $len - length $_;
        $idxLast{$key_bin} = $len;
    }
}
close IN;

open OUT,"> $out" || die $!;
print OUT "#BIN\t$bins\t$size_in\n";
foreach(sort keys %idxFirst){
    my ($chr, $bin_start) = split /\-/, $_;
    next unless($chr !~ /^#/);
    print OUT "$chr\t$bin_start\t$idxFirst{$_}\t$idxLast{$_}\n";
}
close OUT;

对于一般的数据库,就是Chr、Start、End、Score的类型的文件没有问题,可以使用;

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,829评论 1 331
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,603评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,846评论 0 226
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,600评论 0 191
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,780评论 3 272
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,695评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,136评论 2 293
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,862评论 0 182
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,453评论 0 229
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,942评论 2 233
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,347评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,790评论 2 236
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,293评论 3 221
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,839评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,448评论 0 181
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,564评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,623评论 2 249

推荐阅读更多精彩内容