Desde hace ya tiempo hay más secuencias genómicas disponibles que las que podemos recordar y a menudo nos encontramos con nombres de organismos que no podemos ubicar, no sabemos si son bacterias aguas termales o artrópodos, por ejemplo. El siguiente programa puede ayudarnos precisamente a colocar dentro de la clasificación taxonómica vigente del
NCBI un nombre de especie, o en general, cualquier taxón de nombre aceptado. Antes de ejecutarlo debes acceder al servidor FTP en
ftp://ftp.ncbi.nih.gov/pub/taxonomy y descomprimir los archivos de nombres y nodos, como se explica en la cabecera del código Perl, haciendo apuntar la variable global $TAXONOMYPATH al directorio donde se encuentren. Una vez hechos estos ajustes el código ya está listo para usarse desde el terminal. Por ejemplo, si consultamos la taxonomía de
Oryza sativa obtendremos esta respuesta, que podemos comparar con la figura, tomada de
www.biomedcentral.com/1471-2164/8/283:
$ ./prog6.pl Oryza sativa
# input taxon: Oryza sativa
# selected taxonID: 4530 \\
URL: http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=4530
# variants:
> Oryza sativa
> Oryza sativa L.
# lineage : Oryza sativa, Oryza, Oryzeae, Ehrhartoideae, BEP clade, Poaceae,
Poales,commelinids, Liliopsida, Magnoliophyta, Spermatophyta, Euphyllophyta,
Tracheophyta, Embryophyta, Streptophytina, Streptophyta, Viridiplantae,
Eukaryota, cellular organisms
# rank : species, genus, tribe, subfamily, no_rank, family, order, subclass,
class, no_rank, no_rank, no_rank, no_rank, no_rank, no_rank, phylum, kingdom,
superkingdom, no_rank,
Éste es el código:
#!/usr/bin/perl -w
# prog6
# programa que devuelve informacion taxonomica para un taxon de entrada,
# que se apoya en la version texto de la clasificacion Taxonomy del NCBI,
# que se pueden descargar de ftp://ftp.ncbi.nih.gov/pub/taxonomy , en
# particular el fichero taxdump.tar.gz, que debes descomprimir para obtener
# los dos ficheros que se utilizan aqui
# Bruno Contreras Moreira, Junio de 2010
use strict;
$|=1;
if(!$ARGV[0])
{
die "# usage: $0 <taxon name from superkingdom to species, such as:\n".
" 'Arabidopsis thaliana' or 'Viridiplantae' but not 'thaliana'>\n";
}
my $VERBOSE = 0 ;
my $TAXONOMYPATH = '/path/taxonomy/'; # directorio de taxdump.tar.gz
my $NAMESFILE = $TAXONOMYPATH . 'names.dmp';
my $NODESFILE = $TAXONOMYPATH . 'nodes.dmp';
my $NCBIURL = 'http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=';
###################################################################
my ($n_of_taxa,$taxonID,$taxon,$ID,$name,$rank,$parentID) = (0,'','');
my (@names,%node,@lineage,@sorted_lineage,%lineage_name);
my ($lineage_names,$lineage_ranks) = ('','');
## 1) captura el nombre del taxon de entrada, como Arabidopsis thaliana, y
## sustituye caracteres '_' por espacios si los hubiera
if(scalar(@ARGV) > 1){ @ARGV = @ARGV[0,1]; } # corta tercera palabra < especie
else{ @ARGV = split(/_/,$ARGV[0]); @ARGV = @ARGV[0,1]; }
$taxon = join(' ',@ARGV);
$taxon =~ s/^\s+//g; # elimina espacios iniciales y finales
$taxon =~ s/\s+$//g;
print "# input taxon: $taxon\n\n";
## 2) busca el taxonID del taxon de entrada
open(NAMES,$NAMESFILE) || die "# $0 : cannot open $NAMESFILE\n";
while(<NAMES>)
{
next if(!/$taxon/i);
#3702 | Arabidopsis thaliana | | scientific name
if(/^(\d+)\s+\|\s+(.*?)\|/)
{
($ID,$name) = ($1,$2);
$name =~ s/^\s+//g;
$name =~ s/\s+$//g;
# comprueba que el taxon buscado coincide con el principio de $name,
# evita aceptar cadenas que en realidad coinciden con la especie,
# que a veces se pueden encontrar en generos muy distintos
next if(index($name,$taxon) != 0);
# si el taxon buscado contiene una sola palabra evita
# aceptar cadenas de especies, que tendran dos palabras
next if($taxon !~ /\s+/ && $name =~ /\s+/);
if(!$n_of_taxa || $ID eq $taxonID) # guarda primer taxonID
{
$taxonID = $ID;
push(@names,$name);
$n_of_taxa++;
if($VERBOSE){ print }
}
else
{
if($VERBOSE)
{
# muestra otras posibles apariciones del taxon
print "# skip taxon: $ID\n";
}
else{ last } # termina la lectura del fichero
}
}
}
close(NAMES);
if(!$n_of_taxa){ die "# cannot find taxon $taxon in $NAMESFILE , exit\n"; }
else
{
print "# selected taxonID: $taxonID URL: $NCBIURL$taxonID\n";
print "\n# variants:\n";
foreach $name (@names){ print "> $name\n"; }
}
## 3) busca taxon padre del taxonID encontrado
open(NODES,$NODESFILE) || die "# $0 : cannot open $NODESFILE\n";
while(<NODES>)
{
#132596 | 79512 | genus...
#tax_id -- node id in GenBank taxonomy database
#parent tax_id -- parent node id in GenBank taxonomy
#rank -- rank of this node (superkingdom, kingdom, ...)
next if(/forma/ || /varietas/); # evita rangos < especie
if(/^$taxonID\t/)
{
my @data = split(/\t|\t/,$_);
$parentID = $data[2];
$rank = $data[4];
if($VERBOSE){ print }
last;
}
}
close(NODES);
if(!$parentID){ die "\n# cannot find taxon '$taxon' in $NODESFILE , exit\n"; }
## 4) lee arbol taxonomico del nivel de especie hacia arriba
## (son en 2010 35MB, caben en cualquier RAM)
open(NODES,$NODESFILE) || die "# $0 : cannot open $NODESFILE\n";
while(<NODES>)
{
next if(/forma/ || /varietas/);
if(/^(\d+)\t\|\t(\d+)\t\|\t(\S+)/)
{
$node{$1}{'p'} = $2; # p = parent node
$node{$1}{'r'} = $3; # r = rank
}
}
close(NODES);
## 5) reconstruye linaje de $taxonID en terminos de IDs
# 5.1) primero guarda datos de $taxonID
push(@lineage,$taxonID);
$lineage_name{$taxonID} = $rank;
$node{$taxonID}{'r'} = $rank;
# 5.2) ahora guarda los de los taxones superiores
while($parentID ne $taxonID && $parentID > 1)
{
push(@lineage,$parentID);
$parentID = $node{$parentID}{'p'};
}
@sorted_lineage = sort{$a<=>$b}(@lineage);
$taxon = shift(@sorted_lineage);
## 6) recupera los nombres del linaje de $taxonID leyendo NAMES de nuevo,
## que contiene los taxones ordenados de menor a mayor
open(NAMES,$NAMESFILE) || die "# $0 : cannot open $NAMESFILE\n";
while(<NAMES>)
{
next if(!/scientific name/);
if(/^(\d+)\s+\|\s+(.*?)\|/)
{
($ID,$name) = ($1,$2);
next if($ID ne $taxon);
$name =~ s/^\s+//g;
$name =~ s/\s+$//g;
$lineage_name{$ID} = $name;
if(@sorted_lineage){ $taxon = shift(@sorted_lineage) }
else{ last } # termina si no quedan mas IDs
}
}
close(NAMES);
## 7) enumera el linaje de $taxoID
foreach $ID (@lineage)
{
$rank = $node{$ID}{'r'};
if($rank eq 'no'){ $rank = 'no_rank' }
$lineage_names .= "$lineage_name{$ID}, ";
$lineage_ranks .= "$rank, ";
}
print "\n# lineage : $lineage_names\n# rank : $lineage_ranks\n";
