From Raw Data to Clinical Predictions: A Bioinformatics Pipeline for Breast Cancer Bone Relapse Using Python and Machine Learning

Every year, thousands of breast cancer patients undergo surgery with clear lymph nodes seemingly cancer-free. Yet months or years later, some develop bone metastasis that could have been anticipated from the molecular profile of their original tumor. In this article, I’ll show how I built a complete bioinformatics pipeline from downloading raw gene expression data to deploying a live web app that predicts bone relapse risk using 35 genes identified through differential expression analysis.

GSE2034 is one of the most cited breast cancer gene expression datasets in bioinformatics. Published in The Lancet in 2005, it contains Affymetrix microarray data from 286 node-negative breast cancer patients followed for up to 10 years. The key clinical variable: bone relapse (69 patients developed it, 217 did not a 24.1% rate). The biological question: can the gene expression profile of the primary tumor at surgery predict which patients will develop bone relapse years later? I downloaded the dataset directly from NCBI GEO using Python and processed the raw expression matrix through a complete quality control pipeline.

Dataset: GSE2034 (Wang et al., The Lancet, 2005) Patients: 286 | Tool: Python | Live app: Bone-Relapse-Predictor

The Pipeline: 5 Steps

Step 1 — Quality Control

After downloading GSE2034_series_matrix.txt.gz from NCBI, I applied log2 transformation to normalize the expression values (raw values ranged from 0 to 65,000; after log2, the distribution centered around 6–8, as expected for Affymetrix microarray data). I filtered genes with low variability (bottom 25% by standard deviation), reducing 22,283 genes to 16,712 informative features.

Figure 1. Expression distribution across 30 samples from GSE2034 (286 breast cancer patients). Boxes are aligned at similar heights confirming uniform sample quality. No outliers were detected across the full 286-sample cohort. After filtering low-variability genes (bottom 25% by standard deviation), 16,712 informative genes were retained from the original 22,283.

Step 2 — Differential Expression Analysis

I ran a t-test for each of the 16,712 genes comparing the 69 relapse patients vs 217 non-relapse patients, followed by Benjamini-Hochberg FDR correction.

Criteria: FDR < 0.05 AND |log2FC| > 1 Result: 35 differentially expressed genes

• 7 upregulated in bone relapse
• 28 downregulated in bone relapse

The dominance of downregulated genes suggests that bone relapse is associated with silencing of tumor suppressors rather than activation of new oncogenes.

Figure 2. Volcano plot of differential expression analysis between bone relapse (n=69) and no relapse (n=217) patients. Each point represents one of the 16,712 genes tested. Red points (right) are upregulated in bone relapse. Blue points (left) are downregulated. Dashed lines indicate significance thresholds: FDR < 0.05 (horizontal) and |log2FC| > 1 (vertical). 35 genes passed both criteria — 7 up and 28 down — suggesting bone relapse is predominantly associated with gene silencing rather than activation.

Step 3 — Unsupervised Analysis

PCA and K-means clustering (K=3, silhouette=0.165) revealed 3 molecular subgroups with dramatically different relapse rates:

Cluster 1 (n=115): 38.3% relapse — HIGH RISK Cluster 2 (n=57): 8.8% relapse — LOW RISK Cluster 3 (n=114): 17.5% relapse — INTERMEDIATE

This unsupervised result validates that the molecular signal is real. The algorithm found clinically meaningful groups without seeing any clinical labels.

Figure 3. K-means clustering (K=3) of 286 patients projected onto the first two principal components. Colors represent molecular clusters identified without using clinical labels. Circles = patients with bone relapse. Triangles = patients without bone relapse. Despite visual overlap in 2D (PC1+PC2 explain only 12.5% of variance), the three clusters show dramatically different relapse rates: Cluster 1 (38.3%), Cluster 2 (8.8%), and Cluster 3 (17.5%), compared to the dataset average of 24.1%.

Step 4 — Machine Learning Classification

Using the 35 DEGs as features, I trained and compared Logistic Regression, Random Forest, and SVM. Logistic Regression performed best (AUC = 0.697, 5-fold CV), consistent with the literature. Linear models often outperform complex ones when features are carefully pre-selected.

I optimized the classification threshold from 0.50 to 0.30 to maximize recall (detecting actual relapses), increasing detected cases from 4/14 to 6/14 in the test set.

Step 5 — Explainability with SHAP

SHAP values reveal which genes drive each individual prediction, making the model interpretable for clinical researchers.

Figure 4. Mean absolute SHAP values for the 35 DEGs used as model features. Bar length indicates the average impact of each gene on the model’s predictions across all 286 patients. Longer bars represent genes that consistently drive predictions toward high or low risk. SHAP values allow researchers to understand which genes are most relevant to the model’s decision — a critical requirement for clinical interpretability.

From Model to Product

The complete pipeline is deployed as a live web application built with Streamlit. Users can:

• Upload any CSV with Affymetrix HG-U133A gene expression data
• Get bone relapse probability per patient
• See risk classification with optimized threshold
• Explore SHAP explanations gene by gene

Live app: Bone-Relapse-Predictor Code: Bone-relapse-predictor Code

This type of tool represents what bioinformatics is moving toward not just analysis notebooks, but deployable products that researchers and clinicians can use without programming.

Limitations and Next Steps

This model has real limitations that are important to acknowledge:

• Small dataset (286 patients) limits statistical power
• AUC of 0.67 is meaningful but not sufficient for clinical use
• Microarray technology (2005) has lower dynamic range than modern RNA-seq
• External validation on independent cohorts is needed

Next steps I’m working on:

• XGBoost with SMOTE to handle class imbalance
• Testing on additional GEO datasets for validation
• Integrating pathway analysis to interpret the 35 DEGs biologically

The goal is not to replace clinical judgment but to provide an additional layer of molecular information that could help stratify patients for intensive follow-up.

Conclusion

Building this pipeline taught me that the gap between a bioinformatics analysis and a clinical tool is mostly about three things: reproducibility, interpretability, and deployment. The notebooks are reproducible. SHAP makes it interpretable. Streamlit makes it deployable. If you’re working in cancer genomics, computational biology, or clinical AI and want to discuss this work or collaborate, feel free to connect.

References

Wang Y, et al. Gene-expression profiles to predict distant metastasis of lymph-node-negative primary breast cancer. The Lancet. 2005;365(9460):671–679.