Alignment with arbitrary number of mismatches or gaps

Question

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?

Answer 1

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 2

I've ended up using a package called BatMis 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 3

Did you see

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