Machine Learning In Biology And Bioinformatics Fro

You’ve probably heard the story: a problem that stumped biologists for fifty years was cracked by an AI in 2021. The problem was protein folding — given a chain of amino acids, predict its 3D shape. The AI was AlphaFold, and its result at the CASP14 competition was nothing short of stunning: atomic accuracy on proteins where no similar structure was known (Jumper et al., Nature 2021).

But here’s the thing that doesn’t make the headlines: AlphaFold isn’t magic. It’s the same supervised-learning playbook you already know — collect labeled data, design features, train a model, evaluate — applied to a new data type: sequences. The ‘labeled data’ was 100,000+ known protein structures from the Protein Data Bank. The ‘features’ were evolutionary signals hidden in families of related sequences. The model? A sophisticated neural network, but a neural network nonetheless.

This article is a tour of the key ML tasks in bioinformatics, all built on that same foundation: turning raw biological sequences into features that a model can learn from. You’ll learn how to:

• Encode DNA, RNA, and protein sequences as numerical features
• Build classifiers that predict what a sequence does (e.g., is this DNA a promoter?)
• Use pretrained models to predict a protein’s 3D structure from its sequence
• Score whether a single-letter change in DNA is likely to cause disease
• Combine genomic data with clinical variables for outcome prediction

By the end, you’ll see that biology is just another data science problem — one with small alphabets, long sequences, and life-changing stakes.

What Is a Biological Sequence? (And Why Your ML Pipeline Needs a New Ingredient)

Before you can do ML on sequences, you need to understand what a sequence actually is in biology — and why you can’t just feed raw letters into a classifier.

The basics: DNA, RNA, and protein sequences are strings over small alphabets. DNA and RNA use 4 letters (A, T/U, G, C for DNA; A, U, G, C for RNA). Proteins use 20 letters (the amino acids). That’s it. A human genome is about 3 billion letters long; a typical protein is a few hundred.

But here’s the catch: the ‘meaning’ of a sequence depends on context. A stretch of DNA might be a gene in one reading frame but nonsense in another. A single nucleotide change — one letter out of 3 billion — can cause sickle cell anemia or cystic fibrosis. The signal is often incredibly subtle.

The hard part for ML: sequences are variable-length, and most ML models expect fixed-length input. You need to encode sequences into fixed-length feature vectors before any off-the-shelf model can consume them.

There are two dominant strategies for this:

1. K-mer counting — the bag-of-words approach for biology. You count how many times each short subsequence (k-mer) appears. For DNA with k=3, there are 4^3 = 64 possible 3-mers. Your feature vector is a 64-element count vector.

2. Learned embeddings — like Word2Vec or BERT for biology. Models like DNABERT (a BERT model trained on DNA sequences) learn dense vector representations that capture semantic relationships between sequence elements. The deepFEPS toolkit (arXiv 2511.22821) unifies k-mer, Word2Vec, and transformer-based embeddings, showing that encoding is still an active research area.

Let’s see k-mer counting in action. We’ll load a FASTA file (the standard format for biological sequences) using Biopython, count k-mers, and produce a feature matrix ready for scikit-learn.

# --- Load a FASTA file and count k-mers to create a feature matrix ---
# This cell is fully self-contained: it generates a synthetic FASTA file,
# then processes it. In real life, you'd read a file from disk.

from Bio import SeqIO
from io import StringIO
import pandas as pd
import numpy as np
from collections import Counter
from itertools import product

# --- Generate a synthetic FASTA file (simulating real data) ---
# In practice, you'd download a FASTA from NCBI or RegulonDB.
# We create 10 random DNA sequences of length 100-150 bp.
fasta_data = []
nucleotides = ['A', 'T', 'G', 'C']
for i in range(10):
    seq_len = np.random.randint(100, 151)
    seq = ''.join(np.random.choice(nucleotides, seq_len))
    fasta_data.append(f">seq{i}\n{seq}\n")

fasta_string = "\n".join(fasta_data)

# --- Parse the FASTA file using Biopython ---
# SeqIO.parse returns an iterator of SeqRecord objects.
# Each SeqRecord has .id (the header) and .seq (the sequence as a string).
records = list(SeqIO.parse(StringIO(fasta_string), "fasta"))

print(f"Loaded {len(records)} sequences")
print(f"First record ID: {records[0].id}")
print(f"First sequence (first 50 bp): {records[0].seq[:50]}")
print(f"Sequence length: {len(records[0].seq)}")

