This is a test version of Biostars. For the public version, visit https://www.biostars.org.
-taxids option in blast is not working

Hello !

I downloaded the nr database from NCBI :

update_blastdb.pl --decompress nr [*]

and would like to make a seaprate database only for fungi. But when I run

blastdbcmd -db /scratch/caroz/NCBI_nr_database/nr -taxids 4751 -outfmt %a -out /scratch/caroz/NCBI_fungi_nr_database/fungal.fa

I get

Error: [blastdbcmd] Taxonomy ID(s) not found. This could be because the ID(s) provided are not at or below the species level. Please use get_species_taxids.sh to get taxids for nodes higher than species (see https://www.ncbi.nlm.nih.gov/books/NBK546209/).

However when I run

/home/caroz/mambaforge/envs/blast_env/bin/get_species_taxids.sh -n Fungi

Taxid : 4751
 rank : kingdom
 division : fungi
 scientific name : Fungi
 common name :

1 matche(s) found.

I cannot find what I am doing wrong.

blast
• 3,665 views • link
written by caroline.zanchi • caroline.zanchi

1 answer

Did you download the taxID database files (taxdb* from https://ftp.ncbi.nlm.nih.gov/blast/db/). You need those for this to work. They need to be in the same folder as your blast+ index files. Once you get the files this should start working.

Assuming you did, you need to do the following to get all fungal taxID's at species level

$ get_species_taxids.sh -t 4751 > fungi_taxid.txt

BTW: if i run

$ blastdbcmd -db nr -taxids 4751 -outfmt %f

I am able to get fungal sequences (just showing top 10 headers as example below)

$ blastdbcmd -db nr -taxids 4751 -outfmt %f | grep ">" | head -10
>KJX92028.1 hypothetical protein TI39_contig5958g00003 [Zymoseptoria brevis]
>CAL1705590.1 unnamed protein product [Somion occarium]
>GAO13639.1 hypothetical protein UVI_02017470 [Ustilaginoidea virens]
>KAM0157508.1 hypothetical protein ACHAPG_004631 [Botrytis cinerea] >KAM0160633.1 hypothetical protein ACHAQE_004490 [Botrytis cinerea]
>XP_033440462.1 uncharacterized protein GL218_08664 [Daldinia childiae] >KAF3067330.1 hypothetical protein GL218_08664 [Daldinia childiae]
>RHZ78344.1 hypothetical protein Glove_166g22 [Diversispora epigaea]
>KAJ8701588.1 Exocyst complex component S5, variant 2 [Pleurotus ostreatus]
>KAG7829670.1 hypothetical protein KL920_002529 [Ogataea angusta] >KAG7840949.1 hypothetical protein KL942_001937 [Ogataea angusta] >KAG7849900.1 hypothetical protein KL940_002268 [Ogataea angusta] >KAG7860833.1 hypothetical protein KL939_001400 [Ogataea angusta]
>KIJ45469.1 hypothetical protein M422DRAFT_250789 [Sphaerobolus stellatus SS14]
>KAL4786389.1 FAD binding domain protein [Aspergillus varians]
• 0 views • link
written by GenoMax GenoMax

I did not ! And I realize now that I overlooked your previous comment on the matter. How do I know where my blast+ index files are ? Here is an example of what the directory in which I downloaded nr looks like :

ls
nr.000.phd         nr.033.ppd         nr.067.phd         nr.100.ppd
nr.000.phi         nr.033.ppi         nr.067.phi         nr.100.ppi
nr.000.phr         nr.033.psq         nr.067.phr         nr.100.psq
nr.000.pin         nr.033.pxm         nr.067.pin         nr.100.pxm
nr.000.pog         nr.033.tar.gz.md5  nr.067.pog         nr.100.tar.gz.md5
nr.000.ppd         nr.034.phd         nr.067.ppd         nr.101.phd
nr.000.ppi         nr.034.phi         nr.067.ppi         nr.101.phi

...

nr.pal
nr.pdb
nr.pjs
nr.pos
nr.pot
nr.ptf
nr.pto
taxdb.btd
taxdb.bti
taxonomy4blast.sqlite3
• 0 views • link
written by caroline.zanchi • caroline.zanchi

Here is an example of what the directory

That would be the directory to use. I assume you have all pieces of the nr database (132 as of today) in this directory along with taxdb* files. If that is correct then the following should work (replace your directory path below)

$ blastdbcmd -db /path_to_dir_with/nr -taxids 4751 -outfmt %f

I assume you are using the latest blast+ (v.2.17)?

• 0 views • link
written by GenoMax GenoMax

Yes :) Now there is

taxdb.tar.gz
taxdb.tar.gz.md5

in the directory of my database. I suppose I should gunzip the tar.gz ?

• 0 views • link
written by caroline.zanchi • caroline.zanchi

Correct.

• 0 views • link
written by GenoMax GenoMax

When I do, this happens :

tar -xvzf taxdb.tar.gz
taxdb.btd
taxdb.bti
taxonomy4blast.sqlite3

These files were already present in my directory though ?

...and again

Error: [blastdbcmd] Taxonomy ID(s) not found. This could be because the ID(s) provided are not at or below the species level. Please use get_species_taxids.sh to get taxids for nodes higher than species (see https://www.ncbi.nlm.nih.gov/books/NBK546209/).
• 0 views • link
written by caroline.zanchi • caroline.zanchi

This should work. Not sure what is different in your case.

While this should not be necessary, if you are providing the full path, it would not hurt to set the BLASTDB variable to point to the directory you have the files in and see if that helps.

export BLASTDB=/dir_w_blast_db
• 0 views • link
written by GenoMax GenoMax

same outcome...

Error: [blastdbcmd] Taxonomy ID(s) not found. This could be because the ID(s) provided are not at or below the species level. Please use get_species_taxids.sh to get taxids for nodes higher than species (see https://www.ncbi.nlm.nih.gov/$

However, I just realized that I used an other version of BLAST to run this :

module add BLAST+/2.14.1-gompi-2023a
update_blastdb.pl --decompress nr [*]
• 0 views • link
written by caroline.zanchi • caroline.zanchi

As long as you are getting the latest set of blast indexes (v.5) this should work. I can't see of an immediate reason why it is not working.

Try using 162425 (taxID for Aspergillus nidulans) and see if that works. This is a species level taxID.

• 0 views • link
written by GenoMax GenoMax

I could not try because I erased the database in order to re-download it from scratch with the latest version of Blast+. This reminded me why I used an older version before : because I am not managing to download the nr database with the version of Blast+ that I installed in my Conda environment :

mamba activate blast_env

update_blastdb.pl --decompress nr [*]

which: no gsutil in (/home/caroz/mambaforge/envs/blast_env/bin:/home/caroz/mambaforge/condabin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/home/caroz/.local/bin:/home/caroz/bin)
which: no gcloud in (/home/caroz/mambaforge/envs/blast_env/bin:/home/caroz/mambaforge/condabin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/home/caroz/.local/bin:/home/caroz/bin)
Connected to NCBI
Warning: No BLASTDB metadata for [*]
[*] not found, skipping.
Downloading nr.000.tar.gz...Unable to close datastream at /home/caroz/mambaforge/envs/blast_env/bin/update_blastdb.pl line 593.
Undefined subroutine &main::trim_ftp_prefix called at /home/caroz/mambaforge/envs/blast_env/bin/update_blastdb.pl line 606.
• 0 views • link
written by caroline.zanchi • caroline.zanchi

Having old format databases (v.4) can certainly cause this to not work.

If you are not able to use the update script you could manually create the command lines to download the files using curl or wget.

Something like following should do the trick. You could download the files in an interactive session or submit as individual jobs to your cluster scheduler. (remove the word echo when you actually download/submit jobs).

for i in {000..132}; do echo wget https://ftp.ncbi.nlm.nih.gov/blast/db/nr.${i}.tar.gz; done

This will produce the necessary 132 commands (only an example shown).

wget https://ftp.ncbi.nlm.nih.gov/blast/db/nr.000.tar.gz
wget https://ftp.ncbi.nlm.nih.gov/blast/db/nr.001.tar.gz
wget https://ftp.ncbi.nlm.nih.gov/blast/db/nr.002.tar.gz
wget https://ftp.ncbi.nlm.nih.gov/blast/db/nr.003.tar.gz
wget https://ftp.ncbi.nlm.nih.gov/blast/db/nr.004.tar.gz
• 0 views • link
written by GenoMax GenoMax

Sooo I have :

which results in /path/NCBI_nr_database looking like this :

nr.000.phd
nr.000.phi
nr.000.phr
nr.000.pin
nr.000.pog
nr.000.ppd
nr.000.ppi         
nr.000.psq
nr.000.psq    
nr.000.pxm
nr.000.tar.gz 
nr.001.phd
nr.001.phi        
nr.001.phr    
nr.001.pin        
nr.001.pog       
nr.001.ppd    
nr.001.ppi         
nr.001.psq      
nr.001.pxm 

[...]
nr.pal
nr.pdb
nr.pjs
nr.pos
nr.pot
nr.ptf
nr.pto
taxdb.btd
taxdb.bti
taxonomy4blast.sqlite3

I ran

source /home/${USER}/.bashrc mamba activate blast_env export BLASTDB=/path/NCBI_nr_database/nr

blastdbcmd -db path/NCBI_nr_database/nr -taxids 4751 -outfmt %a -out /path/NCBI_fungi_nr_database/fungal.fa

with exactly the same outcome as before.

If I run

source /home/${USER}/.bashrc

mamba activate blast_env

/path/mambaforge/envs/blast_env/bin/get_species_taxids.sh -t 4751 > fungi_taxid.txt

I do get a text file full of taxa numbers.

• 0 views • link
written by caroline.zanchi • caroline.zanchi

Are you certain you are running the blastdbcmd from the latest blast+. It looks like you have multiple versions of blast installed and a previous version may be set as default.

What happens if you run the blastdbcmd command after you activate the mamba environment?

• 0 views • link
written by GenoMax GenoMax 
mamba activate blast_env 
blastdbcmd
BLAST Database error: No alias or index file found for nucleotide database [nr] in search path [/path/NCBI_fungi_nr_database::]

THis is my working directory, it is indeed empty.

• 0 views • link
written by caroline.zanchi • caroline.zanchi

I meant that you should run the complete blastdbcmd command with the full directory path for nr database in it,

• 0 views • link
written by GenoMax GenoMax

I ran it (not as a script), and the outcome is the same.

Edit : I mamba cleaned --all, reinstalled blast+ in a new environment, and so far it seems to be running ! I wonder what happened in the previous one.

• 0 views • link
written by caroline.zanchi • caroline.zanchi

I am stumped. As you can see above I am able to get the command to work with install I have access to.

If species level ID's are working for you then that may be the way to go. Sounds like you at least got those with the accessory script. Use a for loop and go through those taxID.

• 0 views • link
written by GenoMax GenoMax

blastdbcmd is running right now but terribly slowly, and times out after 3 days. I guess I should allocate more resources to it ?

• 0 views • link
written by caroline.zanchi • caroline.zanchi

Are you referring to species level ID;s working? That is likely not an efficient process since the entire index needs to be re-read for each ID.

Using the 4751 taxID should work the best.

• 0 views • link
written by GenoMax GenoMax

I am running this :

blastdbcmd -db /path/NCBI_nr_database/nr -taxids 4751 -outfmt %a -out /path/NCBI_fungi_nr_database/fungal.fa -num_threads ${SLURM_CPUS_PER_TASK}
• 0 views • link
updated by GenoMax GenoMax • written by caroline.zanchi • caroline.zanchi

If the database is stored on external/slow storage then try moving it to the fastest storage you have available.

• 0 views • link
written by GenoMax GenoMax

It is already on the fastest one we have, thought there were problems with it lately. Maybe this could be the cause ? In the meantime I decided to try with the Swissprot database :

update_blastdb.pl --decompress swissprot

then

blastdbcmd -db /scratch/caroz/Swissprot_database/swissprot -taxids 4751 -outfmt %a -out /scratch/caroz/Swissprot_fungi/fungal.fa

which gave me a list of accession numbers. You wrote in an other post that I should not build my database with makeblastdb using these fungal.fa, however this really contains only accession numbers. Do I have to choose an other output format, or is there something that I misunderstood ? Thank you !

• 0 views • link
written by caroline.zanchi • caroline.zanchi

-outfmt %a

You asked to get accessions. You will need to use %f to get fasta formatted sequence with blastdbcmd.

swissprot may be a smaller dataset than nr. As long as you are OK with that.

• 0 views • link
written by GenoMax GenoMax

Thank you ! I would actually love to have a database containing only fungal sequences that I could use for the annotation of a genome I have assembled.

• 0 views • link
written by caroline.zanchi • caroline.zanchi

Log in to answer this question.

More posts like this