# Alignment with arbitrary number of mismatches or gaps

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.

What would be the right way of achieving this using bowtie, or a different tool?

• With 10 mismatches on a 23bp sequence and up to 5bp indels, you'll hit almost the entire genome. I don't know of any mappers that are designed to work with such a high mismatch rate. Could you please add some more context to explain why you want to do this?
– gringer
Feb 18 '19 at 18:10
• 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've adjusted my original question. Feb 18 '19 at 18:12
• 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)
– gringer
Feb 18 '19 at 18:19
• 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. Feb 18 '19 at 18:22
• 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.
– gringer
Feb 18 '19 at 18:27

This is not something mappers will be good at, because you're going to hit most of the genome. 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)]
errors = 0
for rp in 1..23