BLAST

BLAST ("Basic Local Alignment Search Tool") finds regions of similarity between biological sequences. The program compares nucleotide or protein sequences to sequence databases and calculates the statistical significance.

User manual

You can find more information its arguments in the user manual. Or with :

[name@server ~]$ blastn -help

Databases

Some frequently-used sequence databases are installed on Compute Canada clusters. You can find information on the BLAST databases available here : Genomics data.

Accelerating the search

For the examples below, the file ref.fa will be used as the reference database in FASTA format and seq.fa as the queries we are looking for.

makeblastdb

Before running a search, we must build the database. Building the database can be a preprocessing job, where the other job are dependent on the completion of the makeblastdb job. Here is an example of a submission script:

File : makeblastdb.sh

#!/bin/bash

#SBATCH --account=def-<user>  # The account to use
#SBATCH --time=00:02:00       # The duration in HH:MM:SS format
#SBATCH --cpus-per-task=1     # The number of cores
#SBATCH --mem=512M            # Total memory for this task

module load gcc/7.3.0 blast+/2.7.1

# Create the nucleotide database based on `ref.fa`.
makeblastdb -in ref.fa -title reference -dbtype nucl -out ref.fa

Task array

Blast search can greatly benefit from data parallelism by splitting the query file into multiples queries and running these queries against the database.

Preprocess

In order to accelerate the search, the seq.fa file must be split into smaller chunks. The file chunks should be at least 1Mb or greater, but not smaller as it may hurt the parallel file system.

Important: To correctly split a FASTA format file, it must be in its original format and not multiline format. In other words, the sequence must be on one line.

Using the split utility:

[name@server ~]$ split -d -a 1 -l 2 seq.fa seq.fa.

will create 10 files named seq.fa.N where N is in the range of [0..9] for 10 queries (sequences).

Job submission

Once our queries are split, we can create a task for each seq.fa.N file using a job array. The task id from the array will map to the file name containing the query to run.

This solution allows the scheduler to fit the smaller jobs from the array where there's available resources in the cluster.

File : blastn_array.sh

#!/bin/bash

#SBATCH --account=def-<user>  # The account to use
#SBATCH --time=00:02:00       # The duration in HH:MM:SS format of each task in the array
#SBATCH --cpus-per-task=1     # The number of cores for each task in the array
#SBATCH --mem-per-cpu=512M    # The memory per core for each task in the array
#SBATCH --array=0-9           # The number of tasks: 10

module load gcc/7.3.0 blast+/2.7.1

# Using the index of the current task, given by `$SLURM_ARRAY_TASK_ID`, run the corresponding query and write the result
blastn -db ref.fa -query seq.fa.${SLURM_ARRAY_TASK_ID} > seq.ref.${SLURM_ARRAY_TASK_ID}

With the above submission script, we can submit our blast search and it will run after the creation of the database:

[name@server ~]$ sbatch --dependency=afterok:$(sbatch makeblastdb.sh) blastn_array.sh

Once all the tasks from the array are done, the results can be concatenated using:

[name@server ~]$ cat seq.ref.{0..9} > seq.ref

where the 10 files will be catenated into seq.ref file. This could be done from the login node or as a dependent job upon completion of all the tasks from the array.

GNU Parallel

GNU parallel is a great tool to pack many small jobs into one and parallelize it. This solution helps alleviate the issue of too many small files on a parallel file system by querying fixed size chunks from seq.fa and run on one node and multiple cores.

As an example, if your seq.fa file is 3MB, you could read block of 1MB and GNU Parallel will create 3 jobs, thus using 3 cores. If we would have requested 10 cores in our task, we would be wasting 7 cores. Therefore, the block size is important. We can also let GNU Parallel decide, as done below.

See also handling large files with GNU Parallel.

File : blastn_gnu.sh

#!/bin/bash

#SBATCH --account=def-<user>  # The account to use
#SBATCH --time=00:02:00       # The duration in HH:MM:SS format
#SBATCH --cpus-per-task=4     # The number of cores
#SBATCH --mem-per-cpu=512M    # The memory per core

module load gcc/7.3.0 blast+/2.7.1

cmd='blastn -db ref.fa -query - '

# Using the `::::` notation, give the sequences file to GNU parallel
# where
#   --jobs number of core to use, equal $SLURM_CPUS_PER_TASK (the number of cores requested)
#   --keep-order keep same order as given in input
#   --block -1 let GNU Parallel evaluate the block size and adapt
#   --recstart record start, here the sequence identifier `>`
#   --pipepart pipe parts of $cmd together. 
#              `--pipepart` is faster than `--pipe` (which is limited to 500MB/s) as `--pipepart` can easily go to 5GB/s according to Ole Tange.
# and redirect results in `seq.ref`.
parallel --jobs $SLURM_CPUS_PER_TASK --keep-order --block -1 --recstart '>' --pipepart $cmd :::: seq.fa > seq.ref

Note: the file must not be compressed.

Job submission

With the above submission script, we can submit our blast search and it will run after the creation of the database:

[name@server ~]$ sbatch --dependency=afterok:$(sbatch makeblastdb.sh) blastn_gnu.sh

Additional tips

• If it fits into the node-local storage, copy your FASTA database to the localscratch ($SLURM_TMPDIR).
• Lower the number of hits returned (-max_target_seqs, -max_hsps can help), if it is reasonable for your research.
• Limit your hit list using evalue filters to near identical hits (-evalue), if it is reasonable for your research.