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 (&lt;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?
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)]
    for read in reads
        errors = 0
        for rp in 1..23
            if (genomeKmer[rp] != read[rp])
                errors += 1
        if (errors &lt; 7)
            print gp