Removing multi-matched sequences in Bowtie

[gwilson]

Hello,

I'm trying to remove multi-matched sequences from my paired-end solexa dataset and could do with a spot of advice.
I use Bowtie to align my reads and thought that the -m option, if set to -m 1, would remove the multi-matched reads.

I tested this on a subset of 1000 reads, the stderr can be seen below:

# reads processed: 1000
# reads with at least one reported alignment: 561 (56.10%)
# reads that failed to align: 260 (26.00%) # reads with alignments suppressed due to -m: 179 (17.90%) Reported 561 paired-end alignments to 1 output stream(s)

This suggests to me that 260 failed to meet the alignment criteria, whilst 179 met the criteria but were matched more than once and so failed to be reported.

If I was to run the same analysis with the same dataset but remove the -m 1 parameter, I would expect to have 561+179= 740 reported alignments.

However, in fact I get 827:

# reads processed: 1000
# reads with at least one reported alignment: 827 (82.70%)
# reads that failed to align: 173 (17.30%) Reported 827 paired-end alignments to 1 output stream(s)

So my question is - what else is the -m flag affecting that could reduce the reported alignments.

Many Thanks!

[Xi Wang]

Could you provide the command line to show us other parameters when running bowtie?

[gwilson]

Yes, sorry:

bowtie -solexa1.3quals ref_genome -p 4 -S -m 1 -1 read1.fastq -2 read2.fastq

and

bowtie -solexa1.3quals ref_genome -p 4 -S -1 read1.fastq -2 read2.fastq


Thanks

[Xi Wang]

Your results also confused me... you may check the output to tell us why...
for suppressed multiple aligned reads you can save them by specifying the parameter --max, and then do a comparison.

[gwilson]

I printed out the unaligned reads (260) when running with -m 1. I have then tried to align these unaligned reads without the -m tag. 173 failed to align (the same number as previously), however the remaining 87 did align. Some examples are shown below:

HWUSI-EAS492:1:1:7:1425#ATCACG/1 99 7 58333990 255 36M = 58334090 136 ACCGAATCCAAGAACACATTAAAGCAATCATCCATC a`
]Z]_^^__\OW]V]XX\[]ZSJ[ZZY[ZW\SSYU XA:i:0 MD:Z:36 NM:i:0
HWUSI-EAS492:1:1:7:1425#ATCACG/2 147 7 58334090 255 36M = 58333990 -136 AATACATTATATAAACAAACTCAAAGACAAAAACCA ^a
UHP]]ab`aaaba^ab``_[abba`[abba`_aa XA:i:0 MD:Z:36 NM:i:0
HWUSI-EAS492:1:1:7:1818#AACACG/2 163 1 153720395 255 36M = 153720562 203 TGTGTGTGTGTGTGTGGTATGTGTGTATGTGTGTGT ab
abb``baba````_T\L]W``_RVF[`]XWQ_QV XA:i:2 MD:Z:18G7G9 NM:i:2
HWUSI-EAS492:1:1:7:1818#AACACG/1 83 1 153720562 255 36M = 153720395 -203 GGTGTGTGGTGTGTGGTGTGTGGTGTGTGGTGTGTG `]
TW]XGZa\^``\b^^a`[O^aa^aa_a`^aa^`` XA:i:0 MD:Z:36 NM:i:0
HWUSI-EAS492:1:1:7:1783#ATCACG/1 99 9 88890976 255 36M = 88891047 107 AAAATCAGAAATGAAAAGGGAGACATAACAACAGAT aa
a_`]`X^\`__^__XMNSX^Y^_Y]^[[\_XK^] XA:i:1 MD:Z:4A31 NM:i:1
HWUSI-EAS492:1:1:7:1783#ATCACG/2 147 9 88891047 255 36M = 88890976 -107 CAAAAGGGTATACTCAACAAAACTGGAAAACCTGGA B_
a`_NF[[^\]QTO^_Z\aa]ZXUT_a`Z]a\_aa XA:i:0 MD:Z:36 NM:i:0


To me, this output looks fine. So the question remains, why the discrepancy? Does the -m tag affect something in the alignment that I'm not aware of?

Incidentally, I also tried to run the original 1000 reads using the -y tag. This had no great impact.

[Xi Wang]

It seems that some repeat sequences such as "TGTGTGTGTGTGTGTGGTATGTGTGTATGTGTGTGT" cause this problem. As the criteria for align PE sequences may be different from those for SE sequences.
What would the results be if you map only half of the PE sequences?

[lh3]

Another experiment that may help to understand the result is to try a more accurate aligner, for example, novoalign.

[gwilson]

I've put the test data read 1 and read 2 through Bowtie seperately, treating it as single-end. This time the output is as you would expect i.e. the same number fail to align regardless of -m. Hence it only seems to be an issue affecting paired end data.

The single end output is below (-m 1):
bowtie -solexa1.3quals ref_genome -p 4 -m 1 test2
# reads processed: 1000
# reads with at least one reported alignment: 526 (52.60%)
# reads that failed to align: 72 (7.20%)
# reads with alignments suppressed due to -m: 402 (40.20%)
Reported 526 alignments to 1 output stream(s)

and default -m

bowtie -solexa1.3quals ref_genome -p 4 test2
# reads processed: 1000
# reads with at least one reported alignment: 928 (92.80%)
# reads that failed to align: 72 (7.20%)
Reported 928 alignments to 1 output stream(s)


Has anyone else seen this issue with paired end data?

[Xi Wang]

I think we should make the definition of "meads with alignments suppressed due to -m" and the definition of "reads that failed to align" clear.