Alignments and Mappers
Last updated
Last updated
The ‘read mapping’ problem is the alignment of relatively short reads to a reference genome.
Many short read mappers are currently available (Bowtie, BWA, SOAP, TopHat, etc...). They all try to overcome the following challenges presented by the mapping problem:
how do we align efficiently (in time and memory footprint) billions of short reads to a long reference genome?
how do we decide a specific placement when multiple, equally good alignment are available for the same read?
Some mappers are short-read aligners (they do not allow large gaps in their alignments), while others are spliced-read mappers (they allow large gaps corresponding to introns). The choice of the appropriate mapper depends on the particular dataset and application
This tutorial uses a subset of the Illumina IDEA Challenge 2011 dataset (more information can be found ). You can copy the relevant fastq files by typing the following command on Luria:
The dataset consists of paired-end RNA-Seq data of breast carcinoma human cell line BT20. Because the reads are paired-end reads, for each DNA fragment we have sequence data from both ends and the sequence reads are stored in two separate files (one for the data from each end). Additional information about the formats used in this tutorial can be found below:
SAM/BAM format: and .
.
A complete guide on how to run bowtie can be found .
Bowtie can be run in 2 alignment modes (mutually exclusive):
n alignment mode: this is the default mode (specified with flag -n), and determines which alignments are valid according to the following policy:
alignments may have no more than N mismatches in the first L bases of the reads ("seed"). N and L are integers specified through the flags -n and -l.
the sum of the Phred quality values at all mismatched positions (not just in the seed) may not exceed E (set with -e).
v alignment mode:
alignments may have no more than V mismatches (set using the -v option), and quality values are ignored.
Several reporting modes are available, specified through various combinations of the options, such as -k,-a, -m, etc...
When the command has completed, the current directory will contain the index consisting of four new files (ref.1.ebwt, ref.2.ebwt, ref.rev.1.ebwt, and ref.rev.2.ebwt).
How to run Bowtie
If your dataset consists of single-end sequence reads, you can run Bowtie as follows:
If your dataset consists of paired-end sequence reads, you can run Bowtie as follows:
Some additional useful options are:
the flag -k specifies the maximum number of alignments to report if more than one valid alignment is found.
the flag -m suppresses all alignments for a particular read if more than m reportable alignments exist.
the flag -v specifies the maximum number of mismatches allowed in the entire length of the read.
the flag -n specifies the maximum number of mismatches in the high quality seed of the read.
the flag -l specifies the length of the seed to consider for each read.
the flag –-best guarantees that reported singleton alignments are the best in terms of stratum (i.e. number of mismatches or mismatches in the seed, depending on the alignment mode) and in terms of the quality values at the mismatched position(s).
Note that the option --best does not affect which alignments are considered "valid" by bowtie, only which valid alignments are reported. When --best is specified and multiple hits are allowed (via -k or -a), the alignments for a given read are guaranteed to appear in best-to-worst order in bowtie's output. bowtie is somewhat slower when --best is specified.
The last few lines written to the standard error file handle convey summary information about the run (see below). The other lines are written to the standard out file handle and show valid, reportable alignments found by Bowtie.
BWA implements two different algorithms, both based on Burrows-Wheeler Transform (BWT).
ALN:
designed for short reads (up to ~200bp) with low error rate ( less than 3%).
performs gapped global alignment.
supports paired-end reads.
BWA-SW:
designed for longer reads with more errors.
performs heuristic Smith-Waterman-like alignment to find high-scoring local hits.
supports only single-end reads.
The database file (in fasta format) must be first indexed with the index command, which typically takes a few hours. Pre-built indexes are available on rous for the most common model organisms.
If your reference sequence is not among the available pre-built indexes, you can run the following command and create one for it:
For short reads (up to ~200bp) with low error rate (<3%), the algorithm to use is invoked with the command aln, which finds the suffix array (SA) coordinates of good hits of each individual read; the commands samse (for single-end reads) or sampe (for paired-end reads) are then used to convert these SA coordinates into chromosomal coordinates. Note that multi-hits are randomly chosen.
If your dataset consists of single-end sequence reads, you can run BWA as follows:
For example, to run BWA on the query BT20_1.fastq against the reference genome hs18, you must type the following commands:
If your dataset consists of paired-end sequence reads, you must run BWA in two steps as follows:
For example, to run BWA on the query BT20_1.fastq and BT20_2.fastq against the reference genome hs18, you must type the following commands:
Note that by default BWA ALN allows up to K mismatches in the seed (the first L positions of each read) and up to N mismatches in the whole sequence reads. These parameters can be changed using the appropriate options as indicates below.
-l specifies the length of the seed.
-k specifies the maximum edit distance in the seed.
-n specifies the maximum edit distance in the whole sequence read.
-o specifies the maximum number of gaps to open.
-e specifies the maximum number of gaps to extend.
The BWA output in SAM format has the following specific optional fields:
NM Edit distance
MD Mismatching positions/bases
X0 Number of best hits
X1 Number of suboptimal hits found by BWA
XM Number of mismatches in the alignment
XO Number of gap opens
XG Number of gap extensions
XT Type: Unique/Repeat/N/Mate-sw
XA Alternative hits; format: (chr,pos,CIGAR,NM;)*
Bowtie requires an index of the reference sequence against which the reads are mapped. Pre-built indexes for several organisms' assemblies are available on rous or on the Bowtie official (check under left pane section "Pre-built Indexes" ). If your reference sequence is not among the available pre-built indexes, you can run the following command and create one for it:
BWA is another fast aligner based on Burrows-Wheeler Transform (BWT). A complete guide on how to run BWA can be found .