Persephone

• Additionalfunctions
• Seqc_pipeline
  • Seqc_proc
  • Seqc_output
  • Seqc_input

analyseWBMs(mWBMPath, fluxPath, rxnList, varargin)

analyseWBMs predicts the optimal fluxes for a list of user-defined reactions (rxnList). All predicted are further described in analyseWBMsol.m.

USAGE:

analyseWBMs (mWBMPath, fluxPath, rxnList, varargin)

INPUTS:

• mWBMPath – Path (character array) to the folder with the WBMs

• fluxPath – Path to the directory where the results are stored

• rxnList – Cell array of VMH reactions/metabolites to investigate. Example: rxnList = {‘DM_trp_L[bc], DM_met_L[bc],’Brain_trp_L[c], Heart_met_L[x]}. Note that demand reactions are automatically added if they are not present in the models.

OPTIONAL INPUTS rxnSense Character array containing either ‘max’ or ‘min’

to specify the sense of the objective. Option to specify differently for each objective- to do so provide character array the exact length of rxnList. (OPTIONAL, Default = ‘max’).

numWorkersOptimization Number of workers that will perform FBA in parallel. Note

that more workers does not necessarily make the function faster. It is generally not recommended to set numWorkersOptimization equal to the number of available cores (see: feature(‘numCores’)) as linear solvers can already support multi-core linear optimisation, thus resulting in unnessessary overhead. On computers with 8 cores or less, it is recommended to set numWorkersOptimization to 1. On a computer with 36 cores, an optimal configuration of numWorkersOptimization=6 was found.

saveFullRes Boolean (true/false) indicating if all the complete .v, .y.

, and .w vectors are stored in the result. Default = true. It is recommended to set saveFullRes

paramFluxProcessing Structured array with optional parameters:

.NumericalRounding defines how much the predicted flux values are rounded. A defined value of 1e-6 means that a flux value of 2 + 2.3e-8 is rounded to 2. A flux value of 0 + 1e-15 would be rounded to exactly zero. This rounding factor will also be applied to the shadow price values. If microbiome relative abundance data is provided, the same rounding factor will be applied to the relative abundance data.

Default parameterisation: - paramFluxProcessing.NumericalRounding = 1e-6;

Example: - paramFluxProcessing.NumericalRounding = 1e-6;

paramFluxProcessing.NumericalRounding = 1e-6;

.RxnRemovalCutoff defines the minimal number of samples for which a unique reaction flux could be obtained, before removing the reaction for further analysis. This parameter can be expressed as * fraction: the fraction of samples with unique values, * SD: the standard deviation across samples, and * count: the counted number of unique values. If microbiome relative abundance data is provided, the same removal cutoff factor will be applied to the relative abundance data.

Default parameterisation: - paramFluxProcessing.RxnRemovalCutoff = {‘fraction’,0.1};

Examples: - paramFluxProcessing.RxnRemovalCutoff = {‘fraction’,0.1}; - paramFluxProcessing.RxnRemovalCutoff = {‘SD’,1}; - paramFluxProcessing.RxnRemovalCutoff = {‘count’,30};

paramFluxProcessing.RxnRemovalCutoff = {‘fraction’,0.1};

.RxnEquivalenceThreshold defines the minimal threshold of when functionally identical flux values are predicted, and are thus part of the same linear pathways. The threshold for functional equivalence is expressed as the R2 (r-squared) value after performing a simple linear regression between two reactions.

Default parameterisation: - paramFluxProcessing.RxnEquivalenceThreshold = 0.999;

Example: - paramFluxProcessing.RxnEquivalenceThreshold = 0.999;

paramFluxProcessing.RxnEquivalenceThreshold = 0.999;

.fluxMicrobeCorrelationType defines the method for correlating the predicted fluxes with microbial relative abundances. Note that this metric is not used if mWBMs are not present. The available correlation types are: * regression_r2: the R2 (r-squared) value from pairwised linear regression on the predicted fluxes against microbial relative abundances. * spearman_rho: the correlation coefficient, rho obtained from pairwise Spearman nonparametric correlations between predicted fluxes and microbial relative abundances.

Default parameterisation: - paramFluxProcessing.fluxMicrobeCorrelationMetric = ‘regression_r2’;

