<p>Indeed, with suffix arrays, you can map a kmer in ~0.0007s. The way to go is to build a suffix array (for instance using <a href="https://bitbucket.org/alex-warwickvesztrocy/pysais/src" rel="nofollow noreferrer">PySAIS package)</a> and then do a convoluted binary search of the four possible kmers. A full implementation of such mapper can be found <a href="https://github.com/KamilSJaron/smudgeplot/blob/dev/smudgeplot/map.py" rel="nofollow noreferrer">here</a>. </p>

<p>To build suffix array</p>

<pre><code>from Bio import SeqIO
from Bio.Seq import Seq
from PySAIS import sais
from bisect import bisect_right

# load the genome
kmer_genome_file = 'my_genome.fasta'
ffile = SeqIO.parse(kmer_genome_file, "fasta")
sequences = []
scf_names = []
for seq_record in ffile:
    scf_names.append(seq_record.name)
    sequences.append(str(seq_record.seq).upper())

ffile.close()

# paste all scaffolds into a single string and build a suffix array
genome = '$'.join(sequences)
sa = sais(genome)
</code></pre>

<p>To transform genome coordinate to scaffold/position coordinates we define a function <code>evaluated2assembly_position</code></p>

<pre><code>def evaluated2assembly_position(eval_pos):
    scf_index = bisect_right(scf_sizes, eval_pos)
    eval_pos - scf_sizes[scf_index - 1]
    if scf_index &gt; 0:
        return((scf_names[scf_index], eval_pos - scf_sizes[scf_index - 1]))
    else:
        return((scf_names[scf_index], eval_pos))
</code></pre>

<p>and finally the function for a convoluted binary search (assuming all kmers being 21 nt long and having the unknown nucleotide in the middle, on 11th nt). Note that you need to map to both strands of DNA, which is equivalent of mapping kmer and it's reverse complementary sequence</p>

<pre><code>def searchKmer(self, kmer):
    reverse_complementary_kmer = str(Seq(kmer).reverse_complement())
    alignments = []
    for strand in ['+', '-']:
        if strand == '-':
            kmer = reverse_complementary_kmer
        L_kmer = kmer[0:10]
        l = 0
        r = len(sa) - 1
        while l &lt;= r:
            m = (l + r) // 2
            eval_pos = sa[m]
            L_genome_kmer = genome[eval_pos:(eval_pos + 10)]
            if L_kmer == L_genome_kmer:
                r_L_converged = r
                for N in ['A', 'C', 'G', 'T']:
                    kmer = kmer[0:10] + N + kmer[11:21]
                    while l &lt;= r:
                        m = (l + r) // 2
                        eval_pos = sa[m]
                        genome_kmer = genome[eval_pos:(eval_pos + 21)]
                        if kmer == genome_kmer:
                            m_hit = m
                            # following ms
                            while kmer == genome_kmer:
                                scf, pos = evaluated2assembly_position(eval_pos)
                                alignments.append([scf, pos, strand, N])
                                m += 1
                                eval_pos = sa[m]
                                genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            # previous ms
                            m = m_hit - 1
                            eval_pos = sa[m]
                            genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            while kmer == genome_kmer:
                                scf, pos = evaluated2assembly_position(eval_pos)
                                alignments.append([scf, pos, strand, N])
                                m -= 1
                                eval_pos = sa[m]
                                genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            break
                        elif genome_kmer &lt; kmer:
                            l = m + 1
                        else:
                            r = m - 1
                    r = r_L_converged # recover the r after converging on the left sub-kmer
                break
            elif L_genome_kmer &lt; L_kmer:
                l = m + 1
            else:
                r = m - 1
    return(alignments)
</code></pre>

<p>Taking one of the fist kmers in the scaffold one to test it</p>

<pre><code>searchKmer('ACGGCCACCCNGTGTCGTTGT')                                                                                                                                                                                                                                                    
Out[23]: [['scaf000001', 28, 'G']]
</code></pre>

<p>And a softmasked kmer that we can expect more than once</p>

<pre><code>kmer = 'ccgaaccacgNtggaggcctg'.upper()                                                                                                                                                                                                                              

