Introduction to the kmer R package

Introduction

Agglomerative clustering methods that rely on a multiple sequence alignment and a matrix of pairwise distances can be computationally infeasible for large DNA and amino acid datasets. Alternative k-mer based clustering methods involve enumerating all k-letter words in a sequence through a sliding window of length k. The \(n \times 4^k\) matrix of k-mer counts (where \(n\) is the number of sequences) can then be used in place of a multiple sequence alignment to calculate distances and/or build a phylogenetic tree. kmer is an R package for clustering large sequence datasets using fast alignment-free k-mer counting. This can be achieved with or without a multiple sequence alignment, and with or without a matrix of pairwise distances. These functions are detailed below with examples of their utility.

Distance matrix computation

The function kcount is used to enumerate all k-mers within a sequence or set of sequences, by sliding a window of length k along each sequence and counting the number of times each k-mer appears (for example, the \(4^3 = 64\) possible DNA 3-mers: AAA, AAC, AAG, …, TTT). The kdistance function can then compute an alignment-free distance matrix, using a matrix of k-mer counts to derive the pairwise distances. The default distance metric used by kdistance is the k-mer (k-tuple) distance measure outlined in Edgar (2004). For two DNA sequences \(a\) and \(b\), the fractional common k-mer count over the \(4^k\) possible words of length \(k\) is calculated as: \[F = \sum\limits_{\tau}\frac{min (n_a(\tau), n_b (\tau))}{min (L_a , L_b ) - k + 1} \tag{1}\]

where \(\tau\) represents each possible k-mer, \(n_a(\tau)\) and \(n_b(\tau)\) are the number of times \(\tau\) appears in each sequence, \(k\) is the k-mer length and \(L\) is the sequence length. The pairwise distance between \(a\) and \(b\) is then calculated as:

\[d = \frac{log(0.1 + F) - log(1.1)}{log(0.1)} \tag{2}\]

For \(n\) sequences, the kdistance operation has time and memory complexity \(O(n^2)\) and thus can become computationally infeasible when the sequence set is large (e.g. > 10,000 sequences). As such, the kmer package also offers the function mbed, that only computes the distances from each sequence to a smaller (or equal) sized subset of ‘seed’ sequences (Blackshields et al., 2010). The default behavior of the mbed function is to select \(t = (log_2n)^2\) seeds by clustering the sequences (k-means algorithm with \(k = t\)), and selecting one representative sequence from each cluster.

DNA and amino acid sequences can be passed to kcount, kdistance and mbed either as a list of non-aligned sequences or a matrix of aligned sequences, preferably in either the “DNAbin” or “AAbin” raw-byte format (see the ape package documentation for more information on these S3 classes). Character sequences are supported; however ambiguity codes may not be recognized or treated appropriately, since raw ambiguities are counted according to their underlying residue frequencies (e.g. the 5-mer “ACRGT” would contribute 0.5 to the tally for “ACAGT” and 0.5 to that of “ACGGT”). This excludes the ambiguity code “N”, which is ignored.

data(woodmouse, package = "ape")
ape::as.character.DNAbin(woodmouse[1:5, 1:5])
#>         [,1] [,2] [,3] [,4] [,5]
#> No305   "n"  "t"  "t"  "c"  "g" 
#> No304   "a"  "t"  "t"  "c"  "g" 
#> No306   "a"  "t"  "t"  "c"  "g" 
#> No0906S "a"  "t"  "t"  "c"  "g" 
#> No0908S "a"  "t"  "t"  "c"  "g"

woodmouse <- woodmouse[, apply(woodmouse, 2, function(v) !any(v == 0xf0))]

### Compute the full distance matrix and print the first few rows and columns
library(kmer)
woodmouse.kdist <- kdistance(woodmouse, k = 6)
print(as.matrix(woodmouse.kdist)[1:7, 1:7], digits = 2)
#>         No305  No304  No306 No0906S No0908S No0909S No0910S
#> No305   0.000 0.0322 0.0295   0.033   0.036   0.037   0.037
#> No304   0.032 0.0000 0.0051   0.020   0.022   0.032   0.023
#> No306   0.030 0.0051 0.0000   0.016   0.017   0.026   0.018
#> No0906S 0.033 0.0202 0.0162   0.000   0.024   0.033   0.014
#> No0908S 0.036 0.0224 0.0171   0.024   0.000   0.033   0.025
#> No0909S 0.037 0.0322 0.0264   0.033   0.033   0.000   0.034
#> No0910S 0.037 0.0233 0.0176   0.014   0.025   0.034   0.000