Examples: - paramFluxProcessing.fluxMicrobeCorrelationMetric = ‘regression_r2’; - paramFluxProcessing.fluxMicrobeCorrelationMetric = ‘spearman_rho’;

fluxAnalysisPath: Character array with path to directory where all

results will be saved (Default = pwd)

Solver: Validated solvers: ‘cplex_direct’,’ibm_cplex’

‘tomlab_cplex’, ‘gurobi’, ‘mosek’

analyseGF: Boolean indiciating whether or not to investigate

GF models. In the case of personalisation, a germ free iWBM will be created for every sample, this can results in long computation times as there are double the number of models to solve. If you are not interested in solving a germ-free iWBM for each sample, you can set to false. Default is true. When personalisation is skipped, only one germ-free model is made for each sex.

analyseWBMsol(fluxPath, paramFluxProcessing, fluxAnalysisPath, analyseGF)

Load and post-process the FBA solutions produced by analyseWBMs

Loads the FBA solutions produced in analyseWBMs.m, prepares the results for further analysis, and produces summary statistics on the fluxes. It loads the flux solutions and shadow prices, converts the microbial biomass shadow prices to human-readable tables, isolates the microbial component of the fluxes by subtracting the germ-free fluxes in a sex-specific manner, groups reactions with (near-)identical fluxes, and saves all results.

USAGE:

processedFluxResPaths = analyseWBMsol (fluxPath, paramFluxProcessing, fluxAnalysisPath, analyseGF)

INPUTS:

• fluxPath – character array with the path to the .mat files produced in analyseWBMs.m

• paramFluxProcessing – structure with flux-processing parameters (missing fields are filled with defaults). Fields used:

  • .numericalRounding - defines how much the predicted flux and shadow-price values are rounded; e.g. 1e-6 rounds 2 + 2.3e-8 to 2

  • .roundingFactor - number of decimal places used when rounding, derived internally from the numerical rounding value

  • .rxnRemovalCutoff - minimal number of samples for which a unique reaction flux must be obtained before the reaction is removed, given as {‘fraction’, x}, {‘SD’, x} or {‘count’, x}

• fluxAnalysisPath – character array with the path to the directory where all results are saved

• analyseGF – logical, whether the germ-free (microbe-free) flux solutions are also analysed

OUTPUT:

processedFluxResPaths – string array with the paths to the saved processed flux result files

calculateMetrics(data, calculateBrayCurtis)

Calculate diversity and descriptive metrics for a microbiome abundance table

Calculates Pielou’s evenness, the Bray-Curtis dissimilarity, and descriptive statistics for the taxa in a reads or relative-abundance table of a microbiome sample. Developed for MARS for MATLAB and used in runMars.

USAGE:

[bray, pielous, taxonSummary] = calculateMetrics (data, calculateBrayCurtis)

INPUTS:

• data – m x n table with taxa in rows and samples in columns. The first column must be called Taxon and hold (non-duplicated) taxonomic identifiers; the remaining columns hold the per-sample abundances. Fields used:

  • .Properties - table properties; the sample names in .VariableNames label the outputs

  • .Taxon - column of taxonomic identifiers, used as the row names of taxonSummary

• calculateBrayCurtis – logical, true to calculate the Bray-Curtis dissimilarity. Large sample sizes (n > 200) have high computation times

OUTPUTS:

• bray – table with the Bray-Curtis dissimilarity index; rows and columns are samples and the diagonal is 0

• pielous – array with the Pielou evenness (alpha diversity) score

• taxonSummary – table with descriptive statistics on the abundance of the different taxa in the input table

configTemplatePersephone

Configuration template for the Persephone pipeline

Fill in the required and optional variables in this script and pass its path to runPersephone. Running the script builds the paths structure that the pipeline (initPersephone, validatePersephoneInputs, runPersephone) consumes.

createBatchMWBM(mgpipePath, mWBMPath, metadataPath, varargin)

Create a batch of personalised host-microbiome WBM models (mWBMs)

Joins microbiome community models with unpersonalised WBM models in a sex-specific manner, parameterises them on a predefined diet, optionally checks their feasibility, and reports summary statistics on the generated models.

USAGE:

[modelStats, summaryStats, dietInfo, dietGrowthStats] = createBatchMWBM (mgpipePath, mWBMPath, metadataPath, varargin)

INPUTS:

