...

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

grep -c "^+$" SRR030257_1.fastq 

...

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

...

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 bowtie. It is currently the latest and greatest in the eyes of one very picky professor (and his postdoc) in 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, 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.

...

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

...

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

...

module list
bowtie2 --version
module list
which bowtie2

which bowtie2

...

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.Why do so many different mapping programs create an index as a first step you may be wondering?

Finally, map the reads! The command you need is:

bowtie2

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

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

...

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.

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:

...

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

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 performance run time 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 2019 2020 course schedule.