Introduction

Multi-omics integration is now central to mechanistic studies of complex diseases. Combining host transcriptomics (bulk RNA-seq) with gut microbiome profiling (16S rRNA or metagenomics) has revealed key host-microbe interactions in inflammatory bowel disease, tuberculosis, HIV, and metabolic syndrome. Yet the computational toolkit for this specific integration scenario remains fragmented: no single Bioconductor package takes a researcher from raw paired data through joint dimensionality reduction, biomarker discovery, and diagnostic classification.

MultiOmicsBridge closes this gap. The package provides a unified, reproducible workflow with five modules:

Module	Functions	Purpose
1. Harmonization	`loadHostData()`, `loadMicrobiomeData()`, `matchSamples()`	Import, normalize, pair data
2. Joint Dim. Reduction	`jointDimReduction()`	DIABLO sparse multi-block PLS-DA
3. Biomarker Discovery	`biomarkerDiscovery()`	Cross-omics ranked biomarkers
4. Classification	`diagnosticClassifier()`	Host-only vs microbiome-only vs joint
5. Visualization	`plotIntegration()`, `plotBiomarkerNetwork()`, `plotClassifierComparison()`, `plotSankey()`	Publication figures

The one-call wrapper MultiOmicsBridgeAnalysis() executes all five modules and returns a structured MOBResult S4 object.

Installation

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("MultiOmicsBridge")

Simulated Multi-Omics Dataset

We demonstrate the package on a simulated dataset that mimics a real paired study: 20 samples (10 healthy controls, 10 disease cases), 500 host genes, and 80 microbial taxa. To model a biologically realistic signal, we inject upregulation of the first 20 genes and enrichment of the first 5 taxa in the disease group.

library(MultiOmicsBridge)
library(SummarizedExperiment)

set.seed(2026)
n_genes   <- 500
n_taxa    <- 80
n_samples <- 20

# Host RNA-seq count matrix (genes x samples)
host_counts <- matrix(
    rpois(n_genes * n_samples, lambda = 150L),
    nrow = n_genes, ncol = n_samples
)
rownames(host_counts) <- paste0("Gene",   seq_len(n_genes))
colnames(host_counts) <- paste0("Sample", seq_len(n_samples))

# Inject transcriptional signal: genes 1-20 upregulated in disease
host_counts[seq_len(20), seq(11, n_samples)] <-
    host_counts[seq_len(20), seq(11, n_samples)] * 5L

# Microbiome count matrix (taxa x samples)
mb_counts <- matrix(
    rpois(n_taxa * n_samples, lambda = 40L),
    nrow = n_taxa, ncol = n_samples
)
rownames(mb_counts) <- paste0("Taxon",  seq_len(n_taxa))
colnames(mb_counts) <- paste0("Sample", seq_len(n_samples))

# Inject microbial enrichment: taxa 1-5 enriched in disease
mb_counts[seq_len(5), seq(11, n_samples)] <-
    mb_counts[seq_len(5), seq(11, n_samples)] * 4L

# Sample metadata
col_data <- data.frame(
    condition = rep(c("Control", "Disease"), each = 10),
    age       = sample(30:65, n_samples, replace = TRUE),
    sex       = sample(c("M", "F"), n_samples, replace = TRUE),
    row.names = paste0("Sample", seq_len(n_samples))
)

# Simulated data dimensions
dim(host_counts)
#> [1] 500  20
dim(mb_counts)
#> [1] 80 20

Module 1: Data Harmonization

Step 1.1 — Load host RNA-seq data

loadHostData() applies TMM normalization (via edgeR) followed by limma-voom precision weighting. Both raw counts and log2-CPM values are stored in the output SummarizedExperiment.

host_se <- loadHostData(host_counts, col_data = col_data)
host_se
#> class: SummarizedExperiment 
#> dim: 500 20 
#> metadata(0):
#> assays(2): counts voom
#> rownames(500): Gene1 Gene2 ... Gene499 Gene500
#> rowData names(0):
#> colnames(20): Sample1 Sample2 ... Sample19 Sample20
#> colData names(3): condition age sex
assayNames(host_se)
#> [1] "counts" "voom"

Step 1.2 — Load microbiome data

loadMicrobiomeData() applies centered log-ratio (CLR) transformation. This is the compositionally appropriate normalization that removes the unit-sum constraint inherent to relative abundance data.

mb_se <- loadMicrobiomeData(mb_counts, normalization = "CLR",
                             min_prevalence = 0.1)
