¿Necesitas obtener el taxid de una o varias secuencias de NCBI y/o necesitas conocer su linage completo (taxonomía de siete niveles)? este es tu post.

Primero es necesario tener una lista de IDs de las secuencias de NCBI que te interesan. Supondré que el listado lo tienes en un archivo que se llama IDs.txt (típicamente una lista de GIs o Accessions).

Para obtejer los taxID tengo dos formas, la lenta pero segura y la rápida pero aparatosa de implementar.

1.2 TaxID opción1: Usando efetch (leeeeento pero seguro... recomiendo opción2)

Lo primero es crear un script bash con el siguiente contenido... lo llamaré taxID_efetch.sh

for ACC in $@
do
   echo -n -e "$ACC\t"
   curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=protein&id=${ACC}&rettype=fasta&retmode=xml" |\
   grep TSeq_taxid |\
   cut -d '>' -f 2 |\
   cut -d '<' -f 1 |\
   tr -d "\n"
   echo
   sleep 0.25
done

Una vez creado el script lo ejecutamos de la siguiente manera

cat IDs.txt | xargs bash taxID_efetch.sh > IDs_taxID.txt

Listo, en el archivo IDs_taxID.txt tendrás los IDs de taxonomía de cada secuencia.

1.2 TaxID opción2: Usando el archivo prot.accession2taxid de NCBI (rápido algo complejo de implementar)

#descargamos la base de datos y comprobamos la descarga
wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz
wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz.md5
md5sum -c prot.accession2taxid.gz.md5

Crearemos un script de perl que nos permitirá buscar de forma rápida, lo llamaré select_from_first_column.gzipped.pl. Este script asume que el archivo prot.accession2taxid.gz está separado por tabulaciones y que el archivo IDs.txt contiene sólo una columna con los IDs.

#!/usr/bin/perl
use strict;

sub print_uso{
print "Uso:
        select_est archivo_fuente lista_identificadores
        select_est -r archivo_fuente lista_identificadores (para seleccion inversa)\n";
}

my $rev = 0;            # reverso?

if($ARGV[0] eq "-r"){
        $rev = 1;
}

# si los parametros no son la cantidad correcta, salimos
if($#ARGV != (1+$rev)){
        print_uso();
        exit;
}

# leemos lista de identificadores y lo recordamos en un arreglo asociativo
my %ids;
my %ids_estan;
my @array;
open LISTFILE, $ARGV[1+$rev] or die "No pude abrir lista ($ARGV[1+$rev])!\n";
my $s = "";
my $count =0;
my $count2 =0;
while( $s = <LISTFILE> ){
        chomp($s);
        $ids{$s} = 1;
        $ids_estan{$s} = 1;
        $count++;
}

open(SRCFILE, "gunzip -c $ARGV[$rev] |") or die "gunzip $ARGV[$rev]: $!";
while (<SRCFILE>){
        chop;
        @array = split("\t", $_);
        $s = $array[1];
        if(($ids_estan{$s} == 1)){
                $ids_estan{$s} = 2;
        }
        if(($ids{$s} == 1)){            #la buscamos en arr_aso
                print "$array[1]\t$array[2]\n";
                $count2++;
        }
        if($count2 eq $count){
                last;
        }
}

foreach my $key (keys %ids_estan){
        if($ids_estan{$key} == 1){
                print  STDERR "$key no estaba en archivo fasta.\n";
        }
}

Finalmente, ejecutamos el siguiente comando para obtener los IDs

perl select_from_first_column.gzipped.pl prot.accession2taxid.gz IDs.txt > IDs_taxID.txt 2>err.txt

Si el archivo err.txt contiene IDs quiere decir que no se encontraron en el archivo prot.accession2taxid.gz lo cual sucede cuando se eliminan datos desde NCBI por lo cual recomiendo repetir el procedimiento con dichos IDs buscando en el archivo dead_prot.accession2taxid.gz el cual contiene GIs antiguos.

#descargamos la base de datos y comprobamos la descarga
wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/dead_prot.accession2taxid.gz
wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/dead_prot.accession2taxid.gz.md5
md5sum -c dead_prot.accession2taxid.gz.md5

Luego, con este archivo y los IDs faltantes se debe repetir el paso del script perl y el resultado se debe concatenar al primer resultado.

2. Obtener el linaje (taxonomía completa)

Este procedimiento requiere un listado de taxIDs.... supongo que ya lo tienes, de lo contrario lee mas arriba.

Primero necesitamos una base datos que relaciona cada taxID de NCBI con su linaje (si ya la tienes, saltar al paso siguiente). Para ello usaremos el siguiente manual https://github.com/zyxue/ncbitax2lin

Una vez obtenido el archivo ncbi_lineages.csv, lo que queda es hacer magia

awk '{print $2}' IDs_taxID.txt | sort | uniq | sed -e 's/^/^/' -e 's/$/,/' > IDs_taxID2.txt

grep -f IDs_taxID2.txt ncbi_lineages_2021-06-30.csv | awk -F ',' '{print $1"\t"$2";"$3";"$4";"$5";"$6";"$7";"$8}' > IDs_taxID2.txt.lineage

Con esto obtendrás el linaje (siete niveles taxonómicos) para cada taxID.

3. Unir GIs con su linaje

Si quieres hacer un merge entre los GIs y su taxonomía, este es el comando (asumiendo que has seguido mis instrucciones anteriores).

El siguiente comando sólo deja las columnas GI y linaje

join -t$'\t' -i -1 2 -2 1 <(sort -n -t$'\t' -k2 IDs_taxID.txt)  <(sort -n -t$'\t' -k1 IDs_taxID2.txt.lineage) | awk '{print $2"\t"$3}'

El siguiente comando entrega todas las columnas (GI, taxID, linaje).

join -t$'\t' -i -1 2 -2 1 <(sort -n -t$'\t' -k2 IDs_taxID.txt)  <(sort -n -t$'\t' -k1 IDs_taxID2.txt.lineage) | awk '{print $2"\t"$1"\t"$3}'