Bio Applied Machine Learning For Biology

Machine Learning for Biology

When to Use

• Classifying DNA/protein sequences (promoter vs. non-promoter, enzyme vs. non-enzyme, pathogenic vs. benign variant) from sequence-derived features
• Building a baseline classifier or regressor on tabular omics data (expression, mutation, clinical features) and need to pick between logistic regression, random forest, SVM, k-NN
• Comparing model performance with cross-validation, ROC-AUC, or feature importance instead of a single accuracy number
• Reducing high-dimensional biological data (thousands of genes, hundreds of k-mers) with PCA before modeling
• Diagnosing overfitting/underfitting or avoiding data leakage (train/test contamination, subject-level leakage) in a small biological dataset

Version Compatibility

scikit-learn >=1.4, numpy >=1.26, pandas >=2.0, Python >=3.10. APIs shown (Pipeline, StratifiedKFold, GridSearchCV) are stable back to scikit-learn 1.0.

Prerequisites

• pip install scikit-learn pandas numpy matplotlib
• Familiarity with pandas DataFrames and basic statistics (see statistical-analysis skill)
• For sequence I/O before feature extraction, see bio-sequence-io-read-sequences

Feature Engineering for Sequences

Goal: Convert a raw DNA sequence into a fixed-length numeric feature vector suitable for sklearn estimators. Approach: combine k-mer frequency counts with biologically motivated summary statistics (GC content, CpG observed/expected ratio, TATA-box presence, upstream/downstream GC skew around the TSS).

from itertools import product

def kmer_frequencies(sequence, k=2):
    """Normalized k-mer frequencies. k=2 -> 16 features, k=3 -> 64 features (DNA alphabet)."""
    sequence = sequence.upper()
    all_kmers = [''.join(p) for p in product('ACGT', repeat=k)]
    counts = {km: 0 for km in all_kmers}
    total = 0
    for i in range(len(sequence) - k + 1):
        kmer = sequence[i:i + k]
        if kmer in counts:
            counts[kmer] += 1
            total += 1
    if total > 0:
        counts = {km: c / total for km, c in counts.items()}
    return counts


def extract_promoter_features(sequence, tss_pos=None):
    """Sequence-derived features for promoter vs. non-promoter classification."""
    seq = sequence.upper()
    n = len(seq)
    tss = tss_pos if tss_pos is not None else n // 2
    features = {}

    features['gc_content'] = (seq.count('G') + seq.count('C')) / n

    for km, f in kmer_frequencies(seq, k=2).items():
        features[f'di_{km}'] = f

    n_c, n_g = seq.count('C'), seq.count('G')
    n_cpg = seq.count('CG')
    # CpG observed/expected: real promoters are CpG-enriched relative to background
    features['cpg_oe'] = (n_cpg * n) / (n_c * n_g) if n_c > 0 and n_g > 0 else 0

    # Exact-match TATA search is a teaching simplification; MEME/FIMO score a
    # position-weight matrix and tolerate mismatches (consensus TATAWAWR).
    upstream = seq[max(0, tss - 50):tss]
    features['has_tata'] = 1 if 'TATAAA' in upstream else 0

    up = seq[max(0, tss - 500):tss]
    dn = seq[tss:min(n, tss + 500)]
    features['gc_upstream'] = (up.count('G') + up.count('C')) / max(len(up), 1)
    features['gc_downstream'] = (dn.count('G') + dn.count('C')) / max(len(dn), 1)

    return features

Training and Evaluating Classifiers

Goal: Fit and fairly compare several classifiers on the engineered feature matrix. Approach: stratified train/test split, scale features inside a Pipeline (required for logistic regression/SVM/k-NN, unnecessary for trees), evaluate with StratifiedKFold cross-validation and ROC-AUC rather than raw accuracy.

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split, StratifiedKFold, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import SVC
from sklearn.metrics import classification_report, roc_auc_score


def build_feature_matrix(sequences, labels):
    """Extract per-sequence features and return (X, y, feature_names)."""
    feature_list = [extract_promoter_features(seq) for seq in sequences]
    df = pd.DataFrame(feature_list)
    return df.values, np.asarray(labels), list(df.columns)