• mgpipePath – path to the microbiome community models created by the microbiome modelling toolbox

• mWBMPath – path to the directory where the host-microbiome models are saved

• metadataPath – path to the metadata file with sample IDs and sex

OPTIONAL INPUTS (name-value pairs in varargin):

Diet: diet option (default ‘EUAverageDiet’) numWorkersCreation: number of cores used for model creation

(default 4)

numWorkersOptimisation: number of cores used for optimisation

(default 2)

checkFeasibility: logical, run ensureWBMfeasibility to check

that the models grow on the diet (default true)

wbmDirectory: path to a directory with user-defined WBMs;

if empty Harvey/Harvetta are used (default ‘’)

solver: LP solver to use (default ‘gurobi’) maleUnpersonalisedWBMpath: male base WBM to use

(default ‘Harvey_1_03d’)

femaleUnpersonalisedWBMpath: female base WBM to use

(default ‘Harvetta_1_03d’)

OUTPUTS:

• modelStats – table with per-model summary statistics (sex, number of reactions, metabolites, constraints, and taxa)

• summaryStats – table with the mean and standard deviation of the statistics in modelStats

• dietInfo – diet growth information returned by ensureWBMfeasibility (empty when the feasibility check is skipped)

• dietGrowthStats – statistics on which models could grow on the given diet (empty when the feasibility check is skipped)

createForestPlot(estimates, ci, names, pValues, plotTitle, xTitle, hideLegend)

Generate a forest plot displaying confidence intervals for a set of estimates

Points are coloured by significance category (FDR < 0.05, P < 0.05, P > 0.05) and drawn against their confidence intervals.

USAGE:

createForestPlot (estimates, ci, names, pValues, plotTitle, xTitle, hideLegend)

INPUTS:

• estimates – vector of estimates (e.g. effect sizes or log fold changes)

• ci – n x 2 matrix of confidence intervals, lower and upper bounds in the two columns

• names – cell array of labels, one per data point

• pValues – vector of p-values corresponding to each estimate

• plotTitle – char/string, title of the plot

• xTitle – char/string, label for the x-axis (e.g. “Effect Size” or “Log Fold Change”)

OPTIONAL INPUT:

hideLegend – logical, if true the legend is hidden (default false)

createMWBM(microbiota_model, WBM_model, Diet, saveDir)

Create a personalised host-microbiome WBM by joining a microbiome and a WBM model

Combines a microbiome community model with an (un)personalised whole-body metabolic model, parameterises the joined model on a predefined diet, adds setup and microbiota-composition information, sets the whole-body objective, and (optionally) saves the model to disk.

USAGE:

modelHM = createMWBM (microbiota_model, WBM_model, Diet, saveDir)

INPUTS:

• microbiota_model – microbiome community model created by the microbiome modelling toolbox. Field used:

  • .name - identifier used to build the model ID

• WBM_model – Harvey or Harvetta whole-body metabolic model. Fields used:

  • .sex - sex of the model (‘male’ or ‘female’)

  • .SetupInfo - setup information copied into the host-microbiome model

  • .version - WBM version, stored in the setup info

  • .modelAnnotation - model annotation, stored in the setup info

• Diet – diet option to constrain the model, e.g. ‘EUAverageDiet’ (default)

OPTIONAL INPUT:

saveDir – char/string, directory to save the model in; if empty the model is not written to disk (default ‘’)

OUTPUT:

modelHM – personalised host-microbiome WBM model

createVolcanoPlot(estimates, pValues, names, plotTitle, xTitle, yTitle)

Generate a volcano plot of regression estimates against their p-values

Points are coloured by significance category (FDR < 0.05, P < 0.05, P > 0.05); significant points are labelled and reference lines are drawn at estimate = 0 and p = 0.05.

USAGE:

createVolcanoPlot (estimates, pValues, names, plotTitle, xTitle, yTitle)

INPUTS:

• estimates – vector of regression estimates for each data point (e.g. effect sizes or log fold changes)

• pValues – vector of p-values corresponding to each estimate

• names – cell array of labels or names for each data point

• plotTitle – char/string, title of the plot

• xTitle – char/string, label for the x-axis (e.g. “Log Fold Change”)

• yTitle – char/string, label for the y-axis (e.g. “-log10(p-value)”)

Example