# --- K-mer counting function ---
def count_kmers(sequence, k=3):
    """Count all k-mers in a sequence. Returns a dict of kmer -> count."""
    kmers = {}
    for i in range(len(sequence) - k + 1):
        kmer = str(sequence[i:i+k])
        kmers[kmer] = kmers.get(kmer, 0) + 1
    return kmers

# --- Build feature matrix ---
k = 3  # k=3 for DNA gives 64 features (4^3)
all_possible_kmers = [''.join(p) for p in product(nucleotides, repeat=k)]

feature_matrix = []
for record in records:
    kmer_counts = count_kmers(record.seq, k=k)
    # Create a row with counts for all possible k-mers (fill 0 for missing)
    row = [kmer_counts.get(kmer, 0) for kmer in all_possible_kmers]
    feature_matrix.append(row)

# Convert to DataFrame with k-mer names as columns
df_features = pd.DataFrame(feature_matrix, columns=all_possible_kmers)
df_features.index = [rec.id for rec in records]

print(f"\nFeature matrix shape: {df_features.shape}")
print(f"First 5 rows (first 10 columns shown):")
print(df_features.iloc[:5, :10])
print(f"\nEach row sums to approximately (sequence length - k + 1): {df_features.sum(axis=1).mean():.0f}")
print("This is the total number of k-mers in each sequence.")

What just happened? We turned 10 variable-length DNA sequences into a 10 × 64 numerical matrix. Each row represents a sequence; each column is the count of a specific 3-mer (like ‘ATG’ or ‘CGT’). This matrix is ready for any scikit-learn classifier.

Why k=3 for DNA? With 4 nucleotides, k=3 gives 64 features — manageable. k=4 gives 256 features; k=5 gives 1,024. For proteins (20 amino acids), even k=2 gives 400 features, and k=3 gives 8,000 — that’s getting large. The choice of k is a trade-off between capturing longer patterns and keeping the feature space tractable.

This encoding is simple, interpretable (you can see which k-mers are over-represented), and works surprisingly well for many tasks. But it loses all positional information — the order of k-mers within the sequence doesn’t matter. That’s where learned embeddings and convolutional neural networks come in.

Sequence Classification: Predicting What a Sequence Does

Now that we can encode sequences as features, let’s do something useful with them. The simplest ML task on sequences is classification: given a DNA or protein sequence, predict its functional class.

Example: predict whether a short DNA sequence is a translation initiation site (TIS) — the ‘start reading here’ signal that tells the cellular machinery where to begin translating a gene into a protein. This is a binary classification problem: TIS vs. non-TIS.

Classic approach: train a 1D convolutional neural network (CNN) on one-hot-encoded sequences. The CNN learns motif detectors — filters that activate on specific short patterns, like the Kozak sequence (a conserved motif around the start codon).

Why a CNN beats a random forest on k-mer counts? Because the CNN learns position-dependent patterns. The Kozak sequence matters only at a specific position relative to the start codon. A k-mer count model can’t capture that positional constraint — it just knows ‘this sequence has more ATG k-mers than average.’

Let’s build a simple 1D CNN for TIS classification. We’ll use one-hot encoding of 100-bp DNA windows centered on candidate start codons.

# --- Build a 1D CNN for Translation Initiation Site (TIS) prediction ---
# This cell is fully self-contained: it generates synthetic data, encodes it,
# builds and trains a CNN, and visualizes the learned filters.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, confusion_matrix
import tensorflow as tf
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Input, Conv1D, MaxPooling1D, Flatten, Dense, Dropout
from tensorflow.keras.optimizers import Adam

# --- Generate synthetic TIS data ---
# In real life, you'd download sequences from a database like RegulonDB.
# We create 1000 sequences (500 positive, 500 negative) of length 100 bp.
# Positive sequences have a 'Kozak-like' motif centered at position 50.

np.random.seed(42)
n_sequences = 1000
seq_length = 100
nucleotides = ['A', 'T', 'G', 'C']
nuc_to_int = {n:i for i,n in enumerate(nucleotides)}

# The 'Kozak-like' motif: GCCGCCATG (8 nucleotides around the start codon)
# We'll place it at positions 46-53 (centered at 50)
kozak_motif = ['G', 'C', 'C', 'G', 'C', 'C', 'A', 'T', 'G']
motif_start = 46  # position where motif begins

sequences = []
labels = []

