This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract all bacteria sequences from the nr database

Is there a quick way to extract all bacteria sequences from NCBI non-redundant (NR) database using blastdbcmd?

This command applies only to extract all human sequences from the nr database. How about in microbial sequences case?

$ blastdbcmd -db nr -entry all -outfmt "%g %T" | \
   awk ' { if ($2 == 9606) { print $1 } } ' | \
   blastdbcmd -db nr -entry_batch - -out human_sequences.txt
ncbi nr database bacteria
• 14,843 views • link
updated by jgl4 • jgl4 • written by Penny Liu • Penny Liu

thank you very much for replies

• 0 views • link
written by jgl4 • jgl4

BLAST+ 2.8.1 with blastdb v5 allow you to limit your search by taxonomy using information built into the BLAST databases. So you don't need to build blastdb for specific taxids now !

• 0 views • link
written by shenwei356 shenwei356

2 answers

Taxid for human species is 9606. For species of your interest you may find taxid on NCBI website https://www.ncbi.nlm.nih.gov/Taxonomy/TaxIdentifier/tax_identifier.cgi

• 0 views • link
written by Petr Ponomarenko Petr Ponomarenko

BLAST Database error: No alias or index file found for nucleotide database [nr] in search path

Can you explain what vital steps I missed out? Thank you again.

Getting all taxids of bacteria (taxid 2)

./blastdbcmd -db nr -entry all -outfmt "%g %T" | awk ' { if ($2 == 2) { print $1 } } ' | ./blastdbcmd -db nr -entry_batch - -out bacteria_sequences.txt
• 0 views • link
written by Penny Liu • Penny Liu

I guess there's no NR entry with taxid 2, only that of a species has.

• 0 views • link
written by shenwei356 shenwei356

Hi, I got same error even use the ($2 == 9606). 

./blastdbcmd -db nr -entry all -outfmt "%g %T" | awk ' { if ($2 == 9606) { print $1 } } ' | ./blastdbcmd -db nr -entry_batch - -out human_sequences.txt
• 0 views • link
written by Penny Liu • Penny Liu

blastdbcmd uses GI (genbank ID) to taxid table which only has taxids in it on a species level. If you want to use blastdbcm to extract all bacterial NR sequences you first need to have a list of all bacterial species taxids. In order to make such list you can go to the NCBI taxid website that I mentioned https://www.ncbi.nlm.nih.gov/Taxonomy/TaxIdentifier/tax_identifier.cgi and go to the FTP download link there ftp://ftp.ncbi.nih.gov/pub/taxonomy/ There you will see the ftp://ftp.ncbi.nih.gov/pub/taxonomy/gi_taxid.readme which describes GI to taxid tables mentioned above. There you will see also ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump_readme.txt that describes among other things nodes.dmp file that has a child-parent relationship for the taxid tree.

There are tools like TaxonKit to retrieve all taxids on a certain level for a given taxid from that database, but you might be able to do it yourself faster than installing and learning a new tool.

Once you have a list of species taxids species_taxids.list you can use it in your one-liner by replacing if($2==2) to taxid_list[$2] that you populate from a your species_taxids.list file.

The last question is how to read/write files in awk when it is in the middle of the pipe. You can do it in BEGIN part of the awk or anywhere else:

