Analyzing Gene Expression Data with Machine Learning

After Completing This Topic

You'll be able to cluster gene expression data with K-means using scikit-learn, reduce dimensions with PCA, and classify tumor/normal samples.

What Is Machine Learning?

With lab experience, you can look at a gel photo and say "this band pattern is a mutation," right? Machine learning is teaching a computer this kind of "pattern recognition" ability through data.

ML commonly used in bio:

Unsupervised learning: Find patterns in data without labels (clustering, PCA)
Supervised learning: Train a classification model with labeled data (tumor vs normal)

Data Preparation: Simulated Gene Expression Data

We'll create simulated data mimicking actual RNA-seq data. The scenario: measuring expression of 20 genes in 100 samples (50 normal + 50 tumor).

python

import numpy as np
import pandas as pd

np.random.seed(42)

n_samples = 100
n_genes = 20
gene_names = [f"Gene_{i+1:02d}" for i in range(n_genes)]
sample_labels = ["Normal"] * 50 + ["Tumor"] * 50

normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0   # Gene_01~05 upregulated
tumor_expr[:, 15:] -= 2.0  # Gene_16~20 downregulated

expression = np.vstack([normal_expr, tumor_expr])
df = pd.DataFrame(expression, columns=gene_names)
df["Label"] = sample_labels

print(f"Data shape: {df.shape}")
print(f"Samples: Normal {sum(df['Label']=='Normal')}, Tumor {sum(df['Label']=='Tumor')}")
print(f"\nFirst 5 rows:")
print(df.head().to_string())

assert df.shape == (100, 21)
assert sum(df["Label"] == "Normal") == 50
assert sum(df["Label"] == "Tumor") == 50

PCA: Compressing High-Dimensional Data to 2D

We compress expression data from 20 genes into 2 dimensions that humans can visualize. PCA (Principal Component Analysis) finds the directions of greatest variance in the data.

python

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

np.random.seed(42)
n_genes = 20
gene_names = [f"Gene_{i+1:02d}" for i in range(n_genes)]

normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0

expression = np.vstack([normal_expr, tumor_expr])
labels = ["Normal"] * 50 + ["Tumor"] * 50

# 1. Standardize (mean 0, variance 1)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(expression)

# 2. Apply PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

print(f"Original dimensions: {expression.shape[1]} genes")
print(f"Reduced dimensions: {X_pca.shape[1]} principal components")
print(f"Explained variance: PC1={pca.explained_variance_ratio_[0]:.1%}, PC2={pca.explained_variance_ratio_[1]:.1%}")
print(f"Total: {sum(pca.explained_variance_ratio_):.1%}")

assert X_pca.shape == (100, 2)
assert pca.explained_variance_ratio_[0] > pca.explained_variance_ratio_[1]

python

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

np.random.seed(42)
n_genes = 20

normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0

expression = np.vstack([normal_expr, tumor_expr])
labels = np.array(["Normal"] * 50 + ["Tumor"] * 50)

scaler = StandardScaler()
X_scaled = scaler.fit_transform(expression)
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

fig, ax = plt.subplots(figsize=(8, 6))
for label, color in [("Normal", "#2E86AB"), ("Tumor", "#E74C3C")]:
    mask = labels == label
    ax.scatter(X_pca[mask, 0], X_pca[mask, 1], c=color, s=40, alpha=0.7, label=label, edgecolors="white", linewidth=0.5)

ax.set_xlabel(f"PC1 ({pca.explained_variance_ratio_[0]:.1%})", fontsize=12)
ax.set_ylabel(f"PC2 ({pca.explained_variance_ratio_[1]:.1%})", fontsize=12)
ax.set_title("PCA — Normal vs Tumor", fontsize=14, fontweight="bold")
ax.legend(fontsize=11)
ax.grid(True, alpha=0.2)
fig.tight_layout()
fig.savefig("pca_plot.png", dpi=150)
plt.close(fig)

print("pca_plot.png saved")
assert X_pca.shape == (100, 2)

K-means Clustering: Unsupervised Learning

Finding groups from data alone, without labels. K-means is the most intuitive clustering algorithm.

python

import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans

np.random.seed(42)
n_genes = 20

normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0

expression = np.vstack([normal_expr, tumor_expr])
true_labels = [0] * 50 + [1] * 50

scaler = StandardScaler()
X_scaled = scaler.fit_transform(expression)

