Hi,
I have downloaded 2500 genome assembly fasta file and changed them to database formatted files using "makeblastdb". When I blast my query on them using this command (I just show it for one genome file, GCA_000143925.2.fasta as an example) (I run blast for all in parallel):
blastn -query query.fasta -db GCA_000143925.2.fasta -outfmt "6 std qlen slen staxids sscinames" -task dc-megablast -out blast.out
I get this error message:
Warning: [blastn] Taxonomy name lookup from taxid requires installation of taxdb database with ftp://ftp.ncbi.nlm.nih.gov/blast/db/taxdb.tar.gz
(I have downloaded and unzipped taxdb.tar.gz file and added the path to my system .bashrc (export TAXDB=/home/manighanipoorsamami/local/taxdb), but still get the same error.
Also here is the first hit line of blast.out file:
ERV2-1_H.orn#LTR/ERV GL380075.1 81.633 49 7 1 7855 7901 77025 77073 0.018 46.4 8632 90134 0 N/A
as you see it does not give staxids sscinames.
Then, I created the file "taxid_mapping_file.txt" by adding taxid to genome name:
cat taxid_mapping_file.txt
GCA_000143925.2.fasta 135651
The, ran this:
makeblastdb -in GCA_000143925.2.fasta -taxid_map taxid_mapping_file.txt -parse_seqids -dbtype nucl
but got this error:
Building a new DB, current time: 01/25/2024 16:35:44
New DB name: /home/manighanipoorsamami/New_PhD_program/HTT_sea_snake/test/GCA_000143925.2.fasta
New DB title: GCA_000143925.2.fasta
Sequence type: Nucleotide
Keep MBits: T
Maximum file size: 3000000000B
Adding sequences from FASTA; added 3305 sequences in 1.22134 seconds.
Error: [makeblastdb] No sequences matched any of the taxids provided.
Can you please help me resolve this and run a blast command that can add "staxids" and "sscinames" to blast output?
I could not find anything in NCBI blast manuals.
Cheers,
Mani
2 answers
The mapping taxonomy ID file is a bit wrong:
cat taxid_mapping_file.txt
GCA_000143925.2.fasta 135651
In there you need the IDs of each sequence, not the filename itself.
So for a fasta file like this:
>seq1 blablabla
ATCTAGCTAGCTAGCTAGCTAGA
>seq2 blablabla
ATCYAYGTCGACTGATCGA
The taxid_mapping_file would look like this:
seq1 135651
seq2 135651
If each assembly is of a single species, there is no need to use -taxid_map because all sequences in the db have the same, so you can simply use the -taxid parameter:
makeblastdb -in GCA_000143925.2.fasta -taxid 135651 -parse_seqids -dbtype nucl
You can get the taxid from the assembly in different ways, either by issuing an eutils query, or if you have your mapping file already, you can simply use grep to find it in a script, e.g.
TAXID=$(grep -e $ASSEMBLY $mapping_file | cut -f2}
You will still need the taxonomy downloaded to extract names but that should work immediately with these databases. If you want to blast against all databases at once, create a single alias DB with blastdb_aliastool.
Log in to answer this question.