mapping heteryzygous kmers on a genome

Question

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

1. 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?

2. 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.

Answer 1

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 2

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.