Tutorial One: Creating Community Microbiome Models

1. General setup: Core procedures for setting-up the Persephone environment.
2. Bioinformatic processing of sequence data (SeqC): Takes you from raw sequence data to taxonomy assigned read counts.
3. Micorbial Abundances Retrieved from Sequening data (MARS): Takes you from taxonomy assigned read counts to reative abundances mapped to your chosen microbial reconstruction resource e.g., AGORA2.
4. Microbiome model creation with mgPipe: Takes your mapped relative abundance data and creates personalised community micrcobiome models personalised for each sample in your metagenomic data file.

Step 1: Setup

• ibm_cplex
• tomlab_cplex
• gurobi
• mosek

• resultDir - The path where all the results of this tutorial will be stored. This ensure all your results will be in one place and are thus easily accesible. Make sure the the folder is accesible to you.
• metadataPath - The location of the metadata file. This is used in both the generation of human-microbiome models as well as the analysis of results towards the end of the tutorial.
• readsTablePath - The path to the the taxonomic reads file containing the amount of reads per taxonomic assignment. If the taxonomic assignments are in the reads file the first column needs to be called Taxon. If not the first column of the readsFile and the first column of taxaTable needs to have to same name. See the example files.

• initialised - Boolean, indicates if Persephone was succesfully initialised.
• statToolboxInstalled - Boolean, indicates if the statistics toolbox is installed. Parts of Persephone are skipped if false
• updatedMetdataPath - The path to the updated metadata file.

Step 2: Bioinformatic processing of sequencing data

Requirements:

• Docker Desktop or Engine (tested@ 4.37.1 (178610))
• MATLAB (tested@R2024b)
• SeqC Repository (Bioinformatics pipeline for taxonomic classification)
• Apptainer (OPTIONAL; Only if operating in an HPC environment)

Introduction

Section 1: Environment Preparation

Downloading SeqC

Configuring SeqC Paths

Set the paths to all inputs for seqC

Section 2: Running SeqC

Checking Your Operating System

Estimating Disk Usage

Section 3: Running Docker for SeqC

1. Build the Docker image if it does not exist
2. Run the SeqC pipeline

Step 1: Build the Docker Image

Step 2: Run the Pipeline

• Test to confirm image viability
• Establish required databases
• Proccess sequencing files

Section 4: Managing Output Files

Step 3: Running MARS

• First, if the taxonomic assignments and reads tables are separated - it will merge them into a single dataframe
• It removes clade extensions from all taxonomic levels names (e.g. Firmicutes__A & Firmicutes__B will be merged to Firmicutes) to allow for optimal AGORA2 mapping
• It translates certain species based on a pre-defined list to make sure the species names match the ones in the AGORA2 databases
• It removes all reads associated with a taxonomics identification that does not have information up to the specified taxonomic level that is looked at (e.g., looking at species, only reads with information up to the species level will be retained)
• It maps each taxonomic level (kingdom, phylum,species etc.) to the AGORA2 databases. If the taxonomic identification matches it means there is a model present
• The mapped reads are normalised to obtain relative abundances
• All relative abundances under a certain cutoff are removed and the data is re-normalised. The default cutoff is 1e-6.

• readsTable - The path to the the taxonomic reads file containing the amount of reads per taxonomic assignment. If the taxonomic assignments are in the reads file the first column needs to be called Taxon. If not the first column of the readsFile and the first column of taxaTable needs to have to same name. See the example files. (we set this already at the start)
• taxaTablePath - The path to the taxonomic assignments for the taxomomic unit name used in your taxonomic assignment software (e.g., OTU, ASV, OGU). This has to be set to string(missing) if your readsTable already has taxonomic assignment included. Otherwise it need to consists of a column with header similar to the first header in readsTable and the column header "Taxon".
• outputPathMars - The path where the MARS results should be stored. We created this directory in the section 1 and is stored in paths.mars.

• sampleReadCountsCutoff - Numeric value for total read counts per sample under which samples are excluded from analysis. Only applies when readsTable contains absolute read counts (not relative abundances). Defaults to 1, with minimum of 1.
• cutoffMars - The cutoff underwich all relative abundances will be considered 0. If smaller relative abundances are kept, numerical difficulties/infeasibilites could occur later on when solving the microbiome models. Also, relative abundances below this value will most likely not have a noticable impact on the flux results as their contribution to fluxes is minimal due to their low relative abundance. Default is 1e-6.
• flagLoneSpecies - A boolean, flagLoneSpecies, which is true if the species name does NOT have the genus name already there. False if otherwise. Defaults to false.
• taxaDelimiter - A string, taxaSplit, which indicates the delimiter used to separate taxonomic levels in the taxonomic assignment.Defaults to ';'.
• removeClade - A boolean, removeClade, which removes all reads associated with a taxonomics identification that does not have information up to the specified taxonomic level that is looked at. Defaults to true.
• reconstructionDb: A string defining if AGORA2, APOLLO, a combination of both or a user-defined database should be used as model database to check presence in. Allowed Input (case-insensitive): "AGORA2", "APOLLO", "full_db", "user_db". Defaults to "full_db".
• userDbPath: A string containing the full path to the user-defined database, which should be in .csv, .txt, .parquet or .xlsx format and have column names = taxonomic levels. Only required if 'whichModelDatabase' is set to "user_db". Note, that the user database needs to have the same structure as the integrated AGORA2 & APOLLO database to function properly!
• calculateBrayCurtis: Boolean to indicate if Bray Curtis should be calculated for microbiome samples. Defaults to false.
• compoundedDatabase: Boolean; specifies is the reads table is compounded or not. Compounded here means that the reads for a specific taxonomic level are taking into account the taxonomic level above it.

