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.
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.
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!
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']
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 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
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
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 0x11a2e7b00>
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
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!)
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
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>