...
Tool | TACC | Version | Download | Manual | Example | ||
---|---|---|---|---|---|---|---|
Bowtie2 |
You may recall we added this to our .bashrc file yesteray yesterday so it is already loaded | 2.12.06 |
Modules also exist on lonestar5 for bwa.
...
Code Block |
---|
$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.
File Name | Description | Sample |
---|---|---|
| Paired-end Illumina, First of pair, FASTQ format | Re-sequenced E. coli genome |
| Paired-end Illumina, Second of pair, FASTQ format | Re-sequenced E. coli genome |
| Reference Genome in Genbank format | E. coli B strain REL606 |
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.
Code Block |
---|
tacc:~$ ls NC_012967.1.gbk SRR030257_1.fastq SRR030257_2.fastq SRR030257_2.fastq.gz |
Expand | |||||||||
---|---|---|---|---|---|---|---|---|---|
| |||||||||
Remember that to copy an entire folder requires the use of the recursive (-r) option.
|
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:
Warning | ||
---|---|---|
| ||
NGS data can be quite large, a single lane of an Illumina Hi-Seq run generates 2 files each with 100s of millions of lines. Printing all of that can take an enormous amount of time and may crash your terminal long before it finishes. If you find yourself in a seemingly endless scroll of sequence (or anything else for that matter) remember ctrl+c will kill whatever command you just executed |
Reminder about Linux 1 liners
Below are several commands we've already been using, and some new ones put together to improve your skills.
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
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 |
...
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
wc -l * # and then divide by 4 using the your knowledge of fastq files # OR OR 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:
Code Block |
---|
bpmodule load gcc module load bioperl bp_seqconvert.pl |
Exercises
The file NC_012967.1.gbk
is in Genbank
format. The files SRR030257_*.fastq
are in FASTQ
format.
...
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.
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
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):
Expand | |||||||||
---|---|---|---|---|---|---|---|---|---|
| |||||||||
Remember in our earlier tutorial we discussed the use of lonestar's
| |||||||||
Expand | |||||||||
title | OPTIONAL -- How to check what
|
Expand | |||||||
---|---|---|---|---|---|---|---|
| |||||||
Here are a few of the possibilities that will work.
Note that |
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. This command takes a while Remember these are paired-end reads.
Warning | ||
---|---|---|
| ||
This command can take a while (~5 minutes) and is extremely taxing. This is longer than we want to run a job on the head node (especially when all of us are doing it at once). In fact, in previous years, TACC has noticed the spike in usage |
...
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.
Warning | ||||||
---|---|---|---|---|---|---|
| ||||||
Create a Expand | | |||||
| ||||||
Code Block | ||||||
|
Code Block | |||||||||
---|---|---|---|---|---|---|---|---|---|
| |||||||||
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.
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:
- 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
orgrep
or more or using a viewer like IGV, which we will cover later. - 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:
Code Block |
---|
head bowtie2/SRR030257.sam
|
Expand | ||
---|---|---|
| ||
If you thought the answer was the mapping coordinates of the read pairs you were right! |
...
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 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 12 48 processors for the bowtie process.
Expand | |||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| |||||||||||||||||||
You need to use the
Try it out and compare the speed of execution by looking at the log files.and compare the speed of execution by looking at the log files.
|
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.