createVolcanoPlot(estimates, pValues, {‘Metabolite1’, ‘Metabolite2’}, …

‘Volcano Plot’, ‘Log Fold Change’, ‘-log10(p-value)’)

ensureWBMfeasibility(mWBMPath, varargin)

Find a diet that makes all WBMs in a directory feasible and propagate it

Loads the WBMs in mWBMPath and tests which grow on the given diet. If any are infeasible, all diet reactions are opened to test whether growth is possible at all; if so, missing diet components are searched with getMissingDietPersephone (repeated up to 10 iterations). When a feasible diet is found, it is propagated to every model in mWBMPath and each model’s feasibility on the original and updated diet is recorded.

USAGE:

[dietInfo, dietGrowthStats] = ensureWBMfeasibility (mWBMPath, varargin)

INPUT:

mWBMPath – path to the directory where the host-microbiome models are saved

OPTIONAL INPUTS (name-value pairs in varargin):

Diet: diet option (default ‘EUAverageDiet’) solver: LP solver to use (default ‘’, i.e. the current solver) numWorkers: number of workers used for optimisation (default 1)

OUTPUTS:

• dietInfo – structure with the missing diet components (.missingDietComponents) and the updated diet (.updatedDiet); empty if no updated diet was found

• dietGrowthStats – table indicating which WBMs could not grow on the original diet, on any diet, and on the updated diet

fdrBHadjustment(p)

Perform Benjamini-Hochberg FDR correction on a vector of p-values

Returns the FDR-corrected p-values for the input vector, in the same order as the input. The result matches mafdr(p, ‘BHFDR’, true) from the Bioinformatics Toolbox. The procedure follows Benjamini, Y. and Hochberg, Y. (1995), Controlling the false discovery rate: a practical and powerful approach to multiple testing.

USAGE:

q = fdrBHadjustment (p)

INPUT:

p – vector of p-values, assumed to lie in the interval [0, 1]

OUTPUT:

q – vector of FDR-corrected p-values, in the same order as p

Example

pvals = [0.01, 0.04, 0.03, 0.002, 0.07]; qvals = fdrBHadjustment(pvals);

Note

Part of this function was drafted with assistance from an LLM and then manually tested and corrected.

findPvalCategories(pValues)

Categorise a list of p-values into significance groups for plotting

A maximum of three groups are defined: FDR < 0.05, P < 0.05, and P > 0.05. Groups that are not present in the data are removed from the outputs. The FDR values are computed from the p-value list with the Benjamini-Hochberg method.

USAGE:

[groupsToPlot, colours] = findPvalCategories (pValues)

INPUT:

pValues – vector of p-values to categorise

OUTPUTS:

• groupsToPlot – table with one logical column per present category, flagging the p-values that belong to each category

• colours – cell array of RGB triplets, one per present category, for the category-associated plotting colours

generateStackedBarPlot_PhylumMARScoverage(input_relAbundances_preMapping, input_relAbundances_postMapping, saveDir, varargin)

Generate stacked bar plots of phylum relative abundances pre- and post-MARS mapping

Plots the mean relative abundances of phyla before and after mapping to a model database in MARS, and saves the figure to disk.

USAGE:

generateStackedBarPlot_PhylumMARScoverage (input_relAbundances_preMapping, input_relAbundances_postMapping, saveDir, varargin)

INPUTS:

• input_relAbundances_preMapping – char/string, path to the table with phyla and their mean relative abundances pre-mapping (a standard MARS output)

• input_relAbundances_postMapping – char/string, path to the table with phyla and their mean relative abundances post-mapping (a standard MARS output)

• saveDir – char/string, directory where the stacked bar plot is saved

OPTIONAL INPUT (name-value pair in varargin):

mappingDatabase_name: char/string, name of the model database used for

mapping, shown in the plot title (default ‘’)

getMissingDietPersephone(inputModel, missingDietComponents, testInitialFeasibility)

Identify dietary components that are missing for an infeasible WBM or microbiome model

Determines missing dietary compounds in host-microbiome WBMs (or in microbiome community models) that are infeasible, assuming the WBM and microbiome models are individually feasible. The model is first tested for feasibility (this can be skipped). If infeasible, all inactive diet exchange reactions are opened; if still infeasible the function stops. Otherwise diet exchanges are randomly closed in decreasing batch sizes (50, 10, 5, then 1) and reactions not required for feasibility are discarded, leaving the set of components that are needed. Because the batches are drawn at random, repeated runs may return slightly different sets.