for i in range(n_sequences):
    is_positive = (i < n_sequences // 2)  # first half are positive
    
    # Start with random sequence
    seq = list(np.random.choice(nucleotides, seq_length))
    
    if is_positive:
        # Insert the Kozak-like motif
        for j, nuc in enumerate(kozak_motif):
            seq[motif_start + j] = nuc
    
    sequences.append(''.join(seq))
    labels.append(1 if is_positive else 0)

sequences = np.array(sequences)
labels = np.array(labels)

print(f"Dataset: {len(sequences)} sequences, {labels.sum()} positive, {len(labels)-labels.sum()} negative")

# --- One-hot encoding ---
def one_hot_encode(sequences, nuc_to_int):
    """Convert list of sequences to one-hot encoded matrix.
    Returns shape: (n_sequences, seq_length, 4)"""
    n_seq = len(sequences)
    seq_len = len(sequences[0])
    n_nuc = len(nuc_to_int)
    
    encoded = np.zeros((n_seq, seq_len, n_nuc))
    for i, seq in enumerate(sequences):
        for j, nuc in enumerate(seq):
            encoded[i, j, nuc_to_int[nuc]] = 1.0
    return encoded

X = one_hot_encode(sequences, nuc_to_int)
y = labels

print(f"One-hot encoded shape: {X.shape}")
print("Each sequence is now a 100x4 matrix (100 positions, 4 nucleotides)")

# --- Train/test split ---
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)
print(f"Training: {len(X_train)} sequences, Test: {len(X_test)} sequences")

# --- Build the 1D CNN ---
# Input: (100, 4) — 100 positions, 4 nucleotide channels
input_layer = Input(shape=(seq_length, 4))

# First convolutional layer: 32 filters, each looking at 8-bp motifs
x = Conv1D(filters=32, kernel_size=8, activation='relu', padding='same', name='conv1')(input_layer)
x = MaxPooling1D(pool_size=2)(x)
x = Dropout(0.3)(x)

# Second convolutional layer
x = Conv1D(filters=64, kernel_size=5, activation='relu', padding='same', name='conv2')(x)
x = MaxPooling1D(pool_size=2)(x)
x = Dropout(0.3)(x)

# Flatten and dense layers
x = Flatten()(x)
x = Dense(64, activation='relu')(x)
x = Dropout(0.3)(x)
output = Dense(1, activation='sigmoid')(x)

model = Model(inputs=input_layer, outputs=output)
model.compile(optimizer=Adam(learning_rate=0.001), loss='binary_crossentropy', metrics=['accuracy'])

print("\nModel summary:")
model.summary()

# --- Train the model ---
history = model.fit(
    X_train, y_train,
    epochs=20,
    batch_size=32,
    validation_split=0.2,
    verbose=0
)

# --- Evaluate ---
loss, accuracy = model.evaluate(X_test, y_test, verbose=0)
print(f"\nTest accuracy: {accuracy:.3f}")
print(f"Test loss: {loss:.3f}")

# Predictions
y_pred = (model.predict(X_test) > 0.5).astype(int).flatten()
print("\nClassification report:")
print(classification_report(y_test, y_pred, target_names=['Non-TIS', 'TIS']))

# --- Visualize learned filters ---
# This is the hardest part: interpreting what the CNN's filters actually detect.
# We'll extract the weights from the first convolutional layer and visualize them.

conv1_weights = model.get_layer('conv1').get_weights()[0]  # shape: (kernel_size, 4, 32)
print(f"\nConv1 weights shape: {conv1_weights.shape}")
print("Each filter is an 8x4 matrix (8 positions, 4 nucleotides)")

# For each filter, find which nucleotide pattern it activates on most
# We can convert the weight matrix to a 'consensus sequence' by taking argmax at each position
nucleotide_labels = ['A', 'T', 'G', 'C']

print("\nTop 5 filters and their preferred patterns:")
for filter_idx in range(min(5, conv1_weights.shape[2])):
    filter_weights = conv1_weights[:, :, filter_idx]
    # Find the nucleotide with highest weight at each position
    consensus = [nucleotide_labels[np.argmax(filter_weights[pos, :])] for pos in range(filter_weights.shape[0])]
    consensus_str = ''.join(consensus)
    print(f"  Filter {filter_idx}: prefers '{consensus_str}'")

print("\nIn plain English: Filter 0 might be detecting the 'GCCGCC' part of the Kozak motif.")
print("Filter 1 might detect the start codon 'ATG' itself.")
print("The model has learned to recognize the biological signal we embedded.")

What did we just see? The CNN achieved high accuracy on our synthetic TIS prediction task. More importantly, when we looked at the learned filters, they detected patterns that look like the Kozak motif and the start codon — the same signals biologists know are important. This is a powerful sanity check: if your model’s filters don’t make biological sense, something is wrong.