mb_se
#> class: SummarizedExperiment 
#> dim: 80 20 
#> metadata(0):
#> assays(2): counts CLR
#> rownames(80): Taxon1 Taxon2 ... Taxon79 Taxon80
#> rowData names(0):
#> colnames(20): Sample1 Sample2 ... Sample19 Sample20
#> colData names(0):
assayNames(mb_se)
#> [1] "counts" "CLR"

# Verify CLR: column means should be ≈ 0
round(colMeans(assay(mb_se, "CLR")), 10)[1:5]
#> Sample1 Sample2 Sample3 Sample4 Sample5 
#>       0       0       0       0       0

Step 1.3 — Match paired samples

matchSamples() finds the intersection of sample names across both SEs and packages them into a MultiAssayExperiment (MAE). Unpaired samples are transparently reported and excluded.

mae <- matchSamples(host_se, mb_se, min_paired = 10)
mae
#> A MultiAssayExperiment object of 2 listed
#>  experiments with user-defined names and respective classes.
#>  Containing an ExperimentList class object of length 2:
#>  [1] host: SummarizedExperiment with 500 rows and 20 columns
#>  [2] microbiome: SummarizedExperiment with 80 rows and 20 columns
#> Functionality:
#>  experiments() - obtain the ExperimentList instance
#>  colData() - the primary/phenotype DataFrame
#>  sampleMap() - the sample coordination DataFrame
#>  `$`, `[`, `[[` - extract colData columns, subset, or experiment
#>  *Format() - convert into a long or wide DataFrame
#>  assays() - convert ExperimentList to a SimpleList of matrices
#>  exportClass() - save data to flat files

All 20 samples are paired in this dataset. In real studies, matchSamples() will report and exclude any samples missing from either platform.

Module 2: Joint Dimensionality Reduction

jointDimReduction() runs DIABLO (Data Integration Analysis for Biomarker discovery using Latent cOmponents) from the mixOmics package. DIABLO simultaneously maximises the covariance between the host and microbiome latent variates while discriminating between outcome groups, using L1 sparsity to select only the most informative features.

outcome <- col_data$condition   # "Control" or "Disease"

dr_result <- jointDimReduction(
    mae,
    outcome         = outcome,
    n_components    = 2L,
    n_features_host = 40L,
    n_features_mb   = 15L,
    design_off_diag = 0.1
)

# Score matrix dimensions
dim(dr_result$scores)
#> [1] 20  2

# Loading matrix dimensions (host, microbiome)
dim(dr_result$host_loadings)
#> [1] 500   2
dim(dr_result$mb_loadings)
#> [1] 80  2

# Explained variance per component
round(dr_result$explained_variance, 4)
#>  Comp1  Comp2 
#> 0.0978 0.0507

Module 3: Biomarker Discovery

biomarkerDiscovery() ranks host genes and microbial taxa using their DIABLO sparse loading scores (L2 norm across components) and annotates each with its maximum Spearman correlation to a feature in the other omics layer. High-loading features with strong cross-omics correlations represent the most credible multi-omics biomarker candidates.

bm_table <- biomarkerDiscovery(mae, dr_result, n_biomarkers = 30)

# Show top 10 biomarkers
bm_df <- as.data.frame(bm_table)
bm_df_sorted <- bm_df[order(bm_df$loading_score, decreasing = TRUE), ]
head(bm_df_sorted[, c("feature","omics_layer","loading_score",
                        "component","max_cross_cor","top_partner")], 10)
#>         feature omics_layer loading_score component max_cross_cor top_partner
#> Taxon53 Taxon53  microbiome     0.4986963         2     0.5172932     Gene201
#> Taxon39 Taxon39  microbiome     0.4824716         2     0.6375940     Gene201
#> Gene201 Gene201        host     0.4510873         2     0.6375940     Taxon39
#> Gene303 Gene303        host     0.4438466         2     0.6481203     Taxon46
#> Taxon3   Taxon3  microbiome     0.4348072         1     0.8330827      Gene16
#> Taxon5   Taxon5  microbiome     0.4344162         1     0.8601504       Gene5
#> Taxon4   Taxon4  microbiome     0.4295783         1     0.9278195      Gene11
#> Taxon2   Taxon2  microbiome     0.4240883         1     0.8616541      Gene10
#> Taxon1   Taxon1  microbiome     0.4229655         1     0.8300752      Gene19
#> Gene113 Gene113        host     0.3532316         2     0.6887218     Taxon28

# Biomarker counts by omics layer
n_host <- sum(bm_df$omics_layer == "host")
n_mb   <- sum(bm_df$omics_layer == "microbiome")
c(host = n_host, microbiome = n_mb)
#>       host microbiome 
#>         30         30

