MiCARA is designed to be flexible depending on your dataset’s completeness and clinical reporting requirements.
validate_inputs() ensures sample ID alignment across
taxonomic matrices, pathway matrices, and sample metadata. It checks
relative vs. absolute abundances, flags high-missingness covariates
(>30%), and records retained confounders.
# Load example dataset included in MiCARA
data("example_taxa", package = "MiCARA")
data("example_pathways", package = "MiCARA")
data("example_metadata", package = "MiCARA")
# Validate inputs and filter unusable covariates
micara_obj <- validate_inputs(
taxa = example_taxa,
pathways = example_pathways,
metadata = example_metadata,
disease_col = "disease"
)
#>
#> ========== MiCARA Input QC ==========
#> Samples retained: 30
#> Diseases retained: 2
#> Diseases removed: None
#> Taxa scale detected: counts
#> Pathway scale detected: counts
#> Covariates retained: age_category, BMI, gender, study_name
#> Covariates removed: None
#> =====================================Covariate missingness can lead to sample loss during linear model fitting (listwise deletion). MiCARA provides impute_metadata() to address this via MICE: **When to use impute_metadata(): Host covariates (e.g., age_category, BMI) have low-to-moderate missingness (<30%). You want to preserve sample size and statistical power.
**When to skip impute_metadata(): Metadata contains zero missing values (NA). Institutional protocols mandate Complete Case Analysis (CCA).
Note: If you skip imputation, downstream functions will check for NAs in requested covariates and drop or warn on problematic variables before model fitting.
run_diff_abundance() runs bias-corrected log-linear
models across taxa, pathways, or both simultaneously using
ANCOM-BC2.
# Run differential abundance on both feature types
diff_results <- run_diff_abundance(
micara_obj = micara_obj,
feature_type = c("taxa", "pathways"),
reference_level = "healthy",
use_robust = TRUE
)
#> [Info] Setting 'healthy' as the reference level for disease comparisons.
#> [Info] Re-leveling 'disease' with reference level: 'healthy'
#>
#> ========== Running ANCOM-BC2 (taxa) ==========
#> Fixed formula: age_category + BMI + gender + disease
#> Random formula: (1 | study_name)
#> Reference level: healthy
#> Method / Alpha: BH / 0.05
#> ---------------------------------------------------
#> You have inserted a phyloseq object. The phyloseq-based methods will be deprecated. We recommend switching to the (Tree)SummarizedExperiment framework and mia R/Bioconductor package.
#> Warning: The group variable has < 3 categories
#> The multi-group comparisons (global/pairwise/dunnet/trend) will be deactivated
#> Warning: The number of taxa used for estimating sample-specific biases is: 20
#> A large number of taxa (>50) is required for the consistent estimation of biases
#> Loading required package: foreach
#> Loading required package: rngtools
#> Conducting sensitivity analysis for pseudo-count addition to 0s ...
#> For taxa that are significant but do not pass the sensitivity analysis,
#> they are marked in the 'passed_ss' column and will be treated as non-significant in the 'diff_robust' column.
#> For detailed instructions on performing sensitivity analysis, please refer to the package vignette.
#> Significant taxa identified per disease comparison (robust flags):
#> Case: 3
#> ---------------------------------------------------
#> [Info] Re-leveling 'disease' with reference level: 'healthy'
#>
#> ========== Running ANCOM-BC2 (pathways) ==========
#> Fixed formula: age_category + BMI + gender + disease
#> Random formula: (1 | study_name)
#> Reference level: healthy
#> Method / Alpha: BH / 0.05
#> ---------------------------------------------------
#> Warning: The group variable has < 3 categories
#> The multi-group comparisons (global/pairwise/dunnet/trend) will be deactivated
#> Warning: The number of pathways used for estimating sample-specific biases is: 10
#> A large number of pathways (>50) is required for the consistent estimation of biases
#> Conducting sensitivity analysis for pseudo-count addition to 0s ...
#> For taxa that are significant but do not pass the sensitivity analysis,
#> they are marked in the 'passed_ss' column and will be treated as non-significant in the 'diff_robust' column.
#> For detailed instructions on performing sensitivity analysis, please refer to the package vignette.
#> Significant pathways identified per disease comparison (robust flags):
#> Case: 2
#> ---------------------------------------------------
# Inspect significant taxa
head(diff_results$taxa$significant)
#> feature disease lfc se p q direction
#> 1 Taxon_1 Case 2.008631 0.09349275 2.799587e-17 1.866391e-16 enriched
#> 2 Taxon_2 Case 2.189056 0.09271032 3.437782e-18 6.875564e-17 enriched
#> 3 Taxon_3 Case 1.737435 0.08229026 1.877619e-17 1.866391e-16 enriched
# Inspect significant pathways
head(diff_results$pathways$significant)
#> feature disease lfc se p q direction
#> 1 Pathway_1 Case 2.028429 0.07869975 4.762844e-19 2.381422e-18 enriched
#> 2 Pathway_2 Case 1.941612 0.06969169 7.245243e-20 7.245243e-19 enrichedAfter identifying robust differentially abundant features, compute_residual_correlations() evaluates taxon–pathway interactions within each disease cohort: 1. Performs Centered Log-Ratio (CLR) transformations on full feature matrices. 2. Filters to disease-specific significant taxa and pathways. 3. Residualizes feature abundances against host covariates (e.g., age_category, BMI, gender, study_name) using QR decomposition. 4. Computes pairwise Spearman correlations on residuals and applies Benjamini-Hochberg (BH) adjustment within each cohort.
# Compute residual correlations per disease
interactions <- compute_residualCLR_correlations(
micara_obj = micara_obj,
diffab_obj = diff_results,
covariates = c("age_category", "BMI", "gender", "study_name"),
q_cutoff = 0.05,
r_cutoff = 0.20,
verbose = TRUE
)
#> [MiCARA] CLR-transforming complete taxa and pathway matrices.
#> [MiCARA] Case: 3 taxa x 2 pathways = 6 tested pair(s); 1 retained link(s).
#>
#> ========== MiCARA Residual Correlation Summary ==========
#> Cohort: Case (Status: analysed)
#> - Samples: 15 | Residual DF: 10
#> - Features Tested: 3 taxa x 2 pathways (6 pairs)
#> - Retained Links: 1 (1 pos, 0 neg)
#> - Covariates Used: age_category, BMI, gender, study_name
#> ----------------------------------------------------------
# inspect available disease cohorts in the output list
names(interactions)
#> [1] "Case"# View cohort-level analysis summary and covariate usage
attr(interactions, "summary")
#> disease status reason n_samples residual_df n_taxa_tested n_pathways_tested
#> 1 Case analysed <NA> 15 10 3 2
#> n_pairs_tested n_significant_links n_positive_links n_negative_links
#> 1 6 1 1 0
#> covariates_used covariates_dropped
#> 1 age_category, BMI, gender, study_name
# View analysis settings and CLR transformation diagnostics
attr(interactions, "settings")
#> $covariates_requested
#> [1] "age_category" "BMI" "gender" "study_name"
#>
#> $cor_method
#> [1] "spearman"
#>
#> $p_adjust_method
#> [1] "BH_within_disease"
#>
#> $q_cutoff
#> [1] 0.05
#>
#> $r_cutoff
#> [1] 0.2
#>
#> $threshold_rule
#> [1] "q < q_cutoff and abs(rho) > r_cutoff"
#>
#> $taxa_pseudocount
#> [1] 13
#>
#> $pathway_pseudocount
#> [1] 36.5
#>
#> $taxa_all_zero_removed
#> character(0)
#>
#> $pathway_all_zero_removed
#> character(0)
#>
#> $min_samples_per_disease
#> [1] 10
#>
#> $min_residual_df
#> [1] 3
#>
#> $retain_all_pairs
#> [1] TRUE
#>
#> $skipped_no_pathways
#> character(0)
#>
#> $skipped_no_taxa
#> character(0)
# Find the first analyzed cohort dynamically
analyzed_cohorts <- names(interactions)[sapply(interactions, function(x) is.list(x) && x$status == "analysed")]
if (length(analyzed_cohorts) > 0) {
first_cohort <- analyzed_cohorts[1]
message("Displaying links for cohort: ", first_cohort)
head(interactions[[first_cohort]]$links)
}
#> Displaying links for cohort: Case
#> taxon pathway rho p q direction
#> 1 Taxon_1 Pathway_1 0.95 6.086258e-08 3.651755e-07 positiveUse plot_interaction_sankey() to display interactive D3
networks connecting significant taxa to functional pathways:
# Interactive visualization
plot_interaction_sankey(
interactions = interactions,
disease = names(interactions)[1], # Target disease cohort (e.g., "UC")
direction = "both", # Options: "positive", "negative", or "both"
top_n_per_pathway = 5, # Top N correlated taxa per pathway
font_size = 12, # Text label size in pixels
node_width = 35 # Width of Sankey block nodes
)
#> $positive
# Optional: Save interactive HTML and static PNG exports
# Create the results directory first
dir.create("results", showWarnings = FALSE, recursive = TRUE)
plot_interaction_sankey(
interactions = interactions,
disease = names(interactions)[1],
save_path = "results/sankey_network"
)
#> Registered S3 method overwritten by 'webshot':
#> method from
#> print.webshot webshot2
#> `google-chrome`, `chromium-browser` and `chrome` were not found. Try setting the `CHROMOTE_CHROME` environment variable to the executable of a Chromium-based browser, such as Google Chrome, Chromium or Brave or adding one of these executables to your PATH.
#> Warning: Failed to export static PNG image (headless browser unavailable): Invalid path to Chrome
#> Note: The interactive .html file was saved successfully.
#> $positive# Apply global FDR correction across all disease cohorts
interactions_global <- compute_global_fdr(interactions)
#>
#> ========== MiCARA Global FDR Summary ==========
#> Method: BH
#> Global q cutoff: 0.05
#> Correlation cutoff: abs(rho) > 0.2
#> Total tested taxon-pathway pairs: 6
#> Globally significant associations: 1
#> Case: 1 global link(s)
#> -----------------------------------------------
# Plot Sankey network using global significance cutoffs
plot_interaction_sankey(
interactions = interactions_global,
disease = names(interactions_global)[1],
significance = "global",
direction = "both",
top_n_per_pathway = 5
)
#> $positivesessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 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] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] doRNG_1.8.6.3 rngtools_1.5.2 foreach_1.5.2 MiCARA_0.99.2 rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 sys_3.4.3 rstudioapi_0.19.0
#> [4] jsonlite_2.0.0 magrittr_2.0.5 TH.data_1.1-5
#> [7] farver_2.1.2 nloptr_2.2.1 fs_2.1.0
#> [10] vctrs_0.7.3 multtest_2.69.0 minqa_1.2.8
#> [13] ANCOMBC_2.15.1 base64enc_0.1-6 webshot_0.5.5
#> [16] htmltools_0.5.9 forcats_1.0.1 energy_1.7-12
#> [19] haven_2.5.5 cellranger_1.1.0 Formula_1.2-6
#> [22] sass_0.4.10 bslib_0.12.0 htmlwidgets_1.6.4
#> [25] plyr_1.8.9 sandwich_3.1-3 rootSolve_1.8.2.4
#> [28] zoo_1.9-1 cachem_1.1.0 networkD3_0.4.1
#> [31] buildtools_1.0.0 igraph_2.3.4 lifecycle_1.0.5
#> [34] iterators_1.0.14 pkgconfig_2.0.3 webshot2_0.1.2
#> [37] Matrix_1.7-6 R6_2.6.1 fastmap_1.2.0
#> [40] rbibutils_2.4.1 digest_0.6.39 Exact_3.3
#> [43] numDeriv_2016.8-1.1 colorspace_2.1-3 ps_1.9.3
#> [46] S4Vectors_0.51.10 Hmisc_5.3-0 vegan_2.7-6
#> [49] httr_1.4.9 mgcv_1.9-4 compiler_4.6.1
#> [52] proxy_0.4-29 withr_3.0.3 doParallel_1.0.17
#> [55] gsl_2.1-9 htmlTable_2.5.0 S7_0.2.2
#> [58] backports_1.5.1 MASS_7.3-66 biomformat_1.41.0
#> [61] gtools_3.9.5 permute_0.9-10 gld_2.6.8
#> [64] chromote_0.5.1 tools_4.6.1 foreign_0.8-91
#> [67] otel_0.2.0 ape_5.8-1 nnet_7.3-21
#> [70] glue_1.8.1 quadprog_1.5-8 promises_1.5.0
#> [73] nlme_3.1-171 grid_4.6.1 checkmate_2.3.4
#> [76] Rtsne_0.17 cluster_2.1.8.3 reshape2_1.4.5
#> [79] ade4_1.7-24 generics_0.1.4 gtable_0.3.6
#> [82] microbiome_1.35.1 tzdb_0.5.0 class_7.3-24
#> [85] websocket_1.4.4 tidyr_1.3.2 data.table_1.18.6.1
#> [88] lmom_3.3 hms_1.1.4 XVector_0.53.0
#> [91] BiocGenerics_0.59.12 pillar_1.11.1 stringr_1.6.0
#> [94] later_1.4.8 splines_4.6.1 dplyr_1.2.1
#> [97] lattice_0.23-1 survival_3.8-12 tidyselect_1.2.1
#> [100] maketools_1.3.2 Biostrings_2.81.9 knitr_1.52
#> [103] reformulas_0.4.4 gridExtra_2.3.1 phyloseq_1.57.0
#> [106] IRanges_2.47.5 Seqinfo_1.3.2 stats4_4.6.1
#> [109] xfun_0.61 expm_1.0-1 Biobase_2.73.2
#> [112] stringi_1.8.9 yaml_2.3.12 boot_1.3-32
#> [115] evaluate_1.0.5 codetools_0.2-20 data.tree_1.2.0
#> [118] tibble_3.3.1 cli_3.6.6 rpart_4.1.27
#> [121] DescTools_0.99.60 Rdpack_2.6.6 processx_3.9.0
#> [124] jquerylib_0.1.4 Rcpp_1.1.2 readxl_1.5.0.1
#> [127] parallel_4.6.1 ggplot2_4.0.3 readr_2.2.0
#> [130] lme4_2.0-6 mvtnorm_1.4-2 lmerTest_3.2-1
#> [133] scales_1.4.0 e1071_1.7-17 purrr_1.2.2
#> [136] crayon_1.5.3 rlang_1.3.0 multcomp_1.4-32