Overview

Once you know you are working with the best quality data (Evaluating Raw Sequencing data tutorial) possible, the first step in nearly every NGS analysis pipeline is to map sequencing reads to a reference genome. In this tutorial we'll explore these basic principles using bowtie2 on TACC.

The world of read mappers seems to be settling down a bit after being a bioinformatics Wild West where there was a new gun in town every week that promised to be a faster and more accurate shot than the current record holder. Things seem to have reached the point where there is mainly a trade-off between speed, accuracy, and configurability among read mappers that have remained popular. There are over 50 read mapping programs listed here. Each mapper has its own set of limitations (on the lengths of reads it accepts, on how it outputs read alignments, on how many mismatches there can be, on whether it produces gapped alignments). It is possible a different read mapper would be better for your set of experiments. More will be discussed about selecting a good tool on Friday.

Learning Objectives

This tutorial covers the commands necessary to use bowtie2 to map reads to a reference genome, and concepts applicable to many more mappers.

1. Become comfortable with the basic steps of indexing a reference genome, mapping reads, and converting output to SAM/BAM format for downstream analysis.
2. Use bowtie2 to map reads from an E. coli Illumina data set to a reference genome and compare the output.

Theory

Please see the Introduction to mapping presentation on the course outline for more details of the theory behind read mapping algorithms and critical considerations for using these tools and references correctly.

Mapping tools summary

The tutorial currently available on  the Lonestar cluster at TACC is as follows:

module load bowtie/2.3.4

Modules also exist on lonestar5 for bwa.

Tutorial: E. coli genome re-sequencing data

The following DNA sequencing read data files were downloaded from the NCBI Sequence Read Archive via the corresponding European Nucleotide Archive record. They are Illumina Genome Analyzer sequencing of a paired-end library from a (haploid) E. coli clone that was isolated from a population of bacteria that had evolved for 20,000 generations in the laboratory as part of a long-term evolution experiment (Barrick et al, 2009). The reference genome is the ancestor of this E. coli population (strain REL606), so we expect the read sample to have differences from this reference that correspond to mutations that arose during the evolution experiment.

Transferring Data

We have already downloaded data files for this example and put them in the path:

$BI/gva_course/mapping/data

You may recognize this as the same files we used for the fastqc and cutadapt tutorial. If you chose to improve the quality of R2 reads using cutadapt as you did for R1 in the tutorial, you could use the improved reads in this tutorial to see what a difference the improved reads can make for read mapping. 

The easiest way to run the tutorial is to copy this entire directory into a new folder called "GVA_bowtie2_mapping" on your $SCRATCH space and then run all of the commands from inside that directory. See if you can figure out how to do that. When you're in the right place, you should get output like this from the ls command.

tacc:/scratch/<#>/<UserName>/GVA_bowtie2_mapping$ ls
NC_012967.1.gbk  SRR030257_1.fastq  SRR030257_2.fastq  SRR030257_2.fastq.gz

cds
cp -r $BI/gva_course/mapping/data GVA_bowtie2_mapping
cd GVA_bowtie2_mapping
ls

Reminders about working with sequencing files

Remember, from the introduction tutorial, there are multiple ways to look at our sequencing files without using cat:

grep -c "^+$" SRR030257_1.fastq 

sed -n 2p SRR030257_1.fastq | awk -F"[ATCGNatcgn]" '{print NF-1}'

Converting sequence file formats

Occasionally you might download a sequence or have it emailed to you by a collaborator in one format, and then the program that you want to use demands that it be in another format. Why do they have to be so picky? Everybody has own favorite formats and/or those that they are the most familiar with but humans can typically pick the information they need out of comparable formats. Programs can only be written to assume a single type of format (or allow you to specify a format if the author is particularly generous), and can only find things in single locations based on that format. 

The bp_seqconvert.pl script is a common script written in Bioperl that is a helpful utility for converting between many common sequence formats. On TACC, the Bioperl modules are installed, but the helper script isn't. So, we've put it in a place that you can run it from for your convenience. However, remember that any time that you use the script you must have the bioperl module loaded. As it is fairly rare that you need to convert sequence files between different format, bioperl is actually not listed as one of the modules on your .bashrc file in your $HOME directory that you set up yesterday. Additionally it gives an opportunity to have you working with the module system. If you find yourself needing to do lots of sequence conversions, you may want to add a 'module load bioperl/1.007002' line to your .bashrc file.

Run the script without any arguments to get the help message:

module load gcc
module load bioperl
bp_seqconvert.pl

Convert a gbk reference to a embl reference

Convert the Genbank file NC_012967.1.gbk to EMBL format, and name this new file NC_012967.1.embl.

bp_seqconvert.pl --from genbank --to embl < NC_012967.1.gbk > NC_012967.1.embl
head -n 100 NC_012967.1.embl

