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.
Alignment of this prokaryotic data follows the workflow below. Here we will concentrate on steps 1 and 2.
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.
|
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.
First make sure you're in an idev session:
idev -m 180 -p normal -A UT-2015-05-18 -N 1 -n 68 # Answer "y" to the question about using the BIO_DATA_week_1 reservation # But if that hangs ("PD") waiting for a node: # - Ctrl-C to cancel the idev request # - execute the idev command above again # - but answer "n" to the reservation question |
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!
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:
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.
Now we will go back to our scratch area to do the alignment, and set up symbolic links to the index in the work area to simplify the alignment command:
cd $SCRATCH/core_ngs/alignment ln -s -f $WORK/core_ngs/references/bt2/vibCho vibCho |
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.
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 our command would look like this:
bowtie2 -x vibCho/vibCho.O395 -U fastq/cholera_rnaseq.fastq.gz -S cholera_rnaseq.sam |
Parameters are:
|
Create a commands file called bt2_vibCho.cmds with this task definition then generate and submit a batch job for it (time 1 hour, development queue).
Use nano to create the bt2_vibCho.cmds file. Then:
|
When the job is complete you should have a cholera_rnaseq.sam file that you can examine using whatever commands you like.
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.
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:
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):
Before the alignment, of course, we've got to build a mirbase 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:
bowtie2-build <reference_in> <bt2_index_base> |
Following what we did earlier for BWA indexing, namely move our FASTA into place, create the index directory, and establish our symbolic links.
mkdir -p $WORK/core_ngs/references/bt2/mirbase.v20 mv $WORK/core_ngs/references/hairpin_cDNA_hsa.fa $WORK/core_ngs/references/fasta cd $WORK/core_ngs/references/bt2/mirbase.v20 ln -s -f ../../fasta/hairpin_cDNA_hsa.fa ls -la module load perl module load bowtie/2.2.0 |
Now build the mirbase index with bowtie2-build like we did for the V. cholerae index:
bowtie2-build hairpin_cDNA_hsa.fa hairpin_cDNA_hsa.fa |
That was very fast! It's because the mirbase reference genome is so small compared to what programs like this are used to dealing with, which is the human genome (or bigger). You should see the following files:
hairpin_cDNA_hsa.fa hairpin_cDNA_hsa.fa.1.bt2 hairpin_cDNA_hsa.fa.2.bt2 hairpin_cDNA_hsa.fa.3.bt2 hairpin_cDNA_hsa.fa.4.bt2 hairpin_cDNA_hsa.fa.rev.1.bt2 hairpin_cDNA_hsa.fa.rev.2.bt2 |
Now, we're ready to actually try to do the alignment. Remember, unlike BWA, we actually need to set some options depending on what we're after. Some of the important options for bowtie2 are:
Option | Effect |
---|---|
--end-to-end or --local | Controls whether the entire read must align to the reference, or whether soft-clipping the ends is allowed to find internal alignments. Default --end-to-end |
-L | Controls the length of seed substrings generated from each read (default = 22) |
-N | Controls the number of mismatches allowable in the seed of each alignment (default = 0) |
-i | Interval between extracted seeds. Default is a function of read length and alignment mode. |
--score-min | Minimum alignment score for reporting alignments. Default is a function of read length and alignment mode. |
To decide how we want to go about doing our alignment, check out the file we're aligning with less:
cd $SCRATCH/core_ngs/alignment less fastq/human_mirnaseq.fastq.gz |
Lots of reads have long strings of A's, which must be an adapter or protocol artifact. Even though we see how we might be able to fix it using some tools we've talked about, what if we had no idea what the adapter sequence was, or couldn't use cutadapt or other programs to prepare the reads?
In that case, we need a local alignment where the seed length smaller than the expected insert size. Here, we are interested in finding any sections of any reads that align well to a microRNA, which are between 16 and 24 bases long, with most 20-22. So an acceptable alignment should have at least 16 matching bases, but could have more.
If we're also interested in detecting miRNA SNPs, we might want to allow a mismatch in the seed. So, a good set of options might look something like this:
-N 1 -L 16 --local |
Because these are short reads we do not have to adjust parameters like inter-seed distance (-i) or minimum alignment score (--min-score) that are a function of read length. If we were processing longer reads, we might need to use parameters like this, to force bowtie2 to "pretend" the read is a short, constant length: -i C,1,0 Yes, that looks complicated, and it kind of is. It's basically saying "slide the seed down the read one base at a time", and "report alignments as long as they have a minimum alignment score of 32 (16 matching bases x 2 points per match, minimum). See the bowtie2 manual (after you have had a good stiff drink) for a full explanation. |
Let's make a link to the mirbase index directory to make our command line simpler:
cd $SCRATCH/core_ngs/alignment ln -s -f $WORK/core_ngs/references/bt2/mirbase.v20 mb20 |
Putting this all together we have a command line that looks like this.
bowtie2 --local -N 1 -L 16 -x mb20/hairpin_cDNA_hsa.fa -U fastq/human_mirnaseq.fastq.gz -S human_mirnaseq.sam |
Parameters are:
|
Create a commands file called bt2.cmds with this task definition then generate and submit a batch job for it (time 1 hour, development queue).
Use nano to create the bt2.cmds file. Then:
|
When the job is complete you should have a human_mirnaseq.sam file that you can examine using whatever commands you like. An example alignment looks like this.
TUPAC_0037_FC62EE7AAXX:2:1:2607:1430#0/1 0 hsa-mir-302b 50 22 3S20M13S * 0 0 TACGTGCTTCCATGTTTTANTAGAAAAAAAAAAAAG ZZFQV]Z[\IacaWc]RZIBVGSHL_b[XQQcXQcc AS:i:37 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:16G3 YT:Z:UU |
Notes:
Such is the nature of bowtie2 – it it can be a powerful tool to sift out the alignments you want from a messy dataset with limited information, but doing so requires careful tuning of the parameters, which can take quite a few trials to figure out.
Exercise: About how many records in human_mirnaseq.sam represent aligned reads?
We can use our cut / grep trick from Exercise #1, but on the human_mirnaseq.sam file. Since all read names in this file start with TUPAC, we'll use that pattern to select non-header lines.
This expressions returns 221086. |
The mapping quality score is in field 5 of the human_mirnaseq.sam file. We can do this to pull out only that field:
grep -P -v '^@' human_mirnaseq.sam | cut -f 5 | head |
We will use the uniq create a histogram of these values. The first part of the --help for uniq says:
Usage: uniq [OPTION]... [INPUT [OUTPUT]] Filter adjacent matching lines from INPUT (or standard input), writing to OUTPUT (or standard output). With no options, matching lines are merged to the first occurrence. Mandatory arguments to long options are mandatory for short options too. -c, --count prefix lines by the number of occurrences |
To create a histogram, we want to organize all equal mapping quality score lines into an adjacent block, then use uniq -c option to count them. The sort -n command does the sorting into blocks (-n means numerical sort). So putting it all together, and piping the output to a pager just in case, we get:
grep -P -v '^@' human_mirnaseq.sam | cut -f 5 | sort -n | uniq -c | more |
Exercise: What is the flaw in this "program"?
We are looking at mapping quality values for both aligned and un-aligned records, but mapping quality only makes sense for aligned reads. This expression does not distinguish between mapping quality = 0 because the read mapped to multiple locations, and mapping quality = 0 because the sequence did not align. The proper solution will await the use of samtools to filter out unmapped reads. |
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:
Thus, our last exercise will be the alignment of a human long RNA-seq dataset composed (by design) almost exclusively of reads that cross splice junctions.
bwa mem was made available when 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:
Option | Effect |
---|---|
-k | Controls the minimum seed length (default = 19) |
-w | Controls the "gap bandwidth", or the length of a maximum gap. This is particularly relevant for MEM, since it can determine whether a read is split into two separate alignments or is reported as one long alignment with a long gap in the middle (default = 100) |
-r | Controls how long an alignment must be relative to its seed before it is re-seeded to try to find a best-fit local match (default = 1.5, e.g. the value of -k multiplied by 1.5) |
-c | Controls how many matches a MEM must have in the genome before it is discarded (default = 10000) |
-t | Controls the number of threads to use |
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.
|
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!
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?
Because this file was generated exclusively from reads in a larger dataset that cross at least one splice junction. The sequences as they exists in most of the reads do not correspond to a single location in the genome. However subsections of each read do exist somewhere in the genome. So, we need an aligner that is capable aligning different parts of the read to different genomic loci. |
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 |
|
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.
To get an idea of how often each read aligned, and what the 'real' alignment rate is, use the following commands:
This alignment rate is pretty good, but it could get better by playing around with the finer details of bwa mem. |
Be aware that some downstream tools (for example the Picard suite, often used before SNP calling) do not like it when a read name appears more than once in the SAM file. To mark the extra alignment records as secondary, specify the bwa mem -M option. This option leaves the best (longest) alignment for a read as -is but marks additional alignments for the read as secondary (the 0x100 BAM flag). This designation also allows you to easily filter the secondary reads with samtools if desired. |
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:
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!