When One Omic Isn't Enough: Integrating RNA-seq and Proteomics to Find Robust Biomarkers

Every translational scientist knows the frustration: a gene beautifully overexpressed in your RNA-seq data turns out to be undetectable at the protein level in patient plasma. The root cause is a well-documented but underappreciated fact — across the human proteome, mRNA abundance and protein abundance correlate with a Pearson r of only about 0.40. Running a single-omic analysis and treating it as ground truth means you are working with noise masquerading as signal.

The solution is not to choose one layer over the other. It is to run both, integrate the results, and keep only the subset where message and protein agree. That subset — concordant across both molecular layers — is far more likely to reflect genuine biological dysregulation, to survive the journey from discovery cohort to clinical validation, and to be actionable in a diagnostic or therapeutic context.

BioMate executes this as a four-step chained workflow: bulk RNA-seq differential expression, mass spectrometry proteomics, multi-omics integration, and weighted biomarker ranking. Below we walk through each step using a rheumatoid arthritis plasma biomarker discovery scenario grounded in publicly available CPTAC and PRIDE data.

Step 1 — Bulk RNA-seq Differential Expression (DESeq2)

The first step queries synovial tissue transcriptomes from RA patients versus healthy controls. BioMate routes this to the nf-core/rnaseq pipeline, which handles STAR alignment, RSEM quantification, and MultiQC reporting. Counts are passed directly to a DESeq2 R workflow: negative binomial modeling with shrinkage estimators (apeglm), Benjamini-Hochberg correction, and a log2 fold-change threshold of ≥1.0 combined with an adjusted p-value <0.05.

Against ENCODE-curated reference data and GEO Series GSE77298 (RA synovial RNA-seq, n=20 RA, n=10 control), this pass identifies 2,341 differentially expressed genes. The highest-confidence upregulated signals include plasma cell markers (CD38, IGHG1, IGKC), innate immune effectors (S100A8, S100A9), and Fc-receptor family members (FCGR3A, FCGR1A). These are established RA biology — a necessary sanity check before moving to proteomics.

Step 2 — Mass Spectrometry Proteomics (limma on Protein Abundances)

RNA data alone cannot tell you whether a gene product is present in plasma at a detectable and diagnostically meaningful concentration. Step 2 addresses this directly. BioMate ingests label-free quantification (LFQ) intensity matrices from PRIDE dataset PXD013149 — a high-coverage LC-MS/MS study of RA synovial fluid from 18 patients and 10 controls, with 8,547 proteins quantified after filtering — and applies a limma linear model with empirical Bayes variance shrinkage.

limma is the correct choice here rather than a second DESeq2 run: proteomics LFQ data are approximately log-normal after transformation, imputation is required for missing values (MNAR pattern; minimum probability imputation), and the moderate sample sizes make the empirical Bayes shrinkage especially valuable. The output is 892 proteins with significant differential abundance (FDR <0.05, |log2FC| ≥0.75). The top upregulated proteins overlap substantially but imperfectly with the transcriptomic list — which is precisely the point.

"The proteomics layer is not a confirmation of the RNA layer. It is an orthogonal measurement of a different molecular event. Treating them as the same signal is how biomarker projects fail."

Step 3 — Multi-Omics Integration (MOFA+ and Concordance Overlap)

BioMate's integration step runs two complementary approaches and reports both. First, a simple intersection filter: genes present in both the DE gene list and the DE protein list, with concordant direction of change (both up or both down). This is interpretable, fast, and appropriately conservative — it removes the roughly 60% of candidates that are differentially expressed at the transcript level but show no protein-level change, or vice versa.

Second, MOFA+ (Multi-Omics Factor Analysis v2) is run to identify latent factors that jointly explain variance across both data modalities. MOFA+ does not require one-to-one gene-protein correspondence — it discovers coordinated biological programs that manifest across layers even when individual gene-protein correlations are low. Factors with high variance explained in both layers and loadings enriched for immune effector pathways are flagged as high-confidence biological signal.

The combined output after both filters: 487 concordant candidates where mRNA and protein changes are directionally consistent, statistically significant in both assays, and captured by at least one MOFA+ factor with cross-layer variance ≥5%.

Step 4 — Biomarker Candidate Ranking by Combined Evidence Score

The final step converts the 487 concordant candidates into an ordered list. BioMate computes a weighted evidence score combining: the absolute log2FC in RNA (weight 0.35), the absolute log2FC in protein (weight 0.35), the −log10(adjusted p-value) in each assay (weight 0.15 each), and a MOFA+ factor loading bonus (+0.1 if the gene-protein pair appears in the top 20 loadings of any cross-layer factor). The formula is deliberately simple and interpretable — it reflects the evidence structure rather than a black-box ML classifier.

Table 1. Illustrative ranked biomarker candidates (simulated output). Scores are representative of the weighting scheme applied to concordant mRNA+protein signals from GSE77298 and PRIDE PXD013149; actual values depend on the analysis inputs.

Why Concordance Matters Clinically

The canonical mRNA-protein correlation figure from Maier et al. (2009) and the larger-scale analyses from Wilhelm et al. (2014) and Gry et al. (2009) all converge on the same uncomfortable finding: across the human proteome, knowing a gene's mRNA abundance explains only about 16–40% of the variance in its protein abundance. Post-transcriptional regulation — miRNA silencing, RNA-binding protein competition, translational efficiency differences, protein half-life — accounts for most of the gap.

For a biomarker discovery program, this has a direct consequence. If you run only RNA-seq and nominate the top 50 genes, roughly 60 of them will not produce a corresponding protein signal in the biofluid you plan to assay. You will spend months pursuing ghosts. The reverse is equally true: protein markers without transcriptomic support are more likely to reflect post-translational regulation, secretion from a distal tissue, or technical variation in the MS assay — harder to mechanistically interpret and harder to target therapeutically.

Concordant markers — those where both the transcript and the protein move in the same direction with statistical confidence — have passed a two-layer orthogonal filter. They are more likely to be driven by genuine transcriptional dysregulation with functional protein-level consequences. In the RA example, CD38 and S100A8 (calprotectin) are both on this list and both have clinical validation: serum calprotectin correlates with DAS28 disease activity, and CD38 is the target of daratumumab, now in trials for refractory RA.

Running This in BioMate

The entire four-step chain executes as a single conversational workflow request. You describe your data — FASTQ files or count matrices for RNA-seq, LFQ intensity tables for proteomics, sample metadata with group labels — and BioMate sequences the four steps automatically, passing outputs between them without manual file handling. The integration and ranking steps run on AWS Batch workers with full Bioconductor and MOFA+ environments installed; you do not manage packages or compute. The ranked candidate table, scatter plot, Venn diagram, and a downloadable PDF methods report are all generated as outputs of the completed run.

For teams working with existing CPTAC or PRIDE datasets, BioMate can ingest the processed data matrices directly. For prospective cohorts, the nf-core/rnaseq alignment step and the MaxQuant or FragPipe proteomics preprocessing steps are both available as upstream pipeline stages in the same workflow chain.