### Compute and print the embedded distance matrix
suppressWarnings(RNGversion("3.5.0"))
set.seed(999)
seeds <- sample(1:15, size = 3)
woodmouse.mbed <- mbed(woodmouse, seeds = seeds, k = 6)
#> Registered S3 method overwritten by 'openssl':
#>   method      from
#>   print.bytes Rcpp
print(woodmouse.mbed[,], digits = 2)
#>         No0909S No0913S  No304
#> No305    0.0368  0.0391 0.0322
#> No304    0.0322  0.0102 0.0000
#> No306    0.0264  0.0098 0.0051
#> No0906S  0.0332  0.0215 0.0202
#> No0908S  0.0332  0.0273 0.0224
#> No0909S  0.0000  0.0368 0.0322
#> No0910S  0.0341  0.0176 0.0233
#> No0912S  0.0242  0.0322 0.0273
#> No0913S  0.0368  0.0000 0.0102
#> No1103S  0.0171  0.0251 0.0202
#> No1007S  0.0046  0.0368 0.0322
#> No1114S  0.0451  0.0428 0.0373
#> No1202S  0.0345  0.0176 0.0233
#> No1206S  0.0304  0.0251 0.0202
#> No1208S  0.0046  0.0409 0.0359

## compute pairwise distance matrices
dist1 <- ape::dist.dna(woodmouse, model = "K80") 
dist2 <- kdistance(woodmouse, k = 7) 

## build neighbor-joining trees
phy1 <- ape::nj(dist1)
phy2 <- ape::nj(dist2)

## rearrange trees in ladderized fashion
phy1 <- ape::ladderize(phy1)
phy2 <- ape::ladderize(phy2)

## convert phylo objects to dendrograms
dnd1 <- as.dendrogram(phy1)
dnd2 <- as.dendrogram(phy2)

## plot the tanglegram
dndlist <- dendextend::dendlist(dnd1, dnd2)
dendextend::tanglegram(dndlist, fast = TRUE, margin_inner = 5)

Clustering without a distance matrix

To avoid excessive time and memory use when building large trees (e.g. n > 10,000), the kmer package features the function cluster for fast divisive clustering, free of both alignment and distance matrix computation. This function first generates a matrix of k-mer counts, and then recursively partitions the matrix row-wise using successive k-means clustering (k = 2). While this method may not necessarily reconstruct sufficiently accurate phylogenetic trees for taxonomic purposes, it offers a fast and efficient means of producing large trees for a variety of other applications such as tree-based sequence weighting (e.g. Gerstein et al. (1994)), guide trees for progressive multiple sequence alignment (e.g. Sievers et al. (2011)), and other recursive operations such as classification and regression tree (CART) learning.

The package also features the function otu for rapid clustering of sequences into operational taxonomic units based on a genetic distance (k-mer distance) threshold. This function performs a similar operation to cluster in that it recursively partitions a k-mer count matrix to assign sequences to groups. However, the top-down splitting only continues while the highest k-mer distance within each cluster is above a defined threshold value. Rather than returning a dendrogram, otu returns a named integer vector of cluster membership, with asterisks indicating the representative sequences within each cluster.

suppressWarnings(RNGversion("3.5.0"))
set.seed(999)
woodmouse.OTUs <- otu(woodmouse, k = 5, threshold = 0.97, method = "farthest", nstart = 20)
woodmouse.OTUs
#>   No305*    No304   No306*  No0906S  No0908S No0909S*  No0910S  No0912S 
#>        3        1        1        1        1        2        1        2 
#>  No0913S  No1103S  No1007S  No1114S  No1202S  No1206S  No1208S 
#>        1        2        2        3        1        1        2

Concluding remarks

The kmer package is released under the GPL-3 license. Please direct bug reports to the GitHub issues page at http://github.com/shaunpwilkinson/kmer/issues. Any feedback is greatly appreciated.

Acknowledgements

This software was developed with funding from a Rutherford Foundation Postdoctoral Research Fellowship from the Royal Society of New Zealand.