Repository navigation
BLAST and DIAMOND Lab
If you haven't already, log into isaac. Stay in your home directory.
ssh yourusername@login.isaac.utk.edu
Some software packages are available through the module system (more info from OIT here, although they may be out of date.
Lets see what cool software we can play with! This is everything module has available:
module avail
Return to scroll down, q to quit.
Search for something specific (although this is a terrible search, it has to exactly match the first part of the package name)
module avail blast
Can we search the output of the list results in a better way? This is a great example of how passing (piping) the output of one command to another command can be fast and useful. Lets say you are looking for a package with "mpi" in the name
module avail mpi
module avail | grep "mpi"
Grep? Grep is short hand for "globally search a regular expression and print". So what exactly is a regular expression? Well, it's a formal language to define search patterns. Regular expressions (also known as regex) are extremely powerful, and if you are going to be doing a lot of file manipulation and data science, worth learning more about!
Is this the most recent version of blast-plus? Lets find out...
Finally, we actually load the existing software package:
module load blast-plus
Is it working?
blastn -version
blastp -help
BLAST is one of the most common bioinformatics tools out there, and many biology students learn to run BLAST for individual or small sets of sequences via a web interface: you can find a quick tutorial of web BLAST here. The most popular is probably NCBI BLAST. However, you can perform many more BLAST computations on an HPC.
BLAST has a number of possible programs to run depending on whether you have nucleotide or protein sequences:
- nucleotide query and nucleotide db - blastn
- nucleotide query and nucleotide db - tblastx (includes six frame translation of query and db sequences)
- nucleotide query and protein db - blastx (includes six frame translation of query sequences)
- protein query and nucleotide db - tblastn (includes six frame translation of db sequences)
- protein query and protein db - blastp
Make a directory to hold our blast lesson.
mkdir blast_examples
cd blast_examples
Let's start by downloading some data. Grab some cow and human RefSeq proteins from NCBI:
wget ftp://ftp.ncbi.nih.gov/refseq/B_taurus/mRNA_Prot/cow.1.protein.faa.gz
wget ftp://ftp.ncbi.nih.gov/refseq/H_sapiens/mRNA_Prot/human.1.protein.faa.gz
This is only part of the human and cow protein files. check out how big they are:
ls -lh
The database files are both gzipped, so lets unzip them
cd blast_examples
gunzip *gz
ls -lh
How much bigger are they now?
Take a look at the head of each file:
head cow.1.protein.faa
head human.1.protein.faa
Note that the files are in fasta format, even though they end if ".faa" instead of the usual ".fasta". This NCBI's way of denoting that this is a fasta file with amino acids instead of nucleotides.
Learn more about fasta format. This is one of the most common file formats in bioinformatics.
How many sequences are in each one?
grep -c '^>' cow.1.protein.faa
grep -c '^>' human.1.protein.faa
We can practice with wildcards:
grep -c '^>' *faa
This grep command uses the c flag, which reports a count of lines with a match to the pattern. In this case, the pattern is a regular expression, meaning match only lines that begin with a >. In this case, we are asking to search for lines that begin with a greater than sign using the regex convention of the '^'symbol.
This is a bit too big, let's take a smaller set for practice. Let's take the first 5 sequences of the human proteins, which we can see are on the first 19 lines
head -n 20 human.1.protein.faa
head -n 19 human.1.protein.faa > human.small.faa
cat human.small.faa
Now we can blast these five human sequences against the set of cow sequences. First, we need to tell blast about our database. BLAST needs to do some pre-work on the database file prior to searching. This helps to make the software work a lot faster. Use the makeblastdb command:
makeblastdb -in cow.1.protein.faa -dbtype prot
ls
Note that this makes a lot of extra files, with the same name as the database plus new extensions (.pin, .psq, etc). To make blast work, these files, called index files, must be in the same directory as the fasta file.
Now we can run the blast job. We will use blastp, which is appropriate for protein to protein comparisons.
blastp -query human.small.faa -db cow.1.protein.faa -out human_vs_cow_blastp_results.txt
Take a look at the results in nano.
Lets also take a look at the help pages. Unfortunately, there are no man pages (those are usually reserved for shell commands, but some software authors will provide them as well), but there is a text help output
blastp -help
To scroll through slowly
blastp -help | less
To quit the less screen, press the q key.
Parameters of interest: -evalue ( Default is 10?!?) and -outfmt
Lets get only very meaningful matches with a different output format:
blastp \
-query human.small.faa \
-db cow.1.protein.faa \
-out human_vs_cow_blastp_results.tsv \
-evalue .001 \
-outfmt 7
Whoa!? What are those slashes? In a practical sense, these tell the shell "I'm not done with this command yet". In a literal sense, the forward slash is an "escape" character, and it's telling the shell to interpret the newline as just a newline, not a "submit" of the command.
Check out the results with head, cat and/or nano.
Tab-delimited has these default columns (-outfmt 7):
qseqid Query sequence ID
sseqid Subject (ie DB) sequence ID
pident Percent Identity across the alignment
length Alignment length
mismatch # of mismatches
gapopen Number of gap openings
qstart Start of alignment in query
qend End of alignment in query
sstart Start of alignment in subject
send End of alignment in subject
evalue E-value
bitscore Bit score
I find it very useful to add the subject sequence description to the tabular output. If you go through help some more, you will find that you can decide which types of output you would like in the tab-delimited file and what order they should be in. For example:
blastp \
-query human.small.faa \
-db cow.1.protein.faa \
-out human_vs_cow_blastp_results.tsv \
-evalue .0001 \
-outfmt "7 std stitle"
Lets try a medium sized data set next
head -n 200 human.1.protein.faa > human.medium.faa
What size is this db?
grep -c '^>' human.medium.faa
Lets run the blast again, but this time lets not have the comment lines, and lets return only the best hit for each query, and lets filter by e value of 1e-6
blastp \
-query human.medium.faa \
-db cow.1.protein.faa \
-out human_medium_vs_cow_blastp_results.tsv \
-evalue 1e-10 \
-outfmt "6 std stitle" \
-max_target_seqs 1
I picked "max target seqs" for a reason! It's quite a misleading flag - it finds the number of target sequences but they may not be the best matches. See this blog post for more details..
How many results did we get?
wc -l cow_vs_human_blast_results.tsv
This is still very small compared to what we can do. Lets imagine needing to BLAST 50,000 sequences against a database of 10,000,000 proteins (a surprisingly common task!). A single command can do that.
An alternative basic alignment tool to BLAST is diamond. Currently, diamond only works for proteins and translated nucleic acid alignments. See the diamond documentation. Diamond runs much faster than BLAST, but hasn't been around as long.
Module doesnt have diamond, so we'll download it
wget https://github.com/bbuchfink/diamond/releases/download/v2.2.5/diamond-linux64.tar.gz
tar xzf diamond-linux64.tar.gz
We will need to specify its location, because we installed it and the system doesn't know where it is. (To learn more about how to fix this problem go find out about the $PATH environmental variable and how to add to it in your .bashrc)
Test it out
diamond --version
./diamond --version
/nfs/home/mstaton1/blast_examples/diamond --help
Similarly to BLAST, diamond also creates a database:
./diamond makedb --in cow.1.protein.faa --db cow_db
Next run diamond search. Note: here we are using the block flag 'b1' to denote amount of memory for diamond to use. According to diamond's documentation, a block size of 1 (b1) corresponds with 1 billion sequence letters to process at a time, and the ram usage in Gb is approximately 6X this value (e.g. using b2 would use 12Gb, and b3 would use 18gb).
./diamond blastp -q human.medium.faa -d cow_db.dmnd -b1 -o human_vs_cow_diamond_results.tsv
Always read the manual and consider memory usage! Also, whenever using a newer software on GitHub, be sure to swing by the issues page to see if there are any concerns by fellow users than might influence your analysis.
We can check the time savings using the time command
time ./diamond blastp -q human.medium.faa -d cow_db.dmnd -b1 -o human_vs_cow_diamond_results.tsv
time blastp -query human.medium.faa -db cow.1.protein.faa -out human_vs_cow_blastp_results.txt
Why are there three times?
- Real-time (real). The real-life time it takes for the process to run from start to finish. This includes any time taken by other processes and the time spent waiting for them to be complete.
- User time (user). The amount of CPU time spent in user mode during the process. Other processes and blocked time are not included.
- System time (sys). The total CPU time spent in kernel mode during the process. Similar to user time, other processes and time spent blocked by other processes are not counted.
User + Sys is typically used for figuring out the actual resource usage, but for this class and this application we can use real as a perfectly fine proxy.
How much faster is diamond? How about in --very-sensitive mode?
time ./diamond blastp -q human.medium.faa -d cow_db.dmnd -b1 -o human_vs_cow_diamond_results.tsv --very-sensitive