BEGIN{while (getline < "'"$INPUTFILE"'"){split($0,ft,",");id=ft[1];name=ft[2];key=id;data=name;nameArr[key]=data;}close("'"$INPUTFILE"'");
• 0 views • link
written by Petr Ponomarenko Petr Ponomarenko

There are tools like TaxonKit to retrieve all taxids on a certain level for a given taxid from that database

Sorry I did not say it clearly, but TaxonKit retrieves all taxids BELOW a taxid, not ON a certain level.

There maybe some other better solutions. But TaxonKit so easy to install and use, the speed is also very fast. Why not give it try?

Taxonkit provides executable binary files for Linux/OS X/Windows, users just need to to download the tar.gz file, decompress and immediately run. It also provide detailed examples for all the commands, she only need read the example 1 of command taxonkit list, or tutorial.

• 0 views • link
written by shenwei356 shenwei356

How long should this command take to run ? In my case, it is running for a week and it is still not done ? Does that indicate something is wrong ?

• 0 views • link
written by ambarishnag • ambarishnag

Is the file you are capturing the results in growing in size? While it probably should not take that long to run but it may depend on what kind of hardware resources you have access to. If that file is not growing then you should stop the process and re-assess.

• 0 views • link
written by GenoMax GenoMax

Extract all protein sequences of specific taxons from the NCBI nr database.

TaxonKit (can be installed via conda: conda install -c bioconda taxonkit) is used to get all taxids belonging to given taxid(s).

Taxid of Bacteria (superkingdom) is 2. The tutorial use virus as example, which has much few taxids. But I did not wait the extracting proccess finished, because blastdbcmd is so slow.

• 0 views • link
written by shenwei356 shenwei356

The page displaying title is extract all protein sequences from the NCBI nr database. How to get DNA sequences rather than protein sequences?

• 0 views • link
written by Penny Liu • Penny Liu

Sorry, NR is the non-redundant protein database. Use NT instead.

• 0 views • link
written by shenwei356 shenwei356

TaxonKit worked out perfectly. However, there's something strange here. I already got all taxids of bacteria (taxid 2) and extacted accessions also GIs with csvtk (Please refer to this). Subsequently, I extracting nr sequences as the next step $ blastdbcmd -db nr -entry all -outfmt "%g\t%T" | csvtk -t grep -f 2 -P bacteria.taxid.gi.txt | csvtk -t cut -f 1 | blastdbcmd -db nr -entry_batch - -out nr.bacteria.fa. An error occured:

[ERRO] field (2) out of range (1) in file: -

What am I doing wrong?

By the way, which db should I download? (wget ftp://ftp.ncbi.nih.gov/blast/db/nr.**.tar.gz)

• 0 views • link
written by Penny Liu • Penny Liu

It seems that output of blastdbcmd -db nr -entry all -outfmt "%g\t%T" has only one column. It's strange, but I didn't try it.

You can run the commands separately for debug.

blastdbcmd -db nr -entry all -outfmt "%g\t%T" > nr.gi2tax

And check the nr.gi2tax.

Right, nr locates in ftp://ftp.ncbi.nih.gov/blast/db/nr.**.tar.gz

• 0 views • link
written by shenwei356 shenwei356

NCBI has stopped using gi numbers which is probably why you only got one column from the blastdbcmd command.

• 0 views • link
written by GenoMax GenoMax

Right! This may be the point! Use accession instead

blastdbcmd -db nr -entry all -outfmt "%a\t%T" | \
    csvtk -t grep -f 2 -P virus.taxid.acc.txt | \
    csvtk -t cut -f 1 | \
    blastdbcmd -db nr -entry_batch - -out nr.virus.fa
• 0 views • link
written by shenwei356 shenwei356

We can also retrieve bacterial sequences directly from nt FASTA sequences and then makeblastdb.

gzip -c -d nt.gz | seqkit grep -f bacteria.taxid.acc.txt > nt.bacteria.fasta

makeblastdb -dbtype nucl -in nt.bacteria.fasta -out nt.bacteria

seqkit is here.

• 0 views • link
written by shenwei356 shenwei356

I am researching materials for the Biostar Handbook. It appears that most advice in this thread is not quite right anymore - might have become obsolete.

For example

  1. the nt.gz file does not contain the taxid (it does not have a NCBI formatted header which is super weird), hence seqkit grep wouldn't be able to match it.
  2. On my system the the -outfmt parameter for blastdbcmd does not insert tabs for the \t hence that advice won't work either.
• 0 views • link
written by Istvan Albert Istvan Albert

nt blast database does contain taxid's. Instead of the \t a different delimiter could be used (e.g. ,) with the original awk command. I can't get the -entry_batch to work from STDIN. Were you able to?

• 0 views • link
written by GenoMax GenoMax

yes I did get the -entry_batch to work, make sure to put the - there as well.

• 0 views • link
written by Istvan Albert Istvan Albert

nucl_gb.accession2taxid.gz

• 0 views • link
written by shenwei356 shenwei356

Log in to answer this question.

More posts like this