An introduction to k-mers for genome comparison and analysis

k-mers provide sensitive and specific methods for comparing and analyzing genomes.

This notebook provides pure Python implementations of some of the basic k-mer comparison techniques implemented in sourmash, including hash-based subsampling techniques.

Running this notebook.

You can run this notebook interactively via mybinder; click on this button:

A rendered version of this notebook is available at sourmash.readthedocs.io under “Tutorials and notebooks”.

You can also get this notebook from the doc/ subdirectory of the sourmash github repository. See binder/environment.yaml for installation dependencies.

What is this?

This is a Jupyter Notebook using Python 3. If you are running this via binder, you can use Shift-ENTER to run cells, and double click on code cells to edit them.

Contact: C. Titus Brown, ctbrown@ucdavis.edu. Please file issues on GitHub if you have any questions or comments!

Calculating Jaccard similarity and containment

Given any two collections of k-mers, we can calculate similarity and containment using the union and intersection functionality in Python.

def jaccard_similarity(a, b):
    a = set(a)
    b = set(b)

    intersection = len(a.intersection(b))
    union = len(a.union(b))

    return intersection / union

def jaccard_containment(a, b):
    a = set(a)
    b = set(b)

    intersection = len(a.intersection(b))

    return intersection / len(a)

Let’s try these functions out on some simple examples!

a = ['ATGG', 'AACC']
b = ['ATGG', 'CACA']
c = ['ATGC', 'CACA']

jaccard_similarity(a, a)

1.0

jaccard_containment(a, a)

1.0

jaccard_similarity(b, a)

0.3333333333333333

jaccard_similarity(a, c)

0.0

jaccard_containment(b, a)

0.5

%matplotlib inline
from matplotlib_venn import venn2, venn3

venn2([set(a), set(b)])

<matplotlib_venn._common.VennDiagram at 0x132a99340>

venn3([set(a), set(b), set(c)])

<matplotlib_venn._common.VennDiagram at 0x132c42430>

Calculating k-mers from DNA sequences

To extract k-mers from DNA sequences, we walk over the sequence with a sliding window:

def build_kmers(sequence, ksize):
    kmers = []
    n_kmers = len(sequence) - ksize + 1

    for i in range(n_kmers):
        kmer = sequence[i:i + ksize]
        kmers.append(kmer)

    return kmers

build_kmers('ATGGACCAGATATAGGGAGAGCCAGGTAGGACA', 21)

['ATGGACCAGATATAGGGAGAG',
 'TGGACCAGATATAGGGAGAGC',
 'GGACCAGATATAGGGAGAGCC',
 'GACCAGATATAGGGAGAGCCA',
 'ACCAGATATAGGGAGAGCCAG',
 'CCAGATATAGGGAGAGCCAGG',
 'CAGATATAGGGAGAGCCAGGT',
 'AGATATAGGGAGAGCCAGGTA',
 'GATATAGGGAGAGCCAGGTAG',
 'ATATAGGGAGAGCCAGGTAGG',
 'TATAGGGAGAGCCAGGTAGGA',
 'ATAGGGAGAGCCAGGTAGGAC',
 'TAGGGAGAGCCAGGTAGGACA']

In the k-mers that are output, you can see how the sequence shifts to the right - look at the pattern in the middle.

So, now, you can compare two sequences!

seq1 = 'ATGGACCAGATATAGGGAGAGCCAGGTAGGACA'
seq2 = 'ATGGACCAGATATTGGGAGAGCCGGGTAGGACA'
# differences:       ^         ^

K = 10
kmers1 = build_kmers(seq1, K)
kmers2 = build_kmers(seq2, K)

print(K, jaccard_similarity(kmers1, kmers2))

10 0.09090909090909091

Reading k-mers in from a file

In practice, we often need to work with 100s of thousands of k-mers, and this means loading them in from sequences in files.

There are three cut-down genome files in the genomes/ directory that we will use below:

akkermansia.fa
shew_os185.fa
shew_os223.fa

The latter two are two strains of Shewanella baltica, and the first one is an unrelated genome Akkermansia muciniphila.

import screed # a library for reading in FASTA/FASTQ

def read_kmers_from_file(filename, ksize):
    all_kmers = []
    for record in screed.open(filename):
        sequence = record.sequence

        kmers = build_kmers(sequence, ksize)
        all_kmers += kmers

    return all_kmers

akker_kmers = read_kmers_from_file('genomes/akkermansia.fa', 31)

akker_kmers[:5]

['AAATCTTATAAAATAACCACATAACTTAAAA',
 'AATCTTATAAAATAACCACATAACTTAAAAA',
 'ATCTTATAAAATAACCACATAACTTAAAAAG',
 'TCTTATAAAATAACCACATAACTTAAAAAGA',
 'CTTATAAAATAACCACATAACTTAAAAAGAA']

print(len(akker_kmers))

499970

shew1_kmers = read_kmers_from_file('genomes/shew_os185.fa', 31)
shew2_kmers = read_kmers_from_file('genomes/shew_os223.fa', 31)

We can see the relationship between these three like so:

print('akker vs shew1', jaccard_similarity(akker_kmers, shew1_kmers))
print('akker vs shew2', jaccard_similarity(akker_kmers, shew2_kmers))
print('shew1 vs shew2', jaccard_similarity(shew1_kmers, shew2_kmers))

akker vs shew1 0.0
akker vs shew2 0.0
shew1 vs shew2 0.23675152210020398