searchKmer(kmer)                                                                                                                                                                                                                                                    
Out[31]: 
[['CCGAACCACGATGGAGGCCTG', '1_Tps_b3v07_scaf000306', 5196, 'A'],
 ['CCGAACCACGCTGGAGGCCTG', '1_Tps_b3v07_scaf008224', 15781, 'C'],
 ['CCGAACCACGGTGGAGGCCTG', '1_Tps_b3v07_scaf001152', 71261, 'G'],
 ['CCGAACCACGTTGGAGGCCTG', '1_Tps_b3v07_scaf002780', 1721, 'T']]
</code></pre>

<p>So to map all your kmers (assuming a file with a kmer per line)</p>

<pre><code>kmer_file_name = 'my_kmer_file.txt'
with open(kmer_file_name, 'r') as kmer_file:
    kmers = [kmer.rstrip() for kmer in kmer_file]

mapping_list = [searchKmer(kmer) for kmer in kmers]
</code></pre>

<p>Stuff can be probably done in a more elegant way, but this is a working solution that scales to any genome size.</p>

CC BY-SA 4.0
2019-03-21T13:23:32.123Z
2019-04-26T09:04:11.413Z
2019-04-26T09:04:11.413Z
Kamil S Jaron
<p>I don't know if this is the best answer, but one quick-n-dirty approach I commonly use involves tweaking <code>bwa mem</code> 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.</p>
<pre><code>bwa mem -k 21 -T 21 -a -c 5000 assembly.fasta your-21mers.fasta
</code></pre>
<p>The <code>-k 21</code> says to use a minimum seed length of 21, which forces exact matches. The <code>-T 21</code> requires a minimum score of 21, which also enforces an exact match. The <code>-a</code> parameter reports all matches, since a &quot;best&quot; match doesn't make sense in this situation. The <code>-c</code> parameter limits how many matches are reported, which may need to be adjusted depending on how repetitive the 21-mer is.</p>
<p><strong>UPDATE</strong> If you don't want to introduce a dependency like <a href="https://github.com/nanoporetech/bwapy" rel="nofollow noreferrer">bwapy</a>, it's fairly straightforward to implement this in pure Python.</p>
<pre class="lang-py prettyprint-override"><code>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
</code></pre>
<p>Also, one can <a href="https://bioinformatics.stackexchange.com/a/511/96">call bwa mem via a C API</a>, but the mechanics of this are left as an exercise for the reader. </p>

CC BY-SA 4.0
2019-03-20T16:34:34.347Z
2021-09-13T19:08:00.553Z
2021-09-13T19:08:00.553Z
Daniel Standage
<p>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:</p>

<pre><code>AAAAAAAAAANAAAAAAAAAA
</code></pre>

<p>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</p>

<ol>
<li><strong>Using an existing mapper</strong></li>
</ol>

<p>I tried to use python binding of <a href="https://github.com/KamilSJaron/minimap2/tree/master/python#class-mappyaligner" rel="nofollow noreferrer">minimap2</a>, 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 <code>k</code>, <code>w</code> and <code>min_cnt</code>, 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?</p>

<ol start="2">
<li><strong>Hash table</strong></li>
</ol>

<p>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).</p>

<p>I am also suspecting that there will be an elegant solution using a suffix array, but I have not cracked that one yet.</p>

<genome><algorithms><k-mer>
mapping heteryzygous kmers on a genome
# mapping heteryzygous kmers on a genome

(Asked by: username_0 on 2019-03-20)

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.

**Tags:** genome, algorithms, k-mer

---

## 2 answers

