Logistic Regression and SVMs

Logistic regression and support vector machines are the two workhorse linear classifiers you reach for first on sequence data. Both learn a decision boundary of the form w dot x plus b, but they optimize different loss functions and expose different knobs. This lesson covers both conceptually and in code, using k-mer features from DNA sequences, which are the canonical high-dimensional sparse representation in genomics.

Logistic regression

Logistic regression models the log-odds of the positive class as a linear function of the features, then passes it through the sigmoid to produce a class-probability estimate between 0 and 1, which is well calibrated only when the model is correctly specified. The penalty term controls overfitting: L2 (ridge) shrinks coefficients smoothly toward zero, while L1 (lasso) drives many coefficients to exactly zero, giving built-in feature selection that is invaluable when your k-mer vocabulary has thousands of columns. The C hyperparameter is the inverse of regularization strength, so smaller C means stronger shrinkage.

import numpy as np
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score

# seqs: list[str] over {A,C,G,T};  y: 1 = promoter, 0 = background
def kmerize(seq, k=4):
    return " ".join(seq[i:i + k] for i in range(len(seq) - k + 1))

vec = CountVectorizer(analyzer=str.split, lowercase=False)
X = vec.fit_transform(kmerize(s) for s in seqs)          # sparse (n, up to 256)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25,
                                      stratify=y, random_state=0)

l2 = LogisticRegression(penalty="l2", C=1.0, max_iter=1000).fit(Xtr, ytr)
l1 = LogisticRegression(penalty="l1", solver="liblinear", C=0.5,
                        max_iter=1000).fit(Xtr, ytr)

print("L2 test AUC:", round(roc_auc_score(yte, l2.decision_function(Xte)), 3))
print("L1 test AUC:", round(roc_auc_score(yte, l1.decision_function(Xte)), 3))
print("L1 nonzero coefs:", int((l1.coef_ != 0).sum()), "of", l1.coef_.size)

feats = np.array(vec.get_feature_names_out())
top = np.argsort(l2.coef_[0])[-5:][::-1]
print("top promoter 4-mers:", list(feats[top]))

L2 test AUC: 0.94
L1 test AUC: 0.932
L1 nonzero coefs: 48 of 256
top promoter 4-mers: ['TATA', 'ATAA', 'TAAA', 'CAAT', 'GCGC']

Each coefficient is the change in log-odds per unit increase in that k-mer count, holding the others fixed; exponentiating it gives an odds ratio. Positive weights push a sequence toward the promoter class, and here the TATA and CAAT motifs surface exactly as biology predicts. The L1 model kept only 48 of 256 k-mers, yielding a short interpretable motif list, but note that correlated k-mers make individual weights unstable, so interpret groups of related features rather than trusting any single coefficient.

Support vector machines

An SVM finds the maximum-margin separating hyperplane, optimizing hinge loss so only the boundary-defining support vectors matter. The kernel trick lets it compute dot products in an implicit, possibly infinite-dimensional feature space via a kernel function K, without ever materializing that mapping. A linear kernel is just the raw dot product; the RBF kernel is exp of minus gamma times the squared distance, which builds a smooth nonlinear boundary where gamma sets how local each support vector's influence is. Larger gamma means tighter, more wiggly boundaries and a higher overfitting risk.

from sklearn.svm import SVC, LinearSVC
from sklearn.model_selection import cross_val_score

# reuse the sparse k-mer matrix X and labels y from above
# LinearSVC (liblinear) scales to very wide sparse matrices far faster than SVC
models = [
    ("linear SVC", SVC(kernel="linear", C=1.0)),
    ("rbf SVC", SVC(kernel="rbf", C=1.0, gamma="scale")),
    ("LinearSVC", LinearSVC(C=1.0, dual=True, max_iter=5000)),
]
for name, clf in models:
    auc = cross_val_score(clf, X, y, cv=5, scoring="roc_auc").mean()
    print(f"{name:10s} AUC: {auc:.3f}")

rbf = SVC(kernel="rbf", C=1.0, gamma="scale").fit(X, y)
print("rbf support vectors per class:", rbf.n_support_)

linear SVC AUC: 0.951
rbf SVC    AUC: 0.938
LinearSVC  AUC: 0.949
rbf support vectors per class: [141 137]