# Are the injected signal genes recovered?
injected <- paste0("Gene", seq_len(20))
detected <- intersect(bm_df$feature, injected)
length(detected)  # out of 20 injected
#> [1] 20

Module 4: Diagnostic Classification

diagnosticClassifier() trains three Random Forest models using stratified k-fold cross-validation and compares their AUC-ROC values. This comparison directly quantifies the added diagnostic value of multi-omics integration.

clf_results <- diagnosticClassifier(
    mae,
    outcome         = outcome,
    biomarker_table = bm_table,
    cv_folds        = 5L,
    seed            = 42L
)

# Classifier AUC-ROC (mean ± SD across 5-fold CV)
models <- c("host_only", "microbiome_only", "joint")
data.frame(
    Model    = models,
    Mean_AUC = vapply(models, function(m) clf_results[[m]]$mean_auc, numeric(1)),
    SD_AUC   = vapply(models, function(m) clf_results[[m]]$sd_auc,   numeric(1))
)
#>                           Model Mean_AUC SD_AUC
#> host_only             host_only        1      0
#> microbiome_only microbiome_only        1      0
#> joint                     joint        1      0

# AUC gain: multi-omics vs. host-only
delta <- clf_results$joint$mean_auc - clf_results$host_only$mean_auc
round(delta, 3)
#> [1] 0

Full Pipeline with One Call

For standard analyses, MultiOmicsBridgeAnalysis() runs all modules in sequence and returns a MOBResult S4 object:

result <- MultiOmicsBridgeAnalysis(
    mae,
    outcome,
    n_components    = 2L,
    n_features_host = 40L,
    n_features_mb   = 15L,
    n_biomarkers    = 30L,
    cv_folds        = 5L,
    seed            = 42L
)

result
#> MOBResult
#>   Integration    : DIABLO 
#>   Samples        : 20 
#>   Components     : 2 
#>   Biomarkers     : 60 (30 host, 30 microbiome) 
#>   -- Classifier AUC (mean +/- SD) -----------
#>   host_only:       1.000 +/- 0.000
#>   microbiome_only: 1.000 +/- 0.000
#>   joint:           1.000 +/- 0.000
#>   Outcome levels : Control vs Disease

Access individual results with typed accessor methods:

# Integration scores: samples x components
head(integrationScores(result))
#>             comp1      comp2
#> Sample1 -4.577809 -0.2783086
#> Sample2 -4.756835 -1.7897451
#> Sample3 -4.514106 -0.1633972
#> Sample4 -4.714052 -0.1824974
#> Sample5 -4.470973 -1.5173940
#> Sample6 -4.897541  2.4213213

# Top biomarkers
head(as.data.frame(biomarkers(result))[, c("feature","omics_layer",
                                             "loading_score")])
#>         feature omics_layer loading_score
#> Gene201 Gene201        host     0.4510873
#> Gene303 Gene303        host     0.4438466
#> Gene113 Gene113        host     0.3532316
#> Gene147 Gene147        host     0.2376692
#> Gene6     Gene6        host     0.2204654
#> Gene4     Gene4        host     0.2202405

# Classifier performance
perf <- performance(result)
# Joint classifier AUC
round(perf$joint$mean_auc, 3)
#> [1] 1

Module 5: Visualization

Joint integration biplot

The integration biplot projects all samples onto the first two DIABLO latent variates. Points are colored by outcome group, and arrow overlays show the top contributing features from each omics layer. Clear separation between groups along Component 1 indicates that DIABLO has identified a discriminant axis, while the loading arrows reveal which genes and taxa drive that separation.

plotIntegration(result, outcome = outcome, n_loading_arrows = 5)

DIABLO sample scores coloured by outcome group. Loading vectors show the top contributing features per omics layer.

Cross-omics biomarker network

The cross-omics correlation heatmap displays Spearman rank correlations between the top-ranked host genes (rows) and microbial taxa (columns). Rows and columns are reordered by hierarchical clustering to group features with similar correlation profiles. Tiles colored toward red indicate strong positive correlations; blue tiles indicate inverse relationships. Features with consistently strong cross-omics correlations are candidate mediators of host–microbe interactions.

plotBiomarkerNetwork(result, mae, n_host = 15, n_mb = 10)

Clustered heatmap of Spearman correlations between the top host genes and microbial taxa. Strong positive/negative correlations indicate candidate host-microbe interactions.

Classifier comparison (bar chart)

