More Alignment exercises

Exercise #4: Bowtie2 global alignment - Vibrio cholerae RNA-seq

While we have focused on aligning eukaryotic data, the same tools can be used with prokaryotic data. The major differences are less about the underlying data and much more about the external/public databases that store and distribute reference data. If we want to study a prokaryote, the reference data is usually downloaded from a resource like GenBank.  

In this exercise, we will use some RNA-seq data from Vibrio cholerae, published on GEO here, and align it to a reference genome.

Overview of Vibrio cholerae alignment workflow with Bowtie2

Alignment of this prokaryotic data follows the workflow below. Here we will concentrate on steps 1 and 2.

1. Prepare the vibCho reference index for bowtie2 from GenBank records
2. Align reads using bowtie2, producing a SAM file
3. Convert the SAM file to a BAM file (samtools view) 
4. Sort the BAM file by genomic location (samtools sort)
5. Index the BAM file (samtools index)
6. Gather simple alignment statistics (samtools flagstat and samtools idxstat)

Obtaining the GenBank records

First prepare a directory to work in, and change to it:

mkdir -p $SCRATCH/core_ngs/references/vibCho
cd $SCRATCH/core_ngs/references/vibCho

V. cholerae has two chromosomes. We download each separately.

1. Navigate to http://www.ncbi.nlm.nih.gov/nuccore/NC_012582
   • click on the Send to down arrow (top right of page)
     • select Complete Record
     • select File as Destination, and Format FASTA
     • click Create File
       
   • in the Opening File dialog, select Save File then OK
     • Save the file on your local computer as NC_012582.fa
2. Back on the main http://www.ncbi.nlm.nih.gov/nuccore/NC_012582 page
   • click on the Send to down arrow (top right of page)
     • select Complete Record
     • select File as Destination, and Format GFF3
     • click Create File
       
   • in the Opening File dialog, select Save File then OK
     • Save the file on your local computer as NC_012582.gff3
3. Repeat steps 1 and 2 for the 2nd chromosome
   
   • NCBI URL is http://www.ncbi.nlm.nih.gov/nuccore/NC_012583
   • use NC_012583 as the filename prefix for the files you save
   • you should now have 4 files:
     • NC_012582.fa, NC_012582.gff3
     • NC_012583.fa, NC_012583.gff3
       
4. Transfer the files from your local computer to TACC
   • to the ~/scratch/core_ngs/references/vibCho directory created above
     
     • On a Mac or Windows with a WSL shell, use scp from your laptop
       
     • On Windows, use the pscp.exe PuTTy tool
     • See Copying files between TACC and your laptop

mkdir -p $SCRATCH/core_ngs/references/vibCho
cd $SCRATCH/core_ngs/references/vibCho
cp $CORENGS/ncbi/vibCho/NC_* .

Once you have the 4 files locally in your $SCRATCH/core_ngs/references/vibCho directory, combine them using cat:

cd $SCRATCH/core_ngs/references/vibCho
cat NC_01258[23].fa   > vibCho.O395.fa
cat NC_01258[23].gff3 > vibCho.O395.gff3

Now we have a reference sequence file that we can use with the bowtie2 reference builder, and ultimately align sequence data against.

Introducing bowtie2

First make sure you're in an idev session:

idev -m 180 -p normal -A UT-2015-05-18 -N 1 -n 68 

Go ahead and load the bowtie2 module so we can examine some help pages and options.

module biocontainers
module load bowtie

 Now that it's loaded, check out the options. There are a lot of them! In fact for the full range of options and their meaning, Google "Bowtie2 manual" and bring up that page (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml). The Table of Contents is several pages long! Ouch!

This is the key to using bowtie2 - it allows you to control almost everything about its behavior, which make it the go-to aligner for specialized alignment tasks (e.g. aligning miRNA or other small reads). But it also makes it is much more challenging to use than bwa – and it's easier to screw things up too!

Building the bowtie2 vibCho index

Before the alignment, of course, we've got to build a bowtie2- specific index using bowtie2-build. Go ahead and check out its options. Unlike for the aligner itself, we only need to worry about a few things here:

• reference_in file is just the vibCho.O395.fa FASTA we built from GenBank records
• bt2_index_base is the prefix of all the bowtie2-build output file

Here, to build the reference index for alignment, we only need the FASTA file. (This is not always true - extensively spliced transcriptomes requires splice junction annotations to align RNA-seq data properly.)

First create a directory specifically for the bowtie2 index, then build the index using bowtie-build.

mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cd $SCRATCH/core_ngs/references/bt2/vibCho

# Symlink to the fasta file you created
ln -sf $SCRATCH/core_ngs/references/vibCho.O395.fa
# or, to catch up:
ln -sf $CORENGS/references/vibCho.O395.fa

bowtie2-build vibCho.O395.fa vibCho.O395

This should also go pretty fast. You can see the resulting files using ls like before.

Performing the bowtie2 alignment

We'll set up a new directory to perform the V. cholerae data alignment. But first make sure you have the FASTQ file to align and the vibCho bowtie2 index:

# Get the FASTQ to align
mkdir -p $SCRATCH/core_ngs/alignment/fastq
cp $CORENGS/alignment/*fastq.gz $SCRATCH/core_ngs/alignment/fastq/

# Set up the bowtie2 index
mkdir -p $SCRATCH/core_ngs/references/bt2/vibCho
cp $CORENGS/idx/bt2/vibCho/*.*  $SCRATCH/core_ngs/references/bt2/vibCho

Make sure you're in an idev session with the bowtie2 BioContainers module loaded:

idev -p normal -m 120 -A UT-2015-05-18 -N 1 -n 68
module load biocontainers
module load bowtie

Now set up a directlry to do this alignment, with symbolic links to the bowtie2 index directory and the directory containing the FASTQ to align:

mkdir -p $SCRATCH/core_ngs/alignment/vibCho
cd $SCRATCH/core_ngs/alignment/vibCho
ln -sf ../../references/bt2/vibCho
ln -sf ../../alignment/fastq fq

We'll be aligning the V. cholerae reads now in ./fq/cholera_rnaseq.fastq.gz (how many sequences does it contain?)

Note that here the data is from standard mRNA sequencing, meaning that the DNA fragments are typically longer than the reads. There is likely to be very little contamination that would require using a local rather than global alignment, or many other pre-processing steps (e.g. adapter trimming). Thus, we will run bowtie2 with default parameters, omitting options other than the input, output, and reference index. This performs a global alignment.

As you can tell from looking at the bowtie2 help message, the general syntax looks like this:

bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]

So execute this bowtie2 global, single-end alignment command:

cd $SCRATCH/core_ngs/alignment/vibCho
bowtie2 -x vibCho/vibCho.O395 -U fq/cholera_rnaseq.fastq.gz -S cholera_rnaseq.sam 2>&1 | tee aln_global.log

Notes:

• -x  vibCho/vibCho.O395.fa – prefix path of index files
• -U fq/cholera_rnaseq.fastq.gz – FASTQ file for single-end (Unpaired) alignment
• -S cholera_rnaseq.sam – tells bowtie2 to report alignments in SAM format to the specified file
• 2>&1 redirects standard error to standard output
  • while the alignment data is being written to the cholera_rnaseq.sam file, bowtie2 will report its progress to standard error.
• | tee aln.log takes the bowtie2 progress output and pipes it to the tee program
  • tee takes its standard input and writes it to the specified file and also to standard output
  • that way, you can see the progress output now, but also save it to review later (or supply to MultiQC)

Since the FASTQ file is not large, this should not take too long, and you will see progress output like this:

89006 reads; of these:
  89006 (100.00%) were unpaired; of these:
    20675 (23.23%) aligned 0 times
    38226 (42.95%) aligned exactly 1 time
    30105 (33.82%) aligned >1 times
76.77% overall alignment rate

When the job is complete you should have a cholera_rnaseq.sam file that you can examine using whatever commands you like.  Remember, to further process it downstream, you should create a sorted, indexed BAM file from this SAM output.

Exercise: Repeat the alignment performing a local alignment, and write the output in BAM format. How do the alignment statistics compare?

module load samtools
cd $SCRATCH/core_ngs/alignment/vibCho
bowtie2 --local -x vibCho/vibCho.O395 -U fq/cholera_rnaseq.fastq.gz 2>aln_local.log | \
  samtools view -b > cholera_rnaseq.local.bam

89006 reads; of these:
  89006 (100.00%) were unpaired; of these:
    27061 (30.40%) aligned 0 times
    33828 (38.01%) aligned exactly 1 time
    28117 (31.59%) aligned >1 times
69.60% overall alignment rate

Exercise #5: BWA-MEM - Human mRNA-seq

After bowtie2 came out with a local alignment option, it wasn't long before bwa developed its own local alignment algorithm called BWA-MEM (for Maximal Exact Matches), implemented by the bwa mem command. bwa mem has the following advantages:

• It incorporates a lot of the simplicity of using bwa with the complexities of local alignment, enabling straightforward alignment of datasets like the mirbase data we just examined
• It can align different portions of a read to different locations on the genome
  • In a total RNA-seq experiment, reads will (at some frequency) span a splice junction themselves
    • or a pair of reads in a paired-end library will fall on either side of a splice junction.
  • We want to be able to align these splice-adjacent reads for many reasons, from accurate transcript quantification to novel fusion transcript discovery.

Thus, our last exercise will be the alignment of a human total RNA-seq dataset composed (by design) almost exclusively of reads that cross splice junctions.

bwa mem was is made available by loading the bwa BioContainers module

we loaded the bwa module, so take a look at its usage information. The most important parameters, similar to those we've manipulated in the past two sections, are the following:

There are many more parameters to control the scoring scheme and other details, but these are the most essential ones to use to get anything of value at all.

The test file we will be working with is just the R1 file from a paired-end total RNA-seq experiment, meaning it is (for our purposes) single-end. Go ahead and take a look at it, and find out how many reads are in the file.

cd $SCRATCH/core_ngs/alignment
ls fastq
gunzip -c fastq/human_rnaseq.fastq.gz | echo $((`wc -l`/4))

A word about real splice-aware aligners

Using BWA mem for RNA-seq alignment is sort of a "poor man's" RNA-seq alignment method. Real splice-aware aligners like tophat2 or star have more complex algorithms (as shown below) – and take a lot more time!

RNA-seq alignment with bwa aln

Now, try aligning it with bwa aln like we did in Example #1, but first link to the hg19 bwa index directory.  In this case, due to the size of the hg19 index, we are linking to Anna's scratch area INSTEAD of our own work area containing indexes that we built ourselves.

cd $SCRATCH/core_ngs/alignment
ln -s -f /scratch/01063/abattenh/ref_genome/bwa/bwtsw/hg19
ls hg19

You should see a set of files analogous to the yeast files we created earlier, except that their universal prefix is hg19.fa.  

Go ahead and try to do a single-end alignment of the file to the human genome using bwa aln like we did in Exercise #1, saving intermediate files with the prefix human_rnaseq_bwa. Go ahead and just execute on the command line.

bwa aln hg19/hg19.fa fastq/human_rnaseq.fastq.gz > human_rnaseq_bwa.sai
bwa samse hg19/hg19.fa human_rnaseq_bwa.sai fastq/human_rnaseq.fastq.gz > human_rnaseq_bwa.sam

Once this is complete use less to take a look at the contents of the SAM file, using the space bar to leaf through them. You'll notice a lot of alignments look basically like this:

HWI-ST1097:228:C21WMACXX:8:1316:10989:88190     4       *       0       0       *       *       0       0
  AAATTGCTTCCTGTCCTCATCCTTCCTGTCAGCCATCTTCCTTCGTTTGATCTCAGGGAAGTTCAGGTCTTCCAGCCGCTCTTTGCCACTGATCTCCAGCT
  CCCFFFFFHHHHHIJJJJIJJJJIJJJJHJJJJJJJJJJJJJJIIIJJJIGHHIJIJIJIJHBHIJJIIHIEGHIIHGFFDDEEEDDCDDD@CDEDDDCDD

Notice that the contig name (field 3) is just an asterisk ( * ) and the alignment flags value is a 4 (field 2), meaning the read did not align (decimal 4 = hex 0x4 = read did not map).

Essentially, nothing (with a few exceptions) aligned. Why?

RNA-seq alignment with bwa mem

Exercise: use bwa mem to align the same data

Based on the following syntax and the above reference path, use bwa mem to align the same file, saving output files with the prefix human_rnaseq_mem. Go ahead and just execute on the command line.

bwa mem <ref.fa> <reads.fq> > outfile.sam

bwa mem hg19/hg19.fa fastq/human_rnaseq.fastq.gz > human_rnaseq_mem.sam

Check the length of the SAM file you generated with wc -l. Since there is one alignment per line, there must be 586266 alignments (minus no more than 100 header lines), which is more than the number of sequences in the FASTQ file. This is bwa mem can report multiple alignment records for the same read, hopefully on either side of a splice junction. These alignments can still be tied together because they have the same read ID.

# This gives you a view where each read is listed next to the number of entries it has in the SAM file
cut -f 1 human_rnaseq_mem.sam | sort | uniq -c | less

#This gives essentially a histogram of the number of times each read aligned - a plurality of reads aligned twice, which seems reasonable since these are all reads crossing a junction, but plenty aligned more or less
cut -f 1 human_rnaseq_mem.sam | sort | uniq -c | awk '{print $1}' | sort | uniq -c | less	

#This gives a better idea of the alignment rate, which is how many reads aligned at least once to the genome.  Divided by the number of reads in the original file, the real alignment rate is much higher.
cut -f 1 human_rnaseq_mem.sam | sort | uniq | wc -l	

# NOTE: Some of these one-liners are only reasonably fast if the files are relatively small (around a million reads or less). For bigger files, there are better ways to get this information, mostly using samtools.

BWA-MEM vs Tophat

Another approach to aligning long RNA-seq data is to use an aligner that is more explicitly concerned with sensitivity to splice sites, namely a program like Tophat. Tophat uses either bowtie (tophat) or bowtie2 (tophat2) as the actual aligner, but performs the following steps:

• aligns reads to the genome
• reads that do not align to the genome are aligned against a transcriptome, if provided
  • if they align, the transcriptome coordinates are converted back to genomic coordinates, with gaps represented in the CIGAR string, for example as 196N
• reads that do not align to the transcriptome are split into smaller pieces, each of which Tophat attempts to map to the genome

Note that Tophat also reports secondary alignments, but they have a different meaning. Tophat always reports spliced alignments as one alignment records with the N CIGAR string operator indicating the gaps. Secondary alignments for Tophat (marked with the 0x100 BAM flag) represent alternate places in the genome where a  read (spliced or not) may have mapped.

As you can imagine from this series of steps, Tophat is very computationally intensive and takes much longer than bwa mem – very large alignments (hundreds of millions of reads) may not complete in stampede's 48 hour maximum job time!