% Run MARS to perform metagenomic mapping

runMars(readsTablePath, ...

'cutoffMars', cutoffMars, ...

'flagLoneSpecies', flagLoneSpecies, ...

'taxaDelimiter', taxaDelimiter, ...

'removeClade', removeClade, ...

'reconstructionDb', reconstructionDb, ...

'userDbPath', userDbPath, ...

'sampleReadCountCutoff', sampleReadCountsCutoff, ...

'taxaTablePath', taxaTablePath, ...

'outputPathMars', outputPathMars, ...

'calculateBrayCurtis', calculateBrayCurtis, ...

'compoundedDatabase', compoundedDatabase...

);

> Starting MARS > Loading database > Renaming and cleaning up of taxonomic assignments > Mapping and collecting metrics

> Generate metrics visualizations. > MARS finished sucessfully

• metrics - For each taxonomic level, various metric such as alpha and beta diversity are calculated.
• mapped - For each taxonomic level, only the mapped taxa, abundances are not renormalised and still have the values from the pre-mapping normalisation. Columns do not add up to 1
• processed - For each taxonomic level, all the taxa are normalised to the total reads. Columns add up to 1.
• unmapped - For each taxonomic level, only the unmapped taxa, abundances are not renormalised and still have the values from the pre-mapping normalisation. Columns do not add up to 1.
• mapped_forModelling - For each taxonomic level, the mapped taxa are renormalised so that each column adds up to 1.

Step 4: Creating microbiome community models

• It transforms all microbial reconstructions present in the relative abundance file. The external compartment is changed from [e] to [u] and reactions names are adapted. This is done to prepare the models to exchange with the internal environment of the microbiome [u].
• All metabolites that each microbial reconstruction can exchange with its environment are stored
• All reactions present in each reconstruction are stored as well as their associated subsystems. These will be used to calculate reactionsPresence, reactionAbundance and subsystemAbundance and can be used to describe the metabolic capabilities of indiviual microbioem models.
• All microbe reconstructions that are present in a sample are added together in one big model. The reconstructions can exchange metabolites in the [u] compartment created in the beginning of the function. The microbiome model can exchange metabolites with the outside environemt by transporting metabolites from [u] to [e] and vice versa.
• A microbiome biomass reaction is formulated which looks like: a microbe1 + b micobre2 + c microbe3 + ... The letters a,b and c are the relative abundance for the respective microbes as found in the relative abundance file we calculated via MARS.
• A standard diet is added to the microbiome models and checked if the total growth rate of 1 microbiome exchange per day (read as 1 fecal exchange per day) can be reached.
• Statistics on total amount of reactions, metabolites and microbes for each microbiome are collected and plotted in violin plots.

• microbeReconstructionPath - The path to the directory with the pan-species models. These can be downloaded from:

• relAbunFilePath - The path to the file with the AGORA2 mapped relative abundances. This is the mapped_forModellingSpecies.csv file created by MARS in section 2. If you have just run Mars default directory stucture will be assumed. Otherwise you have to define it yourself to point towards your relative abundance file.
• computeProfiles - A boolean (true or false) that tells the pipeline if we want to calculate the exctretion and uptake profiles of the microbiome models. We will not do this during this tutorial. However if you are more interested in the microbiome output, you can put this variable to true. This is a time-consuming step and is sped-up by using ibm_cplex as solver. For more information on the extra output if this is set to true look at this tutorial in the COBRA toolbox https://github.com/opencobra/COBRA.tutorials/tree/0d5becf7171844e53d3d43437035917a2bd9cc61/analysis/microbiomeModelingToolbox

• solver: your chosen solver (see above)
• numWorkersCreation - Number of workers you want to use for parellelisation. More workers means faster computing times. If you are unsure about the value, follow instructions in the code below. That will help make a decision.
• outputPathMgPipe - The path where the mgpipe results and the microbiome community models will be stored. The directory is stored in the paths variable created in section 1

• Individual health status not declared. Analysis will ignore that
• Directory already exists
• The temporary variable 'loadDiet' will be cleared at the beginning of each iteration of the parfor-loop. If 'loadDiet' is used before it is set, a runtime error will occur
• Column headers from the file were modified to make them valid MATLAB identifiers before creating variable names for the table. The original column headers are saved in the VariableDescriptions property. Set 'VariableNamingRule' to 'preserve' to use the original column headers as table variable names.
• Reaction x not in model
• File 'simRes.mat' not found.
• File 'intRes.mat' not found.

• GrowthRates - It gives the maximum growth rate of each microbiome under a rich diet and a pre-defined diet.
• ModelStatistics - for each microbiome model it gives the total amount of reactions, metabolites and microbes. MicrobiomeModel_sizes.png is a visual representation of this information.
• ModelStatsSummary - gives you the mean, median, min and max of the amount reactions metabolites and microbes considering all microbiome models.
• ReactionAbundance - gives the relative abundance of a reaction for each metabolic reaction in each microbiome model. Reaction abundances are based on the relative abundance of the microbe and wether or not the microbe has that specific reaction.
• ReactionPresence - says if a reaction is present in at least 1 microbe in a microbiome model.
• SubsystemAbundance - the same as reactionAbundance, however now all reactions are assigned a subsystem and abundances for the same subsystems are summed.

References (APA style)