Bowtie2 is an ultra-fast program for aligning next generation sequence reads to large genomes written by Ben Langmead and which happens to be my aligner of choice. RTFM here. Download it here or learn how to get it on your Mac OS X with "homebrew" from one of my previous posts, "The Craft of Homebrew-ing for Mac OS". Learn about how it assigns MAPQ scores and related ideas in my previous posts, "How does bowtie2 assign MAPQ scores?" and "The slow death of the term 'uniquely mappable'...". In this post, through mini experiments, I discuss whether or not bowtie2 randomly chooses a site for reads with multiple equal alignment positions. I chose to do this to support or put to rest various concerns raised in forums (such as seqanswers and biostars) that it always puts reads at the same spot or with a heavy bias for 1 or a couple spots.
Bowtie2 searches for as many alignments as it can (with given time and memory constraints) for a given read against the given reference (e.g. human genome) that have alignment scores above the minimum score cutoff. The output contains information on each read's best alignment score (AS) and, if a second alignment existed with a score above the cutoff, it's second best alignment score (XS). If an XS does not exist, this tautologically means there was no second alignment with a score above the minimum cutoff. Though this is a seemingly vacuous interpretation, it is the most robust one being a tautology and all. Some take "no XS" to mean its AS alignment was the only alignment or that it is a "uniread" while simultaneously not considering other reads with AS > XS as "unireads". Instead, this latter group of reads are lumped into the 'multiread' category, which are reads that map to more than one site. Others acknowledge that having no XS is just another way of saying AS > XS and that all reads with AS > XS are "unireads". This latter view acknowledges that having an XS is in part dependent on that minimum score cutoff as well as penalty sizes, whether or not mismatches and/or gaps are allowed, and whether or not you use end-to-end or local mode. Multireads, in this view, are only reads where AS = XS. To be clear, I try to refer to reads with AS = XS as "true multireads" since there is more than one interpretation of "multiread". However, when I do say 'multiread' I do still mean reads with AS = XS. Moreover, since the terms 'uniquely mappable' and 'uniread' has more than one meaning now, I try to refer to reads with AS > XS (including reads with no XS) as "maxireads" since they have a single maximum AS that is greater than all other XS. For practical reasons, some times I refer to reads with no XS as "true unireads", but this is a misnomer and I should stop it. I do it only out of the need to refer to that group of reads with no XS since some people are interested in only them. Perhaps I will just call them "classic unireads". Find out more about the fickleness of the state of being unique in my previous post, "The slow death of the term 'uniquely mappable'…".
This post is interested in "true multireads" where AS = XS. If you are keeping true multireads in your data set, it could be an issue if they all get assigned to the same spot. Of course, this effect would be mitigated if you removed all "redundant reads" (reads that map to the same exact position) except one before moving on to downstream analyses. There also seemed to be some debate out there as to what bowtie2 does with multireads. So I designed a couple very simple experiments to show what it does.
The take away: these experiments support exactly what the bowtie2 manual says about what it does with randomly assigning places for multireads. If you do not know what the manual says than RTFM or read on.
Randomly generated (multi)reads used:
"Chromosomes" were made from the reads:
There is a chromosome with the same sequence as every read (e.g. chr2 for read2) made as tail-to-head repeats of the read in the copy number associated with the read (e.g. chr2 has 2 copies of read 2). This is simply done in python taking advantage of its string multiplication (e.g. string*2 = stringstring). Therefore, every read can map to its respective chromosome >= 2 times with AS = XS and they are all true multireads. The chromosomes were indexed with "bowtie2-build -f chr.fa chr".
Using many copies of the same multireads:
One thousand copies of each mutiread was placed in the fasta file, reads.fa, to study the behavior of mapping the same multiread over multiple trials.
Experiments, results, and conclusions:
First I executed the following bowtie2 command and pipes:
bowtie2 --very-sensitive -N 1 -x chr -U multireads.fa -f | grep -v ^@ | cut -f 1,3,4 | sort | uniq -c
In bowtie2, the output is SAM format and fields 1, 3, and 4 are read name, chr name, and 1-based start position respectively. Thus, this pipeline ends in counting the number of times "read chr pos" were identical. Here are the results:
read_name chr_name pos count
read2 chr2 1 1000
read3 chr3 101 1000
read4 chr4 151 1000
read5 chr5 51 1000
read10 chr10 251 1000
read100 chr100 2901 1000
All one thousand copies of each multiread mapped to the same place. If you RTFM of bowtie2, it explains this behavior. Bowtie2 uses a pseudo-random number generator (for those not familiar with computer science, this is how randomness is simulated in general) to randomize where a read will go given more than 1 choice. So, why do they all go to the same place? The answer is because bowtie2 uses the read name, nucleotide string, and quality string to determine the seed to use for the pseudo-random number generator. Briefly, if you use the same seed every time, then you get the same string of pseudo-random numbers every time. Here all of my multireads have the same name, the same sequence, and since they are fasta files, the same string of quality scores (bowtie2 assumes they are all 40). Thus, this is the expected behavior. In real life, you will not encounter (or at least not often) a situation where you have multiple reads with the same name. Reads can have the same sequence, but it is less likely they will share the same exact quality string. Nonetheless, a different seed would be generated given 2 identical sequences because the read names (and probably quality strings) would differ.
However, for now I am dealing with totally identical reads. If you want fully simulated randomness with identical reads, you need to use the --non-deterministic flag. This turns off the deterministic behavior of picking a seed from the fixed read parameters and instead just picks a seed using the current time (each time it encounters one of the identical reads, it is a different time, so a different seed is used). So, I tested if this randomly distributes the reads with:
bowtie2 --non-deterministic --very-sensitive -N 1 -x chr -U multireads.fa -f | grep -v ^@ | cut -f 1,3,4 | sort | uniq -c
read_name chr_name pos count expectedCountIfuniform
read2 chr2 1 496 500
read2 chr2 51 504 500
read3 chr3 1 338 333.333
read3 chr3 101 326 333.333
read3 chr3 51 336 333.333
read4 chr4 1 249 250
read4 chr4 101 235 250
read4 chr4 151 254 250
read4 chr4 51 262 250
read5 chr5 1 168 200
read5 chr5 101 205 200
read5 chr5 151 201 200
read5 chr5 201 223 200
read5 chr5 51 203 200
It really does not take statistical tests to show that bowtie2 is doing a good job uniformly distributing the reads across all possible sites for a read. The table for all of the read100 datapoints would be too large to show, so instead, here is a plot (made in R):
The only thing left to do is to change the read names to see if that is sufficient for randomized results as well. To do this, I just added a number to the read name such as read2_1. The number added was just the line number as it occurred in the file. For these, since the names were all different and since each read can only map to its respective chromosome, I only counted how many times a read aligned to "chr, pos" with:
bowtie2 --very-sensitive -N 1 -x chr -U diffNameMultiReads.fa -f | grep -v ^@ | cut -f 3,4 | sort | uniq -c
chr_name pos count expectedCountIfuniform
chr2 1 503 500
chr2 51 497 500
chr3 1 362 333.3333
chr3 101 316 333.3333
chr3 51 322 333.3333
chr4 1 277 250
chr4 101 246 250
chr4 151 251 250
chr4 51 226 250
chr5 1 200 200
chr5 101 192 200
chr5 151 216 200
chr5 201 192 200
chr5 51 200 200
Again, it is clear from just the table alone that bowtie2 is doing well at uniformly distributing the multireads as it claims to do.
Bowtie2 definitely assigns 'true multireads' uniformly at random.
A long time ago, I took a very high multiplicity (~12,000) multiread from one of our experiments and tried mapping it with bowtie2 again.. and again… and again. It mapped to the same place every time rather than one of the other 11,999 places. Back then, I was momentarily horrified until I RTFM. This behavior is the expected behavior in deterministic mode, where bowtie2 uses the read name, sequence, and quality string to seed the pseudo-random number generator. The solution would be to use the --non-deterministic flag. The real experiment would have been to change the read's name each time. In our real NGS experiments, none of the reads share a name. So, given the results above, it is safe to say that if 2 true multireads with different names but the same sequence are encountered, they will likely be assigned to different sites uniformly at random. If reads with the same sequence are filtered beforehand to keep only one or if there is just naturally a single instance of a multiread, you can trust that the site it is assigned to is picked at random. Even if it was not, most people filter multireads out. To filter all true multireads out (and some maxireads including classic unireads with poor AS), use MAPQ >= 2. See my previous post for more information on how bowtie2 assigns MAPQ scores. Even if someone does not filter out the multireads, most people for many applications remove all redundant reads except 1 per site (redundant reads are those that map to the same exact site). Thus, even if bowtie2 was not random, the pileup effect would be smashed at that step. Finally, since most people only care about reads of AS > XS for most applications, random assignment does not matter. Those go to the site corresponding to AS.
As a closing remark, if you want to work with multireads, then the a priori uniform distribution over all sites is not the best you can do. One could create a posterior distribution over those sites after looking at the distribution of the high quality maxiread (AS > XS) data. I discuss this more here.