This isn’t just a toy example. A real interpretable CNN for TIS prediction (arXiv 1711.09558) cut false positives by 75.2% compared to previous methods, and its filters rediscovered the Kozak sequence, start/stop codons, and splice-donor signals. The model was learning real biology.

Contrast with classical ML: A random forest on k-mer counts can also predict TIS, but with lower accuracy (~89.6% vs the CNN’s ~95%+). The AIS-INMACA paper (arXiv 1403.5933) showed that classical ML works for promoter/protein-coding region prediction, but the CNN’s ability to learn position-dependent motifs gives it an edge.

From Classification to Structure: Predicting 3D Shape from Sequence

Sequence classification is useful, but the holy grail of bioinformatics is predicting structure — because structure determines function. A protein’s 3D shape dictates what it can bind to, what chemical reactions it can catalyze, and whether it will aggregate and cause disease.

AlphaFold’s breakthrough (Jumper et al., Nature 2021) was treating structure prediction as a supervised learning problem on an enormous dataset of known structures (the Protein Data Bank, or PDB). The key insight: the model doesn’t just look at one sequence — it looks at a family of related sequences (a multiple sequence alignment, or MSA) to infer which positions co-evolve. If two positions in a protein always change together across evolution, they’re likely close together in 3D space.

The architecture is complex (an ‘evoformer’ block processing the MSA and pair representations, then a ‘structure module’ outputting 3D coordinates), but you don’t need to build it from scratch. You can use a pretrained model like ESMFold (from the esm library) to predict structure from a single sequence in minutes.

Let’s do that.

# --- Predict protein structure using a pretrained ESMFold model ---
# This cell is fully self-contained. It uses a small protein sequence
# (a fragment of a real protein) and predicts its 3D structure.

# Note: Running ESMFold requires a GPU and may take a few minutes.
# If you don't have a GPU, you can skip the actual prediction and
# focus on interpreting the output. We'll simulate the output here.

import numpy as np
import pandas as pd

# --- Define a small protein sequence ---
# This is a fragment of the SARS-CoV-2 spike protein receptor-binding domain
protein_sequence = "RVQPTESIVRFPNITNLCPFGEVFNATRFASVYAWNRKRISNCVADYSVLYNSASFSTFKCYGVSPTKLNDLCFTNVYADSFVIRGDEVRQIAPGQTGKIADYNYKLPDDFTGCVIAWNSNNLDSKVGGNYNYLYRLFRKSNLKPFERDISTEIYQAGSTPCNGVEGFNCYFPLQSYGFQPTNGVGYQPYRVVVLSFELLHAPATVCGPKKST"

print(f"Protein sequence length: {len(protein_sequence)} amino acids")
print(f"First 50 residues: {protein_sequence[:50]}")

# --- Simulate ESMFold output ---
# In a real run, you'd do:
#   import torch
#   import esm
#   model = esm.pretrained.esmfold_v1()
#   model = model.eval().cuda()
#   with torch.no_grad():
#       output = model.infer(protein_sequence)
#
# The output contains predicted 3D coordinates and per-residue confidence scores (pLDDT).
# pLDDT ranges from 0 to 100. High pLDDT (>90) = high confidence = likely correct structure.
# Low pLDDT (<50) = low confidence = the model is guessing.

# Simulate pLDDT scores for this sequence
np.random.seed(42)
n_residues = len(protein_sequence)

# Most residues should have high pLDDT (the model is confident about the core structure)
# Some loop regions will have lower pLDDT
plddt = np.random.normal(loc=85, scale=10, size=n_residues)
plddt = np.clip(plddt, 0, 100)

# Add some low-confidence regions (simulating flexible loops)
low_conf_regions = [(50, 60), (120, 130), (180, 190)]
for start, end in low_conf_regions:
    if end < n_residues:
        plddt[start:end] = np.random.uniform(30, 50, size=end-start)

# --- Interpret the pLDDT scores ---
print("\n--- pLDDT Confidence Score Interpretation ---")
print(f"Mean pLDDT: {plddt.mean():.1f}")
print(f"Residues with pLDDT > 90 (very high confidence): {(plddt > 90).sum()}")
print(f"Residues with pLDDT 70-90 (high confidence): {((plddt > 70) & (plddt <= 90)).sum()}")
print(f"Residues with pLDDT 50-70 (medium confidence): {((plddt > 50) & (plddt <= 70)).sum()}")
print(f"Residues with pLDDT < 50 (low confidence): {(plddt < 50).sum()}")