USAGE:

missingDietComponents = getMissingDietPersephone (inputModel, missingDietComponents, testInitialFeasibility)

INPUTS:

• inputModel – host-microbiome or microbiome community model structure. Fields used:

  • .rxns - reaction identifiers, searched for diet and biomass exchange reactions

  • .osenseStr - objective sense, set to ‘max’ before optimisation

• missingDietComponents – list of diet exchange reactions already known to be missing (e.g. from a previous run); pass ‘’ or omit if none are known

OPTIONAL INPUT:

testInitialFeasibility – logical/double, whether the initial feasibility test is performed (default 1)

OUTPUT:

missingDietComponents – list of missing diet components. When a list was supplied as input, the output combines it with the newly discovered components

initPersephone(resultPath, paths)

Initialise the Persephone pipeline: check toolboxes, metadata, and folders

Checks for the required MATLAB toolbox dependencies, reads and validates the metadata file (harmonising the ID and Sex columns), writes a processed metadata file, and creates the output directory structure used by the pipeline.

USAGE:

[initialised, statToolboxInstalled, updatedMetadataPath] = initPersephone (resultPath, paths)

INPUTS:

• resultPath – char/string, path to the main results directory

• paths – structure with the Persephone configuration. Fields used:

  • .General - general settings; .metadataPath gives the metadata file

  • .mgPipe - mgPipe settings (.computeProfiles, .outputPathMgPipe)

  • .seqC - SeqC output settings (.outputPathSeqC)

  • .Mars - MARS output settings (.outputPathMars)

  • .persWBM - personalisation output settings (.outputPathPersonalisation)

  • .mWBM - host-microbiome output settings (.outputPathMWBM)

  • .fba - flux output settings (.outputPathFluxResult, .outputPathFluxAnalysis)

  • .stats - statistics output settings (.outputPathStatistics)

OUTPUTS:

• initialised – logical, true if initialisation was successful

• statToolboxInstalled – logical, true if the Statistics and Machine Learning Toolbox is installed and licensed

• updatedMetadataPath – char, path to the processed metadata file

loadMinimalWBM(modPath)

Load the smallest combination of WBM model fields needed to perform FBA

Loading only the minimal set of fields can reduce model loading times by roughly 6x compared with loading the full model. If any of the required fields are missing from the stored file, the complete model is loaded instead.

USAGE:

model = loadMinimalWBM (modPath)

INPUT:

modPath – char or string, path to the .mat file that stores the whole-body metabolic (WBM) model

OUTPUT:

model – struct, the loaded WBM model. It contains at least the fields ID, S, ub, lb, rxns, mets, c, d, csense, dsense, osenseStr, and C (with sex and SetupInfo added when the minimal load succeeds)

performRegressions(data, metadata, formula, exponentiateLogOdds)

Perform linear or logistic regressions on flux or microbial abundance data

Runs a regression per reaction or taxon, using the supplied formula. The formula must contain either “Flux” or “relative_abundance” as the predicted quantity. Control variables and moderators are supported, and mixed-effect (random-effect) regressions are used when the formula encodes a random effect. Binary responses trigger logistic regression, continuous responses linear regression.

USAGE:

[results, regressions] = performRegressions (data, metadata, formula, exponentiateLogOdds)

INPUTS:

• data – table with flux or microbiome data. Fields used:

  • .Properties - table properties; .VariableNames gives the reaction/taxa columns

  • .ID - sample identifier column (must be the first column)

• metadata – metadata table. Field used:

  • .ID - sample identifier column (must be the first column)

• formula – regression formula; must contain “Flux” or “relative_abundance” as the predicted quantity

OPTIONAL INPUT:

exponentiateLogOdds – logical, whether a logistic-regression estimate is exponentiated to obtain odds ratios (default false)

OUTPUTS:

• results – structure with one results table per predictor, holding the estimates, confidence intervals, p-values, FDR values and adjusted R2 for each reaction/taxon

• regressions – structure with the fitted regression model object for each reaction/taxon (empty entries where no fit was found)

performStatsPersephone(statPath, pathToProcessedFluxes, metadataPath, response, varargin)

