Hi, I'm trying to do an rpsblast using ncbi-blast-2.5.0+. I have a file containing gi numbers as my query and my command line is below.
rpsblast -query t -db results/Cog/Cog -out rps-blast.out -evalue 1e-2 -outfmt 6
Oddly the rpbslast keeps giving me this error:
Warning: [rpsblast] Error initializing remote BLAST database data loader: Protein BLAST database 'Cog/Cog nr' does not exist in the NCBI servers
My query file (t) looks like this. I've tried to remove search using just the numbers (292833481) but it still doesn't solve the problem. gi|292833481 gi|383341230 gi|289693981
However, when I try to the same search but using a fasta file as query, it runs fine and gives the results I need.
rpsblast -query GCF_000005845.2_ASM584v2_protein.faa -db results/Cog/Cog -out rps-blast.out -evalue 1e-2 -outfmt 6
Is there something that I did wrong here? What is the correct way to format query a list of gi's for blast? Thanks!
Hi all, Thank you for the responses. Somehow I am able to use GI numbers as queries for rpsblast using ncbi-blast-2.6.0. This was also tested by Dr. Wayne Matten from NCBI User Services.Here's what I did:
Sort and output only unique GI numbers, else a C++ exception error
will occur (see ncbi error C++ exception ). My query, looks
something like this: gi|545375397 545375577 545375696
Download and untar (tar -xvzf) cdd.tar.gz from ftp://ftp.ncbi.nlm.nih.gov/pub/mmdb/cdd/
Do makeprofiledb on several types of database. (you can read more about this from here: ftp://ftp.ncbi.nlm.nih.gov/pub/mmdb/cdd/README)
makeprofiledb -title SMART.v6.0 -in Smart.pn -out Smart -threshold 9.82 -scale 100.0 -dbtype rps -index true makeprofiledb -title Pfam.v.26.0 -in Pfam.pn -out Pfam -threshold 9.82 -scale 100.0 -dbtype rps -index true makeprofiledb -title COG.v.1.0 -in Cog.pn -out Cog -threshold 9.82 -scale 100.0 -dbtype rps -index true makeprofiledb -title KOG.v.1.0 -in Kog.pn -out Kog -threshold 9.82 -scale 100.0 -dbtype rps -index true makeprofiledb -title CDD.v.3.12 -in Cdd.pn -out Cdd -threshold 9.82 -scale 100.0 -dbtype rps -index true makeprofiledb -title CDD_NCBI.v.3.12 -in Cdd_NCBI.pn -out Cdd_NCBI -threshold 9.82 -scale 100.0 -dbtype rps -index true makeprofiledb -title PRK.v.6.00 -in Prk.pn -out Prk -threshold 9.82 -scale 100.0 -dbtype rps -index true
Before running the job, go to the directory where you save the database
Run rpsblast:
/ncbi-blast-2.6.0+/bin/rpsblast -query /data/test -db Pfam -out /data/test_result -evalue 0.0001 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -max_target_seqs 20 -num_threads 16
As the job runs, I saw this error message:
Warning: [rpsblast] Error initializing remote BLAST database data loader: Protein BLAST database 'Pfam nr' does not exist in the NCBI servers
However, the search results are still being saved to my output file. The results look something like this:
gi|545375397|ref|XP_005651884.1| gnl|CDD|279733 28.856 402 238 13 5 376 64 447 1.44e-64 212 gi|545375577|ref|XP_005651955.1| gnl|CDD|306679 36.207 174 108 2 46 219 1 171 3.07e-37 126 gi|545375696|ref|XP_005652002.1| gnl|CDD|306830 43.590 39 22 0 96 134 1 39 1.97e-08 47.4 gi|545375696|ref|XP_005652002.1| gnl|CDD|306830 48.571 35 18 0 58 92 5 39 6.40e-08 46.2 gi|545375696|ref|XP_005652002.1| gnl|CDD|306830 47.368 38 20 0 185 222 2 39 8.56e-07 42.7 gi|545375696|ref|XP_005652002.1| gnl|CDD|306830 31.818 44 25 2 137 180 1 39 4.33e-06 40.8 gi|545375696|ref|XP_005652002.1| gnl|CDD|306830 36.111 36 21 1 9 44 5 38 8.10e-05 37.3
Woohoo! There you have it. I apologize if the answer is disorganized, not very certain how to use the formatting tools :/ Hope this is helpful!
not very certain how to use the formatting tools
Pretty simple. Highlight code that you want to mark as such and then use the "101" button in the edit window. That should correctly format the text. Note: I am not sure what you did with the text above but I am not able to correctly format it.
I am curious as to how long the gi numbers are going to keep working since NCBI is on record as saying that users need to move away from using them.
I do not have the cog database installed on my machine, but some tests with blastp instead of rpsblast indicate, that it is neccessary to put every GI number in a separate line:
echo 'gi|292833481\ngi|383341230\n' | blastp -task blastp-fast -remote -db nr -outfmt 6
The above command will emit lots of error messages, but after 10 min it has finished successfully.
Another test case to prove the ability of the local blast client to lookup GI numbers and retrieve sequences from NCBI :
echo 'gi|114050348\n' > q.txt
echo 'gi|57284222\n' > s.txt
blastn -query q.txt -subject s.txt -outfmt 6
We specify two sequences by their GI numbers, and then blast the shorter query sequence against the longer subject (like the deprecated bl2seq). The expected outcome is:
AB234058.1 CP000046.1 99.888 889 1 0 1 889 1666415 1665527 0.0 1637
Nevertheless, I regard this automatic sequence retrieval by the local blast client as not really stable and mature.
inloraj, if you already have a list of GI numbers, than you can easily download the sequences with Eutils efetch, and then feed them into rpsblast.
wget 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=protein&rettype=fasta&retmode=txt&id=gi|292833481,gi|383341230,gi|292833481' -O queries.fas
Log in to answer this question.
NCBI has stopped using
ginumbers externally since September 2016. You should substitute theginumbers withaccession numbers.This is probably THE correct answer to this question.
Thank you for your response! Unfortunately, I didn't see any improvement when I converted my query from gi number to accession (eg EFL06024.1). The error regarding "Protein BLAST database 'Cog/Cog nr' does not exist in the NCBI servers" still there....
How old is your
rps-blast? From https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMDWeb&PAGE_TYPE=BlastNews:NCBI is hiding GI numbers from inexperienced kids now, but you can still use them in blast queries and eutils. (see below)
I could reproduce your problem using a valid identifier, maybe it is a bug? Valid identifier worked online.
Fast solution is download the sequences of interest, and use a fasta file.
edit: are you using
rpsblastorrpsblast+?Hi h.mon, I'm using rpsblast+ from blast+ package 2.5.0 and 2.6.0(latest). Yes, I suppose downloading fasta files will be the quickest solution, although my queries are rather large (>700K). NCBI's CD-search accept gi/accession number as query (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) but they only allow 4000 queries each time :/ I've emailed ncbi-help and will update accordingly. Thanks for your help!
Shouldn't you edit your post then? If I try
rpsblastinstead ofrpsblast+, I get a lot of errors due to incorrect parsing of the arguments.