Page History

...

• the yeast genome (sacCer3)
• the human genome (hg19)
• a Vibrio cholerae genome (0395; our name: vibCho)
• the microRNA database mirbase (v20), human subset

NOTE: For the sake of simplicity, these are not necessarily the most recent versions of these references - for example, hg19 is the second most recent human genome, with the most recent called hg38.  Similarly, the most recent mirbase annotation is v21.

Mirbase is a collection of all known microRNAs in all species (and many speculative miRNAs). We will use the human subset of that database as our alignment reference.  This has the advantage of being significantly smaller than the human genome, while likely containing almost all sequences likely to be detected in a miRNA sequencing run.

These are the four reference genomes we will be using today, with some information about them (and here is information about many more genomes):

These are the four reference genomes we will be using today, with some information about them (and here is information about many more genomes):

...

/scratch/01063/abattenh/ref_genome/bwa/bwtsw/hg19

However, we We can index the references for the yeast genome, the human miRNAs, and the V. cholerae genome, because they are all tiny compared to the human genome.  We We will grab the FASTA files for yeast and human miRNAs two references and build each index right before we use them.  We We will also grab the special file that contains the V. cholerae genome sequence and annotations (a ".gbk" file), and generate the reference FASTA and some other interesting information when we get to that exercise. These references are currently at the following locations:

...

First stage all the reference files in your work "$WORK core_ngs" area in a directory called references.  We We will add further structure to this directory later on in specific exercises, but for now the following will suffice:

...

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!

Exercise #2: Bowtie2 and Local Alignment - Human microRNA-seq

Now we're going to switch over to a different aligner that was originally designed for very short reads and is frequently used for RNA-seq data. Accordingly, we have prepared another test microRNA-seq dataset for you to experiment with (not the same one you used cutadapt on). This data is derived from a human H1 embryonic stem cell (H1-hESC) small RNA dataset generated by the ENCODE Consortium – its about a half million reads.

However, there is a problem!  We don't know (or, well, you don't know) what the adapter structure or sequences were. So, you have a bunch of 36 base pair reads, but many of those reads will include extra sequence that can impede alignment – and we don't know where! We need an aligner that can find subsections of the read that do align, and discard (or "soft-clip") the rest away – an aligner with a local alignment mode. Bowtie2 is just such an aligner.

Overview miRNA alignment workflow with bowtie2

If the adapter structure were known, the normal workflow would be to first remove the adapter sequences with cutadapt. Since we can't do that, we will instead perform a local alignment of the single-end miRNA sequences using bowtie2. This workflow has the following steps:

1. Prepare the mirbase v20 reference index for bowtie2 (one time) using bowtie2-build 
2. Perform local alignment of the R1 reads with bowtie2, producing a SAM file directly
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)

This looks so much simpler than bwa – only one alignment step instead of three! We'll see the price for this "simplicity" in a moment...

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!

Exercise #4: Bowtie2 and Local Alignment - Human microRNA-seq

Now we're going to switch over to a different aligner that was originally designed for very short reads and is frequently used for RNA-seq data. Accordingly, we have prepared another test microRNA-seq dataset for you to experiment with (not the same one you used cutadapt on). This data is derived from a human H1 embryonic stem cell (H1-hESC) small RNA dataset generated by the ENCODE Consortium – its about a half million reads.

However, there is a problem!  We don't know (or, well, you don't know) what the adapter structure or sequences were. So, you have a bunch of 36 base pair reads, but many of those reads will include extra sequence that can impede alignment – and we don't know where! We need an aligner that can find subsections of the read that do align, and discard (or "soft-clip") the rest away – an aligner with a local alignment mode. Bowtie2 is just such an aligner.

Overview miRNA alignment workflow with bowtie2

If the adapter structure were known, the normal workflow would be to first remove the adapter sequences with cutadapt. Since we can't do that, we will instead perform a local alignment of the single-end miRNA sequences using bowtie2. This workflow has the following steps:

1. Prepare the mirbase v20 reference index for bowtie2 (one time) using bowtie2-build 
2. Perform local alignment of the R1 reads with bowtie2, producing a SAM file directly
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)

This looks so much simpler than bwa – only one alignment step instead of three! We'll see the price for this "simplicity" in a moment...

As before, we will just do the alignment steps leave samtools for the next section.

Mirbase is a collection of all known microRNAs in all species (and many speculative miRNAs). We will use the human subset of that database as our alignment reference.  This has the advantage of being significantly smaller than the human genome, while likely containing almost all sequences likely to be detected in a miRNA sequencing run.

These are the four reference genomes we will be using today, with some information about them (and here is information about many more genomes):As before, we will just do the alignment steps leave samtools for the next section.

Introducing bowtie2

Go ahead and load the bowtie2 module so we can examine some help pages and options. To do that, you must first load the perl module, and then the a specific version of bowtie2. 

...