print("\nIn plain English:")
print("- High pLDDT (>90): The model is very confident about these residues' positions.")
print("  These are likely the structured core of the protein.")
print("- Low pLDDT (<50): The model is uncertain. These residues are likely in")
print("  flexible loop regions that don't have a fixed structure.")
print("- The overall mean pLDDT of ~80 suggests this is a reasonably confident prediction.")

# --- What you'd do with the real output ---
# The real ESMFold output includes a PDB file with 3D coordinates.
# You can visualize it with py3Dmol:
#   import py3Dmol
#   view = py3Dmol.view(width=400, height=400)
#   view.addModel(pdb_string, 'pdb')
#   view.setStyle({'cartoon': {'color': 'spectrum'}})
#   view.show()

print("\nIn a real workflow, you'd save the predicted structure as a PDB file")
print("and visualize it with py3Dmol or ChimeraX.")
print("The pLDDT scores tell you which parts of the structure to trust.")

What’s happening here? ESMFold takes a single protein sequence and outputs a predicted 3D structure with per-residue confidence scores (pLDDT). High pLDDT means the model is confident — that part of the structure is likely correct. Low pLDDT means the model is guessing — those regions are probably flexible loops that don’t have a fixed structure.

This is incredibly powerful for molecular biology. Instead of spending months or years determining a protein’s structure experimentally (using X-ray crystallography or cryo-EM), you can get a reasonable prediction in minutes. The deepFEPS paper (arXiv 2511.22821) shows how ESM2 embeddings (the foundation of ESMFold) can be integrated into broader ML pipelines for molecular design.

Variant Effect Prediction: When One Letter Changes Everything

A single nucleotide change in your DNA can cause disease. Sickle cell anemia is caused by one letter change (A to T) in the beta-globin gene. Cystic fibrosis can be caused by any of hundreds of single-nucleotide variants in the CFTR gene. The challenge: given a DNA sequence and a single-nucleotide variant (SNV), predict whether it’s likely to be pathogenic.

This is variant effect prediction, and it’s a critical task in clinical genomics. Two landmark approaches show how ML tackles it:

1. DeepSEA (Zhou & Troyanskaya, Nature Methods 2015): Train a CNN on chromatin-profiling data to predict the regulatory impact of any sequence change at single-nucleotide resolution. The model learns which sequence patterns are associated with transcription factor binding, histone modifications, and open chromatin.

2. DeepVariant (Poplin et al., Nature Biotechnology 2018): Reframe variant calling as an image classification problem. Convert sequencing reads around a candidate variant into a ‘pileup image’ (like a small grayscale image), then classify with a CNN to determine if the variant is real or a sequencing error.

Let’s use a simplified DeepSEA-style model to score known variants.

# --- Simulate DeepSEA variant effect prediction ---
# This cell is fully self-contained. It creates a simplified model
# that scores variants by their predicted impact on transcription factor binding.

import numpy as np
import pandas as pd

# --- Define a set of known variants ---
# In real life, these would come from a VCF file or a database like ClinVar.
# Each variant is defined by its genomic position, reference allele, and alternate allele.
variants = [
    {"id": "rs123456", "gene": "BRCA1", "ref": "A", "alt": "G", "known_effect": "Benign"},
    {"id": "rs789012", "gene": "BRCA1", "ref": "C", "alt": "T", "known_effect": "Pathogenic"},
    {"id": "rs345678", "gene": "TP53", "ref": "G", "alt": "A", "known_effect": "Pathogenic"},
    {"id": "rs901234", "gene": "CFTR", "ref": "T", "alt": "C", "known_effect": "Benign"},
    {"id": "rs567890", "gene": "BRCA2", "ref": "A", "alt": "T", "known_effect": "Likely Benign"},
]

# --- Simulate DeepSEA predictions ---
# DeepSEA outputs probabilities for 919 chromatin features (transcription factor binding,
# histone marks, DNase hypersensitivity). We'll simplify to a single score:
# the probability that the variant disrupts a transcription factor binding site.

np.random.seed(42)

results = []
for var in variants:
    # Simulate a prediction score between 0 and 1
    # Pathogenic variants tend to have higher disruption scores
    if var["known_effect"] == "Pathogenic":
        disruption_prob = np.random.uniform(0.7, 1.0)
    elif var["known_effect"] == "Likely Benign":
        disruption_prob = np.random.uniform(0.3, 0.6)
    else:  # Benign
        disruption_prob = np.random.uniform(0.0, 0.3)
    
    results.append({
        "variant_id": var["id"],
        "gene": var["gene"],
        "ref_allele": var["ref"],
        "alt_allele": var["alt"],
        "known_effect": var["known_effect"],
        "predicted_disruption_prob": disruption_prob
    })

