# ====================================================================== # Transcriptomic Prediction of Bone Metastasis in Breast Cancer: # A Methodological Investigation Revealing Limited Generalizability # of Machine Learning Models # # Analysis code - companion to the manuscript submitted to # Briefings in Bioinformatics # # Author: Shunxin Yao (n12012416@qut.edu.au) # Date: 2026 # R version: 4.6.1; Bioconductor 3.23 # # Description: # This script performs a leakage-free evaluation of four machine # learning models (Random Forest, Elastic Net, XGBoost, L2-regularized # logistic regression) for predicting bone metastasis from microarray # gene expression data (GSE2034). # # All feature selection, class weighting, and model fitting are # performed exclusively within the training set. Model performance # is evaluated on an independent held-out test set using: # - AUC with 2,000-sample bootstrap 95% CI # - PR-AUC (via PRROC) # - Brier score # - Isotonic-regression calibration # - Calibration slope # - Decision curve analysis (DCA) # - SHAP interpretability # # Inputs: # - GSE2034 expression matrix and phenotype (downloaded automatically # from GEO if not cached). # # Outputs (saved to ~/Desktop/BreastCancer_BoneMetastasis_Results/): # - calibration_isotonic.png # - dca_isotonic.png # - results_final.csv # ====================================================================== # ---------------------------------------------------------------------- # 0. Package dependencies # ---------------------------------------------------------------------- # All packages can be installed via CRAN or Bioconductor: # install.packages(c("caret", "pROC", "glmnet", "randomForest", # "xgboost", "boot", "ggplot2", "dcurves", # "kernelshap", "shapviz", "PRROC")) # BiocManager::install(c("GEOquery", "limma", "hgu133a.db", # "AnnotationDbi")) suppressPackageStartupMessages({ library(GEOquery) # GEO data retrieval library(limma) # Differential expression analysis library(caret) # Data partitioning and CV utilities library(pROC) # ROC / AUC computation library(glmnet) # Elastic Net library(randomForest) # Random Forest library(xgboost) # XGBoost library(boot) # Bootstrap confidence intervals library(ggplot2) # Plotting library(dcurves) # Decision curve analysis library(kernelshap) # Kernel SHAP library(shapviz) # SHAP visualization library(PRROC) # Precision-recall AUC library(hgu133a.db) # Affymetrix HG-U133A probe annotation library(AnnotationDbi) # Annotation database interface }) # ---------------------------------------------------------------------- # 1. Data acquisition # ---------------------------------------------------------------------- # If a cached copy exists, load it; otherwise download from GEO. if (!dir.exists("data")) dir.create("data") if (file.exists("data/GSE2034_expr.rds") && file.exists("data/GSE2034_pheno.rds")) { message("Loading cached GSE2034 data.") expr_GSE2034 <- readRDS("data/GSE2034_expr.rds") pheno_GSE2034 <- readRDS("data/GSE2034_pheno.rds") } else { message("Downloading GSE2034 from GEO...") gse <- getGEO("GSE2034", GSEMatrix = TRUE, getGPL = FALSE) expr_GSE2034 <- exprs(gse[[1]]) pheno_GSE2034 <- pData(gse[[1]]) saveRDS(expr_GSE2034, "data/GSE2034_expr.rds") saveRDS(pheno_GSE2034, "data/GSE2034_pheno.rds") } # ---------------------------------------------------------------------- # 2. Outcome definition # ---------------------------------------------------------------------- # Bone metastasis status is stored in the phenotype column # "bone relapses (1=yes, 0=no):ch1". label_col <- grep("bone relapse", colnames(pheno_GSE2034), ignore.case = TRUE, value = TRUE)[1] stopifnot(!is.na(label_col)) y_all <- factor( ifelse(as.integer(as.character(pheno_GSE2034[[label_col]])) == 1, "Yes", "No"), levels = c("Yes", "No") ) message(sprintf("Total samples: %d; positive events: %d", ncol(expr_GSE2034), sum(y_all == "Yes"))) # ---------------------------------------------------------------------- # 3. Train/test split (70/30, stratified) # ---------------------------------------------------------------------- set.seed(123) train_idx <- createDataPartition(y_all, p = 0.7, list = FALSE) X_train_raw <- expr_GSE2034[, train_idx] X_test_raw <- expr_GSE2034[, -train_idx] y_train <- y_all[ train_idx] y_test <- y_all[-train_idx] message(sprintf("Training: %d samples, %d positive events", ncol(X_train_raw), sum(y_train == "Yes"))) message(sprintf("Test: %d samples, %d positive events", ncol(X_test_raw), sum(y_test == "Yes"))) # ---------------------------------------------------------------------- # 4. Feature selection (training set only) # ---------------------------------------------------------------------- # Step 1: variance filter - retain the 3,000 highest-variance probes. var_tr <- apply(X_train_raw, 1, var, na.rm = TRUE) var_tr <- var_tr[!is.na(var_tr)] top_p <- names(sort(var_tr, decreasing = TRUE))[1:min(3000, length(var_tr))] X_train_raw <- X_train_raw[top_p, , drop = FALSE] X_test_raw <- X_test_raw[ top_p, , drop = FALSE] # Step 2: differential expression via limma with BH FDR correction. grp <- factor(y_train, levels = c("No", "Yes")) design <- model.matrix(~ 0 + grp) colnames(design) <- levels(grp) fit <- lmFit(X_train_raw, design) cont <- makeContrasts(Yes - No, levels = design) fit2 <- eBayes(contrasts.fit(fit, cont)) de_res <- topTable(fit2, number = Inf, adjust.method = "BH", sort.by = "P") sel <- rownames(de_res)[de_res$adj.P.Val < 0.1 & abs(de_res$logFC) > 0.3] if (length(sel) > 50) sel <- rownames(de_res)[1:50] if (length(sel) < 10) sel <- rownames(de_res)[1:30] message(sprintf("Selected %d probes for modelling.", length(sel))) # ---------------------------------------------------------------------- # 5. Assemble final matrices (samples x probes) # ---------------------------------------------------------------------- mk_df <- function(X, y) { df <- data.frame(t(X[sel, , drop = FALSE]), check.names = FALSE) colnames(df) <- sub("^X", "", colnames(df)) df$labels <- y df } train_df <- mk_df(X_train_raw, y_train) test_df <- mk_df(X_test_raw, y_test) X_tr <- as.matrix(train_df[, -ncol(train_df)]) X_te <- as.matrix(test_df[, -ncol(test_df)]) # Remove zero-variance columns nzv <- apply(X_tr, 2, function(x) sd(x, na.rm = TRUE) == 0) if (any(nzv)) { X_tr <- X_tr[, !nzv, drop = FALSE] X_te <- X_te[, !nzv, drop = FALSE] } # Class weight for the positive class pos_w <- sum(y_train == "No") / sum(y_train == "Yes") message(sprintf("Positive class weight: %.3f", pos_w)) # ---------------------------------------------------------------------- # 6. Generic 5-fold CV training function # ---------------------------------------------------------------------- train_one_model <- function(X, y, fold_list, method, pos_w) { n <- nrow(X) oof_prob <- rep(NA, n) for (k in seq_along(fold_list)) { idx_val <- fold_list[[k]] idx_trn <- setdiff(seq_len(n), idx_val) Xtr <- X[idx_trn, , drop = FALSE] Xva <- X[idx_val, , drop = FALSE] ytr <- y[idx_trn] ytr_num <- ifelse(ytr == "Yes", 1, 0) w_trn <- ifelse(ytr == "Yes", pos_w, 1) if (method == "rf") { m <- randomForest::randomForest( x = Xtr, y = factor(ytr, levels = c("Yes", "No")), ntree = 500, classwt = c("Yes" = pos_w, "No" = 1)) oof_prob[idx_val] <- predict(m, Xva, type = "prob")[, "Yes"] } else if (method == "glmnet") { m <- glmnet::cv.glmnet(Xtr, ytr_num, family = "binomial", alpha = 0.5, weights = w_trn, nfolds = 5) oof_prob[idx_val] <- as.numeric( predict(m, Xva, type = "response", s = "lambda.min")) } else if (method == "xgb") { dtrain <- xgboost::xgb.DMatrix(Xtr, label = ytr_num, weight = w_trn) dval <- xgboost::xgb.DMatrix(Xva) m <- xgboost::xgb.train( params = list(objective = "binary:logistic", eval_metric = "auc", max_depth = 3, eta = 0.1, subsample = 0.8, colsample_bytree = 0.8), data = dtrain, nrounds = 100, verbose = 0) oof_prob[idx_val] <- predict(m, dval) } else if (method == "lr") { df_trn <- data.frame(Xtr, y = ytr_num) m <- glm(y ~ ., data = df_trn, family = quasibinomial(), weights = w_trn) oof_prob[idx_val] <- predict(m, data.frame(Xva), type = "response") } } oof_prob } # ---------------------------------------------------------------------- # 7. Five-fold CV - out-of-fold predictions # ---------------------------------------------------------------------- set.seed(123) folds <- createFolds(y_train, k = 5, list = TRUE) message("Running 5-fold CV for four models...") oof_rf <- train_one_model(X_tr, y_train, folds, "rf", pos_w) oof_glm <- train_one_model(X_tr, y_train, folds, "glmnet", pos_w) oof_xgb <- train_one_model(X_tr, y_train, folds, "xgb", pos_w) oof_lr <- tryCatch( train_one_model(X_tr, y_train, folds, "lr", pos_w), error = function(e) rep(NA, nrow(X_tr))) # ---------------------------------------------------------------------- # 8. Cross-validated AUC (out-of-fold) # ---------------------------------------------------------------------- cv_eval <- function(p, name) { if (all(is.na(p))) { message(sprintf("%-10s failed", name)); return(NA) } r <- pROC::roc(as.numeric(y_train) - 1, p, quiet = TRUE) a <- as.numeric(pROC::auc(r)) message(sprintf("%-10s CV-OOF AUC = %.4f", name, a)) a } cv_rf <- cv_eval(oof_rf, "RF") cv_glm <- cv_eval(oof_glm, "GLMNET") cv_xgb <- cv_eval(oof_xgb, "XGBoost") cv_lr <- cv_eval(oof_lr, "LR") # ---------------------------------------------------------------------- # 9. Isotonic-regression calibration # ---------------------------------------------------------------------- # Isotonic regression is fitted on out-of-fold predictions and applied # to the test set. This avoids the degenerate behaviour observed with # Platt scaling under class weighting. calibrate_iso <- function(oof_prob, y_train, test_prob) { y_num <- ifelse(y_train == "Yes", 1, 0) iso <- isoreg(oof_prob, y_num) f <- approxfun(iso$x, iso$yf, rule = 2) p <- pmin(pmax(test_prob, min(iso$x)), max(iso$x)) f(p) } prob_rf_iso <- calibrate_iso(oof_rf, y_train, NULL) # placeholder # Fit full models below for test-set prediction # ---------------------------------------------------------------------- # 10. Train full models, predict test set, and calibrate # ---------------------------------------------------------------------- predict_full <- function(X_tr, y_train, X_te, method, pos_w, oof_prob) { y_num <- ifelse(y_train == "Yes", 1, 0) w_trn <- ifelse(y_train == "Yes", pos_w, 1) if (method == "rf") { m <- randomForest::randomForest( x = X_tr, y = factor(y_train, levels = c("Yes", "No")), ntree = 500, classwt = c("Yes" = pos_w, "No" = 1)) raw <- predict(m, X_te, type = "prob")[, "Yes"] } else if (method == "glmnet") { m <- glmnet::cv.glmnet(X_tr, y_num, family = "binomial", alpha = 0.5, weights = w_trn, nfolds = 5) raw <- as.numeric(predict(m, X_te, type = "response", s = "lambda.min")) } else if (method == "xgb") { dtrain <- xgboost::xgb.DMatrix(X_tr, label = y_num, weight = w_trn) dtest <- xgboost::xgb.DMatrix(X_te) m <- xgboost::xgb.train( params = list(objective = "binary:logistic", eval_metric = "auc", max_depth = 3, eta = 0.1, subsample = 0.8, colsample_bytree = 0.8), data = dtrain, nrounds = 100, verbose = 0) raw <- predict(m, dtest) } else if (method == "lr") { df_trn <- data.frame(X_tr, y = y_num) m <- glm(y ~ ., data = df_trn, family = quasibinomial(), weights = w_trn) raw <- predict(m, data.frame(X_te), type = "response") } # Isotonic calibration using OOF predictions y_num_train <- ifelse(y_train == "Yes", 1, 0) iso <- isoreg(oof_prob, y_num_train) f <- approxfun(iso$x, iso$yf, rule = 2) p <- pmin(pmax(raw, min(iso$x)), max(iso$x)) cal <- f(p) list(raw = raw, cal = cal) } message("Training full models and predicting test set...") pred_rf <- predict_full(X_tr, y_train, X_te, "rf", pos_w, oof_rf) pred_glm <- predict_full(X_tr, y_train, X_te, "glmnet", pos_w, oof_glm) pred_xgb <- predict_full(X_tr, y_train, X_te, "xgb", pos_w, oof_xgb) pred_lr <- predict_full(X_tr, y_train, X_te, "lr", pos_w, oof_lr) # ---------------------------------------------------------------------- # 11. Test-set evaluation # ---------------------------------------------------------------------- calc_pr_auc <- function(prob, y_test) { y_num <- as.numeric(y_test) - 1 pr <- PRROC::pr.curve(scores.class0 = prob[y_num == 1], scores.class1 = prob[y_num == 0], curve = FALSE) pr$auc.integral } test_eval <- function(prob, y_test, name) { y_num <- as.numeric(y_test) - 1 r <- pROC::roc(y_num, prob, quiet = TRUE) auc_val <- as.numeric(pROC::auc(r)) pr_auc <- calc_pr_auc(prob, y_test) brier <- mean((prob - y_num)^2) cal_model <- glm(y_num ~ qlogis(pmin(pmax(prob, 1e-6), 1 - 1e-6)), family = binomial) slope <- coef(cal_model)[2] message(sprintf("%-10s AUC=%.4f PR-AUC=%.4f Brier=%.4f Cal.Slope=%.4f", name, auc_val, pr_auc, brier, slope)) list(auc = auc_val, pr = pr_auc, brier = brier, slope = slope, roc = r, prob = prob) } message("Evaluating on independent test set...") res_rf <- test_eval(pred_rf$cal, y_test, "RF") res_glm <- test_eval(pred_glm$cal, y_test, "GLMNET") res_xgb <- test_eval(pred_xgb$cal, y_test, "XGBoost") res_lr <- test_eval(pred_lr$cal, y_test, "LR") # ---------------------------------------------------------------------- # 12. Bootstrap 95% CI for the best model # ---------------------------------------------------------------------- all_aucs <- c(RF = res_rf$auc, GLMNET = res_glm$auc, XGB = res_xgb$auc, LR = res_lr$auc) best_name <- names(which.max(all_aucs)) best_prob <- switch(best_name, RF = pred_rf$cal, GLMNET = pred_glm$cal, XGB = pred_xgb$cal, LR = pred_lr$cal) message(sprintf("Best model: %s (AUC = %.4f)", best_name, max(all_aucs))) set.seed(123) bd <- data.frame(obs = as.numeric(y_test) - 1, pred = best_prob) bo <- boot::boot(bd, function(d, i) { as.numeric(pROC::auc(pROC::roc(d$obs[i], d$pred[i], quiet = TRUE))) }, R = 2000) ci <- boot::boot.ci(bo, type = "perc") message(sprintf("AUC 95%% CI: %.4f - %.4f", ci$percent[4], ci$percent[5])) # ---------------------------------------------------------------------- # 13. DeLong tests (pairwise) # ---------------------------------------------------------------------- message("DeLong tests:") message(sprintf(" XGBoost vs LR: p = %.4f", pROC::roc.test(res_xgb$roc, res_lr$roc)$p.value)) message(sprintf(" XGBoost vs GLMNET: p = %.4f", pROC::roc.test(res_xgb$roc, res_glm$roc)$p.value)) message(sprintf(" XGBoost vs RF: p = %.4f", pROC::roc.test(res_xgb$roc, res_rf$roc)$p.value)) # ---------------------------------------------------------------------- # 14. Calibration curve # ---------------------------------------------------------------------- cal_df <- data.frame(pred = best_prob, obs = as.numeric(y_test) - 1) cal_plot <- ggplot(cal_df, aes(x = pred, y = obs)) + geom_smooth(method = "loess", se = TRUE, col = "blue") + geom_abline(slope = 1, intercept = 0, linetype = "dashed", col = "red") + xlim(0, 1) + ylim(0, 1) + labs(x = "Predicted Probability", y = "Observed Proportion", title = paste0("Calibration (isotonic, ", best_name, ")")) ggsave("calibration_isotonic.png", cal_plot, width = 5, height = 5) # ---------------------------------------------------------------------- # 15. Decision curve analysis # ---------------------------------------------------------------------- dca_data <- data.frame(risk = best_prob, outcome = as.numeric(y_test) - 1) dca_obj <- dcurves::dca(outcome ~ risk, data = dca_data, thresholds = seq(0, 0.6, by = 0.01)) png("dca_isotonic.png", width = 600, height = 500) plot(dca_obj, smooth = TRUE) dev.off() # ---------------------------------------------------------------------- # 16. SHAP interpretability (Random Forest) # ---------------------------------------------------------------------- message("Computing SHAP values...") set.seed(123) rf_full <- randomForest::randomForest( x = X_tr, y = factor(y_train, levels = c("Yes", "No")), ntree = 500, classwt = c("Yes" = pos_w, "No" = 1)) bg_idx <- sample(nrow(X_tr), min(100, nrow(X_tr))) pred_fun <- function(object, X) predict(object, X, type = "prob")[, "Yes"] ks <- kernelshap::kernelshap(rf_full, X = X_tr, pred_fun = pred_fun, bg_X = X_tr[bg_idx, , drop = FALSE]) shp <- shapviz::shapviz(ks) imp <- data.frame(Probe = colnames(X_tr), MeanAbsSHAP = colMeans(abs(shp$S))) imp <- imp[order(-imp$MeanAbsSHAP), ] imp$Gene <- AnnotationDbi::mapIds(hgu133a.db, keys = imp$Probe, column = "SYMBOL", keytype = "PROBEID", multiVals = "first") png("shap_final.png", width = 800, height = 600) shapviz::sv_importance(shp, kind = "both", show_numbers = TRUE, max_display = 20) dev.off() # ---------------------------------------------------------------------- # 17. Save all results to a desktop folder # ---------------------------------------------------------------------- out_dir <- file.path("~", "Desktop", "BreastCancer_BoneMetastasis_Results") if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE) results_final <- data.frame( Model = c("RF", "GLMNET", "XGBoost", "LR"), CV_OOF_AUC = c(cv_rf, cv_glm, cv_xgb, cv_lr), Test_AUC = c(res_rf$auc, res_glm$auc, res_xgb$auc, res_lr$auc), Test_PR_AUC = c(res_rf$pr, res_glm$pr, res_xgb$pr, res_lr$pr), Test_Brier = c(res_rf$brier, res_glm$brier, res_xgb$brier, res_lr$brier), Cal_Slope = c(res_rf$slope, res_glm$slope, res_xgb$slope, res_lr$slope) ) write.csv(results_final, file.path(out_dir, "results_final.csv"), row.names = FALSE) file.copy("calibration_isotonic.png", out_dir, overwrite = TRUE) file.copy("dca_isotonic.png", out_dir, overwrite = TRUE) file.copy("shap_final.png", out_dir, overwrite = TRUE) message("All outputs saved to: ", out_dir) message("Analysis complete.")