Hola,
en esta entrada quisiera mostraros una manera de buscar secuencias de proteínas dentro de la base de datos estrictamente no redundante
nr del
NCBI, restringiendo la búsqueda a grupos taxonómicos arbitrarios, sin necesidad de duplicar secuencias en otros archivos. De hecho, por su gran tamaño, 45 millones de secuencias a dia de hoy, desde hace tiempo nosotros ya no trabajamos con el archivo FASTA de
nr, si no con la base de datos ya formateada para blast que descargamos regularmente con
update_blastdb.pl desde
ftp://ftp.ncbi.nlm.nih.gov/blast/db .
La receta para realizar búsquedas sobre subconjuntos de
nr la hemos sacado de
https://www.biostars.org/p/6528, y aquí os mostraré un ejemplo con secuencias de plantas:
1) Si no la tienes ya descárgate con el script
update_blastdb.pl una copia de nr. Ojo, son una serie de archivos que en total ocupan 42G:
$ perl update_blastdb.pl nr
2) Encontrar el identificador taxonómico que mejor represente al subconjunto de
nr que quieres definir, en nuestro caso
Viridiplantae:
http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=33090
3) Buscar todas las proteínas de plantas anotadas en el NCBI y exportar (
Send to File => GI List) sus identificadores de GenBank (
GIs) a un archivo, por ejemplo 'Viridiplantae_protein.gi.txt':
http://www.ncbi.nlm.nih.gov/protein/?term=txid33090[ORGN]
4) Hacer un alias de
nr que contenga sólo proteínas de plantas:
$ ncbi-blast-2.2.29+/bin/blastdb_aliastool -gilist Viridiplantae_protein.gi.txt -db nr -out nr_plants -title nr_plants
En mi ejemplo produce el siguiente output:
Converted 3741334 GIs from Viridiplantae_protein.gi.txt to binary format in nr_plants.p.gil
Created protein BLAST (alias) database nr_plants with 2041649 sequences
Ahora veamos qué diferencias y ventajas tiene usar este subconjunto de
nr en comparación con la base de datos completa:
[ NR ]
$ time ncbi-blast-2.2.29+/bin/blastx -query test.fna -db nr
Database: All non-redundant GenBank CDS translations+PDB+SwissProt+PIR+PRF
excluding environmental samples from WGS projects
45,360,603 sequences; 16,206,973,370 total letters
Query= test
Length=312
Score E
Sequences producing significant alignments: (Bits) Value
dbj|BAK07057.1| predicted protein [Hordeum vulgare subsp. vulgare] 156 8e-47
dbj|BAJ89957.1| predicted protein [Hordeum vulgare subsp. vulgar... 62.0 5e-10
ref|WP_014745933.1| asparagine synthase [Tistrella mobilis] >ref 42.7 0.020
ref|WP_006720063.1| hypothetical protein [Collinsella stercoris]... 35.0 5.0
ref|XP_001817163.2| hypothetical protein AOR_1_2924174 [Aspergil... 35.0 5.3
dbj|BAE55161.1| unnamed protein product [Aspergillus oryzae RIB40] 35.0 6.9
ref|WP_021178774.1| Phosphate transport system permease protein ... 34.7 7.6
ref|WP_021807128.1| Phosphate transport system permease protein ... 34.7 7.9
ref|WP_026166948.1| glucan biosynthesis protein D [Novosphingobi... 34.3 8.8
[...]
real 5m28.923s
user 5m24.564s
sys 0m3.596s
[ NR_plants ]
$ time ncbi-blast-2.2.29+/bin/blastx -query test.fna -db nr_plants
Database: nr_plants
2,041,649 sequences; 772,729,335 total letters
Query= test
Length=312
Score E
Sequences producing significant alignments: (Bits) Value
dbj|BAK07057.1| predicted protein [Hordeum vulgare subsp. vulgare] 156 4e-48
dbj|BAJ89957.1| predicted protein [Hordeum vulgare subsp. vulgar... 62.0 2e-11
gb|EMS57556.1| hypothetical protein TRIUR3_28581 [Triticum urartu] 33.1 0.95
ref|XP_002311179.2| hypothetical protein POPTR_0008s05850g [Popu... 33.5 1.1
ref|XP_004509944.1| PREDICTED: uncharacterized protein LOC101513... 32.0 3.0
ref|XP_004509945.1| PREDICTED: uncharacterized protein LOC101513... 32.0 3.0
gb|EMT15944.1| hypothetical protein F775_17477 [Aegilops tauschii] 32.0 3.3
dbj|BAK00277.1| predicted protein [Hordeum vulgare subsp. vulgare] 32.0 3.7
ref|XP_004234300.1| PREDICTED: uncharacterized protein LOC101244... 31.2 4.7
ref|XP_006358717.1| PREDICTED: LRR receptor-like serine/threonin... 31.2 5.2
ref|XP_002316273.2| hypothetical protein POPTR_0010s20870g [Popu... 31.2 6.6
ref|XP_003571199.1| PREDICTED: probable tocopherol cyclase, chlo... 30.8 8.4
ref|XP_008451998.1| PREDICTED: respiratory burst oxidase homolog... 30.4 9.9
[...]
real 0m16.805s
user 0m16.332s
sys 0m0.372s
Como se puede ver en ambos casos las dos primeras secuencias encontradas son las mismas, con la misma puntuación en bits y distinto E-valor, al ser universos de búsqueda de tamaños distintos.
Además, el tiempo de búsqueda con un sólo procesador es del orden de 20x más rápido con el alias
nr_plants, y a la vez el consumo de RAM en tiempo real mucho menor.
Hasta luego,
Bruno