## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
    collapse  = TRUE,
    comment   = "#>",
    fig.width = 7,
    fig.height = 5,
    message   = FALSE,
    warning   = FALSE
)

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

## ----simulate-----------------------------------------------------------------
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)
dim(mb_counts)

## ----load_host----------------------------------------------------------------
host_se <- loadHostData(host_counts, col_data = col_data)
host_se
assayNames(host_se)

## ----load_mb------------------------------------------------------------------
mb_se <- loadMicrobiomeData(mb_counts, normalization = "CLR",
                             min_prevalence = 0.1)
mb_se
assayNames(mb_se)

# Verify CLR: column means should be ≈ 0
round(colMeans(assay(mb_se, "CLR")), 10)[1:5]

## ----match_samples------------------------------------------------------------
mae <- matchSamples(host_se, mb_se, min_paired = 10)
mae

## ----dim_reduction------------------------------------------------------------
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)

# Loading matrix dimensions (host, microbiome)
dim(dr_result$host_loadings)
dim(dr_result$mb_loadings)

# Explained variance per component
round(dr_result$explained_variance, 4)

## ----biomarker_discovery------------------------------------------------------
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)

## ----biomarker_summary--------------------------------------------------------
# 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)

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

## ----classifier---------------------------------------------------------------
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))
)

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

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

result

## ----accessors----------------------------------------------------------------
# Integration scores: samples x components
head(integrationScores(result))

# Top biomarkers
head(as.data.frame(biomarkers(result))[, c("feature","omics_layer",
                                             "loading_score")])

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

## ----plot_integration, fig.cap="DIABLO sample scores coloured by outcome group. Loading vectors show the top contributing features per omics layer."----
plotIntegration(result, outcome = outcome, n_loading_arrows = 5)

## ----plot_network, fig.cap="Clustered heatmap of Spearman correlations between the top host genes and microbial taxa. Strong positive/negative correlations indicate candidate host-microbe interactions."----
plotBiomarkerNetwork(result, mae, n_host = 15, n_mb = 10)

## ----plot_bar, fig.cap="Mean cross-validated AUC-ROC for host-only, microbiome-only, and joint classifiers. The multi-omics advantage is immediately visible."----
plotClassifierComparison(result, type = "bar")

## ----plot_roc, fig.cap="Overlaid ROC curves from the last cross-validation fold for each classifier configuration."----
plotClassifierComparison(result, type = "roc")

## ----plot_sankey, fig.cap="Feature flow from omics layer through selected biomarkers to outcome classes. Edge width is proportional to loading score."----
plotSankey(result, n_features = 8)

## ----report-------------------------------------------------------------------
generateReport(result, n_top = 8)

## ----real_data, eval=FALSE----------------------------------------------------
# # 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

## ----disease_context, eval=FALSE----------------------------------------------
# # 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")

## ----session------------------------------------------------------------------
sessionInfo()