def compare_classifiers(X, y, random_state=42):
    """Fit LogisticRegression, RandomForest, SVM on a stratified split and
    report held-out classification metrics plus 5-fold CV ROC-AUC for each."""
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=random_state, stratify=y)

    pipelines = {
        'logreg': Pipeline([('scaler', StandardScaler()),
                             ('clf', LogisticRegression(max_iter=1000, random_state=random_state))]),
        'random_forest': Pipeline([('clf', RandomForestClassifier(
            n_estimators=200, max_depth=10, random_state=random_state, n_jobs=-1))]),
        'svm_rbf': Pipeline([('scaler', StandardScaler()),
                              ('clf', SVC(kernel='rbf', probability=True, random_state=random_state))]),
    }

    cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=random_state)
    results = {}
    for name, pipe in pipelines.items():
        pipe.fit(X_train, y_train)
        y_prob = pipe.predict_proba(X_test)[:, 1]
        auc_test = roc_auc_score(y_test, y_prob)
        cv_scores = cross_val_score(pipe, X, y, cv=cv, scoring='roc_auc')
        results[name] = {'test_auc': auc_test, 'cv_auc_mean': cv_scores.mean(), 'cv_auc_std': cv_scores.std()}
        print(f"=== {name} === test AUC={auc_test:.3f}  CV AUC={cv_scores.mean():.3f}+/-{cv_scores.std():.3f}")
        print(classification_report(y_test, pipe.predict(X_test)))
    return results

Feature Importance and Dimensionality Reduction

Goal: Interpret which features drive predictions, and visualize high-dimensional biological data (e.g., expression of hundreds of genes) in 2D. Approach: pull coef_ (logistic regression) or feature_importances_ (random forest) into a sorted pandas.Series; use PCA before clustering or plotting.

import pandas as pd
from sklearn.decomposition import PCA


def top_features(model, feature_names, n=15):
    """Return the top-n most important features for a fitted linear or tree model."""
    if hasattr(model, 'feature_importances_'):
        importance = model.feature_importances_
    elif hasattr(model, 'coef_'):
        importance = abs(model.coef_[0])
    else:
        raise ValueError("model has neither feature_importances_ nor coef_")
    return pd.Series(importance, index=feature_names).sort_values(ascending=False).head(n)


def pca_projection(X, n_components=2, random_state=42):
    """Reduce a feature matrix to n_components for visualization or downstream clustering."""
    pca = PCA(n_components=n_components, random_state=random_state)
    X_reduced = pca.fit_transform(X)
    print(f"Variance explained: {pca.explained_variance_ratio_.sum():.1%}")
    return X_reduced, pca

Algorithm Decision Table

Pitfalls

• Data leakage: fit StandardScaler/imputers only on training data (inside a Pipeline + CV), never on the full dataset before splitting.
• Subject-level leakage: if multiple samples come from the same patient/individual (e.g., repeated recordings), a plain stratify=y split can put that subject's data in both train and test and inflate accuracy — use GroupKFold/GroupShuffleSplit keyed on subject ID instead.
• Class imbalance: 99% negatives → a classifier that predicts all-negative still scores 99% accuracy; use class_weight='balanced', resampling, or report ROC-AUC/F1 instead of accuracy.
• Sequence similarity leakage: train/test sequences from the same gene family or genomic locus share features — use chromosomal hold-out or CD-HIT clustering (proteins) before splitting.
• k-mer curse of dimensionality: k=6 gives 4096 DNA features; apply feature selection or PCA before linear models when samples are few.
• Batch effects: systematic differences between how positive and negative sets were generated/collected confound the model — balance data source, not just labels.
• Overfitting small bio datasets: hundreds of examples is typical; always report cross-validated performance, not a single train/test split.

See Also

• bio-machine-learning-omics-classifiers — classifiers specialized for expression/mutation matrices
• bio-machine-learning-model-validation — deeper cross-validation and calibration strategies
• bio-machine-learning-biomarker-discovery — feature selection for biomarker panels
• bio-machine-learning-prediction-explanation — SHAP/LIME explanation of trained models