mapping heteryzygous kmers on a genome I have a set (a couple of millions) of kmer pairs of known length (usually 21) that I know that are heterozygous in the middle nucleotide. For instance:

AAAAAAAAAANAAAAAAAAAA


Where N has always two states. I would like to map these nucleotides on a fragmented genome assembly that is presumably haploid. I bet that this problem is resolved, but I have hard times to make a practical solution, preferably within python environment. I considered following


Using an existing mapper


I tried to use python binding of minimap2, thinking that I just need to tell it that I search for exact matches only, but I did not manage to set it up. I set up smaller minimisers but it always returned no hits for kmers of length 21 (I tried to tweak parameters k, w and min_cnt, but to be honest I don't have a good understanding of how minimizers work). Is there a mapper that would be easily set up to do these searches?


Hash table


The second thing I tried was to build a hash table of all k - 1 / 2 mers of the genome to map the left and right sides of the kmers. Where sequences were mapped to lists of corresponding genomic positions, but this have blown my memory very soon (~50m of kmers were indexed using ~8GB of memory; i.e. this approach does not seems to be viable on a genomic scale).

I am also suspecting that there will be an elegant solution using a suffix array, but I have not cracked that one yet.
I don't know if this is the best answer, but one quick-n-dirty approach I commonly use involves tweaking bwa mem to report only perfect matches. So in your case, the idea is to make two copies of each 21-mer, replace the N with the appropriate allele, and then map the resulting sequences to your bwa-indexed assembly like so.

bwa mem -k 21 -T 21 -a -c 5000 assembly.fasta your-21mers.fasta


The -k 21 says to use a minimum seed length of 21, which forces exact matches. The -T 21 requires a minimum score of 21, which also enforces an exact match. The -a parameter reports all matches, since a "best" match doesn't make sense in this situation. The -c parameter limits how many matches are reported, which may need to be adjusted depending on how repetitive the 21-mer is.

UPDATE If you don't want to introduce a dependency like bwapy, it's fairly straightforward to implement this in pure Python.

with tempfile.TemporaryFile() as kmerfile:
    for kmer in kmers:
        print('&gt;{k}\n{k}'.format(k=kmer), file=kmerfile)
    kmerfile.flush()             # Often works without these two lines
    os.fsync(kmerfile.fileno())  # But best to be safe
    bwacmd = ['bwa', 'mem', '-k', ksize, '-t', ksize, '-a', '-c', '5000', refrfile, kmerfile.name]
    bwaproc = subprocess.Popen(bwacmd, stdout=subprocess.PIPE, universal_newlines=True)
    for line in bwaproc.stdout:
        # process SAM records line by line


Also, one can call bwa mem via a C API, but the mechanics of this are left as an exercise for the reader.