Perform regression or non-parametric statistics on processed fluxes and abundances

Runs statistical analyses on processed flux and gut-microbiome relative abundance data for large cohorts (minimum ~50 samples). The analysis depends on the inputs: binary response without confounders uses Wilcoxon tests, binary response with confounders uses logistic regression, continuous response without confounders uses simple linear regression, and continuous response with confounders uses multiple linear regression. Fluxes are z-transformed and relative abundances log-transformed before testing.

USAGE:

results = performStatsPersephone (statPath, pathToProcessedFluxes, metadataPath, response, varargin)

INPUTS:

• statPath – path (character array) to the working directory

• pathToProcessedFluxes – path to the processed flux data

• metadataPath – path to the metadata file

• response – char/string, name of the response variable

OPTIONAL INPUTS (name-value pairs in varargin):

pathToWbmRelAbundances: path to the microbial relative abundances

(default [])

confounders: cell array with the names of confounding

variables to include (default ‘’)

microbeCutoff: minimal fraction of samples in which a

microbe must be present to be analysed (default 0.1)

alternativevolcanoPlotTitle: char/string, alternative title for the

volcano plot (default ‘’)

OUTPUT:

results – structure with the statistical result tables produced for the fluxes and relative abundances

persWBM(metadata, varargin)

Create physiologically personalised whole-body models (WBMs) from metadata

Takes physiological parameters for one or several individuals and adjusts the parameters of a provided WBM (or of Harvey/Harvetta) to create a personalised WBM. Physiological parameters are computed from whatever data are available (e.g. cardiac output is used directly if given, otherwise derived from stroke volume). Details of every calculation are stored in the personalised model output and in the accompanying Excel file.

USAGE:

[iWBM, iWBMcontrol_female, iWBMcontrol_male, persParams] = persWBM (metadata, varargin)

INPUT:

metadata – physiological parameters, either as a struct (to personalise one model) or a path to an Excel file (batch personalisation). When a struct, field used:

• .sex - string, “male” or “female”

OPTIONAL INPUTS (name-value pairs in varargin):

persPhysiology: cell array of the non-metabolite parameters in the

metadata to personalise the model with (default {})

femaleWBM: a female WBM, or a path to it; if empty Harvetta is

loaded (default ‘’)

maleWBM: a male WBM, or a path to it; if empty Harvey is loaded

(default ‘’)

resPath: path in which to store the personalised model and

outputs (default the current directory)

persMetabolites: cell array of metabolites to also personalise from

metabolomic data; if empty this step is skipped (default {})

Diet: diet option (default ‘EUAverageDiet’) solver: LP solver to use (default ‘’)

OUTPUTS:

• iWBM – personalised model with updated physiological parameters (described in model.IndividualisedParameters)

• iWBMcontrol_female – unpersonalised control WBM for female subjects (a copy of the supplied model or of Harvetta)

• iWBMcontrol_male – unpersonalised control WBM for male subjects (a copy of the supplied model or of Harvey)

• persParams – details of the updated parameters and how they were calculated (also written to an Excel file)

persWBMmetabolomics(sex, metabolomicParams, varargin)

Create metabolomically personalised whole-body models (WBMs) from metabolite data

Takes metabolomic parameters for one or several individuals and adjusts the metabolite constraints of a provided WBM (or of Harvey/Harvetta) to create a personalised WBM. When supplied as a table, the sample and metabolite names are read from metabolomicParams.Properties.VariableNames.

USAGE:

[iWBM, iWBMcontrol, personalisationOverview] = persWBMetabolomics (sex, metabolomicParams, varargin)

INPUTS:

• sex – string, sex of the subject(s) (“male” or “female”)

• metabolomicParams – metabolomic parameters as a cell array, a table, or a path to an Excel file. Each column holds the metabolite concentrations for one individual, with the compartment and unit given in the leading rows

OPTIONAL INPUTS (name-value pairs in varargin):

iWBM: an already physiologically personalised WBM to further

personalise (default ‘’)

iWBMcontrol: a control WBM to carry through unchanged (default ‘’) resPath: path in which to store the personalised model and outputs

(default the current directory)

Diet: diet in the form of a text file or named .mat file from

the COBRA toolbox (default ‘EUAverageDietNew’)

OUTPUTS:

• iWBM – model with updated metabolite constraints (updated parameters described in model.IndividualisedParameters)

• iWBMcontrol – control WBM with no personalised adjustments

• personalisationOverview – table summarising, per individual and biofluid, the metabolite min/max concentrations applied

readMetadataForPersephone(metadataPath)

DESCRIPTION:

This function loads a metadata file into a MATLAB table. It accounts for: - Preserving the full variable names by disabling automatic truncation. - Capturing variable units if provided in the second row of the file. - Storing the captured units in the ‘VariableUnits’ property of the table.

USAGE:

metadata = readMetadataForPersephone (metadataPath)

INPUTS:

metadataPath – A string specifying the path to the metadata file (in CSV or XLSX format).

OUTPUTS:

metadata – A MATLAB table containing the processed metadata. Variable units, if present in the second row of the file, are stored in the table’s ‘VariableUnits’ property.

Notes

• If variable names exceed 64 characters, they are truncated, but the true variable names are preserved in the ‘VariableDescriptions’ property.

• Warnings related to variable name truncation are suppressed to avoid unnecessary console output.

AUTHORS

Tim Hensen, January 2025

runMars(readsTablePath, varargin)

This function processes microbiome taxonomy and read abundance data and maps microbial species on a microbial reconstruction database, such as AGORA2 and APOLLO.

USAGE:

runMars (readsTablePath, varargin)

INPUTS:

readsTablePath – String; path to the reads abundance file. If

taxonomic assignment is not present in this file, provide taxonomy in taxaTablePath.

cutoffMars: Numeric; value under which individual taxa

relative abundances are considered to be zero. Optional, defaults to 1e-6

flagLoneSpecies: Boolean; indicates if the genus name is in the

name of the species. E.g., if the species name is Prevotella copri, set to false. If the species name is copri set to true. Optional, defaults to true.

taxaDelimiter: String; delimiter used to separate taxonomic

levels. Optional, defaults to ;

removeClade: Boolean; specifies to remove clade name

extensions from all taxonomic levels of microbiome taxa. If set to false, MARS might find significantly less models in AGORA2 and APOLLO databases, as clade extensions are not included there. Optional, defaults to true.

reconstructionDb: String; defining if AGORA2, APOLLO, a

combination of both (full) or a user-defined database should be used as the model database to check presence in. Allowed Input (case-insensitive): “AGORA2”, “APOLLO”, “full_db”, “user_db”. Optional, defaults to full_db.

userDbPath: String; The path to the user database if

reconstructionDb is set to user_db. Optional, defaults to ‘’.

sampleReadCountCutoff: Numeric; value for total read counts per sample

under which samples are excluded from analysis. If the reads table is already normalised, set this value to 0 or 0.1 to ensure that samples are not removed.

taxaTablePath: String; path to the file where taxonomies are

matched to taxonomic unit (OTU/ASV etc.) This requires that the taxonomic unit is present in the reads table as well in order to match the taxonomy with the reads. Ensure that the column with the taxonomic unit has the same header in both the reads table and taxonomy table and that the column with the taxonomies in the taxonomy table is called ‘taxon’. Optional, defaults to ‘’.

outputPathMars: String; path to the directory where the output

of MARS is stored. Optional, defaults to [pwd, filesep, ‘resultMars’].

calculateBrayCurtis: Boolean; Specifies if the Bray-Curtis

dissimilarity index is calculated. Caution! Putting this to true can greatly increase the time MARS needs to run. Defaults to false.

compoundedDatabase: Boolean, indicates if the database reads are

compounded or not. E.g., if the total number of reads for an order o__x is the same as the genus g__x. This means the reads assigned to g__x are in the number of reads of o__x. This requires a different data handling as compared where o__x only has reads that could be assocated up to that level. Defaults to false.

OUTPUTS:

• The function does not return variables but writes processed results

• to the specified output directory in the MARS pipeline.

AUTHOR: Tim Hensen, July 2025

Bram Nap, August 2025

runNonparametricTests(data, metadata, predictor, response)

Perform non-parametric statistical tests (Wilcoxon or Kruskal-Wallis) on data

The test is selected automatically from the number of response groups: a Wilcoxon rank-sum test for two groups, a Kruskal-Wallis test for more. Data are reshaped to long format, missing values removed, and the results FDR-corrected with the Benjamini-Hochberg method.

