Repository accompanying the manuscript:
Brown S, Xiao B, Klein K, Dupuis J, Zhang Q and Greenwood G. "MOLA: a novel topological data analysis framework for analyzing multiomic loops in precision medicine." bioRxiv, 2026.
A previous multiomic analysis framework (Gurnari, 2025) demonstrated the utility of topological data analysis (TDA) for identifying subtle non-linear structure in high-dimensional omic data, i.e., coordinated loop patterns of gene omic profiles represented as gene-valued vectors. However, methods for interpreting and relating these loop features to sample-level characteristics remain limited.
To address this gap, we developed MOLA (MultiOmic Loop Analysis), a framework for the characterization and interpretation of multiomic loops. MOLA extends previous workflows with tools for:
- Loop visualization for exploratory analysis and interpretation
- Loop filtering to prioritize well-sampled structures
- Epigenetic alterations analysis to quantify associations between loops and sample-level characteristics
- Loop Fourier enrichment analysis to identify pathways that cluster around loops
MOLA is designed for the analysis of multi-sample multiomic datasets that include sample-level characteristics. By linking topological structure to biological and clinical variables, MOLA enables the discovery of complex signals that may be overlooked by conventional linear approaches, with anticipated applications in precision medicine, biomarker discovery, and risk modeling.
The initial steps of all MOLA analyses are:
(A) Identification and visualization of loop structures. Left: a single multiomic sample, projected onto three (RNA-seq, promoter methylation, and enhancer methylation) of the total 13 omic dimensions in the dataset, with each point representing a gene. Using the full 13-dimensional dataset, a gene–gene correlation distance matrix is computed (center-left) and provided as input to the (Gurnari, 2025) framework, which generates persistence diagrams and associated (harmonic) loop weight vectors. In the persistence diagram (center), each point represents a loop: the x-axis denotes the birth radius, the correlation distance at which genes connect around the loop circumference, while the y-axis denotes the death radius, the distance at which genes connect across the loop and fill it in. For a highlighted loop, the multiomic sample can be coloured by the (Gurnari, 2025) loop gene vector (center-right), and genes can be projected onto a 2D embedding of the loop (right).
(B) Loop filtering procedure and downstream analysis approaches. A transformed loop statistic based on the death-to-birth ratio is computed for all loops across all samples (left), the distribution of which follows a standard Gumbel distribution (left-center, see Bobrowski, 2023). Only loops with unusually large death/birth ratios, including the highlighted example, are retained for downstream modeling analysis across samples (right). We provide two options for such analyses, indicated by the double arrow: evaluating the association between loop gene weights and sample-level characteristics, either by predicting the loop weight of gene
The methods we applied for studying both directions of associations were:
-
Epigenetic alteration analysis. For each gene, a linear model is fit to the retained loop patterns predicting that gene's logit-transformed loop weight from the sample-level characteristics; coefficient t-values for each characteristic are aggregated (e.g., all "age" t-values across loops) and the ordered gene statistics are input to GSEAs such as GO, KEGG or Reactome.
-
Survival modelling. Cox regression models of clinical characteristics and gene weights identify genes significantly associated with clinical outcomes and prognosis.
In datasets where samples are dominated by a single loop (as was the case in our breast cancer study) our Fourier enrichment approach can detect which pathways/terms cluster in the dataset:
-
Left: A breast cancer dataset loop in a two-dimensional embedding, where each point represents a gene assigned an angular coordinate (
$\theta_{g_i}$ ). Genes in the endometrial cancer pathway (highlighted in red) cluster near$\theta_{g_i} = 180^\circ$ . -
Center-left and center-right: Pathway membership and
$\theta_{g_i}$ values are used to fit a Fourier model (intercept,$\cos\theta_{g_i}$ , and$\sin\theta_{g_i}$ terms) and a baseline intercept-only model. A likelihood ratio test (LRT) compares the models, testing the hypothesis that the sin and cos coefficients are both 0, and resulting pathway p-values are the outcomes of interest. -
Right: Pathway p-values are corrected for multiple testing to identify Fourier-enriched pathways. For this loop, the endometrial cancer pathway is significantly enriched.
git clone https://github.com/shaelebrown/MOLA.git
cd MOLA- R >= 4.4.2
- Tested on Mac OSX and Linux systems - the maTilDA dependency doesn't seem to build on Windows.
- Significant RAM memory for larger numbers (> 500) of genes.
At startup, mola.R installs any missing R package dependencies and configures the required reticulate Conda environment. Before running mola.R, ensure that Python 3.13 is installed on your system.
You must also clone the maTilDA repository into the same parent directory as the MOLA repository.
Here is a minimal example of how to identify loops, and their corresponding gene weight vectors, in several multiomic datasets. This code assumes that the directory ../multiomic_datasets contains files sample1.csv, sample2.csv, ..., sample50.csv, where all rows correspond to the same set (and order) of genes and the columns correspond to the same omic measures (e.g., "RNAseq", "ATACseq", etc.), of which there must be 3 or more.
source("mola.R")
# process all samples
lapply(X = 1:50, FUN = function(X){
# read in that sample's data
# the genes can be the rownames, but should not be in a column "X"
df <- read.csv(paste0("../multiomic_datasets/sample", X, ".csv"))
df$X <- NULL
# assume directory "../mola_results" already exists
results_dir <- paste0("../mola_results/sample", X)
# run MOLA and save results
# keep all loops to jointly filter across all samples
# if data contains missing values a message will be printed with the
# minimum number of shared non-missing entries in all pairs of rows
# in the dataset
run_mola(df = df, results_dir = results_dir, alpha = 1.0, verbose = T)
})
Results for sample X will be saved in the directory ../mola_results/sampleX, with the following files:
-
gene_weights.RData- the gene weights for each identified (and retained if$\alpha < 1$ ) loop. Loading this file will create a variable calledpoint_weights_mat, in which rows are the genes (in the same row order as all the samples), columns are the loop (named by loop key) and entries are the weights. -
gene_interaction_weights.RData- the gene interaction (i.e., edge) weights for each loop. -
persistence_diagram.csv- the computed persistence diagram with loop key as a column. -
params.RData- the alpha and enclosing radius used in calculations. -
log.txt- informative messages and diagnostics in the case of errors.
If no loops were found in the dataset, only the files log.txt and params.txt will be generated.
To filter loops jointly across all samples, you can use the filter_loops function, like so:
# read in the diagrams for all samples
joint_diagram <- do.call(rbind, lapply(X = 1:50, FUN = function(X){
fname <- paste0('../mola_results/sample', X, '/persistence_diagram.csv')
if(file.exists(fname))
{
diag <- read.csv(paste0('../mola_results/sample', X, '/persistence_diagram.csv'))
diag$samp <- X
}else
{
diag <- data.frame(birth = numeric(), death = numeric(), p_value = numeric(), key = numeric(), samp = numeric())
}
return(diag)
}))
# filter loops - subsets joint_diagram for loop p-values < 0.05
filtered_diagram <- filter_loops(diagram = joint_diagram, alpha = 0.05)
# the columns "samp" and "key" will give the necessary identifiers for
# using retained loops
weights <- do.call(cbind, lapply(X = 1:nrow(filtered_diagram), FUN = function(X){
samp <- filtered_diagram$samp[[X]]
key <- filtered_diagram$key[[X]]
# read in that sample's point_weights_mat and subset for the key
load(paste0("../mola_results/sample", samp, "/gene_weights.RData"))
harmonic_weights <- point_weights_mat[, as.character(key)]
harmonic_weights <- harmonic_weights/max(harmonic_weights)
return(matrix(data = harmonic_weights, ncol = 1))
}))
The plot_loop function can aid in the visualization of loop patterns:
# visualize the first retained loop in filtered_diagram
# get sample number and loop key
samp <- filtered_diagram$samp[[1]]
key <- filtered_diagram$key[[1]]
# read in the dataset for that sample
df <- read.csv(paste0('../multiomic_datasets/sample', samp, '.csv'))
df$X <- NULL
# load the weight vector
loop_weight_vec <- weights[, 1L]
# plot visualization
visualize_loop(df = df, weight_vector = loop_weight_vec)
To carry out an epigenetic alteration analysis you will need to modify the epigenetic_alteration_analysis function to adjust for your own sample-level characteristics and modelling procedures - outputs will be a file tvals.RData stored in a user-supplied results directory:
# read in sample-level characteristics file
pheno <- read.csv('../pheno.csv')
pheno$ancestry <- factor(pheno$ancestry, levels = c("European American", "African American"))
pheno$disease <- factor(pheno$disease, levels = c("Non-infected", "Flu-infected"))
# match with retained loops, assuming "samp" is the sample ID column in
# both dataframes
pheno <- merge(filtered_diagram, pheno, by = c("samp"))
# read in or define gene vector, corresponding to the rows of "weights"
genes <- c("ERLIN2", "CLIC6", "SCGB2A2", "GSTM1", "ZNF703",
"CDC42BPA", "ANKRD30A", "HSPA12B", "ASH2L", "IKZF1",
"VSTM2A", "BAG4", "PSMG3", "TACC1", "TAT",
"EPO", "FGF3", "TFAP2B", "SHH", "CCND1",
"ADIPOQ", "PCDHGB6", "ANO1", "CEACAM5", "MUC5B",
"HOXB2", "MT1G", "CALML5", "STARD3", "PGAP3",
"TCAP", "HOXB13", "HSPB6", "GRIA2", "IKZF3",
"NRG1", "GRB7", "CRABP1", "BMPR1B", "WSCD2",
"MYEOV", "BEST4", "GFRA1", "CDK12", "PSMD3",
"PCBP1", "A2ML1", "CLCA2", "GIMAP7", "FBXL20")
# define results directory to save the output file (results_dir/tvals.RData)
results_dir <- "../mola_results"
# run the epigenetic alteration analysis
epigenetic_alteration_analysis(weights = weights, pheno = pheno, genes = genes, results_dir = results_dir)
# load the results
load(paste0(results_dir, '/tvals.RData'))
# variables are named lists of gene t-values, each in decreasing order
# these lists can be input to GSEA analyses
age_tvals[1:5] # these variables depend on your modifications
For a survival analysis you will need to modify the survival_analysis function to adjust for your own variables and modelling procedures - outputs will be a list, either empty if no genes are significant after FDR correction (at the
# assuming pheno contains survival data columns PFI and PFI.time:
# (but adjust for your own variables)
significant_genes <- survival_analysis(surv_data = pheno, weights = weights, genes = genes)
# inspect
names(significant_genes)
significant_genes[["RET"]]
The fourier_enrichment function takes as input a GSEA term-gene binary membership matrix (generated by the generate_term_map function), and produces a list of FDR adjusted p-values for whether the genes belonging to each pathway form one strong cluster somewhere around the loop (lower values indicate stronger evidence of clustering):
# first generate a term map (this can take a while..)
# i.e., the binary term-gene membership matrix
# note that this function may need to be modified if
# some of the genes do not map to entrez IDs (error may be thrown)
term_map <- generate_term_map(genes = genes)
# filter for terms that have at least a few genes
term_map <- term_map[which(rowSums(term_map) >= 3),]
# run Fourier enrichment for one loop
samp <- filtered_diagram$samp[[1]]
key <- filtered_diagram$key[[1]]
df <- read.csv(paste0('../multiomic_datasets/sample', samp, '.csv'))
df$X <- NULL
# the rows of df should correspond exactly to the genes vector
load(paste0('../mola_results/sample', samp, '/gene_weights.RData'))
w <- point_weights_mat[, as.character(key)]
w <- w/max(w)
loop1_fourier <- fourier_enrichment(df = df, genes = genes, weight_vector = w, term_map = term_map)
# inspect the results
loop1_fourier$pvals_adjusted
# run Fourier enrichment for another loop, using the precomputed
# id2entrez vector
samp2 <- filtered_diagram$samp[[2]]
key2 <- filtered_diagram$key[[2]]
df2 <- read.csv(paste0('../multiomic_datasets/sample', samp2, '.csv'))
df2$X <- NULL
# the rows of df should correspond exactly to the genes vector
load(paste0('../mola_results/sample', samp2, '/gene_weights.RData'))
w2 <- point_weights_mat[, as.character(key2)]
w2 <- w/max(w2)
loop2_fourier <- fourier_enrichment(df = df2, genes = genes, weight_vector = w2, term_map = term_map, id2entrez_vec = loop1_fourier$id2entrez_vec)
Runs the MOLA pipeline on a gene-level multiomic dataset.
Arguments
| Parameter | Description |
|---|---|
df |
Multiomic data frame with one row per gene. |
results_dir |
Directory where all output files will be written. Defaults to the parent directory, /... |
alpha |
p-value threshold for loop filtering. Defaults to 1, retaining all loops so filtering can be performed after joint-sample processing. |
Returns
No object is returned. Results are written to files in results_dir.
Filters a joint-sample persistence diagram to retain statistically significant loops using the universal null procedure (Bobrowski, 2023).
Arguments
| Parameter | Description |
|---|---|
diagram |
Persistence diagram containing at least the birth and death columns. |
alpha |
p-value threshold for significance. Default value is 0.05. |
Returns
A filtered persistence diagram (i.e., dataframe) containing only loops with p-values less than alpha.
Visualizes a single loop.
Arguments
| Parameter | Description |
|---|---|
df |
Multiomic dataset used to compute the loop. Must not contain missing values. |
weight_vector |
Loop weight vector with values between 0 and 1. |
Returns
No object is returned; a plot is generated.
Performs an epigenetic alteration analysis on MOLA results.
Arguments
| Parameter | Description |
|---|---|
weights |
Matrix of loop weights, with one row per loop and one column per gene. Typically loaded from a gene_weights.RData file generated by run_mola(). |
pheno |
Data frame containing sample-level characteristics, with one row per loop. |
genes |
Vector of shared gene names corresponding to the columns of weights. |
results_dir |
Directory where the output file (tvals.RData) will be saved. |
Returns
No object is returned; results are saved as tvals.RData.
Performs survival analysis on MOLA-derived loop weights.
Arguments
| Parameter | Description |
|---|---|
surv_data |
Data frame containing sample-level characteristics and survival information. Each row should correspond to a sample. |
weights |
Matrix of loop weights, with one row per loop and one column per gene. |
genes |
Vector of shared gene names corresponding to the columns of weights. |
Returns
A named list containing Cox proportional hazards regression results for each significant gene.
Constructs a gene set membership matrix for Gene Set Enrichment Analysis (GSEA) using GO, KEGG, and Reactome annotations.
Arguments
| Parameter | Description |
|---|---|
genes |
Vector of gene names. |
Returns
A binary membership matrix with pathways/terms as rows and genes as columns. Matrix entries indicate pathway membership (0/1).
Performs Fourier-based pathway enrichment analysis for a single loop.
Arguments
| Parameter | Description |
|---|---|
df |
Multiomic data frame used to compute the loop. Row names must correspond to genes. |
genes |
Vector of gene names. |
weight_vector |
Weight vector for the loop of interest. |
term_map |
Binary pathway membership matrix with pathways as rows and genes as columns. |
id2entrez_vec |
Optional named vector mapping gene identifiers to Entrez IDs. Providing this avoids recomputing the mapping. |
Returns
A named list containing:
- FDR-adjusted enrichment p-values for each pathway.
id2entrez_vec, allowing the mapping to be reused in subsequent analyses.
The MOLA manuscript provided two case study analyses of an Influenza A virus (IAV) dataset and a TCGA breast cancer dataset.
The raw data from the IAV dataset are available from the European Genome-phenome Archive (EGA), under accession numbers EGAD00001008422 and EGAD00001008359. The dataset was introduced in (Aracena et al, 2023).
The TCGA breast cancer dataset analyzed during the current study were accessed via the R package UCSCXenaTools version 1.7.0 (Wang, 2019) on March 12, 2026.
If you use MOLA, please cite:
@article{MOLA2026,
title = {MOLA: a Novel Topological Data Analysis Framework for Analyzing Multiomic Loops in Precision Medicine},
author = {Brown, Shael and Xiao, Bowei and Oros Klein, Kathleen and Tonin, Patricia N. and Dupuis, Josée and Zhang, Qihuang and Greenwood, Celia M.T.},
journal = {bioRxiv},
year = {2026},
doi = {10.1101/2026.xx.xx.xxxxxx},
note = {Preprint}
}This project is licensed under the MIT License. See the LICENSE file for details.
For questions please contact:
- Shael Brown
- [email protected]
The authors acknowledge funding from the DNA-to-RNA (D2R) Foundational Projects (Cycle 1) grant: Modelling individual differences in patterns of epigenetic regulation.
- Initial release accompanying manuscript publication.
Gurnari, D., Guzmán-Sáenz, A., Utro, F. et al. Probing omics data via harmonic persistent homology. Sci Rep 15, 38836 (2025).
Bobrowski, O., Skraba, P. A universal null-distribution for topological data analysis. Sci Rep 13, 12274 (2023).
Aracena KA, Lin YL, Luo K, Pacis A, Gona S, Mu Z, Yotova V, Sindeaux R, Pramatarova A, Simon MM, Chen X, Groza C, Lougheed D, Gregoire R, Brownlee D, Boye C, Pique-Regi R, Li Y, He X, Bujold D, Pastinen T, Bourque G, Barreiro LB. Epigenetic variation impacts individual differences in the transcriptional response to influenza infection. Nat Genet. 2024 Mar;56(3):408-419. doi: 10.1038/s41588-024-01668-z. Epub 2024 Feb 29. PMID: 38424460.
Wang et al., (2019). The UCSCXenaTools R package: a toolkit for accessing genomics data from UCSC Xena platform, from cancer multi-omics to single-cell RNA-seq. Journal of Open Source Software, 4(40), 1627, https://doi.org/10.21105/joss.01627