根據相似性原理,序列相似,功能相似,所有功能注釋無非是用比對工具將輸入序列比對到數據庫序列,再將輸入ID對應數據庫ID,進一步對應到功能條目的關系。
數據庫要么建到本地,要么聯網調用API,一般的軟件或包做注釋都是通過聯網來獲得,或者調用依賴的一些專門注釋的包(文件較大)。工業生產中,一般需要構建本地數據庫。
如果不對原始數據庫按物種或其他分類來進行拆分的話,整個數據庫會很大,比對和注釋消耗的時間和資源都會很大,顯然不經濟,而且也會有一些假陽性的結果。比如將人特有的功能比對到了小鼠上,客戶無法結果。所以分庫是很有必要的,只是怎么分以及分多大的問題。
KEGG本地庫文件
- 序列數據庫文件
如kegg_all_clean.fa
- ko系列文件(ko與其他ID的對應關系),ko與不同類型數據庫
我們這里要用到的是ko和geneID的對應關系,其他數據庫類似
-
物種文件
misc下的taxonomy文件,按物種分庫的依據。
-
map目錄,通路圖。
每條通路有三個文件:png是通路圖,html是網頁通路,conf是通路的配置
conf文件內容
-
map_title.tab文件,是通路的三個層級
-
ko_map.tab文件,是K與通路的全部物種對應文件
是聯系注釋結果之間對應關系的必需文件。
-
komap目錄下,是各個物種(三個字母縮寫)的通路圖(png)及其配置(conf),以及該物種對應的通路。
如人的komap/hsa目錄:
當然也可以不細分到單物種,可以划分物種大類,如動物、植物等,相對應地文件animal_ko_map.tab、plant_ko_map.tab
利用上面的這些文件,其實我們就可以進行KEGG Pathway功能注釋了,即存在這樣的關系:蛋白——序列ID(基因)——K號——ko(pathway)——Level1-3——通路圖。這樣得到的通路圖,都是map開頭,即reference pathway;如果是物種特異通路,即ko開頭,則用komap目錄結果。KEGG的5種通路類型等基礎知識這里不講,不懂可去查。
如果要按物種進行拆庫,則需要將上面的文件都按物種進行分類,使用是一樣的。
KEGG數據庫非常龐大,除了Pathway,genes等數據庫外,還有很多其他的文件,比如:
-
links目錄,pathway與其他ID的對應關系
如pathway_ko.list
-
ligand目錄,即配體數據庫,不做介紹。
按物種拆分KEGG數據庫
1.獲得物種分類信息
按物種拆分可大可小:最大就是原始庫,最小就是單一物種,中間可以按不同分類來拆。關鍵取決於你的輸入序列是什么成分,當然不大不小恰好能全部包含是最理想的分庫結果。比如taxonomy文件(misc目錄下)格式是:
我們可以按Eukaryotes、Animals、Vertebrates、Mammals中的任一個層級來分。也可以自定義分類,將不同物種添加到一起進行歸類。
這里寫一個簡單腳本來用上面文件中的第二層級物種來進行數據庫分庫:
#!urs/bin/perl
open F , $ARGV[0];
while (<F>){
chomp;
if (/^## (.+)/){
$spec=$1;
open OUT, ">>$spec.specie.xls";
}
elsif (!/^#/){
@aa=split/\t/,$_;
print OUT "$spec\t$aa[1]\t$aa[3]\n";
}
}
拆分后得到Animals.specie.xls、 Archaea.specie.xls、 Bacteria.specie.xls、 Fungi.specie.xls、 Plants.specie.xls、 Protists.specie.xls
等系列分類文件。如Animals.specie.xls
文件如下:
2.獲得物種分類的序列信息並建庫
從全部物種的原始序列數據庫中拆分出以上分類物種的序列,編寫如下get_fasta.pl
腳本:
#!/usr/bin/perl
use strict;
use warnings;
use diagnostics;
unless(@ARGV>=2){
print "perl $0 list.txt db.fa specie.fa\n";
exit(1);
}
my %hash;
open L,"$ARGV[0]" or die "$!\n";
my %list;
while(<L>){
chomp;
my @aa=split/\t/,$_;
$list{$aa[1]}='';
}
close L;
my $num_need = scalar keys %list;
print("begin fetch $num_need sequence...\n");
open F,"$ARGV[1]" or die "$!\n";
my %seq;
while(<F>){
LINE: #if(m/^>([^\|]+)/){
if(/^>([^:]+):([^\s]+)/){
chomp;
my $acc = $1;
my $line=$_;
my $idd=$1.':'.$2;
$hash{$idd}=$line;
next unless exists $list{$acc};
while(<F>){
goto LINE if /^>/;
s/[^a-z]//gsi;
$seq{$idd} .= $_;
}
last;
}
elsif (/^>(12122[^\s]+)(.*)/){
chomp;
my $acc = $1;
$hash{$1}=$_;
next unless exists $list{$acc};
while(<F>){
goto LINE if /^>/;
s/[^a-z]//gsi;
$seq{$acc} .= $_;
}
last;
}
}
close F;
open O,">$ARGV[2]" or die "$!\n";
foreach my $acc (keys %seq){
print O "$hash{$acc}\n$seq{$acc}\n";
}
my $yesn = scalar keys %seq;
my $non = $num_need - $yesn;
if($non>=1){
print "have $non sequence not found!\n";
}else{
print "you have successfully got $yesn sequence!\n";
}
close O;
獲得分類后的序列:
perl get_fasta.pl Animals.specie.xls kegg_all_clean.fa animals.fa
獲得分類后的數據庫后,可用blast/diamond等軟件進行建庫,以便進行輸入序列的比對。
3.獲得物種分類的K-ko對應文件
從全部物種的ko_map.tab
文件(K與通路的對應關系文件)中獲取物種分類后的子文件。編寫腳本get_species_komap.pl
:
#!/usr/bin/perl
=pod
this script is subsplit species komap
perl $0 species.xls ko_genes.list ko_map.tab species_ko_map.tab
=cut
my $spe=shift;
my $ko_gene=shift;
my $ko_map=shift;
my $map=shift;
my (%spe,%ko);
open F,"<$spe";
while(<F>){
chomp;
my @F=split/\t/,$_;
$spe{$F[1]}=1;
}
close F;
open F,"<$ko_gene";
while(<F>){
chomp;
my @F=split/\t/,$_;
my $ko=(split/:/,$F[0])[1];
my $spe=(split/:/,$F[1])[0];
if(exists $spe{$spe}){
$ko{$ko}=1;
}
}
close F;
open F,"<$ko_map";
open O,">$map";
while(<F>){
chomp;
my @F=split/\t/,$_;
if(exists ($ko{$F[0]})){
print O "$_\n";
}
}
close F;
close O;
得到animal_ko_map.tab
文件,其他分類物種也是類似。
拆分子庫后比對獲得該子庫中的功能信息,后續注釋的數據處理其實和不分庫時是一樣的,都是一些文本的格式轉換以及可視化。
比如我們可將KEGG數據庫拆分:動物animal.fa、植物plant.fa、真菌fungi.fa、真核eukaryotes.fa、原核prokaryote.fa、原生生物microorganism.fa,以及包含原核、真菌和原生生物三種組合的微生物庫other.fa,除動植物、真菌、原核、原生之外但在KEGG數據庫中的其他生物unknow.fa。
后面比對注釋時只需設置物種參數即可。