USAGE:

resultTable = runNonparametricTests (data, metadata, predictor, response)

INPUTS:

• data – m x n table of flux or abundance values (m samples, n reactions/taxa). Fields used:

  • .Properties - table properties; .VariableNames must include an ID column and the measurement columns

  • .ID - sample identifier column (first column)

• metadata – table of sample metadata. Fields used:

  • .Properties - table properties; .VariableNames must include an ID column and the response column

  • .Response - response column derived from response and converted to group indices

  • .ID - sample identifier column (first column)

• predictor – char/string, name of the predictor variable (‘Flux’ or ‘relative_abundance’)

• response – char/string, name of the response variable in metadata

OUTPUT:

resultTable – table of statistical results (one row per reaction/taxon) with the test method, group sizes, test statistic, p-value, FDR-adjusted p-value and effect size

runPersephone(configPath)

Orchestrate the end-to-end Persephone human-microbiome WBM pipeline

Runs the full pipeline for constructing and analysing human-microbiome whole-body models: sequencing data processing (SeqC), metagenomic mapping (MARS), microbiome community model generation (mgPipe), WBM personalisation, host-microbiome model creation, flux balance analysis, and statistical evaluation. Which sections run is controlled by the config file, and progress is recorded so the pipeline can resume.

USAGE:

progress = runPersephone (configPath)

INPUT:

configPath – path to a .mat file holding all inputs needed to run some or all sections of the pipeline. A template is at src/analysis/persephone (configTemplatePersephone.m) or at https://vmh2.life/persephone

OUTPUT:

progress – structure recording which pipeline steps have completed and their run times; it is also saved to progress.mat in the results directory so a run can be resumed

runSeqC(repoPathSeqC, outputPathSeqC, fileIDSeqC, procKeepSeqC, maxMemSeqC, maxCpuSeqC, maxProcSeqC, debugSeqC, runApptainer)

Build and run the SeqC pipeline from MATLAB as a wrapper around its container

Builds and runs the SeqC Docker (or Apptainer) image from a MATLAB environment, ensures the necessary databases exist, and coordinates execution with the MARS pipeline. Requires Docker (and, optionally, Apptainer) to be installed and on the system path.

USAGE:

status = runSeqC (repoPathSeqC, outputPathSeqC, fileIDSeqC, procKeepSeqC, maxMemSeqC, maxCpuSeqC, maxProcSeqC, debugSeqC, runApptainer)

INPUTS:

• repoPathSeqC – char, path to the SeqC repository

• outputPathSeqC – char, path for the SeqC output

• fileIDSeqC – char, unique identifier for file processing

• procKeepSeqC – logical, keep all intermediate files (true/false)

• maxMemSeqC – integer, maximum memory allocation for SeqC

• maxCpuSeqC – integer, maximum CPU allocation for SeqC

• maxProcSeqC – integer, maximum number of processes for SeqC

• debugSeqC – logical, enable debug mode (true/false)

• runApptainer – logical, wrap execution with Apptainer (true/false)

OUTPUT:

status – double, 0 if SeqC completed successfully, 1 otherwise

validatePersephoneInputs(paths, resultPath)

Validate the inputs supplied in the paths structure for the Persephone pipeline

Ensures that the required fields, flags, and parameters for each pipeline component (SeqC, MARS, mgPipe, WBM personalisation, host-microbiome creation, FBA, and statistics) are present, correctly typed, and within their expected constraints. Validation errors are raised by the underlying validateattributes/validatestring calls.

USAGE:

validated = validatePersephoneInputs (paths, resultPath)

INPUTS:

• paths – structure with all Persephone input fields. Fields used:

  • .General - shared settings (.solver, .diet, .metadataPath, worker counts)

  • .seqC - SeqC flag and parameters

  • .Mars - MARS flag and parameters

  • .mgPipe - mgPipe flag and parameters

  • .persWBM - WBM personalisation flag and parameters

  • .mWBM - host-microbiome creation flag and parameters

  • .fba - FBA flag and parameters (including .rxnList and .paramFluxProcessing)

  • .stats - statistics flag and parameters

• resultPath – char/string, path to the results directory (validated as non-empty)

OUTPUT:

validated – logical scalar, true if validation succeeded