DNA Sequence Classification With Machine Learning

Almost every classical DNA classifier is the same two-part machine: a fixed-length feature vector built from substring counts, and an ordinary classifier on top of it. The first part is where all the modelling decisions live, and it is usually the part explained least.

The representation problem

A logistic regression, a gradient-boosted tree ensemble and a linear support vector machine all want the same thing: a vector of fixed length, the same length for every example. DNA does not arrive that way. One read is 150 bases, one gene is 2,400, one contig is 80,000. You cannot pad to the longest, because the longest is orders of magnitude longer than the median, and you cannot truncate, because the discriminating region may sit anywhere.

The k-mer approach solves this by throwing away position entirely. Count how often each length-k substring occurs, and you have a vector whose length depends only on k and the alphabet, never on the input. A sequence of any length becomes a point in the same space. That is the trade: you get a usable feature vector, and you lose the order the substrings appeared in.

A k-mer count vector, worked

Take a 15-base sequence and k equal to 3. Slide a window of width three one base at a time. A sequence of length L yields L minus k plus one windows, so 15 minus 3 plus 1 is 13:

sequence:  A T G C G A T A C G C T T G A
window 1:  ATG
window 2:   TGC
window 3:    GCG
window 4:     CGA
window 5:      GAT
window 6:       ATA
window 7:        TAC
window 8:         ACG
window 9:          CGC
window 10:          GCT
window 11:           CTT
window 12:            TTG
window 13:             TGA

With a four-letter alphabet there are 4 to the power of k possible 3-mers, which is 64. Thirteen of the 64 slots get a count of one and the other 51 stay at zero. To make sequences of different lengths comparable, divide each count by the number of windows: every observed 3-mer here has frequency 1 divided by 13, about 0.077. That normalisation matters more than it looks. Without it, a classifier trained on a mix of 150-base reads and 5,000-base contigs learns length, because raw counts scale with length and length correlates with almost everything in a real dataset.

There is one more step that is skipped surprisingly often. DNA is double-stranded, and a read can come off either strand, so ATG on one strand is CAT read on the other. If your labels do not depend on strand, then treating those as two different features doubles your feature count and halves your effective counts for no gain. The fix is canonicalisation: for each k-mer, compute its reverse complement and keep whichever of the two is lexicographically smaller. In the example above, GCG has reverse complement CGC, so both windows 3 and 9 collapse onto the canonical form CGC, which now has a count of two. For odd k the canonical alphabet is exactly half the size, so canonical 3-mers number 32 rather than 64.

Choosing k, and what each step costs

The feature space grows as 4 to the power of k, and it grows fast. k equal to 3 gives 64 features. k equal to 6 gives 4,096. k equal to 8 gives 65,536. k equal to 12 gives 16,777,216, which is more features than most people have training examples, and at that point almost every feature is zero in almost every row.

• Small k is dense and generic. With k equal to 3 or 4, essentially every k-mer occurs in essentially every sequence, and what the classifier sees is close to a compositional signature — GC content and codon bias in a more detailed form. This is enough to separate coarse categories such as coding from non-coding, or broadly divergent organisms.
• Large k is sparse and specific. A 31-mer is almost unique in a bacterial genome, which is exactly why exact-match taxonomic classifiers use lengths in that range — see taxonomic assignment from short reads. But a sparse binary feature is fragile: one sequencing error destroys the k k-mers that overlap it.
• A single error damages k features, not one. This is the arithmetic that decides your upper bound on k for noisy data. A substitution at position p falsifies every window containing p, which is k windows in the interior of the sequence. At k equal to 31, one miscalled base costs you 31 k-mers.

In practice people either fix k somewhere between 4 and 8 for a classifier over whole genes, or use several values of k and concatenate the vectors, which lets the model use coarse composition and specific motifs at once at the price of a wider matrix.

What classifier goes on top

Once you have the matrix the modelling is ordinary. A multinomial naive Bayes over raw counts is the historical baseline and is fast enough to be worth running first. Logistic regression with L2 regularisation on frequency-normalised features is the usual strong baseline, and it has the practical virtue that the coefficients are per-k-mer, so you can read off which substrings the model is using and check whether they correspond to something known. Gradient-boosted trees do well when the signal is a threshold on a few k-mers rather than a linear combination of many.

The count matrix is sparse and high-dimensional, which is the shape text classifiers were designed for, so the tooling transfers directly: TF-IDF weighting down-weights k-mers that appear in every sequence and therefore separate nothing, and truncated SVD reduces the dimension without materialising a dense matrix. If you would rather have a learned representation than a counted one, that is a different family of methods entirely; see what a sequence model learns instead of counting.

The four ways this goes wrong

• Homology leakage across the train/test split. This is the big one and it is almost always present in a first attempt. Biological sequences are related by descent, so a random split puts near-identical sequences on both sides and your held-out accuracy measures memorisation. The correction is to cluster sequences by identity first and split by cluster, so no cluster spans the boundary. Report the identity threshold you clustered at, because the number changes a lot between 90 percent and 30 percent.
• GC content as a confound. k-mer frequencies encode GC content strongly. If your two classes were sequenced from organisms with different GC, or prepared with different protocols, a classifier can reach high accuracy on that alone and generalise to nothing. Check by fitting on GC content by itself and seeing how much of your accuracy it already explains.
• Position is gone. A k-mer vector cannot represent “this motif appears 30 bases upstream of that one”. Where spacing is the signal — promoter architecture, splice site context — a bag of k-mers is the wrong representation and a convolutional or attention-based model over the one-hot sequence is the right one.
• Class imbalance in the labels, not the sequences.Reference databases are built from what people sequenced, which is heavily skewed toward clinically and agriculturally important organisms. Accuracy on such a set is dominated by the majority class; report per-class recall and a confusion matrix instead.