Blast locally on multiple fasta files with multiple database

Question

[fec2]

Hi all,

I need to run blast locally on multiple fasta files contain in a directory. Previously, I have used command below:

ls *.fasta | parallel -a - blastp -query {} -db my_database -evalue 0.00001 -qcov_hsp_perc 50 -outfmt 6 -max_target_seqs 1 -out {.}.tsv

Since I need to do pairwise comparison and need to blast all fasta files against each other, the command above will need to performe multiple time because I have to specify my database. Is it possible to have a command that is able to blast all fasta file against each other? And generate output file with file name that combined the name of my database plus the query file name?

Thank you.

[GenoMax]

Small nitpick but

Since I need to do pairwise comparison and need to blast all fasta files against each other

does not match the title of this post.

Blast locally on multiple fasta files with multiple database

May want to edit one or other.

[AK]

• the command above will need to performed multiple time because I have to specify my database
• And generate output file name that combined of the name of my database
• -db mydatabase

[fec2]

Hi,

Thank you. Any suggestion for me to change the title of the post?

[fec2]

Thanks all, the command work well!

Answer 1

[AK]

Something like this:

#!/bin/bash
set -euo pipefail

for q in *.fasta; do
  for s in *.fasta; do
    # Ignore self-comparison
    if [[ "${q}" != "${s}" ]]; then
      # Before blastp, check if exist or create database for ${s} (I leave this to you)
      blastp -query ${q} -db ${s} -evalue 0.00001 -qcov_hsp_perc 50 -outfmt 6 -max_target_seqs 1 -out ${q}_2_${s}.tsv
    fi
  done
done

[GenoMax]

This is NOT doing pairwise blast comparisons. One would need to use

blastp -query ${q} -subject ${s} -evalue 0.00001 -qcov_hsp_perc 50 -outfmt 6 -max_target_seqs 1 -out ${q}_2_${s}.tsv

with fasta format files.

[AK]

Thanks, genomx for introducing the -subject argument. Will blastp -query ${q} -subject ${s} build a db for ${s} each time? Or it's without the use of db (how's tradeoff for big fasta file?)

Because as I commented in the response:

Before blastp, check if exist or create database for ${s} (I leave this to you)

which assumed there is a db with prefix ${s}, or build the db if not exist, only the first time.

[h.mon]

Other people already asked about possible differences between -subject and -db:

Different results with BLAST depending on if subject is formatted database or FASTA-file

blast - difference between database and subject

The first post indicates it might result in very different output due to how blast parses the subject in each case, though I would like to see confirmation before believing (and I am feeling lazy to test this today, maybe another day).

[GenoMax]

-subject allows direct comparisons of two sequence files against each other. No database needed.

Answer 2

[h.mon]

I will borrow from SMK loop above to provide a GNU parallel solution.

First save the all vs all filenames into a file:

for q in *.fasta; do
  for s in *.fasta; do
    # Ignore self-comparison
    if [[ "${q}" != "${s}" ]]; then
      echo "${q} ${s}" >> allpairs.txt
    fi
  done
done

Then cat the filenames into parallel to build the multiple blast commands:

cat pairs.txt | parallel --colsep ' ' -j 1 \
    'blastp -query {1} -subject {2} -evalue 0.00001 -qcov_hsp_perc 50 -outfmt 6 -max_target_seqs 1 -out {1.}_vs_{2.}.tsv'

If you have multiple processors cpu cores, you can use -j n to perform n blast searches simultaneously. You don't need to build blast databases, blast can perform fasta vs fasta comparison with -query file1 -subject file2 , but in case each fasta file is large, then it would be best if you build the databases beforehand, and replace -subject {2} by -db {2.}.

edits / comments:

• in view of my comment above, you should check if indeed -subject and -db behave the same way or not.
• with older blast versions, scaling wasn't very good and it was faster to perform n parallel searches than using -num_threads n. This may have improved in recent blast versions (current is BLAST+ 2.9.0), but I didn't test.
• keep in mind if you run several blast searches simultaneously with GNU parallel, memory usage will increase correspondingly compared to just one blast search.

[GenoMax]

Tried a query against the human genome. Was pretty quick. So does not look like OP will need to build databases.

[fec2]

Hi, is it possible to do reverse search using the to match with same BLAST cut-off?

[GenoMax]

If you switch

echo "${q} ${s}" >> allpairs.txt

to

echo "${s} ${q}" >> allpairs_rev.txt

that should get you the reverse comparisons.

[fec2]

Actually I need to do an analysis on 422 genome, but system in my lab can only allow the job to run for 3 days. Therefore, I can't use the other tool online (all of them using phyton and I can't understand the script). Therefore, I wish to use script similar to above (I edited the allpairs.txt file and continue for another 3 days).

Basically they perform BLAST from the query genome against the reference genome with cut-offs of at least 30% identity and at least 70% coverage. Then they took the top match and performed the reverse search using BLAST with the same cut-offs. I don't understand how they use the top match and performed the reverse search.

[GenoMax]

If I understand your question, one does not take just the top match and perform the reverse search. Idea behind a reciprocal blast search is that homologs will show up as (one of) top hits, no matter which direction one does the search in. That confirms their evolutionary relatedness.

[fec2]

Yes, in the website of the online tool, they did mention about reciprocal best hits between two genomic datasets of proteins. Is that mean that I can create a allpairs_rev.txt, and average 2 values?

[GenoMax]

I am not sure what two values you are referring to. What exactly are you trying to do?

[fec2]

I have to do an analysis called average nucleotide identity (AAI). There are many tool available, for example compareM that I previously used. However, this time I need to analyse 422 genome, but my system in the lab can only allows a job to run for 3 days. I have tried a few tools but non of these tools can finish the job in 3 days. So I am thinking to do it manually using some commands like in this post. The command in this post is for another analysis called percentage of conserved protein. Because I can edit the allpairs.txt file and continue the BLAST after 3 days, so I able to do it manually and continuously until the BLAST job done for all genome. But for AAI, I dont really understand the meaning of reciprocal best hits between two genomic datasets of proteins.

[fec2]

From I checked online, the AAI analysis actually only using a representative of paralog genes (with highest identity), so the total number of genes calculated will also deduct the the paralogous gene that not included in the analysis.