This is a test version of Biostars. For the public version, visit https://www.biostars.org.
BLAST, setting maximum number of hits

Hello,

I am trying to set the number of maximum hits to 5, so that the procedure can finish sooner, but I still get 100s of hits found.

# TBLASTX 2.2.29+
# Query: Locus_40_Transcript_185/186_Confidence_0.224_Length_4778
# Database: ../../Genome/Genome
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 714 hits found

I am running:

tblastx -db ../../Genome/Genome -query all_merged_k125.fa -evalue 1e-10 -outfmt 7 -out tblastx/all_merged_k125.fmt7 -num_threads 16 -max_target_seqs 5

Any idea why it's still reporting so many hits?

Adrian

blast tblastx parallel
• 42,966 views • link
updated by Ram Ram • written by apelin20 apelin20

3 answers

You can get e.g. 10 hits from one long target sequence (-max_target_seqs 1), i.e. max_target_seqs doesn't specify the maximum number of hits per query, but the maximum number of target sequences for hits per query.

• 0 views • link
updated by Ram Ram • written by 5heikki 5heikki

Ok that makes sense. How do I limit the amount of hits?

• 0 views • link
updated by Ram Ram • written by apelin20 apelin20

I guess you can always give a relatively stringent e-value and filter the resulting hits later.

• 0 views • link
written by Ram Ram

What I wanted is to speed up the blasting.

• 0 views • link
updated by Ram Ram • written by apelin20 apelin20

I doubt limiting the number of hits like that would speed up your blasting significantly. It still has to go through the whole db for every query, so the only difference would be in how long it takes to write 5 or 10 lines (or whatever) to the output file. Instead, if your db is small (or you have a ton of RAM), you should parallelize blast (e.g. with GNU Parallel) by running multiple single-threaded blasts on split input instead of using -num_threads X..

• 0 views • link
updated by Ram Ram • written by 5heikki 5heikki

^True. You will benefit from multi-threading, and trying both tblastx and blastall -p tbalstx before choosing one of them. For shorter query sequences, I've seen the latter be significantly faster than the former.

• 0 views • link
written by Ram Ram

I would do it post blast (with -outfmt 6 output):

Make sure the file is sorted based on query and best hits (here bitscore > evalue > perc identity):

export LC_ALL=C LC_LANG=C; sort -k1,1 -k12,12gr -k11,11g -k3,3gr outputFile > sortedFile

Then get the top 5 hits for every query:

for next in $(cut -f1 sortedFile | sort -u); do grep -w -m 5 "$next" sortedFile; done > topFivePerQuery
• 0 views • link
updated by Ram Ram • written by 5heikki 5heikki

Yes, the above command made my day!!!!!!!!!

Thanks

• 0 views • link
written by archana.bioinfo87 archana.bioinfo87

Setting -max_target_seqns to 1 will give only 1 subject/hit but several HSPs if they are present.

Setting -max_hsps to 1 will give only 1 HSP per subject but for all subject/hits in the database.

If you really want only 5 HSPs per subject, set the -max_target_seqns to 1 and -max_hsps to 5.

• 0 views • link
updated by Ram Ram • written by Siva Siva

Makes sense, thank you!

• 0 views • link
updated by Ram Ram • written by Adrian Pelin Adrian Pelin

(typo: should be -max_target_seqs instead of -max_target_seqns)

• 0 views • link
updated by Ram Ram • written by al-ash al-ash

There seems to a problem with -outfmt 7. Can you check if this problem persists if you use the default output format?

• 0 views • link
written by Ram Ram

I am trying outfmt 6. Same thing, gives out many results, beyond 5.

• 0 views • link
written by apelin20 apelin20

Try not giving the outfmt param. If it works, then it could be a out format group specific problem.

• 0 views • link
written by Ram Ram

If I am not wrong with -outfmt > 4, -max_target_seqs is ignored. At least this is true for psiblast

• 0 views • link
updated by Ram Ram • written by dago dago

max_target_seqs - Number of aligned sequences to keep. Use with report formats that do not have separate definition line and alignment sections such as tabular (all outfmt > 4). Not compatible with num_descriptions or num_alignments.

• 0 views • link
written by apelin20 apelin20

You are right. I checked my code and it refers to `-num_descriptions`

• 0 views • link
written by dago dago

Log in to answer this question.

More posts like this