Use of uninitialized value in substitution (s///) at /opt/apps/bioperl/1.6.901/Bio/SeqIO/embl.pm line 777, <STDIN> line 164674.
Use of uninitialized value in concatenation (.) or string at /opt/apps/bioperl/1.6.901/Bio/SeqIO/embl.pm line 779, <STDIN> line 164674.

head -n 100 NC_012967.1.embl

less NC_012967.1.embl

more NC_012967.1.embl

Converting from fastq to fasta format

Sometimes you only want to work with a subset of a full data file to check for overall trends, or to try out a new piece of code. Convert only the first 10,000 lines of SRR030257_1.fastq to FASTA format.

head -n 10000 SRR030257_1.fastq | bp_seqconvert.pl --from fastq --to fasta > SRR030257_1.fasta

head SRR030257_1.fastq
head SRR030257_1.fasta

Mapping with bowtie2

Bowtie2 is a complete rewrite of an older program bowtie. In terms of configurability, sensitivity, and speed it is useful for a wide range of projects. After years of teaching bwa mapping along with bowtie2, bowtie2 alone is now taught as I never recommend anyone use bwa, and based on positive feedback we continue with this set up. For some more details about how read mappers work see the bonus presentation about read mapping details and file formats on the course home page, and if you find a compelling reason to use bwa (or any other read mapper) rather than bowtie2 after the course is over, I'd love to hear from you.

Create a fresh output directory named bowtie2. We are going to create a specific output directory for the bowtie2 mapper within the directory that has the input files so that you can compare the results of other mappers if you choose to do the other tutorials.

mkdir bowtie2

First you need to convert the reference file from GenBank to FASTA using what you learned above. Name the new output file NC_012967.1.fasta and put it in the same directory as NC_012967.1.gbk.

bp_seqconvert.pl --from genbank --to fasta < NC_012967.1.gbk > NC_012967.1.fasta

Next, we want to make sure the bowtie2 module is loaded (we use module spider to get the current name, which may not be bowtie/2.3.4 if you re-run this tutorial later):

module list bowtie

# Currently Loaded Modules Matching: bowtie
#  None found.
# Inactive Modules Matching: bowtie
#   1) bowtie/2.3.4

Further, when we try to load bowtie/2.3.4 we get an error message.

Lmod has detected the following error:  These module(s) or extension(s) exist but cannot be loaded as requested: "bowtie/2.3.4"
   Try: "module spider bowtie/2.3.4" to see how to load the module(s).

See if you can figure out how to load bowtie using the information above.

module load intel/18.0.2
module load bowtie/2.3.4

module list
bowtie2 --version
which bowtie2

For many read mappers, the first step is quite often indexing the reference file regardless of what mapping program is used. Put the output of this command into the bowtie directory we created a minute ago. The command you need is:

bowtie2-build

Try typing this alone in the terminal and figuring out what to do from the help show just from typing the command by itself.

bowtie2-build NC_012967.1.fasta bowtie2/NC_012967.1

Take a look at your output directory using ls bowtie2 to see what new files have appeared. These files are binary files, so looking at them with head or tail isn't instructive and can cause issues with your terminal. If you insist on looking at them (or accidentally do so before you read this) and your terminal begins behaving oddly, simply close it and log back into lonestar with a new terminal, and start a new idev session.

Finally, it's time to map the reads! 

bowtie2

bowtie2 -t -x bowtie2/NC_012967.1 -1 SRR030257_1.fastq -2 SRR030257_2.fastq -S bowtie2/SRR030257.sam  # the -t command is not required for the mapping, but it can be particularly informative when you begin comparing different mappers

Your final output file is in SAM format. It's just a text file, so you can peek at it and see what it's like inside. Two warnings though:

1. SAM files can be enormously humongous text files (potentially >1 gigabytes). Attempting to open the entire file at once can cause your computer to lock up or your text editor to crash. You are generally safer only looking at a portion at a time using linux commands like head or grep or more or using a viewer like IGV, which we will cover in a later tutorial.
2. SAM files have some rather complicated information encoded as text, like a binary encoded FLAGS field and CIGAR strings. We'll take a look at some of these later, if we have time, or they are covered in the bonus presentation about read mapping and file formats which you can find on the home page.

Still, you should recognize some of the information on a line in a SAM file from the input FASTQ, and some of the other information is relatively straightforward to understand, like the position where the read mapped. Give this a try:

head bowtie2/SRR030257.sam

More reading about SAM files

• The official SAM file specification
• http://genome.sph.umich.edu/wiki/SAM

Multithreaded execution

We have actually massively under-utilized Lonestar in this example. We ran the command using only a single processor (a single "thread") rather than the 48 we have available. For programs that support multithreaded execution (and most mappers do because they are obsessed with speed) we could have sped things up by using all 48 processors for the bowtie process.

bowtie2 -t -p 48 -x bowtie2/NC_012967.1 -1 SRR030257_1.fastq -2 SRR030257_2.fastq -S bowtie2/SRR030257.sam

If you want to launch many processes as part of one job, so that they are distributed one per node and use the maximum number of processors available, then you need to think about the "wayness" of how you request nodes on Lonestar (we'll go over this more on this on Friday), or make better use of running commands in the background using the & symbol at the end of the command line.

One consequence of using multithreading that might be confusing is that the aligned reads might appear in your output SAM file in a different order than they were in the input FASTQ. This happens because small sets of reads get continuously packaged, "sent" to the different processors, and whichever set "returns" fastest is written first. You can force them to appear in the same order (at a slight cost in speed) by adding the --reorder flag to your command, but is typically only necessary if the reads are already ordered or you intend to do some comparison between the input and output. 

What comes after mapping?

The next steps are often to view the output using a specific viewer on your local machine, or to begin identifying variant locations where the reads differ from the reference sequence. These will be the next things we cover in the course.

Optional exercises 

• In the bowtie2 example, we mapped in --local mode. Try mapping in --end-to-end mode (aka global mode).

• Do the BWA tutorial so you can compare their outputs.
  • Did bowtie2 or BWA map more reads?
  • In our examples, we mapped in paired-end mode. Try to figure out how to map the reads in single-end mode and create this output.
  • Which aligner took less time to run? Are there any options you can change that:
    • Lead to a larger percentage of the reads being mapped? (increase sensitivity)
    • Speed up run time without causing many fewer reads to be mapped? (increase performance)

Here is a link to help you return to the GVA 2020 course schedule.