1
$\begingroup$

I have been having some trouble with running command line BLASTn remotely. Therefore, upon recommendations, and considering I have >100 searches planned, I am considering just downloading the nt database (thankfully possible).

My issue is that one of the settings that proved valuable with my data in WebBlast is running an Entrez query. I get more useful results for each query when instead of one search, I run three, with each of the sub-searches delimited by subject sequence length (e.g. -entrez_query 700:1000[slen] to only get the results which are between 700 and 1000 nucleotides long).

WebBlast and remote Blast have that, local Blast does not.

Since I only really need 2 Entrez queries, I guess it would not be conceptually impossible to get my copy of nt database as FASTA, and then run a script to split it into three directories (700-1000nt, 1000-1400nt, the rest while discarding all<700 nt), then make three blastable databases out of it. However it does feel rather inelegant. I was wondering if anyone might have a better solution!

$\endgroup$
1
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – gringer
    Oct 15, 2022 at 1:19

2 Answers 2

1
$\begingroup$

I looked back at the question and have a more comprehensive solution.

  1. There is no option for -entrez_query on a local blast search. The information is here

entrez_query Restrict search with the given Entrez query. Remote searches only.

There is an option without -outfmt for length and that is worth investigating.

length means Alignment length

  1. Perl one liner solution to NCBI's many parameter problem

a) Preprocess your fasta file into separate databases of the required length ... by turning a fasta flatfile db into one where all the subject sequences are 700 - 1000.

perl -0777 -nle 'while (/^(>.+\n^[AGCT-NRYMKWS]{700,1000})\n/gm){print $1,"\n"}' example.fa > example700-1000.fa

b) Now make the database on the restricted length fasta flatfile db,

makeblastdb -in example700-1000.fa -parse_seqids -blastdb_version 5  -out laura700-1000.db -title "Laura 700 - 1000" -dbtype nucl

This db will now contain an equivalent of the -entrez_query 700:1000[slen]. Obviously, you would need to make a new db for each change in the sequence length query. Note the option -taxid_map.

c) Finally, run the blast search on the restricted length blast db.

blastn -query laura.fas -db laura700-1000.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

Code testing

Input is test.fa

>here
AGAGACGTAGAGAAAAAAAAAA
>here2
TAGAGA
>here3
AGAGACGTAGAGA

Modified perl oneliner because sequence length is different,

perl -0777 -nle 'while (/^(>.+\n^[AGC-T]{6,13})\n/gm){print $1,"\n"}' test.fa

Output

>here2
TAGAGA
>here3
AGAGACGTAGAGA

Thus the sequence that was outside the range 6 to 13 (sequence here2 and here3) is excluded (here). This seems alot easier to me that all those blastn options.

Note NRYMKWS refer to the degenerate nucleotide code

$\endgroup$
0
$\begingroup$

"Often, one needs to search ... a specific subset of sequences within an existing database. For these type of searches a convenient way to conduct them is by creating a virtual BLAST database"

From "Use blastdb_aliastool to manage the BLAST databases" https://www.ncbi.nlm.nih.gov/books/NBK569848/

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.