# Bio Applied Machine Learning For Biology

> Engineer k-mer/GC/CpG DNA features, train scikit-learn classifiers (LogisticRegression, RandomForest, SVC), evaluate with CV/ROC-AUC. Use for promoter/variant classifiers or model comparison on omics features.

- Skill: `pavel-kravchenko/bio-applied-machine-learning-for-biology` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-machine-learning-for-biology`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-machine-learning-for-biology/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: AI & ML
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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).

```python
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.

```python
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.

```python
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

| Algorithm | Strengths | Weaknesses | When |
|-----------|-----------|------------|------|
| Logistic Regression | Interpretable coefficients, fast | Linear boundary only | Fit first as baseline |
| Random Forest | Non-linear, robust, feature importance | Less interpretable | Good default |
| SVM (RBF) | High-dimensional, small datasets | Slow to tune, needs scaling | Try when RF underperforms |
| Gradient Boosting (XGBoost) | Best tabular accuracy | More hyperparameters | When RF plateaus |
| Neural net / CNN | Learns motifs directly from sequence | Data-hungry, needs GPU | Large sequence/expression datasets |

## 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