# K-means: classify into 2 clusters
kmeans = KMeans(n_clusters=2, random_state=42, n_init=10)
clusters = kmeans.fit_predict(X_scaled)

# Check how well clusters match actual labels
from sklearn.metrics import adjusted_rand_score
ari = adjusted_rand_score(true_labels, clusters)
print(f"Cluster 0: {sum(clusters == 0)} samples")
print(f"Cluster 1: {sum(clusters == 1)} samples")
print(f"Adjusted Rand Index: {ari:.3f} (1.0 = perfect match)")

assert len(set(clusters)) == 2
assert ari > 0.5
print("K-means clustering successfully distinguished normal/tumor")

Classification: Supervised Learning

This time, we train a model with labeled data (Normal/Tumor) and classify new samples.

python

import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, classification_report

np.random.seed(42)
n_genes = 20

normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0

X = np.vstack([normal_expr, tumor_expr])
y = np.array([0] * 50 + [1] * 50)  # 0=Normal, 1=Tumor

# 1. Split data (80% train, 20% test)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)

# 2. Standardize
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 3. Train Random Forest
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train_scaled, y_train)

# 4. Predict & evaluate
y_pred = clf.predict(X_test_scaled)
accuracy = accuracy_score(y_test, y_pred)

print(f"Accuracy: {accuracy:.1%}")
print(f"\nTraining data: {len(X_train)} samples, Test data: {len(X_test)} samples")
print(f"\nClassification Report:")
print(classification_report(y_test, y_pred, target_names=["Normal", "Tumor"]))

assert accuracy > 0.8
assert len(X_train) == 80
assert len(X_test) == 20

Feature Importance: Which Genes Matter?

Random Forest reveals which genes (features) contributed most to the classification.

python

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier

np.random.seed(42)
n_genes = 20
gene_names = [f"Gene_{i+1:02d}" for i in range(n_genes)]

normal_expr = np.random.randn(50, n_genes) * 1.0 + 5.0
tumor_expr = np.random.randn(50, n_genes) * 1.5 + 5.0
tumor_expr[:, :5] += 3.0
tumor_expr[:, 15:] -= 2.0

X = np.vstack([normal_expr, tumor_expr])
y = np.array([0] * 50 + [1] * 50)

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_scaled, y)

importances = clf.feature_importances_
sorted_idx = np.argsort(importances)[-10:]  # Top 10

fig, ax = plt.subplots(figsize=(8, 5))
ax.barh(range(len(sorted_idx)), importances[sorted_idx], color="#2E86AB")
ax.set_yticks(range(len(sorted_idx)))
ax.set_yticklabels([gene_names[i] for i in sorted_idx], fontsize=10)
ax.set_xlabel("Feature Importance", fontsize=12)
ax.set_title("Top 10 Important Genes", fontsize=14, fontweight="bold")
fig.tight_layout()
fig.savefig("feature_importance.png", dpi=150)
plt.close(fig)

top_gene = gene_names[sorted_idx[-1]]
print(f"Most important gene: {top_gene}")
print(f"Top 5 genes: {[gene_names[i] for i in sorted_idx[-5:]]}")

assert len(importances) == 20
assert sum(importances) > 0.99

Try It Yourself (Faded Example)

Fill in the blanks to complete scikit-learn's fit-predict pattern.

Fill in the Blankspython

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = scaler.fit_transform(data)
kmeans = KMeans(n_clusters=)
clusters = kmeans.(X_scaled)

Common Errors & Solutions

Q: ValueError: could not convert string to float

String columns (e.g., gene names, labels) are included in the data. Use df.select_dtypes(include=[np.number]) to select only numeric columns.

Q: Accuracy is only 50%

This means the data is no different from random. Check whether you standardized (StandardScaler) and whether your features (genes) have sufficiently different patterns.

Q: ConvergenceWarning: Number of distinct clusters found

K-means failed to converge. Try increasing n_init=10 or max_iter=300, or check your data scaling.

Congratulations! You've completed all 5 MVP topics in the Common + A-1 tracks of DevBench.

BioPlayground

Analyzing Gene Expression Data with Machine Learning

Analyzing Gene Expression Data with Machine Learning

After Completing This Topic

What Is Machine Learning?

Data Preparation: Simulated Gene Expression Data

PCA: Compressing High-Dimensional Data to 2D

K-means Clustering: Unsupervised Learning

Classification: Supervised Learning

Feature Importance: Which Genes Matter?

Try It Yourself (Faded Example)

Common Errors & Solutions

💬 Questions & Comments