The Statistics of k-mers from a Sequence Undergoing a Simple Mutation Process Without Spurious Matches

Abstract

k-mer-based methods are widely used in bioinformatics, but there are many gaps in our understanding of their statistical properties. Here, we consider the simple model where a sequence S (e.g., a genome or a read) undergoes a simple mutation process through which each nucleotide is mutated independently with some probability r, under the assumption that there are no spurious k-mer matches. How does this process affect the k-mers of S? We derive the expectation and variance of the number of mutated k-mers and of the number of islands (a maximal interval of mutated k-mers) and oceans (a maximal interval of nonmutated k-mers). We then derive hypothesis tests and confidence intervals (CIs) for r given an observed number of mutated k-mers, or, alternatively, given the Jaccard similarity (with or without MinHash). We demonstrate the usefulness of our results using a few select applications: obtaining a CI to supplement the Mash distance point estimate, filtering out reads during alignment by Minimap2, and rating long-read alignments to a de Bruijn graph by Jabba.

Keywords: Jaccard similarity; MinHash; confidence intervals; k-mers; mutation process; sketching.

FIG. 1.

An example of the simple mutation process, with $L=36$ and $k=5$. There are five nucleotides that are mutated (marked with an x). For example, the mutation at position 10 mutates the k-spans $K_6,\dots ,K_{10}$ (marked in red). Note that an isolated nucleotide mutation (e.g., at position 10) can affect up to k k-spans (e.g., $K_6,\dots ,K_{10}$), but nearby nucleotide mutations can affect the same k-span (e.g., mutation of nucleotides at positions 23 and 27 both affect $K_{23}$.) There are two islands (marked in blue) and three oceans, and $N_{mut}=21$. For example, $K_{19},\dots ,K_{34}$ is an island, and $K_{35}$ is an ocean.

FIG. 2.

Estimates of sequence divergence as done by mimimap2 ($\widehat{\delta }$) and by us (${\widehat{r}}_1$). Reads are simulated from a random 10 kbp sequence introducing mutations at the given r₁ rate. For each r₁ value, 100 reads are used. As in Li (2018), we use $k=15$ and, using a random hash function, identify as seeds the k-mer minimizers, one for every window of 25 k-mers. In the case when $\widehat{\delta }$ is undefined, we set $\widehat{\delta }=1$.

FIG. 3.

Box and whisker plot of $C_{ber}$ scores for 5000 replicates of random strings of length 10,000nt, with mutations introduced at a rate of $r_1=0.1$. The solid black line corresponds to the empirical median of $C_{ber}$, while the dashed top line corresponds to $E\left[C_{ber}\right]+z_{0.05}\sqrt{Var\left(C_{ber}\right)}$ and the bottom dot-dashed line corresponds to $E\left[C_{ber}\right]-z_{0.05}\sqrt{Var\left(C_{ber}\right)}$, both computed from Theorem 11.