$ ./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";
No hay comentarios:
Publicar un comentario