### Answer by username_1 on 2019-03-20. Score: 5 | Accepted answer

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('>{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.


---

### Answer by username_0 on 2019-03-21. Score: 4

Indeed, with suffix arrays, you can map a kmer in ~0.0007s. The way to go is to build a suffix array (for instance using PySAIS package) and then do a convoluted binary search of the four possible kmers. A full implementation of such mapper can be found here. 
To build suffix array
from Bio import SeqIO
from Bio.Seq import Seq
from PySAIS import sais
from bisect import bisect_right

# load the genome
kmer_genome_file = 'my_genome.fasta'
ffile = SeqIO.parse(kmer_genome_file, "fasta")
sequences = []
scf_names = []
for seq_record in ffile:
    scf_names.append(seq_record.name)
    sequences.append(str(seq_record.seq).upper())

ffile.close()

# paste all scaffolds into a single string and build a suffix array
genome = '$'.join(sequences)
sa = sais(genome)

To transform genome coordinate to scaffold/position coordinates we define a function evaluated2assembly_position
def evaluated2assembly_position(eval_pos):
    scf_index = bisect_right(scf_sizes, eval_pos)
    eval_pos - scf_sizes[scf_index - 1]
    if scf_index > 0:
        return((scf_names[scf_index], eval_pos - scf_sizes[scf_index - 1]))
    else:
        return((scf_names[scf_index], eval_pos))

and finally the function for a convoluted binary search (assuming all kmers being 21 nt long and having the unknown nucleotide in the middle, on 11th nt). Note that you need to map to both strands of DNA, which is equivalent of mapping kmer and it's reverse complementary sequence
def searchKmer(self, kmer):
    reverse_complementary_kmer = str(Seq(kmer).reverse_complement())
    alignments = []
    for strand in ['+', '-']:
        if strand == '-':
            kmer = reverse_complementary_kmer
        L_kmer = kmer[0:10]
        l = 0
        r = len(sa) - 1
        while l <= r:
            m = (l + r) // 2
            eval_pos = sa[m]
            L_genome_kmer = genome[eval_pos:(eval_pos + 10)]
            if L_kmer == L_genome_kmer:
                r_L_converged = r
                for N in ['A', 'C', 'G', 'T']:
                    kmer = kmer[0:10] + N + kmer[11:21]
                    while l <= r:
                        m = (l + r) // 2
                        eval_pos = sa[m]
                        genome_kmer = genome[eval_pos:(eval_pos + 21)]
                        if kmer == genome_kmer:
                            m_hit = m
                            # following ms
                            while kmer == genome_kmer:
                                scf, pos = evaluated2assembly_position(eval_pos)
                                alignments.append([scf, pos, strand, N])
                                m += 1
                                eval_pos = sa[m]
                                genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            # previous ms
                            m = m_hit - 1
                            eval_pos = sa[m]
                            genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            while kmer == genome_kmer:
                                scf, pos = evaluated2assembly_position(eval_pos)
                                alignments.append([scf, pos, strand, N])
                                m -= 1
                                eval_pos = sa[m]
                                genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            break
                        elif genome_kmer < kmer:
                            l = m + 1
                        else:
                            r = m - 1
                    r = r_L_converged # recover the r after converging on the left sub-kmer
                break
            elif L_genome_kmer < L_kmer:
                l = m + 1
            else:
                r = m - 1
    return(alignments)

Taking one of the fist kmers in the scaffold one to test it
searchKmer('ACGGCCACCCNGTGTCGTTGT')                                                                                                                                                                                                                                                    
Out[23]: [['scaf000001', 28, 'G']]

And a softmasked kmer that we can expect more than once
kmer = 'ccgaaccacgNtggaggcctg'.upper()                                                                                                                                                                                                                              

searchKmer(kmer)                                                                                                                                                                                                                                                    
Out[31]: 
[['CCGAACCACGATGGAGGCCTG', '1_Tps_b3v07_scaf000306', 5196, 'A'],
 ['CCGAACCACGCTGGAGGCCTG', '1_Tps_b3v07_scaf008224', 15781, 'C'],
 ['CCGAACCACGGTGGAGGCCTG', '1_Tps_b3v07_scaf001152', 71261, 'G'],
 ['CCGAACCACGTTGGAGGCCTG', '1_Tps_b3v07_scaf002780', 1721, 'T']]

So to map all your kmers (assuming a file with a kmer per line)
kmer_file_name = 'my_kmer_file.txt'
with open(kmer_file_name, 'r') as kmer_file:
    kmers = [kmer.rstrip() for kmer in kmer_file]

mapping_list = [searchKmer(kmer) for kmer in kmers]

Stuff can be probably done in a more elegant way, but this is a working solution that scales to any genome size.


---

