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.
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.
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.
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)
a = ['ATGG', 'AACC'] b = ['ATGG', 'CACA'] c = ['ATGC', 'CACA']
%matplotlib inline from matplotlib_venn import venn2, venn3 venn2([set(a), set(b)])
<matplotlib_venn._common.VennDiagram at 0x10cb677b8>
venn3([set(a), set(b), set(c)])
<matplotlib_venn._common.VennDiagram at 0x111e74c18>
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
['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))
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)
['AAATCTTATAAAATAACCACATAACTTAAAA', 'AATCTTATAAAATAACCACATAACTTAAAAA', 'ATCTTATAAAATAACCACATAACTTAAAAAG', 'TCTTATAAAATAACCACATAACTTAAAAAGA', 'CTTATAAAATAACCACATAACTTAAAAAGAA']
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 0x11a2e7b00>
We need to pick a hash function that takes DNA k-mers and converts them into numbers.
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
should be equivalent to
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) 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:
The same input word always returns the same number:
as does its reverse complement:
and nearby words return very different numbers:
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)
(ok, it changes it a little, because of the canonical k-mer calculation!)
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)
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
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.
venn3([set(akker_kmers), set(shew1_kmers), set(shew2_kmers)])
<matplotlib_venn._common.VennDiagram at 0x123ede358>
venn3([set(akker_sub), set(shew1_sub), set(shew2_sub)])
<matplotlib_venn._common.VennDiagram at 0x1252f0e48>