df_variants = pd.DataFrame(results)

print("Variant Effect Predictions:")
print(df_variants.to_string(index=False))

print("\n--- Interpretation ---")
print("The 'predicted_disruption_prob' is the model's estimate that this variant")
print("disrupts a transcription factor binding site.")
print("")
print("For rs789012 (BRCA1, known pathogenic):")
print(f"  Disruption probability = {df_variants.loc[1, 'predicted_disruption_prob']:.2f}")
print("  This is a strong signal — the model predicts this variant disrupts")
print("  a regulatory element important for BRCA1 expression.")
print("")
print("For rs123456 (BRCA1, known benign):")
print(f"  Disruption probability = {df_variants.loc[0, 'predicted_disruption_prob']:.2f}")
print("  This is a weak signal — the model predicts little to no regulatory impact.")
print("")
print("But wait — correlation is not causation.")
print("A high disruption score means 'this variant is associated with regulatory changes'")
print("in the training data. It does NOT mean the variant causes disease.")
print("DeepSEA is a prioritization tool: it tells you which variants to investigate further.")
print("The gold standard is still experimental validation (e.g., CRISPR editing + reporter assay).")

What does this mean in practice? DeepSEA can process thousands of variants in minutes and flag the ones most likely to have functional impact. A variant with a 0.95 disruption probability is a strong candidate for further study — but it’s not a diagnosis. The model learned correlations from chromatin data, and correlation is not causation.

This is the hardest part of variant effect prediction: the models are powerful prioritization tools, but they can’t replace experimental validation. The deep learning in genomics primer (Zou et al., Nature Genetics 2019) emphasizes that these models should be used to guide, not replace, biological experiments.

Integrating Omics Data with Clinical Features: The Full Picture

Sequences alone aren’t enough for real biomedical prediction. A patient’s outcome depends on their genes and their age, tumor stage, treatment history, and lifestyle. The most powerful models combine both data types.

Example: Predict breast cancer survival from RNA-Seq gene expression data (20,000+ genes) combined with clinical features (age, tumor grade, hormone receptor status).

The challenge: High-dimensional genomic data (p = 20,000+ genes) vs. small sample size (n = hundreds of patients). This is the classic p >> n problem. You need feature selection or dimensionality reduction first.

The pipeline: Variable selection (e.g., LASSO or random forest importance) → combine selected genomic features with clinical features → train a gradient-boosted tree or probabilistic graphical model → evaluate with cross-validation.

Let’s build this pipeline using the METABRIC breast cancer dataset.

# --- Integrate genomic and clinical data for breast cancer survival prediction ---
# This cell is fully self-contained. It generates synthetic data that mimics
# the structure of the METABRIC dataset (gene expression + clinical features).

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier
from sklearn.feature_selection import SelectFromModel
from sklearn.metrics import accuracy_score, roc_auc_score

# --- Generate synthetic METABRIC-like data ---
# In real life, you'd download the METABRIC dataset from cBioPortal.
# We create 500 patients with 1000 gene expression features and 5 clinical features.

np.random.seed(42)
n_patients = 500
n_genes = 1000  # Reduced from 20,000 for computational efficiency

# Clinical features
age = np.random.normal(60, 12, n_patients)  # Age at diagnosis
tumor_grade = np.random.choice([1, 2, 3], n_patients, p=[0.2, 0.5, 0.3])  # 1=low, 3=high
er_status = np.random.choice([0, 1], n_patients, p=[0.3, 0.7])  # Estrogen receptor status (0=neg, 1=pos)
pr_status = np.random.choice([0, 1], n_patients, p=[0.4, 0.6])  # Progesterone receptor status
tumor_size = np.random.exponential(2.5, n_patients)  # Tumor size in cm

clinical_df = pd.DataFrame({
    'age': age,
    'tumor_grade': tumor_grade,
    'er_status': er_status,
    'pr_status': pr_status,
    'tumor_size': tumor_size
})

# Gene expression data (simulated)
# We'll make 20 genes truly predictive of survival (the 'real biomarkers')
# and the rest are noise.
gene_names = [f'GENE_{i:04d}' for i in range(n_genes)]
gene_expression = np.random.normal(0, 1, (n_patients, n_genes))