print('akker vs shew1', jaccard_containment(akker_kmers, shew1_kmers))
print('akker vs shew2', jaccard_containment(akker_kmers, shew2_kmers))
print('shew1 vs shew2', jaccard_containment(shew1_kmers, shew2_kmers))

akker vs shew1 0.0
akker vs shew2 0.0
shew1 vs shew2 0.38397187523995907

venn3([set(akker_kmers), set(shew1_kmers), set(shew2_kmers)])

<matplotlib_venn._common.VennDiagram at 0x15aa8d8b0>

Let’s hash!

Choose a hash function!

We need to pick a hash function that takes DNA k-mers and converts them into numbers.

Both the mash software for MinHash, and the sourmash software for modulo and MinHash, use MurmurHash:

https://en.wikipedia.org/wiki/MurmurHash

this is implemented in the ‘mmh3’ library in Python.

The other thing we need to do here is take into account the fact that DNA is double stranded, and so

hash_kmer('ATGG')

should be equivalent to

hash_kmer('CCAT')

Following mash’s lead, for every input k-mer we will choose a canonical k-mer that is the lesser of the k-mer and its reverse complement.

import mmh3

def hash_kmer(kmer):
    # calculate the reverse complement
    rc_kmer = screed.rc(kmer)

    # determine whether original k-mer or reverse complement is lesser
    if kmer < rc_kmer:
        canonical_kmer = kmer
    else:
        canonical_kmer = rc_kmer

    # calculate murmurhash using a hash seed of 42
    hash = mmh3.hash64(canonical_kmer, 42)[0]
    if hash < 0: hash += 2**64

    # done
    return hash

This is now a function that we can use to turn any DNA “word” into a number:

hash_kmer('ATGGC')

13663093258475204077

The same input word always returns the same number:

hash_kmer('ATGGC')

13663093258475204077

as does its reverse complement:

hash_kmer('GCCAT')

13663093258475204077

and nearby words return very different numbers:

hash_kmer('GCCAA')

1777382721305265773

Note that hashing collections of k-mers doesn’t change Jaccard calculations:

def hash_kmers(kmers):
    hashes = []
    for kmer in kmers:
        hashes.append(hash_kmer(kmer))
    return hashes

shew1_hashes = hash_kmers(shew1_kmers)
shew2_hashes = hash_kmers(shew2_kmers)

print(jaccard_similarity(shew1_kmers, shew2_kmers))

0.23675152210020398

print(jaccard_similarity(shew1_hashes, shew2_hashes))

0.2371520123045373

(ok, it changes it a little, because of the canonical k-mer calculation!)

Implementing subsampling with modulo hashing

We are now ready to implement k-mer subsampling with modulo hash.

We need to pick a sampling rate, and know the maximum possible hash value.

For a sampling rate, let’s start with 1000.

The MurmurHash function turns k-mers into numbers between 0 and 2**64 - 1 (the maximum 64-bit number).

Let’s define these as variables:

scaled = 1000
MAX_HASH = 2**64

Now, choose the range of hash values that we’ll keep.

keep_below = MAX_HASH / scaled
print(keep_below)

1.844674407370955e+16

and write a filter function:

def subsample_modulo(kmers):
    keep = []
    for kmer in kmers:
        if hash_kmer(kmer) < keep_below:
            keep.append(kmer)
        # otherwise, discard

    return keep

Now let’s apply this to our big collections of k-mers!

akker_sub = subsample_modulo(akker_kmers)
shew1_sub = subsample_modulo(shew1_kmers)
shew2_sub = subsample_modulo(shew2_kmers)

print(len(akker_kmers), len(akker_sub))
print(len(shew1_kmers), len(shew1_sub))
print(len(shew2_kmers), len(shew2_sub))

499970 502
499970 513
499970 503

So we go from ~500,000 k-mers to ~500 hashes! Do the Jaccard calculations change??

print('akker vs akker, total', jaccard_similarity(akker_kmers, akker_kmers))
print('akker vs akker, sub', jaccard_similarity(akker_sub, akker_sub))

akker vs akker, total 1.0
akker vs akker, sub 1.0

print('akker vs shew1, total', jaccard_similarity(akker_kmers, shew1_kmers))
print('akker vs shew1, sub', jaccard_similarity(akker_sub, shew1_sub))

akker vs shew1, total 0.0
akker vs shew1, sub 0.0

print('shew1 vs shew2, total', jaccard_similarity(shew1_kmers, shew2_kmers))
print('shew1 vs shew2, sub', jaccard_similarity(shew1_sub, shew2_sub))

shew1 vs shew2, total 0.23675152210020398
shew1 vs shew2, sub 0.2281795511221945

And you can see that the numbers are different, but not very much - the Jaccard similarity is being estimated, so it is not exact but it is close.

Let’s visualize –

venn3([set(akker_kmers), set(shew1_kmers), set(shew2_kmers)])

<matplotlib_venn._common.VennDiagram at 0x163888190>

venn3([set(akker_sub), set(shew1_sub), set(shew2_sub)])

<matplotlib_venn._common.VennDiagram at 0x1638e11c0>

Other pointers

Sourmash: a practical guide

Classifying signatures taxonomically

Pre-built search databases

A full list of notebooks

An introduction to k-mers for genome comparison and analysis

Some sourmash command line examples!

Working with private collections of signatures.

Using the LCA_Database API.