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?