<p>I've ended up using a package called <a href="https://www.ncbi.nlm.nih.gov/pubmed/22689389" rel="nofollow noreferrer">BatMis</a> for this, which can detect alignments with an arbitrary number of mismatches (not sure about indels). It is quite performant and yields about the order of magnitude of alignments that I was expecting (several thousands for alignment of a 20bp read to the whole human genome with up to 6 mismatches).</p>

CC BY-SA 4.0
2019-03-06T13:20:16.503Z
2019-03-06T13:20:16.503Z
Flagon13
<p>Did you see</p>

<blockquote>
  <p>-k  Report up to  valid alignments per read or pair (default: 1).</p>
</blockquote>

CC BY-SA 4.0
2019-02-18T17:17:24.517Z
Did see that, but I thought --all (or -a in short) means it outputs all found alignments, ie supersedes -k? Using -k 5 instead of --all still only yields one result.
Flagon13
CC BY-SA 4.0
2019-02-18T17:16:09.523Z
2019-02-18T17:16:09.523Z
swbarnes2
<p>This is not something mappers will be good at, because you're going to hit most of the genome. Bowtie won't do it, and I don't think LAST will either. The majority of fast mapping programs rely on seed-based mapping where the seed is essentially free of errors.</p>
<p>Doing only mismatches makes it much easier to determine, but it's still too relaxed. Consider that the example sequence you have provided is 10 mismatches away from <code>GGGGGGGGGGGGGGGGAGAGAGG</code>, and also 10 mismatches away from <code>AAAAAAAAAAAAAAAGAGAGAGG</code> (with the last 3 bp matching exactly)</p>
<p>Depending on how many reads you have it may be better to work in reverse, iterating over all the kmers in the genome and working out whether or not they are fewer than 7 mismatches (plus indels) away from any reads:</p>
<pre><code>for gp in genomePositions
    genomeKmer = genome[gp:(gp+23-1)]
    for read in reads
        errors = 0
        for rp in 1..23
            if (genomeKmer[rp] != read[rp])
                errors += 1
        if (errors &lt; 7)
            print gp 
</code></pre>

CC BY-SA 4.0
2019-02-18T18:24:33.357Z
2022-10-21T18:33:21.030Z
2022-10-21T18:33:21.030Z
gringer
<p>I have 23bp long reads and want to find all possible alignments of them to the human genome (hg19, hg38) for an arbitrary number of mismatches (&lt;7), possibly also small indels. I've read in literature that people use bowtie2 for this, so I've tried the following:</p>
<pre><code>echo &quot;GGAGAACAAGTGGCAGATAGAGG&quot; | bowtie2 --all -f -U- -x hg19 -r --no-hd --quiet --score-min C,-36 -p 3 -L 4
</code></pre>
<p>I'm not very experienced with bowtie2, and I hoped that this would give me all alignments of the piped-in sequence with up to 6 mismatches (6*-6 = -36) to the hg19 genome which I previously indexed from fasta files. It however only gives me the correct full alignment:</p>
<pre><code>0   0   chr2    173327389   255 23M *   0   0   GGAGAACAAGTGGCAGATAGAGG IIIIIIIIIIIIIIIIIIIIIII AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:23 YT:Z:UU
</code></pre>
<p>I've read that bowtie2 by design only allows a maximum of one mismatch per seed. Even if I implement the smallest possible seed length and interval (-L 4 -i C,1 -N 1), I only get this one result which I think is highly unlikely, given that I'm aligning to the whole genome.</p>
<p>The last three bp will have to match exactly - I was going to enforce this manually later on, maybe there's a better option for doing this. Taking the indels into account is also a secondary aim.</p>
<p>I'm fine with ~10k alignments per query for &lt;=6 mismatches plus the last 3bp matching exactly. Haven't done the combinatorics but that sounds realistic to me - question is how to make bowtie do it, or a different tool.</p>
<p>What would be the right way of achieving this using bowtie, or a different tool?</p>

<phylogenetics><bowtie2>
Alignment with arbitrary number of mismatches or gaps
# Alignment with arbitrary number of mismatches or gaps

(Asked by: username_0 on 2019-02-18)

I have 23bp long reads and want to find all possible alignments of them to the human genome (hg19, hg38) for an arbitrary number of mismatches (<7), possibly also small indels. I've read in literature that people use bowtie2 for this, so I've tried the following:
echo "GGAGAACAAGTGGCAGATAGAGG" | bowtie2 --all -f -U- -x hg19 -r --no-hd --quiet --score-min C,-36 -p 3 -L 4

I'm not very experienced with bowtie2, and I hoped that this would give me all alignments of the piped-in sequence with up to 6 mismatches (6*-6 = -36) to the hg19 genome which I previously indexed from fasta files. It however only gives me the correct full alignment:
0   0   chr2    173327389   255 23M *   0   0   GGAGAACAAGTGGCAGATAGAGG IIIIIIIIIIIIIIIIIIIIIII AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:23 YT:Z:UU

I've read that bowtie2 by design only allows a maximum of one mismatch per seed. Even if I implement the smallest possible seed length and interval (-L 4 -i C,1 -N 1), I only get this one result which I think is highly unlikely, given that I'm aligning to the whole genome.
The last three bp will have to match exactly - I was going to enforce this manually later on, maybe there's a better option for doing this. Taking the indels into account is also a secondary aim.
I'm fine with ~10k alignments per query for <=6 mismatches plus the last 3bp matching exactly. Haven't done the combinatorics but that sounds realistic to me - question is how to make bowtie do it, or a different tool.
What would be the right way of achieving this using bowtie, or a different tool?

**Tags:** phylogenetics, bowtie2

---

## 3 answers

### Answer by username_2 on 2019-02-18. Score: 7 | Accepted answer

This is not something mappers will be good at, because you're going to hit most of the genome. Bowtie won't do it, and I don't think LAST will either. The majority of fast mapping programs rely on seed-based mapping where the seed is essentially free of errors.
Doing only mismatches makes it much easier to determine, but it's still too relaxed. Consider that the example sequence you have provided is 10 mismatches away from GGGGGGGGGGGGGGGGAGAGAGG, and also 10 mismatches away from AAAAAAAAAAAAAAAGAGAGAGG (with the last 3 bp matching exactly)
Depending on how many reads you have it may be better to work in reverse, iterating over all the kmers in the genome and working out whether or not they are fewer than 7 mismatches (plus indels) away from any reads:
for gp in genomePositions
    genomeKmer = genome[gp:(gp+23-1)]
    for read in reads
        errors = 0
        for rp in 1..23
            if (genomeKmer[rp] != read[rp])
                errors += 1
        if (errors < 7)
            print gp


---

### Answer by username_0 on 2019-03-06. Score: 2

I've ended up using a package called [BatMis](https://www.ncbi.nlm.nih.gov/pubmed/22689389) for this, which can detect alignments with an arbitrary number of mismatches (not sure about indels). It is quite performant and yields about the order of magnitude of alignments that I was expecting (several thousands for alignment of a 20bp read to the whole human genome with up to 6 mismatches).


---

### Answer by username_1 on 2019-02-18. Score: 0

Did you see

> -k Report up to valid alignments per read or pair (default: 1).

> **Comments:**
> **username_0**: Did see that, but I thought --all (or -a in short) means it outputs all found alignments, ie supersedes -k? Using -k 5 instead of --all still only yields one result.  

---