The classifier comparison bar chart displays the mean AUC-ROC ± standard deviation from 5-fold cross-validation for host-only, microbiome-only, and joint classifiers. Comparing the three models side by side directly quantifies the diagnostic gain from multi-omics integration. Higher bars indicate stronger discriminative power; overlapping error bars suggest no significant difference.

plotClassifierComparison(result, type = "bar")

Mean cross-validated AUC-ROC for host-only, microbiome-only, and joint classifiers. The multi-omics advantage is immediately visible.

Classifier comparison (ROC curves)

Overlaid ROC curves from the last cross-validation fold illustrate the sensitivity–specificity trade-off for each classifier configuration. A curve closer to the top-left corner indicates better classification performance. Comparing curves visually reveals whether the joint model consistently outperforms single-omics classifiers across all operating points.

plotClassifierComparison(result, type = "roc")

Overlaid ROC curves from the last cross-validation fold for each classifier configuration.

Feature flow (Sankey) diagram

The Sankey (alluvial) diagram traces individual features from their omics layer of origin, through the biomarker selection step, to the outcome classes they help distinguish. Edge width is proportional to the DIABLO loading score, so wider flows represent more influential features. This visualization illustrates how MultiOmicsBridge “bridges” information across host and microbiome data types.

plotSankey(result, n_features = 8)

Feature flow from omics layer through selected biomarkers to outcome classes. Edge width is proportional to loading score.

Structured text report

The generateReport() function produces a structured text summary of all analysis results, suitable for inclusion in supplementary materials or quick review of the integrated analysis outcomes.

generateReport(result, n_top = 8)

Working with Real Data

HMP2 / IBDMDB dataset

The IBDMDB dataset (Franzosa et al. 2019) is the recommended primary validation dataset. To use it with MultiOmicsBridge:

# The HMP2 data is available from the IBDMDB portal.
# After downloading, load as:
library(MultiOmicsBridge)

# Host RNA-seq
host_se <- loadHostData(
    counts   = host_count_matrix,   # genes x samples
    col_data = sample_metadata
)

# Microbiome 16S
mb_se <- loadMicrobiomeData(
    taxa_table    = taxa_count_matrix,   # taxa x samples
    normalization = "CLR"
)

# Match and run
mae    <- matchSamples(host_se, mb_se)
result <- MultiOmicsBridgeAnalysis(mae, outcome = sample_metadata$diagnosis)
result

Complex disease contexts

MultiOmicsBridge is designed to generalize to complex disease contexts. For tuberculosis studies (e.g. GEO: GSE79362), the same workflow applies:

# TB study: blood RNA-seq + sputum microbiome
host_se <- loadHostData(blood_counts, col_data = patient_metadata)
mb_se   <- loadMicrobiomeData(sputum_taxa, normalization = "CLR")
mae     <- matchSamples(host_se, mb_se)

result  <- MultiOmicsBridgeAnalysis(
    mae,
    outcome         = patient_metadata$tb_status,
    n_features_host = 100,
    n_features_mb   = 30,
    cv_folds        = 5
)
plotClassifierComparison(result, type = "bar")

Design Principles

Why CLR for microbiome data?

Microbiome count data is compositional: only relative abundances are observed. Analyzing composites with Euclidean distances or Pearson correlations leads to spurious results (the Aitchison problem). The CLR transformation maps compositional data to real space by:

\[clr(x_j) = \log\left(\frac{x_j + \delta}{g_\delta(x)}\right)\]

where \(g_\delta(x) = \exp\left(\frac{1}{D}\sum_{k}\log(x_k + \delta)\right)\) is the geometric mean with pseudocount \(\delta\). This removes the unit-sum constraint and enables standard Euclidean geometry, making correlations between host genes and microbial taxa statistically valid.

Why DIABLO for integration?

Unlike concatenation (simply merging feature matrices) or independent analysis of each omics layer, DIABLO enforces cross-block correlation: the latent variates from the host and microbiome blocks are constrained to covary. This means DIABLO identifies features that change together across conditions, which is biologically more meaningful than features that are simply marginally associated with the outcome in each dataset independently.

Quantifying the multi-omics advantage

A key scientific contribution of every multi-omics study is demonstrating that integrating two data types improves diagnostic performance over either alone. diagnosticClassifier() makes this comparison automatic and transparent, following best practices from clinical prediction modelling: stratified k-fold cross-validation to estimate predictive performance.

References

