# Integrate RNA-seq Microarray Machine Learning: Step-by-Step Guide By the end of this guide you will be able to take a microarray dataset and an RNA-seq dataset for the same disease, harmonize them into a single analysis-ready matrix, select candidate biomarker genes with statistical and network methods, train machine learning classifiers on the combined data, and validate the resulting signature on an independent cohort from a different platform. Every step includes the actual R commands, the parameters that matter, the errors you are most likely to hit, and how to confirm each stage worked before moving on. You need R (4.2 or newer) and RStudio installed, about 16 GB of RAM for typical GEO-scale datasets, and access to the internet for downloading data. The packages you will use are GEOquery, limma, edgeR or DESeq2, sva, WGCNA, glmnet, randomForest, e1071, xgboost, Boruta, pROC, and ggplot2. This article uses human disease examples because that is where the published cross-platform biomarker literature is densest, but the same pipeline works for mouse, rat, and plant [transcriptomics](/knowledge/bioinformatics/modern-transcriptomics-bulk-single-cell-spatial) with minor identifier changes. ## Why Integrate RNA-seq and Microarray Data With Machine Learning RNA-seq ([RNA sequencing](/blog/guides/rna-sequencing)) measures transcript abundance by counting sequencing reads. Microarray measures it by hybridizing labeled cDNA to fixed probes on a chip. The two technologies produce different scales, different dynamic ranges, and different systematic biases, which is why naïve merging fails. When the platforms are harmonized correctly, integration gives you three things a single platform cannot provide. First, sample size. A disease with 30 microarray samples in one study and 40 RNA-seq samples in another gives you a discovery cohort of 70 when integrated, which materially improves machine learning stability. Second, platform independence as evidence. A gene signature that classifies disease when trained on microarray data and still performs when applied to RNA-seq data is far more credible than one that only works on its training platform. Published work in mucoepidermoid carcinoma, dilated cardiomyopathy, and glioblastoma all followed this logic, training on microarray cohorts and validating on independent RNA-seq cohorts [1][2][3]. Third, legacy data. Microarray studies are steadily becoming a minority of new publications, but decades of well-annotated microarray cohorts remain in GEO (the [Gene Expression Omnibus](/blog/guides/gene-expression-omnibus), NCBI's public repository of [gene expression](/blog/guides/gene-expression) data). A recent review of diabetes gene expression studies found microarrays were still the most frequent data acquisition technique, appearing in 70 of 96 included studies, while LASSO (43 uses) and PCA (47 uses) dominated the machine learning methods [4]. Integrating microarray archives with new RNA-seq data lets you reuse that investment. The core technical challenge is real, not cosmetic. Differences in measurement scales, probe chemistries, and signal distributions create systematic discrepancies that make naive large-scale integration nearly impossible. Probabilistic cross-platform translation frameworks such as PXN have been developed specifically to bridge the microarray and RNA-seq divide, and they outperform simpler normalization approaches in cross-platform benchmarks [5]. This guide covers the classical harmonization route in detail and notes where learned translation tools fit. ## What You Need Installed Run this block once. It installs every package used in this guide. ```r install.packages("BiocManager") BiocManager::install(c( "GEOquery", "limma", "edgeR", "DESeq2", "sva", "WGCNA", "glmnet", "randomForest", "e1071", "xgboost", "Boruta", "pROC", "caret", "ggplot2", "ComplexHeatmap" )) install.packages(c("pROC", "caret", "ggplot2", "xgboost")) ``` Set a working directory that has enough disk space, because raw GEO downloads for a single series can exceed 1 GB. ```r setwd("~/projects/cross_platform_ml") dir.create("data", showWarnings = FALSE) dir.create("results", showWarnings = FALSE) ``` ## Step 1: Download Data From GEO Identify two cohorts measuring the same phenotype on different platforms. A useful pattern, drawn from published integrative studies, is a microarray discovery cohort plus an RNA-seq validation cohort. For example, a polycystic ovary syndrome study integrated two microarray datasets (GSE34526 and GSE137684) as the discovery cohort and validated findings in the RNA-seq dataset GSE168404 [6]. A dilated cardiomyopathy study used GSE57338 as the discovery cohort, three microarray series (GSE26887, GSE42955, GSE79962) as external microarray validation, and GSE116250 for independent RNA-seq validation [2]. Decide this architecture before you write any code. Download and load a microarray series with GEOquery: ```r library(GEOquery) gse <- getGEO("GSE34526", GSEMatrix = TRUE, destdir = "data") eset <- gse[[7]] expr_raw <- exprs(eset) dim(expr_raw) ``` Expected output shape (your numbers will differ): ``` dim(expr_raw) [7] 54675 10 ``` That output means 54,675 probe rows and 10 samples. Check the sample metadata immediately, because most integration mistakes trace back to mismatched phenotype labels: ```r pData(eset)[, c("title", "characteristics_ch1")] table(pData(eset)$`group:ch1`) ``` For RNA-seq data deposited as a counted matrix rather than raw reads, download the supplementary file directly: ```r getGEOSuppFiles("GSE168404", base_dir = "data") counts <- read.delim("data/GSE168404/GSE168404_counts.txt", row.names = 1) ``` If the series only has raw FASTQ files, you must align and count them yourself with a standard pipeline (STAR for alignment, then featureCounts, or a quasi-aligner such as Salmon) before normalization [8]. Budget several hours of compute per sample for alignment on a workstation. ## Step 2: Normalize Each Platform Separately This is the step people skip and regret. Normalize each dataset within its own platform first. Cross-platform correction only works after each platform's internal preprocessing is sound. ### Microarray preprocessing with limma Most Affymetrix and Illumina microarray data on GEO is already background-corrected and quantile-normalized when you receive it as a processed series matrix, but you must confirm this from the series header. If you downloaded raw CEL files instead, process them with the oligo or affy packages using the RMA algorithm (Robust Multi-array Average, which performs background correction, quantile normalization, and log2 transformation). For already-processed matrices, log2-transform if the values are on the linear scale: ```r library(limma) # Log2 transform if the data is not already log-scale if (max(expr_raw, na.rm = TRUE) > 100) { expr_ma <- log2(expr_raw + 1) } else { expr_ma <- expr_raw } # Quantile normalization across samples (harmless if already done, # mandatory if arrays were processed independently) expr_ma <- normalizeBetweenArrays(expr_ma, method = "quantile") boxplot(expr_ma, outline = FALSE) ``` How to check the result: every sample's boxplot should sit at roughly the same median, and the value distribution should be centered near a log2 scale (most values between roughly 4 and 16 for typical expression arrays). ### RNA-seq preprocessing with DESeq2 ![Worked example of the DESeq2 filter that keeps genes with at least 10 counts in at least 3 samples, plus the 5-sample version.](/article-images/illustrations/integrate-rna-seq-microarray-machine-learning-vst-filter-threshold.webp) For raw counts, use variance-stabilizing transformation (VST), which produces log-like values with variance stabilized across the count range. This is the transformation shown to work well for microarray and RNA-seq harmonization in a published cross-platform plant stress pipeline [9], and it is the standard recommendation from the DESeq2 documentation because VST output lives on a scale comparable to log2 microarray intensities. ```r library(DESeq2) count_mat <- as.matrix(counts) group <- factor(c("case", "case", "control", "control")) # match your metadata dds <- DESeq2::DESeqDataSetFromMatrix( countData = round(count_mat), colData = data.frame(row.names = colnames(count_mat), group = group), design = ~ group ) # Filter very low-count genes: this reduces both noise and runtime keep <- rowSums(DESeq2::counts(dds) >= 10) >= 3 dds <- dds[keep, ] dds <- DESeq(dds) vsd <- vst(dds, blind = FALSE) expr_seq <- assay(vsd) dim(expr_seq) ``` The `blind = FALSE` argument keeps the design in the transformation, which you want here because you are not using the transform for quality control discovery. Filtering genes with fewer than 10 counts in at least 3 samples is a common convention, not a law. The point is to remove rows that carry almost no information, because dead rows inflate multiple-testing burden and destabilize downstream models. How to check the result: plot mean versus variance of the transformed data. After VST, variance should flatten across the expression range instead of rising with the mean. ## Step 3: Map Probes and Counts to a Shared Gene Identifier You cannot merge a matrix keyed by Affymetrix probe IDs with a matrix keyed by Ensembl gene IDs. Collapse everything to official gene symbols. For microarray, map probes using the platform annotation: ```r # Extract probe-level annotation fdat <- fData(eset) head(fdat[, grep("Gene", colnames(fdat))]) # Aggregate multiple probes per gene by taking the probe with the # highest average expression (a widely used, defensible rule) probe_map <- fdat[["Gene Symbol"]] keep_idx <- !is.na(probe_map) & probe_map != "" expr_ma <- expr_ma[keep_idx, ] colnames(expr_ma) <- probe_map[keep_idx] # Collapse probes to genes expr_ma_gene <- avereps(expr_ma, ID = rownames(expr_ma)) ``` For RNA-seq, Ensembl IDs convert cleanly with the org.Hs.eg.db annotation package or a GTF-derived mapping table. ```r library(AnnotationDbi) library(org.Hs.eg.db) ens <- gsub("\\..*$", "", rownames(expr_seq)) # strip Ensembl version suffix sym <- mapIds(org.Hs.eg.db, keys = ens, keytype = "ENSEMBL", column = "SYMBOL", multiVals = "first") expr_seq <- expr_seq[!is.na(sym), ] rownames(expr_seq) <- sym[!is.na(sym)] # Collapse any duplicate symbols by mean expr_seq_gene <- avereps(expr_seq, ID = rownames(expr_seq)) ``` Two decisions matter here. When several probes map to one gene on a microarray, taking the probe with maximum average expression or the mean of all probes are both standard. The maximum-expression rule is common because it reduces noise from poorly hybridizing probes. When collapsing duplicates, document your rule in the methods section, because different rules change results at the margin. ## Step 4: Intersect Genes and Assess Batch Structure Take the intersection of genes present on both platforms. Microarray probes cover a smaller, biased subset of the genome, typically around 20,000 genes, so your merged matrix will usually be limited by the microarray's gene space. A mucoepidermoid carcinoma study working with two microarray cohorts ended up with 19,565 genes after harmonization, which is a realistic scale [1]. ```r common_genes <- intersect(rownames(expr_ma_gene), rownames(expr_seq_gene)) expr_merged <- cbind(expr_ma_gene[common_genes, ], expr_seq_gene[common_genes, ]) dim(expr_merged) meta <- data.frame( row.names = c(colnames(expr_ma_gene), colnames(expr_seq_gene)), platform = c(rep("microarray", ncol(expr_ma_gene)), rep("rnaseq", ncol(expr_seq_gene))), group = c(group_ma, group_seq) ) ``` Before any correction, look at where the variation sits. Principal component analysis (PCA) will almost always show samples separating by platform rather than by disease, and that is expected: ```r pca_res <- prcomp(t(expr_merged), center = TRUE, scale. = TRUE) plot(pca_res$x[, 1], pca_res$x[, 2], col = factor(meta$platform), pch = factor(meta$group), xlab = "PC1", ylab = "PC2") legend("topright", legend = c("microarray", "rnaseq"), col = 1:2, pch = 16) ``` If samples separate by platform on PC1, you have a batch effect (systematic non-biological variation between processing groups). Correct it with ComBat from the sva package, protecting the biological variable of interest in the model matrix: ```r library(sva) mod <- model.matrix(~ group, data = meta) expr_corrected <- ComBat(dat = expr_merged, batch = meta$platform, mod = mod) pca_res2 <- prcomp(t(expr_corrected), center = TRUE, scale. = TRUE) plot(pca_res2$x[, 1], pca_res2$x[, 2], col = factor(meta$platform), pch = factor(meta$group)) ``` The `mod` argument is the whole trick. It tells ComBat to remove platform-driven variation while preserving disease-control differences. Omitting `mod` risks removing your biology along with your batch. A pooled, batch-corrected framework built on variance-stabilized RNA-seq values was used to harmonize microarray and RNA-seq data in an Arabidopsis stress study, and that pipeline detected 102 consensus differentially expressed genes that neither platform alone would have supported as strongly [9]. An alternative worth knowing: learned cross-platform translation models. PXN is a probabilistic machine learning framework trained to represent biological signal jointly across expression technologies, and once trained it can translate data between platforms, including between microarray and RNA-seq, while removing technology-specific bias [5]. If your integration involves many cohorts and platforms, a trained translator can outperform per-project ComBat, at the cost of depending on an external model. ## Step 5: Define the Candidate Feature Set Do not feed 19,000 genes into a classifier with 80 samples. You need to shrink the feature space to biologically defensible candidates first. Three routes dominate the published literature, and strong studies use two or more in combination. **Route 1, differential expression.** Fit a linear model with limma on the harmonized matrix, using the empirical Bayes moderation that limma provides: ```r design <- model.matrix(~ group, data = meta) fit <- lmFit(expr_corrected, design) fit <- eBayes(fit) deg_tab <- topTable(fit, coef = "groupcase_vs_control", number = Inf) deg_sig <- subset(deg_tab, adj.P.Val < 0.05 & abs(logFC) > 1) ``` Note that log fold changes computed after ComBat are shrunken toward zero. Use the corrected matrix for modeling, but consider computing fold changes on the original per-platform data if effect sizes matter for your narrative. **Route 2, co-expression modules.** Weighted gene co-expression network analysis (WGCNA) groups genes into modules correlated with a trait, then ranks genes within modules by connectivity (hub status). This is used heavily in the integration literature. A PCOS study intersected differentially expressed genes with WGCNA module genes and reported a module of 143 genes most correlated with disease [6]. A dilated cardiomyopathy study intersected 309 differentially expressed genes with 2,093 WGCNA hub genes to define candidates [2]. A mucoepidermoid carcinoma study intersected WGCNA hub genes with differentially expressed genes and then refined candidates through 100 resampling iterations, keeping genes selected in at least 60 percent of iterations [1]. ```r library(WGCNA) allowWGCNAThreads() datExpr <- t(expr_corrected) # samples in rows for WGCNA datExpr <- goodGenesSamples(datExpr, verbose = 3) powers <- c(1:10, 12:20) sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 3) soft_power <- sft$powerEstimate # use printed suggestion net <- blockwiseModules(datExpr, power = soft_power, TOMType = "unsigned", minModuleSize = 30, mergeCutHeight = 0.25, numericLabels = TRUE) module_labels <- net$colors ``` The soft threshold power is the single most important WGCNA parameter. `pickSoftThreshold` prints a table of scale-free topology fit against power. Choose the lowest power that reaches a fit index above 0.8. If no power achieves it, your data may have too few samples or too much residual batch structure. **Route 3, curated gene sets.** If your disease has a mechanistic hypothesis, restrict candidates to a curated list. A sepsis study scored glycolysis activity with AUCell using the HALLMARK_GLYCOLYSIS set and derived candidates from that program [7]. A major depressive disorder study intersected differentially expressed genes with mitochondrial quality control genes [10]. A tuberculosis study started from mitochondria-related genes and identified ten differentially expressed candidates enriched in oxidative phosphorylation [11]. Curated sets shrink the search space and give reviewers a biological story, not just a classifier. The candidate set you carry into machine learning is the intersection or union of these routes, typically a few hundred genes. ## Step 6: Split Your Data Before You Touch a Model ![Decision path asking whether genes were chosen on all samples before cross-validation, leading to a leakage fix or an honest nested resampling design.](/article-images/illustrations/integrate-rna-seq-microarray-machine-learning-leakage-triage.webp) This ordering is where most published pipelines quietly cheat. If you select features on all samples and then cross-validate the model on the same samples, performance estimates are inflated because the test folds already influenced feature choice. Split first. ```r set.seed(42) library(caret) train_idx <- createDataPartition(meta$group, p = 0.7, list = FALSE) X <- t(expr_corrected[candidate_genes, ]) # samples in rows, genes in columns y <- factor(meta$group, levels = c("control", "case")) X_train <- X[train_idx, ] y_train <- y[train_idx] X_test <- X[-train_idx, ] y_test <- y[-train_idx] ``` With small cohorts, a better design is leave-one-out cross-validation (LOOCV), where each sample takes one turn as the test set. A study of oral submucosal fibrosis used LOOCV for model development and held out entirely independent microarray datasets for validation [12]. The critical rule: run the entire pipeline, including any feature re-selection inside each fold, within resampling. Nested cross-validation is the honest way to do this when sample size permits. ## Step 7: Train Machine Learning Models Now the models. The published cross-platform biomarker literature converges on a small set of algorithms, and the right move is to run several and compare rather than pick one in advance. A renal allograft delayed graft function study used machine learning to pinpoint hub genes and build predictive models from integrated microarray and RNA-seq datasets, arriving at a seven-gene hub set including ADM, SERPINH1, SLC2A3, BAG3, NFKBIA, KLF6, and CEBPD [13]. A BK polyomavirus nephropathy study compared nine algorithms and found a support vector machine model performed best, with an AUC of 0.923 for the final diagnostic model built from macrophage-associated features [14]. | Algorithm | R package | What it does best | Watch out for | |--|--|--|--| | LASSO (L1-regularized regression) | glmnet | Shrinks coefficients to zero, giving a compact gene signature | Coefficient choice depends on lambda, always report it | | Elastic net | glmnet | LASSO with correlated features handled better via alpha between 0 and 1 | Needs alpha tuning in cross-validation | | Random forest | randomForest | Captures nonlinear effects and interactions, gives importance scores | Unstable variable importance when features are correlated | | SVM (linear or RBF kernel) | e1071 | Strong performance on small, high-dimensional data | Cost and gamma parameters need tuning | | XGBoost | xgboost | Often the top classifier on structured data | Overfits small cohorts without careful early stopping | | SVM-RFE | caret or custom | Recursive feature elimination, direct feature count control | Slow on large feature sets | | Boruta | Boruta | Wrapper around random forest that keeps all features confirmed as relevant | Conservative, returns more features than a final classifier needs | | Cox-LASSO | glmnet(type = "cox") | Prognostic signatures with survival endpoints | Requires correct survival format and time scaling | A sepsis study combined LASSO, random forest, and Boruta to prioritize five glycolysis-associated candidates (GLRX, MDH1, MDH2, TGFBI, COPB2) that discriminated disease in bulk blood data [7]. A glioblastoma study used LASSO regression alongside seven classifiers and found SVM and k-nearest neighbors gave the highest accuracy on unseen datasets, with PLOD1, COL4A2, and COL5A2 consistently top-ranked across models [3]. The lesson is consistent: agreement between algorithm families, not the performance of any single algorithm, is what identifies robust genes. Run the core three: ```r library(glmnet) library(pROC) # LASSO with internal cross-validation on the training set only cv_fit <- cv.glmnet(x = X_train, y = y_train, family = "binomial", alpha = 1, type.measure = "auc", nfolds = 10) lasso_genes <- coef(cv_fit, s = "lambda.min") |> as.matrix() |> as.data.frame() lasso_genes <- subset(lasso_genes, `1` != 0) rownames(lasso_genes) # your LASSO-selected genes lasso_prob <- predict(cv_fit, newx = X_test, s = "lambda.min", type = "response") lasso_auc <- auc(y_test, as.vector(lasso_prob)) lasso_auc ``` Expected output shape: ``` lasso_auc Area under the curve: 0.88 ``` ```r library(randomForest) rf_fit <- randomForest(x = X_train, y = y_train, ntree = 1000, importance = TRUE) rf_prob <- predict(rf_fit, newdata = X_test, type = "prob")[, "case"] auc(y_test, rf_prob) library(e1071) svm_fit <- svm(x = X_train, y = y_train, probability = TRUE, kernel = "linear", cost = 1) svm_prob <- attr(predict(svm_fit, X_test, probability = TRUE), "probabilities")[, "case"] auc(y_test, svm_prob) ``` For an even-handed comparison across many algorithm families, consider an ensemble framework. An IBD study ran an ensemble of 113 machine learning algorithm combinations across its pipeline to prioritize seven hub genes [15]. A glioblastoma prognostic study used an automated machine learning pipeline that tested over 100 algorithmic combinations and arrived at a 22-gene signature [16]. You do not need 113 combinations for a first analysis. A practical approach is the 10-algorithm panel used in several studies, typically LASSO, ridge, elastic net, random forest, SVM variants, XGBoost, and a few distance- or tree-based classifiers, scored by cross-validated AUC [10][17][18]. ## Step 8: Evaluate Performance Honestly Receiver operating characteristic (ROC) analysis summarizes classifier performance as the area under the curve (AUC), where 0.5 is chance and 1.0 is perfect separation. Plot the curves and compute confidence intervals: ```r roc_lasso <- roc(y_test, as.vector(lasso_prob)) ci(roc_lasso) plot(roc_lasso) ``` AUC alone is not enough for a diagnostic model. The dilated cardiomyopathy framework explicitly evaluated tissue-level model performance with ROC analysis, calibration assessment, nomogram visualization, and decision curve analysis, and it deliberately separated diagnostic performance from localization evidence and exploratory mechanistic analysis [2]. Calibration asks whether predicted probabilities match observed event rates. Decision curve analysis asks whether acting on the model beats acting on everyone or no one across a range of clinical thresholds. Calibration curves are straightforward with the rms package: ```r library(rms) dat_eval <- data.frame(prob = as.vector(lasso_prob), y = as.numeric(y_test) - 1) fit_cal <- lrm(y ~ prob, data = dat_eval, x = TRUE, y = TRUE) cal <- calibrate(fit_cal, method = "boot", B = 200) plot(cal) ``` A well-calibrated model tracks the diagonal. A model with AUC 0.9 and poor calibration can still misclassify the borderline patients you care most about. ## Step 9: External and Cross-Platform Validation Internal resampling is necessary but never sufficient. Apply the frozen signature (the exact genes and the fitted model) to an independent cohort, ideally one measured on the other platform. This is the pattern across the field: the mucoepidermoid carcinoma eleven-gene signature trained on harmonized microarray data and validated in an independent RNA-seq cohort, demonstrating cross-platform reproducibility [1]. The dilated cardiomyopathy models were tested in three external microarray cohorts and one RNA-seq cohort [2]. The nonsmall cell lung cancer immunogenic cell death signature was trained in TCGA and validated in three independent GEO cohorts, one of them RNA-seq based [17]. Validation requires two conversions. First, map the validation dataset's identifiers to your signature's gene symbols. Second, apply the same transformation logic. If you trained on ComBat-corrected data, you cannot apply ComBat fitted on your training set to a new dataset without care. The cleanest options are to refit batch correction with the validation cohort included as a separate batch, or, better, to z-score each gene within each dataset before applying a classifier trained on z-scored data: ```r zscore <- function(x) (x - mean(x)) / sd(x) X_valid <- t(apply(expr_seq_gene[signature_genes, ], 1, zscore)) valid_prob <- predict(cv_fit, newx = X_valid, s = "lambda.min", type = "response") auc(y_valid, as.vector(valid_prob)) ``` Report the validation AUC with its confidence interval. If the AUC drops from 0.9 internally to 0.65 externally, that is a finding about your signature, not a failure of the pipeline, and it is far better that you discover it than a reviewer. ## Step 10: Interpret the Signature and Validate It Biologically A gene list from a classifier is a hypothesis, not a result. Three layers of interpretation turn it into a paper. **Cell-type context.** Bulk signals are mixtures of cell types, and knowing which cell population drives your signature changes the story entirely. The sepsis study scored glycolysis with AUCell at single-cell resolution and found the strongest signal in monocytes and plasma cells, with TGFBI monocyte-enriched and centrally positioned in a ligand-receptor communication network [7]. The major depressive disorder study inferred gene set activity across 18 annotated cell types in a single-nucleus dataset to examine cell-type-specific alterations [10]. If you have a relevant single-cell dataset, score your signature genes per cell type with AUCell or AddModuleScore and report which cell types carry the signal. ```r # Conceptual AUCell call on a single-cell matrix library(AUCell) ranks <- AUCell_buildRankings(sc_expr) auc_scores <- AUCell_calcAUC(list(Signature = signature_genes), ranks) summary(getAuc(auc_scores)["Signature", ]) ``` **Explainability.** For tree-based or ensemble models, SHAP (SHapley Additive exPlanations) values attribute each gene's contribution to each prediction. The glioblastoma prognostic study applied SHAP analysis to explain each gene's contribution to risk in its 22-gene model, which achieved a concordance index of approximately 0.72 and separated patients into survival-distinct high and low risk groups [16]. Explainability also functions as a sanity check. If a SHAP analysis says your model runs almost entirely on one immunoglobulin gene, you may be detecting batch or immunoglobulin contamination rather than disease biology. **Experimental confirmation.** The gold standard endpoint in this literature is independent expression measurement. The sepsis TGFBI finding was confirmed by qRT-PCR in a cecal ligation and puncture mouse model [7]. The tuberculosis mitochondrial signature (STAT2, CASP1, COX7B prioritized by LASSO and SVM-RFE) was validated by qPCR in an H37Rv-infected THP-1 cell model [11]. The glioblastoma PLOD1 finding was confirmed in patient samples by qPCR [3]. Plan a small confirmation experiment, even a cell line stimulation model, because computational validation alone increasingly fails review. ## Troubleshooting: Common Errors and Exact Fixes Most cross-platform projects stall on a handful of recurring failures. Here are the ones you will actually see, with their real messages and fixes. **Error:** `Error in estimateSizeFactorsForMatrix(counts(object)) : every gene contains at least one zero, cannot compute log geometric means` This DESeq2 error means your gene filter was too permissive and every remaining row has at least one zero count. Fix by filtering more aggressively, for example `rowSums(counts(dds) >= 10) >= 5`, before calling `DESeq()`. **Error:** `Error in checkFullRank(modelMatrix) : the model matrix is not full rank` DESeq2 detected a design where one term is perfectly confounded with another, for example batch identical to condition. You cannot estimate both. Either merge confounded levels or drop the redundant term. This is the single most common RNA-seq error in integrative designs because platform and disease often correlate in public data. **Error:** `Error in solve(...) : system is computationally singular` ComBat or a covariance operation hit a singular matrix, usually because samples outnumber independent conditions or because a batch contains a single sample. Ensure every batch level has multiple samples in every condition, and drop genes with zero variance before correction. **Error:** `Error: cannot allocate vector of size 4.1 Gb` Memory exhaustion, typically during WGCNA's topological overlap matrix on thousands of genes. Reduce `maxBlockSize` in `blockwiseModules`, pre-filter to your candidate set, or run on a machine with more RAM. **Error:** `Error in download.file(...) : cannot open URL` GEO occasionally blocks parallel rapid downloads or the FTP path has changed. Retry with a fresh R session, confirm the accession exists on the GEO website, and use `destdir` with a writable path. Downloading through the browser as a fallback is acceptable. **Error:** `Error in cv.glmnet(...) : x should be a matrix with 2 or more columns` Your feature selection left zero or one gene. Either your intersection was empty (check identifier mapping, the most common cause) or your filters were chained too aggressively. Print `length(common_genes)` and `length(candidate_genes)` at each stage. **Error:** `NA/NaN/Inf in foreign function call (arg 0)` A [random forest](/blog/data-analysis/what-is-random-forest) or SVM hit non-finite values. After ComBat and log transforms, run `sum(is.na(expr_corrected))` and `sum(!is.finite(expr_corrected))` and impute or drop offending rows. Microarray matrices frequently carry NA rows for probes with no annotation. **Error:** `Found 2 genes with zero variance` warnings in limma Zero-variance genes break linear model fitting. Run `expr_corrected <- expr_corrected[apply(expr_corrected, 1, var) > 0, ]` before `lmFit`. The general bioinformatics troubleshooting principle: print dimensions and check for NA values after every transformation. A pipeline that prints `dim()` at each stage tells you within one run which step broke. ## Common Mistakes and Limitations These are the errors that produce confident-looking papers that do not replicate. **Feature leakage.** Selecting genes on all samples before cross-validation is the most damaging mistake in this workflow, and it inflates AUC by a wide margin in small cohorts. Every feature choice must happen inside the training folds. **Correcting away the biology.** Aggressive batch correction applied without a preserved design matrix can remove the disease signal you are hunting. Always run ComBat with `mod` including your group variable, and always compare PCA before and after correction to confirm the group separation you expect survives. **Treating scales as interchangeable.** Raw RNA-seq counts and microarray intensities are not comparable, and even log2 counts per million differ in distribution from log2 microarray values. VST-transformed RNA-seq values, possibly followed by ComBat or z-scoring, is the minimum defensible harmonization [9]. Direct concatenation without transformation is never acceptable. **Probe mapping carelessness.** Probes that cross-hybridize to multiple genes, genes measured by a single probe, and version-suffixed Ensembl IDs all silently corrupt the intersection. A mismatch in identifier handling is a leading cause of "empty intersection" and "AUC collapses in validation" failures. **Composition artifacts.** Whole blood, PBMCs, and tissue biopsies have different cellular compositions, and disease changes that composition. A signature that mostly tracks an expanded monocyte fraction will not translate between blood compartments. Single-cell overlay, as used in the sepsis and tuberculosis studies, is the practical check [7][11]. **Small samples with many features.** With 20 cases and 20 controls, any model can achieve perfect internal separation. Report [confidence intervals](/blog/data-analysis/confidence-interval-formula-examples), prefer external validation over internal resampling, and resist claims built on a single held-out set of ten samples. **Overclaiming mechanism from correlation.** Co-expression, differential expression, and classifier importance are associative. Some studies appropriately layer mechanistic tools, for example [transcription factor](/knowledge/molecular-biology/transcription-factor) network inference or Mendelian randomization in the cardiomyopathy work, and frame them explicitly as hypothesis generating [2][8]. Match your claims to that standard. **Clinical overreach.** An AUC from retrospective public data is a research result. Individual patient care decisions require qualified clinicians, prospective cohorts, and regulatory pathways that no public-data pipeline substitutes for. ## Practical Tips That Save Days Standardize on official gene symbols early and store a mapping log, because re-deriving mappings six months later is worse than doing it carefully once. Write the version numbers of R, Bioconductor, and every package into a `sessionInfo()` file saved with your results, since cross-platform harmonization results can shift with package versions. Freeze the signature gene list to a plain text file before validation and never edit it after seeing validation results, because post hoc signature adjustment is indistinguishable from overfitting. When comparing many algorithms, compare them on identical folds using caret or a fixed seed, otherwise differences in resampling noise masquerade as algorithm differences. And when a result looks too good, an AUC above 0.99 from a small public cohort usually means leakage, label errors, or a trivially confounded variable such as batch or sample source, so investigate before celebrating. ## Frequently Asked Questions ### Can I merge RNA-seq counts and microarray intensities directly into one matrix? No. The two technologies produce different scales and signal distributions, and direct merging produces a matrix in which platform dominates every axis of variation. Transform RNA-seq counts with variance-stabilizing transformation, log2-transform microarray data if needed, intersect on gene symbols, and then apply batch correction such as ComBat [9]. ### How many samples do I need for a credible machine learning biomarker study? There is no hard threshold, but the published integrative studies cited here typically work with tens to hundreds of samples per group and rely on independent external validation rather than sample volume alone. With fewer than 20 samples per group, treat any internal performance estimate as unstable and put your effort into external validation. ### Which machine learning algorithm should I start with? Start with LASSO [logistic regression](/blog/data-analysis/logistic-regression-definition-formula), because it produces a compact, interpretable gene signature and is the most widely used method in this literature [4]. Add a random forest and a support vector machine as a comparison panel, and consider XGBoost if you have enough samples. Agreement across algorithm families is more informative than any single model. ### Do I need batch correction if my datasets come from different laboratories? Yes. Different laboratories, reagent lots, and platforms all introduce systematic non-biological variation. Even datasets on the same platform from different labs usually require correction. Protect the biological contrast in the ComBat design matrix and verify with PCA that disease separation survives correction. ### How do I validate a signature discovered on microarray data in an RNA-seq cohort? Map the RNA-seq dataset to the same gene symbols, transform it comparably (VST followed by z-scoring or batch correction), apply the frozen model without retraining, and report the AUC with a confidence interval. This cross-platform validation step is now standard practice in published biomarker studies [1][2][17]. ### What is data leakage in this workflow and how do I avoid it? Data leakage is any information from the test set influencing model construction. The classic form here is differential expression testing or WGCNA module identification on all samples before cross-validation. Perform every feature selection step inside training folds, and keep the final signature completely frozen before it touches any external cohort. ### Should I include single-cell RNA-seq data in this pipeline? Use it for interpretation, not for classifier training. Single-cell data tells you which cell types express your signature genes and whether the bulk signal reflects a genuine per-cell program or a compositional shift. Several published pipelines used single-cell or single-nucleus datasets specifically for this localization step [10]. ### What should I do when a gene exists on one platform but not the other? Drop it from the integrated analysis and record the drop. Restricting to the intersection of platforms is the safe default. If a candidate gene is absent from the microarray platform entirely, you can still report it from the RNA-seq data alone, but it cannot participate in the integrated model or the cross-platform validation claim. ## Sources 1. [A cross-platform validated biomarker signature for mucoepidermoid carcinoma identified by integrated WGCNA and machine learning.](https://pubmed.ncbi.nlm.nih.gov/42664892/) 2. [Identification and multi-layered validation of seven diagnostic biomarkers for dilated cardiomyopathy via integrative machine learning, single-cell transcriptomics, and Mendelian randomization.](https://pubmed.ncbi.nlm.nih.gov/42343898/) 3. [Integrative multi-platform gene expression and machine learning analysis for glioblastoma biomarker discovery with experimental validation.](https://pubmed.ncbi.nlm.nih.gov/41705725/) 4. [Pathway to precision: machine learning and bioinformatics in diabetes gene expression studies.](https://pubmed.ncbi.nlm.nih.gov/41424950/) 5. [PXN Unlocks the Power of Public Gene Expression Data Through Cross-Technology Integration.](https://pubmed.ncbi.nlm.nih.gov/42182239/) 6. [Integrated weighted gene co-expression network analysis and machine learning analysis identifies SEC14L5 as a potential biomarker for polycystic ovary syndrome.](https://pubmed.ncbi.nlm.nih.gov/41792859/) 7. [Immunometabolic reprogramming and glycolysis-associated signatures in sepsis: insights from single-cell RNA sequencing and machine learning.](https://pubmed.ncbi.nlm.nih.gov/42220506/) 8. [Machine learning, whole-transcriptome and integrative omics analysis reveals key regulatory networks governing human spermatogonial stem cells.](https://pubmed.ncbi.nlm.nih.gov/41665764/) 9. [Integrative cross-platform transcriptomic analysis coupled with machine learning identifies candidate biomarkers of abiotic stress resilience in Arabidopsis thaliana.](https://pubmed.ncbi.nlm.nih.gov/41764430/) 10. [Integrative multi-omics and machine learning analysis identifies candidate biomarkers associated with mitochondrial quality control in major depressive disorder.](https://pubmed.ncbi.nlm.nih.gov/42656552/) 11. [Identifying mitochondria-related signatures for tuberculosis diagnosis through machine learning on single-cell transcriptomics data and experimental verification.](https://pubmed.ncbi.nlm.nih.gov/41967724/) 12. [Delineating novel diagnostic biomarkers and therapeutic targets for oral submucosal fibrosis: an integrative multi-omics and machine learning approach.](https://pubmed.ncbi.nlm.nih.gov/42158597/) 13. [Immune Infiltration-Related Genes as Potential Biomarkers and Predicted Targets for Renal Allograft Delayed Graft Function and Survival Outcome: An Integrated Machine Learning Approach and Drugs Analysis.](https://pubmed.ncbi.nlm.nih.gov/42283098/) 14. [BKPyVN-associated macrophage subpopulations and diagnostic biomarker identification - combine machine learning with single-cell analysis.](https://pubmed.ncbi.nlm.nih.gov/42953158/) 15. [Integrative Multi-Omics and Machine Learning Reveal Senescent Fibroblast-Associated ADAMTS2 as a Critical Indicator of Intestinal Inflammation in IBD.](https://pubmed.ncbi.nlm.nih.gov/42773852/) 16. [Interpretable prognostic modeling of glioblastoma using cross-cohort transcriptome integration and machine learning approaches.](https://pubmed.ncbi.nlm.nih.gov/41713199/) 17. [Integrated transcriptomic analysis and machine learning identify immunogenic cell death genes as prognostic markers and therapeutic targets in non-small cell lung cancer.](https://pubmed.ncbi.nlm.nih.gov/42180906/) 18. [Machine learning integration with multi-omics data constructs a robust prognostic model and identifies PTGES3 as a therapeutic target for precision oncology in lung adenocarcinoma.](https://pubmed.ncbi.nlm.nih.gov/41089696/)