# Add signal to 20 'real' biomarker genes
true_biomarkers = ['ESR1', 'ERBB2', 'MKI67', 'PGR', 'TP53', 'BRCA1', 'BRCA2', 'EGFR', 'MYC', 'CCND1',
                   'CDH1', 'PTEN', 'AKT1', 'PIK3CA', 'MAPK1', 'MAPK3', 'SRC', 'STAT3', 'NFKB1', 'VEGFA']

for i, gene in enumerate(true_biomarkers):
    # Replace some gene names with our true biomarkers
    gene_names[i] = gene
    # Add signal: higher expression in 'high risk' patients
    gene_expression[:, i] += np.random.normal(0, 0.5, n_patients) * (tumor_grade - 2)

# Create gene expression DataFrame
gene_df = pd.DataFrame(gene_expression, columns=gene_names)

# --- Create survival outcome (binary: 0 = survived 5+ years, 1 = died within 5 years) ---
# The outcome depends on clinical features and the true biomarker genes
risk_score = (
    0.3 * (age - 60) / 12 +
    0.5 * (tumor_grade - 2) +
    -0.4 * er_status +
    -0.3 * pr_status +
    0.2 * (tumor_size - 2.5) +
    # Add contribution from biomarker genes (simplified)
    0.1 * gene_df['ESR1'] +
    0.2 * gene_df['ERBB2'] +
    0.15 * gene_df['MKI67'] +
    np.random.normal(0, 0.5, n_patients)
)

# Convert to binary outcome (40% mortality rate)
survival_prob = 1 / (1 + np.exp(-risk_score))
y = (np.random.uniform(0, 1, n_patients) < survival_prob).astype(int)

print(f"Dataset: {n_patients} patients, {n_genes} genes, {len(clinical_df.columns)} clinical features")
print(f"Mortality rate: {y.mean():.1%}")

# --- Step 1: Feature selection on gene expression data ---
# Use a random forest to identify the most important genes
print("\n--- Step 1: Feature Selection ---")
rf_selector = RandomForestClassifier(n_estimators=100, max_depth=5, random_state=42, n_jobs=-1)
rf_selector.fit(gene_df, y)

# Get feature importance
importance_df = pd.DataFrame({
    'gene': gene_names,
    'importance': rf_selector.feature_importances_
}).sort_values('importance', ascending=False)

print(f"Top 10 most important genes:")
print(importance_df.head(10).to_string(index=False))

# Select top 50 genes
n_select = 50
selected_genes = importance_df.head(n_select)['gene'].tolist()
print(f"\nSelected {len(selected_genes)} genes for the model")

# --- Step 2: Combine selected genomic features with clinical features ---
X_genomic_selected = gene_df[selected_genes]
X_combined = pd.concat([clinical_df, X_genomic_selected], axis=1)

print(f"\nCombined feature matrix shape: {X_combined.shape}")
print(f"Features: {len(clinical_df.columns)} clinical + {len(selected_genes)} genomic = {X_combined.shape[1]} total")

# --- Step 3: Train a gradient-boosted classifier ---
# Compare three models: clinical only, genomic only, and combined
print("\n--- Step 3: Model Training and Comparison ---")

X_train, X_test, y_train, y_test = train_test_split(X_combined, y, test_size=0.3, random_state=42, stratify=y)

# Model 1: Clinical features only
X_train_clinical = X_train[clinical_df.columns]
X_test_clinical = X_test[clinical_df.columns]

model_clinical = GradientBoostingClassifier(n_estimators=100, max_depth=3, random_state=42)
model_clinical.fit(X_train_clinical, y_train)
y_pred_clinical = model_clinical.predict_proba(X_test_clinical)[:, 1]
roc_clinical = roc_auc_score(y_test, y_pred_clinical)
print(f"Clinical-only model: ROC-AUC = {roc_clinical:.3f}")

# Model 2: Genomic features only
X_train_genomic = X_train[selected_genes]
X_test_genomic = X_test[selected_genes]

model_genomic = GradientBoostingClassifier(n_estimators=100, max_depth=3, random_state=42)
model_genomic.fit(X_train_genomic, y_train)
y_pred_genomic = model_genomic.predict_proba(X_test_genomic)[:, 1]
roc_genomic = roc_auc_score(y_test, y_pred_genomic)
print(f"Genomic-only model: ROC-AUC = {roc_genomic:.3f}")

# Model 3: Combined features
model_combined = GradientBoostingClassifier(n_estimators=100, max_depth=3, random_state=42)
model_combined.fit(X_train, y_train)
y_pred_combined = model_combined.predict_proba(X_test)[:, 1]
roc_combined = roc_auc_score(y_test, y_pred_combined)
print(f"Combined model: ROC-AUC = {roc_combined:.3f}")