Aitchison J (1986). The Statistical Analysis of Compositional Data. Chapman & Hall.
Breiman L (2001). “Random Forests.” Machine Learning, 45, 5–32.
Franzosa EA, Sirota-Madi A, Avila-Pacheco J et al. (2019). “Gut microbiome structure and metabolic activity in inflammatory bowel disease.” Nature Microbiology, 4, 293–305.
Law CW, Chen Y, Shi W, Smyth GK (2014). “voom: Precision weights unlock linear model analysis tools for RNA-seq read counts.” Genome Biology, 15, R29.
Ramos M, Schiffer L, Re A et al. (2017). “Software for the Integration of Multiomics Experiments in Bioconductor.” Cancer Research, 77(21), e39–e42.
Robinson MD, Oshlack A (2010). “A scaling normalization method for differential expression analysis of RNA-seq data.” Genome Biology, 11, R25.
Robin X, Turck N, Hainard A et al. (2011). “pROC: an open-source package for R and S+ to analyze and compare ROC curves.” BMC Bioinformatics, 12, 77.
Rohart F, Gautier B, Singh A, Le Cao KA (2017). “mixOmics: An R package for omics feature selection and multiple data integration.” PLoS Computational Biology, 13(11), e1005752.
Singh A, Shannon CP, Gautier B et al. (2019). “DIABLO: an integrative approach for identifying key molecular drivers from multi-omics assays.” Bioinformatics, 35(17), 3055–3062.
Wright MN, Ziegler A (2017). “ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R.” Journal of Statistical Software, 77(1), 1–17.

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] SummarizedExperiment_1.43.0 Biobase_2.73.1             
#>  [3] GenomicRanges_1.65.1        Seqinfo_1.3.0              
#>  [5] IRanges_2.47.2              S4Vectors_0.51.5           
#>  [7] BiocGenerics_0.59.10        generics_0.1.4             
#>  [9] MatrixGenerics_1.25.0       matrixStats_1.5.0          
#> [11] MultiOmicsBridge_0.99.1     BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] ellipse_0.5.0               gtable_0.3.6               
#>  [3] xfun_0.60                   bslib_0.11.0               
#>  [5] ggplot2_4.0.3               ggrepel_0.9.8              
#>  [7] lattice_0.22-9              vctrs_0.7.3                
#>  [9] tools_4.6.1                 parallel_4.6.1             
#> [11] tibble_3.3.1                rARPACK_0.11-0             
#> [13] pkgconfig_2.0.3             Matrix_1.7-5               
#> [15] RColorBrewer_1.1-3          S7_0.2.2                   
#> [17] mixOmics_6.37.0             lifecycle_1.0.5            
#> [19] stringr_1.6.0               compiler_4.6.1             
#> [21] farver_2.1.2                statmod_1.5.2              
#> [23] codetools_0.2-20            htmltools_0.5.9            
#> [25] sys_3.4.3                   buildtools_1.0.0           
#> [27] sass_0.4.10                 yaml_2.3.12                
#> [29] tidyr_1.3.2                 pillar_1.11.1              
#> [31] jquerylib_0.1.4             MASS_7.3-65                
#> [33] BiocParallel_1.47.0         DelayedArray_0.39.3        
#> [35] cachem_1.1.0                limma_3.69.2               
#> [37] abind_1.4-8                 RSpectra_0.16-2            
#> [39] tidyselect_1.2.1            locfit_1.5-9.12            
#> [41] digest_0.6.39               stringi_1.8.7              
#> [43] purrr_1.2.2                 reshape2_1.4.5             
#> [45] dplyr_1.2.1                 labeling_0.4.3             
#> [47] maketools_1.3.2             fastmap_1.2.0              
#> [49] grid_4.6.1                  cli_3.6.6                  
#> [51] SparseArray_1.13.2          magrittr_2.0.5             
#> [53] S4Arrays_1.13.0             withr_3.0.3                
#> [55] edgeR_4.11.4                corpcor_1.6.10             
#> [57] scales_1.4.0                rmarkdown_2.31             
#> [59] XVector_0.53.0              igraph_2.3.3               
#> [61] otel_0.2.0                  gridExtra_2.3.1            
#> [63] ranger_0.18.0               evaluate_1.0.5             
#> [65] knitr_1.51                  MultiAssayExperiment_1.39.0
#> [67] rlang_1.3.0                 Rcpp_1.1.2                 
#> [69] glue_1.8.1                  BiocManager_1.30.27        
#> [71] pROC_1.19.0.1               jsonlite_2.0.0             
#> [73] plyr_1.8.9                  R6_2.6.1

Integrative Multi-Omics Analysis with MultiOmicsBridge