A bioinformatics example: BLAST


Teaching: 0 min
Exercises: 15 min
  • Run a real-world bioinformatics application in a Docker container

Running BLAST from a container with Docker

We’ll be running a BLAST (Basic Local Alignment Search Tool) example with a container from BioContainers. BLAST is a tool bioinformaticians use to compare a sample genetic sequence to a database of known seqeuences; it’s one of the most widely used bioinformatics tools.

To begin, we’ll pull the BLAST container (this will take a little bit):

$ docker pull biocontainers/blast:v2.2.31_cv2
v2.2.31_cv2: Pulling from biocontainers/blast
22dc81ace0ea: Pull complete 
1a8b3c87dba3: Pull complete 
Digest: sha256:238717ec69830ec62a19fc05c6f70183f218a13f7678864060f0157dc63dc54f
Status: Downloaded newer image for biocontainers/blast:v2.2.31_cv2

We can run a simple command to verify the container works:

$ docker run biocontainers/blast:v2.2.31_cv2 blastp -help
  blastp [-h] [-help] [-import_search_strategy filename]
   Compute locally optimal Smith-Waterman alignments?

Let’s download some data to start blasting:

$ mkdir blast_example
$ cd blast_example
$ wget http://www.uniprot.org/uniprot/P04156.fasta

This is a human prion FASTA sequence. We’ll also need a reference database to blast against:

$ curl -O ftp://ftp.ncbi.nih.gov/refseq/D_rerio/mRNA_Prot/zebrafish.1.protein.faa.gz
$ gunzip zebrafish.1.protein.faa.gz

We need to prepare the zebrafish database with makeblastdb for the search, so we’ll run the following (see a previous episode for details on -v):

$ docker run -v `pwd`:/data/ biocontainers/blast:v2.2.31_cv2 makeblastdb -in zebrafish.1.protein.faa -dbtype prot

After the container has terminated, you should see several new files in the current directory. We can now do the final alignment step using blastp:

$ docker run -v `pwd`:/data/ biocontainers/blast:v2.2.31_cv2 blastp -query P04156.fasta -db zebrafish.1.protein.faa -out results.txt

The final results are stored in results.txt;

$ less results.txt
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

  XP_017207509.1 protein piccolo isoform X2 [Danio rerio]             43.9    2e-04
  XP_017207511.1 mucin-16 isoform X4 [Danio rerio]                    43.9    2e-04
  XP_021323434.1 protein piccolo isoform X5 [Danio rerio]             43.5    3e-04
  XP_017207510.1 protein piccolo isoform X3 [Danio rerio]             43.5    3e-04
  XP_021323433.1 protein piccolo isoform X1 [Danio rerio]             43.5    3e-04
  XP_009291733.1 protein piccolo isoform X1 [Danio rerio]             43.5    3e-04
  NP_001268391.1 chromodomain-helicase-DNA-binding protein 2 [Dan...  35.8    0.072

We can see that several proteins in the zebrafish genome match those in the human prion (interesting?).

Running BLAST on HPC with Shifter

First, let us pull the BLAST container:

$ module load shifter
$ sg $PAWSEY_PROJECT -c 'shifter pull biocontainers/blast:v2.2.31_cv2'

The following script permits to execute the same bioinformatics example using Shifter and the SLURM scheduler:

#!/bin/bash -l

#SBATCH --account=<your-pawsey-project>
#SBATCH --partition=workq
#SBATCH --ntasks=1
#SBATCH --time=00:05:00
#SBATCH --export=NONE
#SBATCH --job-name=blast

module load shifter

# download sample inputs
wget http://www.uniprot.org/uniprot/P04156.fasta
curl -O ftp://ftp.ncbi.nih.gov/refseq/D_rerio/mRNA_Prot/zebrafish.1.protein.faa.gz
gunzip zebrafish.1.protein.faa.gz

# prepare database
srun --export=all shifter run biocontainers/blast:v2.2.31_cv2 makeblastdb -in zebrafish.1.protein.faa -dbtype prot

# align with BLAST
srun --export=all shifter run biocontainers/blast:v2.2.31_cv2 blastp -query P04156.fasta -db zebrafish.1.protein.faa -out results.txt

Now you can create a test directory somewhere under $MYSCRATCH or $MYGROUP, e.g.

$ mkdir blast_example
$ cd blast_example

use your favourite text editor to copy paste the script above in a file called blast.sh (remember to specify your Pawsey project ID in the script!),

and then submit this script using SLURM:

$ sbatch --reservation <your-pawsey-reservation> blast.sh

Best practices

  • Don’t re-invent the wheel, it’s worth looking to see who’s done what you want already

Key Points

  • There are a lot of applications (not just bioinformatics) already wrapped up in container images

  • Here’s a small list of some of the registries we use at Pawsey:

  • Docker Hub

  • Biocontainers

  • Quay^

  • Nvidia GPU Cloud (NGC)^

  • ^The last two require you to create an account and login to access containers