Page Comparison

...

module load bowtie/2.2.6

...

$BI/gva_course/mapping/data

...

File Name

...

Description

...

Sample

...

SRR030257_1.fastq

...

Paired-end Illumina, First of pair, FASTQ format

...

You may recognize this as the same files we used for the fastqc and fastx_toolkit tutorial. If you chose to improve the quality of R2 reads using fastx_toolkit 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 "bowtie2mappingBDIB_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:~$ ls
NC_012967.1.gbk  SRR030257_1.fastq  SRR030257_2.fastq

  SRR030257_2.fastq.gz

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

Useful commands

Often you will have general questions about your sequencing files that you want to answer before or after starting your actual analysis. Here we show you some very handy commands after a warning:

Reminder about Linux 1 liners

Below are several commands we've already been using, and some new ones put together to improve your skills.

head * # will show the top 10 lines to give you an idea of the file type/structure
tail * # will show the last 10 lines to determine if the file is a repetitive structure

...

wc -l *  # and then divide by 4 using the your knowledge of fastq files
 
# OROR 
 
grep ^@SRR030257 SRR030257_1.fastq | wc -l
 
# OR
grep --count ^@SRR030257 SRR030257_1.fastq

...

The bp_seqconvert.pl script that is installed as part of Bioperl is one 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. We also took care of this for you when we edited your ~/.profile_userbashrc file in the Linux introduction.

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

module load gcc
module load bioperl
bp_seqconvert.pl

Exercises

...

Bowtie2 is a complete rewrite of bowtie. It is currently the latest and greatest in the eyes of one very picky instructor (and his postdoc/gradstudent) in terms of configurability, sensitivity, and speed.

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.

...

After years of teaching bwa mapping along with bowtie2, we've decided that you will be the first class to use only bowtie2 since we never recommend anyone use bwa. For some more details about the differences between them see the bonus presentation, and if you find a compelling reason to use bwa rather than bowtie2, we'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

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

module spiderlist bowtie
 
module# load bowtie/2.1.0

bowtie2 --version
module list
which bowtie2

This should return the following 2 lines:
# Currently Loaded Modules Matching: bowtie
#  1) bowtie/2.2.6
 
# If you do not see the above output, run these commands, and repeat this block
module spider bowtie
module load bowtie/2.2.6

bowtie2 --version
module list
which bowtie2

Now 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.

...

Try reading the help to figure out how to run the command yourself. Remember these are paired-end reads. 

...

...

So, you will want to submit the full job to the cluster like you learned in the introduction.

But first, try to figure out the command and start it in interactive mode. Remember these are paired-end reads. Use control-c to stop the job once you are sure it is running without an immediate error! Then, submit your command that is working to the TACC queue.

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

launcher_creator.py -n "bowtie2" -t 00:10:00

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 (maybe >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 using a viewer like IGV, which we will cover later.
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.

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 submitted a job that reserved a single node on the cluster, but that node has 12 processors. Bowtie was only using one of those processors (a single "thread")! 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 12 processors for the bowtie process.

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

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 (maybe >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 later.
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 tutorial.

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 learn about the "wayness" of how you request nodes on Lonestar and possibly edit your *.sge scriptslurm script in the future (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.

Optional Exercises for your free time

• 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 performance without causing many fewer reads to be mapped? (increase performance)

...

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. Here is a link to help you return to the GVA 2015 2016 course schedule.