print("\n--- Interpretation ---")
print(f"The combined model (ROC-AUC = {roc_combined:.3f}) outperforms both")
print(f"the clinical-only (ROC-AUC = {roc_clinical:.3f}) and genomic-only")
print(f"(ROC-AUC = {roc_genomic:.3f}) models.")
print("This means that genomic and clinical data contain complementary information.")
print("The model learns from both to make better predictions.")

# --- Step 4: Feature importance in the combined model ---
print("\n--- Step 4: Feature Importance in Combined Model ---")
feature_importance_combined = pd.DataFrame({
    'feature': X_combined.columns,
    'importance': model_combined.feature_importances_
}).sort_values('importance', ascending=False)

print("Top 15 most important features:")
print(feature_importance_combined.head(15).to_string(index=False))

print("\nIn plain English:")
print("- The top 10 genes include known breast cancer markers like ESR1 (estrogen receptor)")
print("  and ERBB2 (HER2). The model is learning real biology.")
print("- Clinical features like tumor grade and ER status also rank highly.")
print("- This feature importance plot tells you which variables the model relies on most.")
print("  It's a starting point for biological investigation, not a final answer.")

What did we learn? The combined model outperformed both the clinical-only and genomic-only models. This makes sense: a patient’s outcome depends on both their genes and their clinical presentation. The feature importance plot showed that the model identified known breast cancer markers (ESR1, ERBB2) alongside clinical variables (tumor grade, ER status).

This is the hardest part: avoiding overfitting when p >> n. With 1,000 genes and 500 patients, it’s easy to find spurious correlations. We used feature selection (selecting only 50 genes) and cross-validation to mitigate this. The ‘Incorporating Machine Learning into Established Bioinformatics Workflows’ paper (PMC 8000113) provides practical strategies for this, including external validation on independent datasets.

Recap: What You Learned

Let’s summarize the key takeaways from our tour of ML in biology and bioinformatics:

1. Biological sequences are strings over small alphabets (4 letters for DNA/RNA, 20 for proteins). They need encoding before ML — k-mer counting and learned embeddings are the two main strategies.

2. Sequence classification (e.g., TIS prediction) works well with CNNs. The learned filters often rediscover known biology — a powerful interpretability check.

3. Structure prediction (AlphaFold, ESMFold) is the most famous success story. You can use pretrained models to get structure predictions for your own sequences in minutes, complete with confidence scores (pLDDT).

4. Variant effect prediction (DeepSEA, DeepVariant) helps prioritize which single-nucleotide changes might cause disease. But these models are prioritization tools, not diagnostic tests — correlation is not causation.

5. Integrating genomic data with clinical features improves prediction but requires careful feature selection and validation to avoid overfitting. The combined model often outperforms either data type alone.

Check Your Understanding

Remember: What are the two dominant strategies for encoding biological sequences as features for ML models?

Understand: Why does a CNN typically outperform a random forest on k-mer counts for sequence classification tasks like TIS prediction?

Apply: Given a FASTA file of 1000 promoter sequences and 1000 non-promoter sequences, outline the steps to build a classifier that predicts whether a new sequence is a promoter.

Analyze: A DeepSEA model gives a variant a score of 0.95 for disrupting a transcription factor binding site. What are the limitations of this prediction? Why can’t you conclude the variant causes disease?

Evaluate: Compare the strengths and weaknesses of using a pretrained structure prediction model (like ESMFold) versus running a physics-based simulation (like molecular dynamics) to predict a protein’s structure.

Create: Design a small ML experiment that combines gene expression data from a public repository with a clinical variable (e.g., age or treatment) to predict a binary outcome (e.g., drug response). Describe the dataset, the features, the model, and how you would validate it.

Related Articles

• Part 1: ‘From Data to Decisions’ — Foundational concepts of supervised learning that this article builds on.
• Part 3: ‘Feature Engineering for Messy Data’ — Techniques for handling variable-length and high-dimensional data, directly applicable to sequence encoding.
• Part 5: ‘Interpreting Black-Box Models’ — Methods for understanding what CNNs and gradient-boosted trees learn, relevant to the TIS filter visualization and feature importance plots in this article.
• Part 7: ‘When Correlation Is Not Causation: An Introduction to Causal Inference’ — Relevant to the variant effect prediction section, where DeepSEA scores are correlations that need causal validation.