0
$\begingroup$

I had created databases of different sets of metagenome datasets - one with 6 runs, other with 48 runs, another with 100 runs, etc using the accession list for each of these datasets and makeblastdb command. Since I wanted to check the abundance of one specific 16S rRNA gene sequence, I ran a blastn command using this 16S sequence (F7) as the query against each of the created databases with a txt file (F7_droughtS.txt, F7_Dy.txt, so on) as the output. Ideally I should have got differently sized files with different number of OTU matches. I did get different files since the SRR ids were different for each database, but the results did not go beyond 500 lines for any database search. Can anyone explain why this is so?

Following is the command I used for blastn of: sequence.fasta against database firstdb :

blastn -query sequence.fasta -db firstdb -outfmt 6 -out output1.txt

Also ran the following to know the default number of lines : blastn -help -max_target_seqs

Number of database sequences to show one-line descriptions for

Not applicable for outfmt > 4 Default = 500' * Incompatible with: max_target_seqs -num_alignments <Integer, >=0> Number of database sequences to show alignments for Default = 250' * Incompatible with: max_target_seqs

$\endgroup$
6
  • 2
    $\begingroup$ I suggest you have a look at the settings for the blast version you are using. These can usually be accessed with 'blastn -help'. In particular, '-max_target_seqs' may be of interest. $\endgroup$
    – Laura
    May 9, 2023 at 8:54
  • 1
    $\begingroup$ Please edit your question and i) remove the images of the files, that isn't telling us anything. You can show us the output of wc -l file to show the number of lines. ii) show us the exact commands you used. You are probably using a version of blast that shows the top 500 hits by default, but we need to know what blast you use and what exact commands you ran in order to help. $\endgroup$
    – terdon
    May 9, 2023 at 10:43
  • 1
    $\begingroup$ @Laura, Thank you for your suggestions. When I did run the code provided by you, this is what I got : Number of database sequences to show one-line descriptions for Not applicable for outfmt > 4 Default = 500' * Incompatible with: max_target_seqs -num_alignments <Integer, >=0> Number of database sequences to show alignments for Default = 250' * Incompatible with: max_target_seqs How can I change these settings? $\endgroup$
    – K_081
    May 14, 2023 at 14:55
  • 1
    $\begingroup$ Not an expert on 16S, but is BLAST really the correct tool for this(?) Maybe one of these tool suites like QIME $\endgroup$ May 16, 2023 at 17:40
  • $\begingroup$ QIIME is an interesting suggestion. Here I'd say BLAST - although there are performance issues - offset parallelising. BLAST is a known approach to scoop all the genetic diversity for a given locus and 16S is a taxonomic goldstandard. This isn't the analysis - the analytics are in the next step, but how you've prepped the data for the next step needs to be understood and be solid. There are a lot of other approaches but blast has flexibility e.g. to shift tblastn .. etc ... and its why I like it. $\endgroup$
    – M__
    May 16, 2023 at 19:26

1 Answer 1

1
$\begingroup$

I've reproduced the issue. The solution is to simply set -max_target_seqs to whatever number of hits you want. E.g. -max_target_seqs 1000 gives 1000 and so on ... The problem is specific to the -outfmt option and not default searches.

blastn version I'm using v2.12.0, Feb 8 2022


Proof My database contains around 260000 fasta sequences.

makeblastdb -in mystuff.aln.0.fa   -parse_seqids -blastdb_version 5  -out test.db -title "test" -dbtype nucl 
blastn -query query.fas -db test.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

Then wc blasthits.txt

>5000 hits

If we now use -outfmt 6, but remove the stringency of -evalue 1e-6 then ...

blastn -query query.fa -db test.db -outfmt 6 -num_threads 4 -out blasthits.txt

Now perform wc blasthits.txt

500

So the OP is correct, I didn't know -outfmt 6 restricts the output to 500.

If we now introduce the max_target_seqs 1000 parameter, suggested by @Laura then,

blastn -query query.fa -db test.db -outfmt 6 -max_target_seqs 1000 -num_threads 4 -out blasthits.txt

Then wc blasthits.txt

1000

If its 550 wc ... is 550 etc...

Thats it. Note in my blasts I usually parallelise across 4 processes i.e. -num_threads 4.

$\endgroup$
2
  • 1
    $\begingroup$ Thank you very much, it worked! $\endgroup$
    – K_081
    May 17, 2023 at 15:13
  • $\begingroup$ Thanks @KirtiKulanthaivel glad it worked and much appreciated. $\endgroup$
    – M__
    May 19, 2023 at 8:30

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.