A reproducible workflow for Füssel et al. 2025

If you have any questions, notice an issue, and/or are unable to find an important piece of information here, please feel free to leave a comment down below, send an e-mail to us, or get in touch with us through Discord:

Questions? Concerns? Find us on

Study description

Background

Microorganisms in the surface ocean remineralize the majority (~84%) of photosynthetically fixed carbon within minutes to days, with ~16% persisting for weeks to years, and <1% entering the long-lived reservoir of dissolved organic matter (DOM) that is comparable in size to the atmospheric carbon dioxide pool. The activity of microbial communities shapes the molecular composition of the marine DOM pool and drives the sequential transformtion of labile to persistent DOM. It is challenging to gain mechanistic insights into these microbially mediated processes in complex natural environments. Methodological and technological limitations such as the incomplete functional annotation of genes and the selective and incomplete recovery of dissolved organic compounds from seawater complicate the effective, direct integration of biological and chemical data. To explore the nature of microbial interactions driving DOM transformation, we used four bacterial isolates from the same North Sea water sample in a fully factorial setup, allowing us to study and compare the fates of several thousand largely unknown molecular formulas across eleven microbial co-cultures.

Cultures

We cultured the four isolate strains belonging to the Roseobacter clade individually and in co-cultures of two, three, and all four strains. We set up each treatment in triplicates. The artificial seawater minimal medium used in the experiment contained 1 g/L of glucose, trace elements, vitamins, and a bicarbonate buffer.

The four strains were found to have divergent metabolic capabilities from their genomes and substrate utilization preferences in culture. Pelagimonas varians SH4-1 (SH4) has a more extensive set of sugar metabolism genes than the other three strains, and grew on a variety of organic acids and monosaccharides, as well as a few polysaccharides. In the glucose minimal medium, SH4 had a negligible lag phase and grew to a higher optical density compared to the other strains. Phaeobacter sp. SH40 and Sulfitobacter sp. SH22-1 (SH22) grew well on organic acids and relatively poorly on sugars, and Sulfitobacter sp. SH24-1b (SH24) exhibited limited growth on all tested substrates. After a longer lag phase compared to SH4, the three other strains also grew more slowly to stationary phase.

Growth in co-culture contrasted with growth in monoculture. The observed growth exceeded modeled growth based on competitive glucose consumption, especially in co-cultures with SH4. The discrepancy between co-culture growth curves and models generally increased with the addition of strains, suggesting metabolic cross-feeding across strains.

Untargeted metabolomics

We extracted DOM from the filtered (0.2 µm) culture supernatant of each biological replicate at the beginning of the experiment and after 255 hours. DOM was extracted via Priority PolLutant (PPL) SPE cartridges, which preferentially retain hydrophobic organic compounds. Analytes were then measured by FT-ICR-MS in negative ion mode using electrospray ionization. The mass error was <0.1 ppm for all samples following calibration to endogenous peaks. Only masses detected in all replicates of a culture and not present in blanks were retained. Molecular formulas were assigned to spectra by ICBM-OCEAN software.

Molecular formula fates were tracked from pure culture to co-cultures. The fastest growing strain, SH4, yielded 2,216 formulas, or 89% of unique formulas found in pure cultures of the four strains. Across all co-cultures, 2,066 formulas were also found in pure cultures, while 2,508 were not. A majority of formulas novel to the co-cultures were unique to a single co-culture.

Compound prediction

A formula can represent various isomers, so we used the metabolic networks predicted for each of the strains and groups of strains in co-cultures to propose molecular identifications of the formulas, as described in this workflow. This approach involves the anvi’o reaction-network. Reaction networks are constructed from KEGG Ortholog (KO) annotations of genes (see anvi-run-kegg-kofams) and associated reaction and compound entries from the ModelSEED Biochemistry Database (see anvi-setup-modelseed-database). KOs are often annotated with KEGG reactions and EC numbers, indicating potential reactions that may be catalyzed by a gene protein product. Genomic reaction networks of co-cultured strains were merged to produce networks representing the combined metabolic potential of the community.

For each molecular feature in a culture, we matched its neutral formula, formula with one subtracted proton and charge of -1, and formula with two subtracted protons and charge of -2 to the formulas of compounds in the culture reaction network. The network often contains compounds in the protonation state that would exist in aqueous solution, so it is necessary to also search for -1 and -2 variants of the neutral formula to capture metabolites such as mono- and dicarboxylates.

Criteria

Formula matches to reaction network compounds were screened using a set of filters. Some of these criteria are implemented automatically while others require interpretation.

Multiple compound matches

A formula can match multiple compounds in a reaction network, and the strength of the evidence supporting each match can vary. We chose to retain formulas that match multiple closely related metabolites, such as isomers occurring in the same KEGG Pathway. Otherwise, we ignored formulas that match compounds with different metabolic roles. The search for deprotonated versions of each formula increases the likelihood of discarded uncertain matches to multiple compounds.

Compound consistency across cultures

If a formula is found in multiple cultures, a compound match must occur in all of the cultures’ reaction networks. If the formula is from cultures A and B, but the matching compound is only in the culture A network and not the culture B network, then the compound match would be ignored.

KO annotation specificity

Matching compounds must be strongly associated with KO annotations. We evaluated whether there was sufficient evidence to include compounds in the reaction network given the specificity of KO and associated reaction annotations involving the compounds. We ignored compounds included in the network via KOs associated with higher EC categories, such as 1.1.1.- and 2.3.-.-, or broad EC categories, such as 1.1.1.1 (alcohol dehydrogenase), that are linked to numerous ModelSEED reactions. Likewise, we ignored compounds included in the network via KOs with unconstrained catalytic capabilities, such as K00128 (aldehyde dehydrogenase), which is associated with a variety of reactions not necessarily catalyzed by the particular gene product. Many KOs, however, are associated with a single reaction, reducing the uncertainty that participating compounds are actually involved in the organism’s metabolism.

Production pathway

Matching compounds must be produced by reactions in a network, not just consumed. Furthermore, reactions are more likely to occur in the organism when they are well-connected to other reactions encoded by the network rather than isolated from other parts of the network, particularly where reaction substrates and products do not arise from and feed into other reactions in the network. We checked KEGG pathway maps for reaction connectivity. Gene KO annotations are also occasionally wrong, with a lower-ranking KO hit to the gene sequence rather than the top hit representing the true protein product, and co-occurrence of a KO with others in a KEGG pathway bolsters confidence in the KO. Additionally, annotation of multiple genes with the same KO bolsters annotation confidence.

Compound chemistry

Chemical considerations support the existence of a matching compound. Predicted compounds are more likely to exist in the sample when they have properties consistent with sample extraction and ionization. The SPE cartridges used in our study are more likely to retain hydrophobic compounds, and the negative ion mode in which the mass spectrometer was run favors ionization of compounds that can attain a -1 charge, such as carboxylic and phenolic acids.

Known biological isomers

There is the possibility that the true compound represented by a formula is not encoded in the reaction network. It is therefore sensible to compare the number of compounds with the formula in the network to the number in a large database of metabolites. We find the number of isomeric compounds in the ModelSEED Biochemistry compound database. This database includes pesticides and other synthetic compounds, many of which are not represented in the KEGG compound database, one of the databases incorporated into the ModelSEED database. Thus we also subset the isomeric ModelSEED compounds to those in the KEGG database. Furthermore, we subset isomeric KEGG compounds that participate in KEGG reactions, as these tend to be more common biological substrates. All else equal, matching compounds with fewer “potential false negative” isomers in the reference databases are more likely to actually be in the culture than compounds with more isomers.

Reproducing this workflow

Computational environment

This workflow uses the development version of anvi’o (8-dev), which you can install and activate following anvi’o installation instructions. Any more recent version of anvi’o should also work successfully. Load the anvi’o conda environment before running the workflow. The ModelSEED database should be installed in the default location for the anvi’o environment by anvi-setup-modelseed-database.

The computational demands of reproducing the workflow are minimal. All commands below should run within a few minutes or less on a modest laptop.

The data pack

Below you will find brief descriptions of individual files used in our downstream analyses. If you would like to follow this workflow, you can download the following data pack that includes the four Roseobacter genomes and the metabolomics table associated with each culture experiment. For this, please open a terminal, create a work directory, and type the following commands (or replace directory names manually):

# make sure there is a Downloads directory at your home
mkdir -p ~/Downloads

# change your current directory
cd ~/Downloads

# download the data pack
curl -o roseobacter-metabolomics.tar.gz https://merenlab.org/data/roseobacter-metabolomics/files/roseobacter-metabolomics.tar.gz

# unpack the data pack
tar -zxvf roseobacter-metabolomics.tar.gz

# go into the resulting data directory:
cd roseobacter-metabolomics

If you are here, you should be looking at a directory structure like this:

.
├── SH4-CONTIGS.db
├── SH40-CONTIGS.db
├── SH24-CONTIGS.db
├── SH22-CONTIGS.db
├── roseobacter-metabolomics-data.tsv

Genomes

The files with the extension .db represent the four isolate genomes sequenced with PacBio Hifi long reads. To include them in our computational workflows we used the anvi’o program anvi-gen-contigs-database to turn the FASTA files into so-called contigs-db files for downstream analyses. This file format contains much more information than a FASTA file, including gene coordinates, function annotations, and metabolic module membership of individual genes that will be essential to have in this workflow.

You can use the program anvi-db-info to learn more about the contents of a given contigs-db:

anvi-db-info SH22-CONTIGS.db

DB Info (no touch)
===============================================
Database Path ................................: SH22-CONTIGS.db
description ..................................: [Not found, but it's OK]
db_type ......................................: contigs (variant: unknown)
version ......................................: 24


DB Info (no touch also)
===============================================
project_name .................................: S_marinus_SH22
contigs_db_hash ..............................: hash52f2e51b
split_length .................................: 20000
kmer_size ....................................: 4
num_contigs ..................................: 4
total_length .................................: 4087537
num_splits ...................................: 201
gene_level_taxonomy_source ...................: None
genes_are_called .............................: 1
external_gene_calls ..........................: 0
external_gene_amino_acid_seqs ................: 0
skip_predict_frame ...........................: 0
splits_consider_gene_calls ...................: 1
trna_taxonomy_was_run ........................: 0
trna_taxonomy_database_version ...............: None
creation_date ................................: 1717747570.84562
modules_db_hash ..............................: a2b5bde358bb
scg_taxonomy_was_run .........................: 1
scg_taxonomy_database_version ................: GTDB: v214.1; Anvi'o: v1
gene_function_sources ........................: COG20_FUNCTION,Transfer_RNAs,CAZyme,KOfam,KEGG_Module,COG20_CATEGORY,KEGG_BRITE,KEGG_Class,COG20_PATHWAY
reaction_network_ko_annotations_hash .........: 1e5748cd73acfd2c24692de4d2c488044059aa32
reaction_network_kegg_database_release .......: 5a9644d40061
reaction_network_modelseed_database_sha ......: 194ac8afe48f8a606c0dd07ba3c7af10c02ba2fd

* Please remember that it is never a good idea to change these values. But in some
  cases it may be absolutely necessary to update something here, and a
  programmer may ask you to run this program and do it. But even then, you
  should be extremely careful.


AVAILABLE GENE CALLERS
===============================================
* 'prodigal' (3,851 gene calls)
* 'Transfer_RNAs' (45 gene calls)
* 'Ribosomal_RNA_23S' (3 gene calls)
* 'Ribosomal_RNA_16S' (3 gene calls)


AVAILABLE FUNCTIONAL ANNOTATION SOURCES
===============================================
* CAZyme (129 annotations)
* COG20_CATEGORY (3,208 annotations)
* COG20_FUNCTION (3,208 annotations)
* COG20_PATHWAY (830 annotations)
* KEGG_BRITE (2,396 annotations)
* KEGG_Class (511 annotations)
* KEGG_Module (511 annotations)
* KOfam (2,400 annotations)
* Transfer_RNAs (45 annotations)


AVAILABLE HMM SOURCES
===============================================
* 'Archaea_76' (76 models with 35 hits)
* 'Bacteria_71' (71 models with 72 hits)
* 'Protista_83' (83 models with 3 hits)
* 'Ribosomal_RNA_12S' (1 model with 0 hits)
* 'Ribosomal_RNA_16S' (3 models with 3 hits)
* 'Ribosomal_RNA_18S' (1 model with 0 hits)
* 'Ribosomal_RNA_23S' (2 models with 3 hits)
* 'Ribosomal_RNA_28S' (1 model with 0 hits)
* 'Ribosomal_RNA_5S' (5 models with 0 hits)
* 'Transfer_RNAs' (61 models with 45 hits)

You can get a standard FASTA file for a given genome using the program anvi-export-contigs:

anvi-export-contigs -c SH22-CONTIGS.db -o SH22.fa

Metabolomics table

The other file in this data pack, roseobacter-metabolomics-data.tsv, contains the processed spectral data, including monoisotopic molecular formulas and sample abundances. This is the same file that appears in our Füssel et al. publication as SI Table 2b.

Here are the first few lines of this table, so you can browse the individual columns that are included:

Creating the reaction networks

The workflow was originally implemented in a Jupyter notebook that used anvi’o libraries to process the anvi’o reaction-network.

Python cells of the Jupyter notebook are split up into sections of this workflow document with accompanying explanations and analyses of the output. If you wish to follow along interactively in Python, you can enter python3 in your terminal and run the code blocks sequentially. The following package imports are first required in the Python interactive shell.

import os
import sys
import rdkit
import itertools
import numpy as np
import pandas as pd

from rdkit import Chem
from copy import deepcopy
from typing import Iterable
from collections import defaultdict

Imports from anvi’o will fail if the anvi’o package isn’t in the Python module search path, a problem that can arise in Jupyter notebooks. If you have used the standard installation instructions on the anvi’o installation page, running this command should solve that issue:

sys.path.append('~/github/anvio')

Then you should be able to run these two lines without any errors:

import anvio
import anvio.reactionnetwork as rn

Genomic networks

The four contigs-db files in our data pack for the four strains we have worked with contain gene calls with KO annotations and reaction networks based on the KOs. Genes were annotated with KOs using anvi-run-kegg-kofams, and networks were constructed with anvi-reaction-network.

List the strains and their database files. Load reaction networks into memory. The dictionary of reaction networks is keyed by tuples, as co-culture “metagenomic” reaction networks keyed by tuples of strain IDs will be added to the dictionary.

all_strains = ['SH22', 'SH24', 'SH4', 'SH40']
strain_names = {
    'SH22': 'Sulfitobacter sp. SH22-1',
    'SH24': 'Sulfitobacter sp. SH24-1b',
    'SH4': 'Pelagimonas varians SH4-1',
    'SH40': 'Phaeobacter sp. SH40'
}
all_contigs_dbs = [f'{strain}-CONTIGS.db' for strain in all_strains]

con = rn.Constructor()
all_networks: dict[tuple[str], rn.GenomicNetwork] = {}
for contigs_db in all_contigs_dbs:
    strain = contigs_db[: contigs_db.index('-CONTIGS.db')]
    all_networks[(strain, )] = con.load_contigs_database_network(contigs_db, quiet=True)

Remove EC categories from networks

Avoid the inclusion of reactions on the basis of higher EC categories, such as 1.1.1.- or 2.3.-.-, that annotate KOs. Higher categories encompass a range of ModelSEED reactions that cannot be confidently attributed to the particular enzyme. Inclusion of these reactions increases the likelihood of false positive formula matches to compounds that are not actually produced by the organism. Networks filtered to removed EC categories are called “refined” networks. These network are used in formula matching.

all_refined_networks: dict[tuple[str], rn.GenomicNetwork] = {}
for strain_combo, unrefined_network in all_networks.items():
    modelseed_reaction_ids_to_retain = []
    for ko in unrefined_network.kos.values():
        modelseed_reaction_ids_to_check = []
        for modelseed_reaction_id, ec_numbers in ko.ec_number_aliases.items():
            for ec_number in ec_numbers:
                if '-' not in ec_number:
                    modelseed_reaction_ids_to_retain.append(modelseed_reaction_id)
                    break
            else:
                modelseed_reaction_ids_to_check.append(modelseed_reaction_id)
        for modelseed_reaction_id in modelseed_reaction_ids_to_check:
            if modelseed_reaction_id in ko.kegg_reaction_aliases:
                modelseed_reaction_ids_to_retain.append(modelseed_reaction_id)
    modelseed_reaction_ids_to_retain = set(modelseed_reaction_ids_to_retain)
    refined_network = unrefined_network.subset_network(reactions_to_subset=modelseed_reaction_ids_to_retain)
    all_refined_networks[strain_combo] = refined_network

Networks based on KEGG reactions

Compare the sizes of reaction networks constructed in two ways: first, using the default of both KEGG reaction and EC number annotations of KOs, and second, using just KEGG reaction annotations. KEGG reactions are more specific than EC numbers, which often map to a larger group of reactions in the ModelSEED database, as explained above in KO annotation specificity. The “EC+KEGG” network is prone to more false positive formula-compound matches that must be evaluated and fewer false negatives, or missing formula-compound matches, than the “just KEGG” network.

Although it would be useful to design a flag in anvi-reaction-network that allows a network to be constructed from KEGG reactions excluding EC numbers, for now we will remove the parts of the “EC+KEGG” networks that are based solely on EC numbers. This is achieved using the function from the anvi’o library that subsets networks by select items.

all_kegg_networks: dict[tuple[str], rn.GenomicNetwork] = {}
for strain_combo, ec_kegg_network in all_networks.items():
    modelseed_reaction_ids_to_retain = []
    for ko in ec_kegg_network.kos.values():
        for modelseed_reaction_id in ko.kegg_reaction_aliases:
            modelseed_reaction_ids_to_retain.append(modelseed_reaction_id)
    modelseed_reaction_ids_to_retain = set(modelseed_reaction_ids_to_retain)
    kegg_network = ec_kegg_network.subset_network(reactions_to_subset=modelseed_reaction_ids_to_retain)
    all_kegg_networks[strain_combo] = kegg_network

Co-culture “metagenomic” networks

Merge genomic reaction networks to represent co-culture “metagenomic” reaction networks. The network merge function avoids duplicate entries, such as KOs or reactions shared by both networks. Genes with identical anvi’o gene caller IDs (GCIDs) in different genomes would be considered the same in merging, so the identity of the genes must be maintained by adjusting integer GCIDs to be non-overlapping. Since the number of genes in these genomes is less than 10,000, add 10,000 to SH22 genome GCIDs, 20,000 to SH24 GCIDs, 30,000 to SH4 GCIDs, and 40,000 to SH40 GCIDs. Each gene in the network can thereby be traced back to the source genome, with SH22 genes, for example, having GCIDs between 10,000 and 20,000.

def make_gcids_nonoverlapping(networks: dict[tuple[str], rn.GenomicNetwork], increment: int = 10000) -> None:
    i = increment
    for network in networks.values():
        gcids_to_remove = []
        for gcid, gene in network.genes.items():
            assert gcid < increment
            new_gcid = i + gcid
            gene.gcid = new_gcid
            gcids_to_remove.append(gcid)
        for gcid in gcids_to_remove:
            gene = network.genes.pop(gcid)
            network.genes[gene.gcid] = gene
        i += increment

def merge_networks(networks: dict[tuple[str], rn.GenomicNetwork]) -> None:
    merged_networks = {}
    for r in range(2, len(networks) + 1):
        for combo in itertools.combinations(networks.items(), r):
            merged_strains = tuple()
            merged_network = None
            for strains, network in combo:
                merged_strains += strains
                if merged_network is None:
                    merged_network = network
                else:
                    merged_network = merged_network.merge_network(network)
            merged_networks[merged_strains] = merged_network
    networks.update(merged_networks)

make_gcids_nonoverlapping(all_networks)
merge_networks(all_networks)

make_gcids_nonoverlapping(all_refined_networks)
merge_networks(all_refined_networks)

make_gcids_nonoverlapping(all_kegg_networks)
merge_networks(all_kegg_networks)

List the strain combination tuples identifying the co-culture networks.

all_strain_combos = list(all_networks)

Compare networks constructed with different KO annotations

Compare the three types of networks constructed on the basis of varying KO annotations: KEGG reactions and all EC numbers (“default networks”), KEGG reactions and EC numbers but not higher EC categories (“refined networks”), and just KEGG reactions (“KEGG networks”). How many compounds are removed from the default networks excluding higher EC categories and EC numbers altogether?

header = ['strains', 'EC+KEGG_network_compounds', 'refined_network_compounds', 'KEGG_network_compounds']
rows = []
for strain_combo, ec_kegg_network in all_networks.items():
    refined_network = all_refined_networks[strain_combo]
    kegg_network = all_kegg_networks[strain_combo]
    row = []
    row.append('_'.join(strain_combo))
    row.append(len(ec_kegg_network.metabolites))
    row.append(len(refined_network.metabolites))
    row.append(len(kegg_network.metabolites))
    rows.append(row)
network_compound_counts = pd.DataFrame(rows, columns=header).set_index('strains')
network_compound_counts['refined_compound_fraction'] = network_compound_counts['refined_network_compounds'] / network_compound_counts['EC+KEGG_network_compounds']
network_compound_counts['KEGG_compound_fraction'] = network_compound_counts['KEGG_network_compounds'] / network_compound_counts['EC+KEGG_network_compounds']
print(network_compound_counts.to_string())

mean_refined_compound_fraction = network_compound_counts['refined_compound_fraction'].mean()
mean_kegg_compound_fraction = network_compound_counts['KEGG_compound_fraction'].mean()
print(f"An average of {round((1 - mean_refined_compound_fraction) * 100, 1)}% of compounds in the \"EC+KEGG\" network are removed ignoring higher EC categories in the \"refined\" network")
print(f"{round((1 - mean_kegg_compound_fraction) * 100, 1)}% of compounds in the \"EC+KEGG\" network are removed ignoring EC numbers and only considering KEGG reactions in the \"KEGG\" network")

On average 40.8% of compounds in the default “EC+KEGG” network are removed ignoring higher EC categories in the “refined” network. On average 73.2% of compounds in the default “EC+KEGG” network are removed ignoring EC numbers and only considering KEGG reactions in the “KEGG” network.

Prepare metabolomics data

Load the metabolomics data table, SI Table 2b from the paper. Each row represents a monoisotopic molecular feature.

roseobacter_metabolomics_df = pd.read_csv('roseobacter-metabolomics-data.tsv', sep='\t', header=0)

Confirm that a unique molecular formula was assigned to each feature.

len(roseobacter_metabolomics_df) == roseobacter_metabolomics_df['formula_isotopefree'].nunique()

Add deprotonated formulas

Add formulas for deprotonated versions of compounds as they may exist in the aqueous solution of cultures and the ModelSEED database used to populate compounds in reaction networks. Allow up to 2 hydrogens, 1 per oxygen, to be removed from each neutral formula. It does not make sense to remove 3 hydrogens in searching for common metabolites, since there are few with a -3 charge – primarily the tricarboxylic acids citrate, isocitrate, and aconitate in the TCA cycle.

formula_data = roseobacter_metabolomics_df[['formula', 'formula_isotopefree', 'O', 'H']]

deprot_rows = []
for _, row in formula_data.iterrows():
    formula_isotopefree = row.formula_isotopefree

    atom_count = {}
    for atomic_entry in row.formula.split():
        atom, count = atomic_entry.split('_')
        count = int(count)
        atom_count[atom] = count

    deprot_row = []
    for num_protons_subtracted in range(1, 3):
        if num_protons_subtracted > row.O:
            deprot_row.append('')
            continue

        new_atom_count = atom_count.copy()
        new_atom_count['H'] = atom_count['H'] - num_protons_subtracted

        new_formula_isotopefree = ''
        for atom, count in new_atom_count.items():
            new_formula_isotopefree += f'{atom}{count}' if count > 1 else atom
        deprot_row.append(new_formula_isotopefree)
    deprot_rows.append(deprot_row)

header = [f'formula_isotopefree_minus_{num_protons_subtracted}_H' for num_protons_subtracted in range(1, 3)]
deprot_table = pd.DataFrame(deprot_rows, columns=header)

cols = roseobacter_metabolomics_df.columns.tolist()
col_idx = cols.index('formula_isotopefree')
before = roseobacter_metabolomics_df[cols[: col_idx + 1]]
after = roseobacter_metabolomics_df[cols[col_idx + 1: ]]
feature_table = pd.concat([before, deprot_table, after], axis=1)

Make a new version of the table with a row per formula protonation state.

new_rows = []
new_idx = 0
for _, row in feature_table.iterrows():
    new_row = row.drop(['formula_isotopefree_minus_1_H', 'formula_isotopefree_minus_2_H'])
    new_row.name = new_idx
    new_row['search_formula'] = row['formula_isotopefree']
    new_row['search_charge'] = 0
    new_rows.append(new_row)
    new_idx += 1
    break

    if not row['formula_isotopefree_minus_1_H']:
        continue
    new_row = row.drop(['formula_isotopefree_minus_1_H', 'formula_isotopefree_minus_2_H'])
    new_row.name = new_idx
    new_row['search_formula'] = row['formula_isotopefree_minus_1_H']
    new_row['search_charge'] = -1
    new_rows.append(new_row)
    new_idx += 1

    if not row['formula_isotopefree_minus_2_H']:
        continue
    new_row = row.drop(['formula_isotopefree_minus_1_H', 'formula_isotopefree_minus_2_H'])
    new_row.name = new_idx
    new_row['search_formula'] = row['formula_isotopefree_minus_2_H']
    new_row['search_charge'] = -2
    new_rows.append(new_row)
    new_idx += 1

feature_table = pd.DataFrame(new_rows)
last_col_names = ['search_formula', 'search_charge']
first_col_names = feature_table.columns.tolist()[: -2]
feature_table = feature_table[last_col_names + first_col_names]

Match formulas to compounds

Find database isomers

Find compounds in the ModelSEED Biochemistry database with molecular formulas, including deprotonated formulas. To help evaluate the number of possible biomolecular isomers that could exist as part of controlling false positive compound matches (see the section, Known biological isomers), subset isomeric compounds in the KEGG compound database, and those that participate in KEGG reactions.

# Keys are (<formula>, <charge>), values are {<source of isomers>: [(<ModelSEED compound ID>, <ModelSEED compound name>)]}.
compound_isomers: dict[tuple[str, int], dict[str, list[tuple[str, str]]]] = {}
# Load the ModelSEED database from the default anvi'o installation location.
modelseed_db = rn.ModelSEEDDatabase()
compounds_table = modelseed_db.compounds_table

# Subset compounds with KEGG aliases.
compounds_with_kegg_alias_table = compounds_table[compounds_table['KEGG'].notna()]

# Subset compounds that participate in KEGG reactions.
kegg_reactions_table = modelseed_db.kegg_reactions_table
kegg_reaction_compound_ids = []
for compound_ids in kegg_reactions_table['compound_ids']:
    if not isinstance(compound_ids, str):
        continue
    compound_ids: str
    if compound_ids.strip() == '':
        continue
    for compound_id in compound_ids.split(';'):
        kegg_reaction_compound_ids.append(compound_id)
kegg_reaction_compound_ids = sorted(set(kegg_reaction_compound_ids))
select_rows = []
for row in compounds_with_kegg_alias_table.itertuples():
    if row.Index in kegg_reaction_compound_ids:
        select_rows.append(row)
compounds_with_kegg_reaction_table = pd.DataFrame(select_rows).set_index('Index')

for feature_row in feature_table.itertuples():
    formula = feature_row.search_formula
    charge = feature_row.search_charge
    compound_isomers[(formula, charge)] = isomers = {'modelseed_isomers': [], 'kegg_isomers': [], 'kegg_isomers_with_reaction': []}
    for compound_row in compounds_table[(compounds_table['formula'] == formula) & (compounds_table['charge'] == charge)].itertuples():
        isomers['modelseed_isomers'].append((compound_row.Index, compound_row.name))
    for compound_row in compounds_with_kegg_alias_table[(compounds_with_kegg_alias_table['formula'] == formula) & (compounds_with_kegg_alias_table['charge'] == charge)].itertuples():
        isomers['kegg_isomers'].append((compound_row.Index, compound_row.name))
    for compound_row in compounds_with_kegg_reaction_table[(compounds_with_kegg_reaction_table['formula'] == formula) & (compounds_with_kegg_reaction_table['charge'] == charge)].itertuples():
        isomers['kegg_isomers_with_reaction'].append((row.Index, row.name))

charge_isomer_stats: dict[int, dict[str, list[int]]] = {}
for charge in [0, -1, -2]:
    charge_isomer_stats[charge] = {'modelseed_isomers': [], 'kegg_isomers': [], 'kegg_isomers_with_reaction': []}
for (formula, charge), isomers in compound_isomers.items():
    isomer_stats = charge_isomer_stats[charge]
    for db_source, entries in isomers.items():
        if len(entries):
            isomer_stats[db_source].append(len(entries))

for charge, isomer_stats in charge_isomer_stats.items():
    formula_count = len(feature_table[feature_table['search_charge'] == charge])
    print(f"{formula_count} formulas with a charge of {charge} will be searched against reaction networks")
    for db_source, entry_counts in isomer_stats.items():
        print(f"- {len(entry_counts)} match {db_source.replace('_', ' ')}, {round(np.mean(entry_counts), 1)} isomers per formula on average")

A minority of molecular formulas match database compounds. A greater proportion of neutral formulas match database compounds than speculative deprotonated formulas with a -1 charge, and more -1 formulas match database compounds than -2 formulas.

• 4522 formulas with a charge of 0 will be searched against reaction networks
  • 667 match ModelSEED compounds, 3.4 isomers per formula on average
  • 554 match KEGG compounds, 2.8 isomers per formula on average
  • 280 match KEGG compounds in a reaction, 2.2 isomers per formula on average
• 4414 formulas with a charge of -1 will be searched against reaction networks
  • 243 match ModelSEED compounds, 2.7 isomers per formula on average
  • 189 match KEGG compounds, 2.2 isomers per formula on average
  • 130 match KEGG compounds in a reaction, 1.9 isomers per formula on average
• 4181 formulas with a charge of -2 will be searched against reaction networks
  • 73 match ModelSEED compounds, 2.1 isomers per formula on average
  • 50 match KEGG compounds, 1.8 isomers per formula on average
  • 40 match KEGG compounds in a reaction, 2.0 isomers per formula on average

Since reaction network compounds must be in ModelSEED, these statistics also show the upper bound of the number of formulas that may be identified in the genomes.

Search reaction network compounds

Match molecular formulas to compounds predicted in the reaction networks. If a feature is observed in a particular culture, match it to that culture’s network. Match to the “refined” network, which ignores higher EC category annotations of KOs, and match to the “KEGG” network, which ignores EC number annotations altogether.

def match_formulas(networks: dict[tuple[str], rn.GenomicNetwork]) -> tuple[
    dict[tuple[str, int], dict[tuple[str], list[rn.ModelSEEDCompound]]],
    dict[tuple[str, int], dict[tuple[str], rn.GenomicNetwork]],
    dict[tuple[str], dict[tuple[str, int], list[rn.ModelSEEDCompound]]],
    dict[tuple[str], dict[tuple[str, int], rn.GenomicNetwork]]
]:
    formula_culture_compounds: dict[tuple[str, int], dict[tuple[str], list[rn.ModelSEEDCompound]]] = {}
    formula_culture_subnetwork: dict[tuple[str, int], dict[tuple[str], rn.GenomicNetwork]] = {}
    culture_formula_compounds: dict[tuple[str], dict[tuple[str, int], list[rn.ModelSEEDCompound]]] = {}
    culture_formula_subnetwork: dict[tuple[str], dict[tuple[str, int], rn.GenomicNetwork]] = {}
    for row in feature_table.itertuples():
        formula = row.search_formula
        charge = row.search_charge
        formula_culture_compounds[(formula, charge)] = culture_compounds = {}
        formula_culture_subnetwork[(formula, charge)] = culture_subnetwork = {}
        for strain_combo, network in networks.items():
            final_abund = getattr(row, f"{'_'.join(strain_combo)}_Final")
            start_abund = getattr(row, f"{'_'.join(strain_combo)}_Start")
            if final_abund - start_abund == 0:
                continue

            matcher = rn.FormulaMatcher(network)
            compounds, subnetwork = matcher.match_metabolites_network(formula, charge=charge)
            culture_compounds[strain_combo] = compounds
            culture_subnetwork[strain_combo] = subnetwork

            try:
                formula_compounds = culture_formula_compounds[strain_combo]
            except KeyError:
                culture_formula_compounds[strain_combo] = formula_compounds = {}
            formula_compounds[(formula, charge)] = compounds

            try:
                formula_subnetwork = culture_formula_subnetwork[strain_combo]
            except KeyError:
                culture_formula_subnetwork[strain_combo] = formula_subnetwork = {}
            formula_subnetwork[(formula, charge)] = subnetwork
    return formula_culture_compounds, formula_culture_subnetwork, culture_formula_compounds, culture_formula_subnetwork

formula_culture_refined_compounds, formula_culture_refined_subnetwork, culture_formula_refined_compounds, culture_formula_refined_subnetwork = match_formulas(all_refined_networks)
formula_culture_kegg_compounds, formula_culture_kegg_subnetwork, culture_formula_kegg_compounds, culture_formula_kegg_subnetwork = match_formulas(all_kegg_networks)

formula_count = len(feature_table[feature_table['search_charge'] == charge])
charge_match_stats: dict[int, dict[str, int]] = {}
for charge in [0, -1, -2]:
    charge_match_stats[charge] = {'search_formulas': 0, 'refined': 0, 'kegg': 0}
for (formula, charge), culture_compounds in formula_culture_refined_compounds.items():
    match_stats = charge_match_stats[charge]
    match_stats['search_formulas'] += 1
    for compounds in culture_compounds.values():
        if compounds:
            match_stats['refined'] += 1
            break
for (formula, charge), culture_compounds in formula_culture_kegg_compounds.items():
    match_stats = charge_match_stats[charge]
    for compounds in culture_compounds.values():
        if compounds:
            match_stats['kegg'] += 1
            break

for charge, match_stats in charge_match_stats.items():
    print(f"{match_stats['search_formulas']} formulas with a charge of {charge} were searched against reaction networks")
    print(f"- {match_stats['refined']} match refined network compounds")
    print(f"- {match_stats['kegg']} match KEGG network compounds")

Here are the numbers of formulas of different charges that match reaction network compounds.

• 4522 formulas with a charge of 0 were searched against reaction networks
  • 71 match refined network compounds
  • 53 match KEGG network compounds
• 4414 formulas with a charge of -1 were searched against reaction networks
  • 29 match refined network compounds
  • 20 match KEGG network compounds
• 4181 formulas with a charge of -2 were searched against reaction networks
  • 15 match refined network compounds
  • 15 match KEGG network compounds

Evaluate compound matches

Evaluate the strength of compound matches using the criteria given earlier.

Detailed compound match report

We generated a report for each compound match to evaluate the criteria. The report is structured like the reaction network, showing the genes, then KO annotations, then reaction annotations that are the basis of the inclusion of the compound in the network.

# Filter the feature table to rows representing search formulas that match network compounds.
matching_formulas: list[tuple[str, int]] = []
for (search_formula, search_charge), culture_compounds in formula_culture_refined_compounds.items():
    for compounds in culture_compounds.values():
        if compounds:
            matching_formulas.append((search_formula, search_charge))
matching_feature_table = feature_table[feature_table[['search_formula', 'search_charge']].apply(tuple, axis=1).isin(matching_formulas)]

indent_increment = 4
for formula_isotopefree, group_table in matching_feature_table.groupby('formula_isotopefree'):
    # Print neutral formulas that have at least one search formula match network compounds.
    print(f"Feature neutral formula: {formula_isotopefree}")

    group_matching_formulas: list[tuple[str, int]] = []
    for group_row in group_table.itertuples():
        group_matching_formulas.append((group_row.search_formula, group_row.search_charge))

    for strain_combo in all_strain_combos:
        formula_compounds = culture_formula_refined_compounds[strain_combo]
        for group_matching_formula in group_matching_formulas:
            if group_matching_formula in formula_compounds:
                break
        else:
            continue
        # Print cultures with a network compound matching a search formula.
        print(f"{' ' * indent_increment}Culture: {'_'.join(strain_combo)}")

        for group_matching_formula in group_matching_formulas:
            try:
                matching_compounds = formula_compounds[group_matching_formula]
            except KeyError:
                continue
            # Print search formulas that match the culture network.
            print(f"{' ' * indent_increment * 2}Search formula: {group_matching_formula[0]} [{group_matching_formula[1]}]")

            isomers = compound_isomers[group_matching_formula]
            modelseed_isomer_count = len(isomers['modelseed_isomers'])
            kegg_isomer_count = len(isomers['kegg_isomers'])
            kegg_with_reaction_isomer_count = len(isomers['kegg_isomers_with_reaction'])
            # Print database isomer counts.
            print(f"{' ' * indent_increment * 2}- ModelSEED database isomer count: {modelseed_isomer_count}")
            print(f"{' ' * indent_increment * 2}- ModelSEED database KEGG compound isomer count: {kegg_isomer_count}")
            print(f"{' ' * indent_increment * 2}- ModelSEED database KEGG compound in KEGG reaction isomer count: {kegg_with_reaction_isomer_count}")

            formula_subnetwork = formula_culture_refined_subnetwork[group_matching_formula][strain_combo]
            for compound in matching_compounds:
                # Print compound matches.
                print(f"{' ' * indent_increment * 3}ModelSEED {compound.modelseed_id} {compound.modelseed_name}")
                print(f"{' ' * indent_increment * 3}- KEGG compound aliases: {' '.join(compound.kegg_aliases)}")

                compound_subnetwork = formula_subnetwork.subset_network(metabolites_to_subset=[compound.modelseed_id])
                for gcid, gene in compound_subnetwork.genes.items():
                    # Print genes linked to the compound.
                    print(f"{' ' * indent_increment * 4}Gene {gcid}")

                    for ko_id in gene.ko_ids:
                        # Print KO annotations of the gene. Print ModelSEED reactions associated
                        # (via EC numbers and KEGG reactions) with the KO. Print all EC numbers and
                        # KEGG reactions associated with the KO.
                        ko = compound_subnetwork.kos[ko_id]
                        print(f"{' ' * indent_increment * 5}KO {ko_id} {ko.name}")
                        print(f"{' ' * indent_increment * 5}- KO-associated ModelSEED reaction IDs: {' '.join(ko.reaction_ids)}")
                        message = ""
                        for modelseed_reaction_id, kegg_reaction_ids in ko.kegg_reaction_aliases.items():
                            message += f"{modelseed_reaction_id}: {' '.join([kegg_reaction_id for kegg_reaction_id in kegg_reaction_ids])} ; "
                        message = message[: -3]
                        print(f"{' ' * indent_increment * 5}- KO-associated ModelSEED reaction KEGG reaction aliases: {message}")
                        message = ""
                        for modelseed_reaction_id, ec_numbers in ko.ec_number_aliases.items():
                            message += f"{modelseed_reaction_id}: {' '.join([ec_number for ec_number in ec_numbers])} ; "
                        message = message[: -3]
                        print(f"{' ' * indent_increment * 5}- KO-associated ModelSEED reaction EC number aliases: {message}")

                        for reaction_id in ko.reaction_ids:
                            # Print ModelSEED reactions involving the compound. Print all KEGG and
                            # EC number aliases of the reaction.
                            reaction = compound_subnetwork.reactions[reaction_id]
                            equation = rn.get_chemical_equation(
                                reaction,
                                use_compound_names=[compound_subnetwork.metabolites[compound_id].modelseed_name for compound_id in reaction.compound_ids],
                                ignore_compartments=True
                            )
                            print(f"{' ' * indent_increment * 6}ModelSEED reaction {reaction_id}")
                            print(f"{' ' * indent_increment * 6}{equation}")
                            print(f"{' ' * indent_increment * 6}- KEGG reaction aliases: {' '.join(reaction.kegg_aliases)}")
                            print(f"{' ' * indent_increment * 6}- EC number aliases: {' '.join(reaction.ec_number_aliases)}")
                            try:
                                kegg_reaction_aliases = " ".join(ko.kegg_reaction_aliases[reaction_id])
                            except KeyError:
                                kegg_reaction_aliases = ""
                            print(f"{' ' * indent_increment * 6}- KO KEGG reaction associations: {kegg_reaction_aliases}")
                            try:
                                ec_number_aliases = " ".join(ko.ec_number_aliases[reaction_id])
                            except KeyError:
                                ec_number_aliases = ""
                            print(f"{' ' * indent_increment * 6}- KO EC number associations: {ec_number_aliases}")

Part of the report for the first match in the output is shown below. A feature was assigned the neutral molecular formula of C₁₀H₁₀O₆. The deprotonated variants of the formula with a -1 and -2 charge were also searched against the reaction networks of the culture with this feature. As shown, one network was from the pure culture of SH4, while another was from the co-culture of SH22 and SH4. The feature matched compounds in the SH4 network and all SH4 co-culture networks, suggesting that it was produced by SH4 but not fully consumed by SH22, SH24, and SH40 in co-culture. C₁₀H₈O₆^-2 was the only deprotonated formula that matched compounds in the networks. To evaluate the potential breadth of the match, we found all isomers with the formula in the ModelSEED Biochemistry compound database, and two subsets of the database in KEGG. There were three isomers with the formula in each of the three sets of reference compounds. These three compounds – prephenate, chorismate, and isochorismate – are also in the reaction networks (note that since they are in the SH4 network, they must be in the co-culture networks which are supersets of the SH4 network). The absence of other isomeric compounds in the ModelSEED database besides those in the reaction network reduces the likelihood of missing biological compounds that may actually represent the molecular feature. Isochorismate is ignored because it is only included on the basis of an enzyme which consumes it, isochorismate pyruvate lyase (K04782), not any enzymes that produce it.

Prephenate and chorismate are related compounds in the shikimate pathway for biosynthesis of aromatic amino acids and other compounds. The report presents genomic evidence for production of these compounds. The SH4 genome encodes chorismate mutase, the key enzyme responsible for prephenate biosynthesis from chorismate. The genome also encodes cyclohexadienyl dehydratase and prephenate dehydrogenase, enzymes which react prephenate to form the precursors of phenylalanine and tyrosine, respectively. Chorismate mutase (K04092) has KEGG reactions and an EC number linked to three ModelSEED reactions, redundant entries with different IDs for the same chorismate mutase reaction of chorismate to prephenate.

Feature neutral formula: C10H10O6
    Culture: SH4
        Search formula: C10H8O6 [-2]
        - ModelSEED database isomer count: 3
        - ModelSEED database KEGG compound isomer count: 3
        - ModelSEED database KEGG compound in KEGG reaction isomer count: 3
            ModelSEED cpd00219 Prephenate
            - KEGG compound aliases: C00254
                Gene 30271
                    KO K04092 chorismate mutase [EC:5.4.99.5]
                    - KO-associated ModelSEED reaction IDs: rxn01256 rxn19309 rxn33299
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn01256: R01715 ; rxn19309: R01715 ; rxn33299: R01715
                    - KO-associated ModelSEED reaction EC number aliases: rxn01256: 5.4.99.5 ; rxn19309: 5.4.99.5 ; rxn33299: 5.4.99.5
                        ModelSEED reaction rxn01256
                        1 Chorismate -> 1 Prephenate
                        - KEGG reaction aliases: R01715
                        - EC number aliases: 5.4.99.5
                        - KO KEGG reaction associations: R01715
                        - KO EC number associations: 5.4.99.5
                        ModelSEED reaction rxn19309
                        1 Chorismate -> 1 Prephenate
                        - KEGG reaction aliases: R01715
                        - EC number aliases: 5.4.99.5
                        - KO KEGG reaction associations: R01715
                        - KO EC number associations: 5.4.99.5
                        ModelSEED reaction rxn33299
                        1 Chorismate -> 1 Prephenate
                        - KEGG reaction aliases: R01715
                        - EC number aliases: 5.4.99.5
                        - KO KEGG reaction associations: R01715
                        - KO EC number associations: 5.4.99.5
                Gene 34329
                    KO K00220 cyclohexadieny/prephenate dehydrogenase [EC:1.3.1.43 1.3.1.12]
                    - KO-associated ModelSEED reaction IDs: rxn01268 rxn28086 rxn33078
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn01268: R01728 ; rxn28086: R01728 ; rxn33078: R01728
                    - KO-associated ModelSEED reaction EC number aliases: rxn01268: 1.3.1.12 ; rxn28086: 1.3.1.12 ; rxn33078: 1.3.1.12
                        ModelSEED reaction rxn01268
                        1 NAD + 1 Prephenate -> 1 NADH + 1 CO2 + 1 p-hydroxyphenylpyruvate
                        - KEGG reaction aliases: R01728
                        - EC number aliases: 1.3.1.12
                        - KO KEGG reaction associations: R01728
                        - KO EC number associations: 1.3.1.12
                        ModelSEED reaction rxn28086
                        1 NAD + 1 Prephenate -> 1 NADH + 1 CO2 + 1 p-hydroxyphenylpyruvate
                        - KEGG reaction aliases: R01728
                        - EC number aliases: 1.3.1.12
                        - KO KEGG reaction associations: R01728
                        - KO EC number associations: 1.3.1.12
                        ModelSEED reaction rxn33078
                        1 NAD + 1 Prephenate -> 1 NADH + 1 CO2 + 1 p-hydroxyphenylpyruvate
                        - KEGG reaction aliases: R01728
                        - EC number aliases: 1.3.1.12
                        - KO KEGG reaction associations: R01728
                        - KO EC number associations: 1.3.1.12
                    KO K04517 prephenate dehydrogenase [EC:1.3.1.12]
                    - KO-associated ModelSEED reaction IDs: rxn01268 rxn28086 rxn33078
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn01268: R01728 ; rxn28086: R01728 ; rxn33078: R01728
                    - KO-associated ModelSEED reaction EC number aliases: rxn01268: 1.3.1.12 ; rxn28086: 1.3.1.12 ; rxn33078: 1.3.1.12
                        ModelSEED reaction rxn01268
                        1 NAD + 1 Prephenate -> 1 NADH + 1 CO2 + 1 p-hydroxyphenylpyruvate
                        - KEGG reaction aliases: R01728
                        - EC number aliases: 1.3.1.12
                        - KO KEGG reaction associations: R01728
                        - KO EC number associations: 1.3.1.12
                        ModelSEED reaction rxn28086
                        1 NAD + 1 Prephenate -> 1 NADH + 1 CO2 + 1 p-hydroxyphenylpyruvate
                        - KEGG reaction aliases: R01728
                        - EC number aliases: 1.3.1.12
                        - KO KEGG reaction associations: R01728
                        - KO EC number associations: 1.3.1.12
                        ModelSEED reaction rxn33078
                        1 NAD + 1 Prephenate -> 1 NADH + 1 CO2 + 1 p-hydroxyphenylpyruvate
                        - KEGG reaction aliases: R01728
                        - EC number aliases: 1.3.1.12
                        - KO KEGG reaction associations: R01728
                        - KO EC number associations: 1.3.1.12
                Gene 32286
                    KO K01713 cyclohexadienyl dehydratase [EC:4.2.1.51 4.2.1.91]
                    - KO-associated ModelSEED reaction IDs: rxn01000 rxn28085 rxn33346 rxn33962
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn01000: R01373 ; rxn28085: R01373 ; rxn33346: R01373
                    - KO-associated ModelSEED reaction EC number aliases: rxn01000: 4.2.1.51 4.2.1.91 ; rxn28085: 4.2.1.51 4.2.1.91 ; rxn33346: 4.2.1.51 4.2.1.91 ; rxn33962: 4.2.1.51
                        ModelSEED reaction rxn01000
                        1 H+ + 1 Prephenate -> 1 H2O + 1 CO2 + 1 Phenylpyruvate
                        - KEGG reaction aliases: R01373
                        - EC number aliases: 4.2.1.51 4.2.1.91
                        - KO KEGG reaction associations: R01373
                        - KO EC number associations: 4.2.1.51 4.2.1.91
                        ModelSEED reaction rxn28085
                        1 H+ + 1 Prephenate -> 1 H2O + 1 CO2 + 1 Phenylpyruvate
                        - KEGG reaction aliases: R01373
                        - EC number aliases: 4.2.1.51 4.2.1.91
                        - KO KEGG reaction associations: R01373
                        - KO EC number associations: 4.2.1.51 4.2.1.91
                        ModelSEED reaction rxn33346
                        1 H+ + 1 Prephenate -> 1 H2O + 1 CO2 + 1 Phenylpyruvate
                        - KEGG reaction aliases: R01373
                        - EC number aliases: 4.2.1.51 4.2.1.91
                        - KO KEGG reaction associations: R01373
                        - KO EC number associations: 4.2.1.51 4.2.1.91
                        ModelSEED reaction rxn33962
                        1 Prephenate <-> 1 H2O + 1 CO2 + 1 Chloroplast Phenylpyruvate
                        - KEGG reaction aliases:
                        - EC number aliases: 4.2.1.51
                        - KO KEGG reaction associations:
                        - KO EC number associations: 4.2.1.51
                Gene 30738
                    KO K04518 prephenate dehydratase [EC:4.2.1.51]
                    - KO-associated ModelSEED reaction IDs: rxn01000 rxn28085 rxn33346 rxn33962
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn01000: R01373 ; rxn28085: R01373 ; rxn33346: R01373
                    - KO-associated ModelSEED reaction EC number aliases: rxn01000: 4.2.1.51 ; rxn28085: 4.2.1.51 ; rxn33346: 4.2.1.51 ; rxn33962: 4.2.1.51
                        ModelSEED reaction rxn01000
                        1 H+ + 1 Prephenate -> 1 H2O + 1 CO2 + 1 Phenylpyruvate
                        - KEGG reaction aliases: R01373
                        - EC number aliases: 4.2.1.51 4.2.1.91
                        - KO KEGG reaction associations: R01373
                        - KO EC number associations: 4.2.1.51
                        ModelSEED reaction rxn28085
                        1 H+ + 1 Prephenate -> 1 H2O + 1 CO2 + 1 Phenylpyruvate
                        - KEGG reaction aliases: R01373
                        - EC number aliases: 4.2.1.51 4.2.1.91
                        - KO KEGG reaction associations: R01373
                        - KO EC number associations: 4.2.1.51
                        ModelSEED reaction rxn33346
                        1 H+ + 1 Prephenate -> 1 H2O + 1 CO2 + 1 Phenylpyruvate
                        - KEGG reaction aliases: R01373
                        - EC number aliases: 4.2.1.51 4.2.1.91
                        - KO KEGG reaction associations: R01373
                        - KO EC number associations: 4.2.1.51
                        ModelSEED reaction rxn33962
                        1 Prephenate <-> 1 H2O + 1 CO2 + 1 Chloroplast Phenylpyruvate
                        - KEGG reaction aliases:
                        - EC number aliases: 4.2.1.51
                        - KO KEGG reaction associations:
                        - KO EC number associations: 4.2.1.51
            ModelSEED cpd00216 Chorismate
            - KEGG compound aliases: C00251
                Gene 32359
                    KO K01657 anthranilate synthase component I [EC:4.1.3.27]
                    - KO-associated ModelSEED reaction IDs: rxn00726 rxn00727 rxn27709 rxn32242 rxn33991 rxn35359 rxn38042 rxn38043
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn00726: R00985 ; rxn00727: R00986 ; rxn27709: R00986 ; rxn32242: R00985 ; rxn35359: R00985 ; rxn38042: R00986 ; rxn38043: R00986
                    - KO-associated ModelSEED reaction EC number aliases: rxn00726: 4.1.3.27 ; rxn00727: 4.1.3.27 ; rxn27709: 4.1.3.27 ; rxn32242: 4.1.3.27 ; rxn33991: 4.1.3.27 ; rxn35359: 4.1.3.27 ; rxn38042: 4.1.3.27 ; rxn38043: 4.1.3.27
                        ModelSEED reaction rxn00726
                        1 NH3 + 1 Chorismate -> 1 H2O + 1 Pyruvate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00985
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00985
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn00727
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn27709
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn32242
                        1 NH3 + 1 Chorismate -> 1 H2O + 1 Pyruvate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00985
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00985
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn33991
                        1 PRPP + 1 Chorismate + 1 Glutamine -> 1 PPi + 1 Pyruvate + 1 L-Glutamate + 2 H+ + 1 N-5-phosphoribosyl-anthranilate
                        - KEGG reaction aliases:
                        - EC number aliases: 2.4.2.18 4.1.3.27
                        - KO KEGG reaction associations:
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn35359
                        1 NH3 + 1 Chorismate -> 1 H2O + 1 Pyruvate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00985
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00985
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn38042
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn38043
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                Gene 30271
                    KO K04092 chorismate mutase [EC:5.4.99.5]
                    - KO-associated ModelSEED reaction IDs: rxn01256 rxn19309 rxn33299
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn01256: R01715 ; rxn19309: R01715 ; rxn33299: R01715
                    - KO-associated ModelSEED reaction EC number aliases: rxn01256: 5.4.99.5 ; rxn19309: 5.4.99.5 ; rxn33299: 5.4.99.5
                        ModelSEED reaction rxn01256
                        1 Chorismate -> 1 Prephenate
                        - KEGG reaction aliases: R01715
                        - EC number aliases: 5.4.99.5
                        - KO KEGG reaction associations: R01715
                        - KO EC number associations: 5.4.99.5
                        ModelSEED reaction rxn19309
                        1 Chorismate -> 1 Prephenate
                        - KEGG reaction aliases: R01715
                        - EC number aliases: 5.4.99.5
                        - KO KEGG reaction associations: R01715
                        - KO EC number associations: 5.4.99.5
                        ModelSEED reaction rxn33299
                        1 Chorismate -> 1 Prephenate
                        - KEGG reaction aliases: R01715
                        - EC number aliases: 5.4.99.5
                        - KO KEGG reaction associations: R01715
                        - KO EC number associations: 5.4.99.5
                Gene 31901
                    KO K00766 anthranilate phosphoribosyltransferase [EC:2.4.2.18]
                    - KO-associated ModelSEED reaction IDs: rxn33991
                    - KO-associated ModelSEED reaction KEGG reaction aliases:
                    - KO-associated ModelSEED reaction EC number aliases: rxn33991: 2.4.2.18
                        ModelSEED reaction rxn33991
                        1 PRPP + 1 Chorismate + 1 Glutamine -> 1 PPi + 1 Pyruvate + 1 L-Glutamate + 2 H+ + 1 N-5-phosphoribosyl-anthranilate
                        - KEGG reaction aliases:
                        - EC number aliases: 2.4.2.18 4.1.3.27
                        - KO KEGG reaction associations:
                        - KO EC number associations: 2.4.2.18
                Gene 31902
                    KO K01658 anthranilate synthase component II [EC:4.1.3.27]
                    - KO-associated ModelSEED reaction IDs: rxn00726 rxn00727 rxn27709 rxn32242 rxn33991 rxn35359 rxn38042 rxn38043
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn00726: R00985 ; rxn00727: R00986 ; rxn27709: R00986 ; rxn32242: R00985 ; rxn35359: R00985 ; rxn38042: R00986 ; rxn38043: R00986
                    - KO-associated ModelSEED reaction EC number aliases: rxn00726: 4.1.3.27 ; rxn00727: 4.1.3.27 ; rxn27709: 4.1.3.27 ; rxn32242: 4.1.3.27 ; rxn33991: 4.1.3.27 ; rxn35359: 4.1.3.27 ; rxn38042: 4.1.3.27 ; rxn38043: 4.1.3.27
                        ModelSEED reaction rxn00726
                        1 NH3 + 1 Chorismate -> 1 H2O + 1 Pyruvate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00985
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00985
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn00727
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn27709
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn32242
                        1 NH3 + 1 Chorismate -> 1 H2O + 1 Pyruvate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00985
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00985
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn33991
                        1 PRPP + 1 Chorismate + 1 Glutamine -> 1 PPi + 1 Pyruvate + 1 L-Glutamate + 2 H+ + 1 N-5-phosphoribosyl-anthranilate
                        - KEGG reaction aliases:
                        - EC number aliases: 2.4.2.18 4.1.3.27
                        - KO KEGG reaction associations:
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn35359
                        1 NH3 + 1 Chorismate -> 1 H2O + 1 Pyruvate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00985
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00985
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn38042
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn38043
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                Gene 33798
                    KO K01736 chorismate synthase [EC:4.2.3.5]
                    - KO-associated ModelSEED reaction IDs: rxn01255 rxn32460
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn01255: R01714 ; rxn32460: R01714
                    - KO-associated ModelSEED reaction EC number aliases: rxn01255: 4.2.3.5 ; rxn32460: 4.2.3.5
                        ModelSEED reaction rxn01255
                        1 5-O--1-Carboxyvinyl-3-phosphoshikimate -> 1 Phosphate + 1 Chorismate
                        - KEGG reaction aliases: R01714
                        - EC number aliases: 4.2.3.5
                        - KO KEGG reaction associations: R01714
                        - KO EC number associations: 4.2.3.5
                        ModelSEED reaction rxn32460
                        1 5-O--1-Carboxyvinyl-3-phosphoshikimate -> 1 Phosphate + 1 Chorismate
                        - KEGG reaction aliases: R01714
                        - EC number aliases: 4.2.3.5
                        - KO KEGG reaction associations: R01714
                        - KO EC number associations: 4.2.3.5
            ModelSEED cpd00658 Isochorismate
            - KEGG compound aliases: C00885
                Gene 33754
                    KO K04782 isochorismate pyruvate lyase [EC:4.2.99.21]
                    - KO-associated ModelSEED reaction IDs: rxn04454
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn04454: R06602
                    - KO-associated ModelSEED reaction EC number aliases: rxn04454: 4.2.99.21
                        ModelSEED reaction rxn04454
                        1 Isochorismate -> 1 Pyruvate + 1 SALC
                        - KEGG reaction aliases: R06602
                        - EC number aliases: 4.2.99.21
                        - KO KEGG reaction associations: R06602
                        - KO EC number associations: 4.2.99.21
    Culture: SH22_SH4
        Search formula: C10H8O6 [-2]
        - ModelSEED database isomer count: 3
        - ModelSEED database KEGG compound isomer count: 3
        - ModelSEED database KEGG compound in KEGG reaction isomer count: 3
            ModelSEED cpd00216 Chorismate
            - KEGG compound aliases: C00251
                Gene 11205
                    KO K01657 anthranilate synthase component I [EC:4.1.3.27]
                    - KO-associated ModelSEED reaction IDs: rxn00726 rxn00727 rxn27709 rxn32242 rxn33991 rxn35359 rxn38042 rxn38043
                    - KO-associated ModelSEED reaction KEGG reaction aliases: rxn00726: R00985 ; rxn00727: R00986 ; rxn27709: R00986 ; rxn32242: R00985 ; rxn35359: R00985 ; rxn38042: R00986 ; rxn38043: R00986
                    - KO-associated ModelSEED reaction EC number aliases: rxn00726: 4.1.3.27 ; rxn00727: 4.1.3.27 ; rxn27709: 4.1.3.27 ; rxn32242: 4.1.3.27 ; rxn33991: 4.1.3.27 ; rxn35359: 4.1.3.27 ; rxn38042: 4.1.3.27 ; rxn38043: 4.1.3.27
                        ModelSEED reaction rxn00726
                        1 NH3 + 1 Chorismate -> 1 H2O + 1 Pyruvate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00985
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00985
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn00727
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn27709
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn32242
                        1 NH3 + 1 Chorismate -> 1 H2O + 1 Pyruvate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00985
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00985
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn33991
                        1 PRPP + 1 Chorismate + 1 Glutamine -> 1 PPi + 1 Pyruvate + 1 L-Glutamate + 2 H+ + 1 N-5-phosphoribosyl-anthranilate
                        - KEGG reaction aliases:
                        - EC number aliases: 2.4.2.18 4.1.3.27
                        - KO KEGG reaction associations:
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn35359
                        1 NH3 + 1 Chorismate -> 1 H2O + 1 Pyruvate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00985
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00985
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn38042
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
                        ModelSEED reaction rxn38043
                        1 L-Glutamine + 1 Chorismate -> 1 Pyruvate + 1 L-Glutamate + 1 H+ + 1 Anthranilate
                        - KEGG reaction aliases: R00986
                        - EC number aliases: 4.1.3.27
                        - KO KEGG reaction associations: R00986
                        - KO EC number associations: 4.1.3.27
...

Pathway integration

It is important to assess the connectivity of putative compounds in the metabolic network of an organism to evaluate the likelihood that compounds are actually in the culture. Pathway-level analysis using maps was performed for each matching compound. The genomic capacity of the four strains to cycle prephenate and chorismate in the shikimate pathway is displayed in a KEGG pathway map produced by anvi-draw-kegg-pathways. The possibility of prephenate and chorismate being erroneous matches is reduced by the extensive genomic evidence for production of these compounds by SH4.

Extraction and ionization chemistry

Putative exometabolites are more likely when their structures are consistent with retention in the sample extract and the conditions of electrospray ionization. The PPL SPE cartridges used for DOM extraction in this study retain aromatic compounds well. The mass spectrometer was run in negative ion mode in this study, so compounds such as carboxylic and phenolic acids that can assume a -1 charge through a ready change in protonation state are easily ionized. Prephenate and chorismate are aromatic carboxylates that are well-suited for extraction and ionization, increasing the likelihood that they are indeed exometabolites.

# Here is the source on how to render molecules as a grid of SVG images in RDKit:
# https://iwatobipen.wordpress.com/2020/05/01/draw-molecules-as-svg-in-horizontal-layout-drawing-rdkit-memo/

class HorizontalDisplay:
    def __init__(self, *args):
        self.args = args

    def _repr_html_(self):
        template = '<div style="float:left;padding:10px;">{0}</div>'
        return "\n".join(template.format(arg) for arg in self.args)

def draw_structures(
    compounds: list[rn.ModelSEEDCompound] = None,
    panel_width: int = 250,
    panel_height: int = 250
) -> HorizontalDisplay:
    drawing_texts = []

    for compound in compounds:
        mol: Chem.Mol = Chem.MolFromSmiles(compound.smiles)
        mol.SetProp('_Name', compound.modelseed_name)
        d2d = rdMolDraw2D.MolDraw2DSVG(panel_width, panel_height)
        d2d.DrawMolecule(mol, legend=mol.GetProp('_Name'))
        d2d.FinishDrawing()
        text = d2d.GetDrawingText()
        drawing_texts.append(text)

    if not drawing_texts:
        return
    print(len(drawing_texts))

    return HorizontalDisplay(*drawing_texts)

chorismate = all_refined_networks[('SH4', )].metabolites['cpd00216']
prephenate = all_refined_networks[('SH4', )].metabolites['cpd00219']
draw_structures([prephenate, chorismate, isochorismate])

Evaluation table

The evaluation of formula matches to network compounds is summarized in the following table. Compound matches that were not found in all of the reaction networks of cultures containing the particular formula were first screened out and not included in the table (see Compound consistency across cultures). A molecular feature’s “Neutral formula” in the table has an “Ionized formula” that was searched against the reaction networks to yield a “Compound match”. “Database isomers” of the ionized formula are the numbers of isomeric ModelSEED Biochemistry database compounds, the subset aliasing KEGG compounds, and the subset aliasing KEGG compounds in reactions (see Known biological isomers).

Subsequent columns after “Compound match” evaluate criteria for compound assignment confidence. A value of 1 indicates that the compound passes the filter, and a value of 0 indicates removal by the filter. The last column, “Passes filters”, has a value of 1 if the compound match has a value of 1 in each of the individual filter columns. To expedite the process of evaluating compounds, the series of filters was applied from left to right, and evaluation was stopped if a value of 0 was recorded: this is the reason for absent values in the individual filter columns to the right of a 0.

Examples of compound matches that are retained and discarded illustrate the filters. As discussed above, C₁₀H₁₀O₆ matched three compounds in the culture reaction networks containing the formula: the deprotonated variant of the formula, C₁₀H₈O₆^-2, matched prephenate, chorismate, and isochorismate. These represented the only three isomeric compounds recorded in the ModelSEED Biochemistry database and the two KEGG subsets (the value of “3,3,3” in the “Database isomers” column). We judged that this relatively low number allowed the compounds to pass the “Database isomer specificity” filter, with a value of 1. The “Annotation specificity” filter was passed because the compounds were populated in the reaction networks via KO enzyme annotations with KEGG reactions and EC numbers that always specifically involve the compounds. Isochorismate did not pass the next “Metabolic integration” filter since the compound was only included in reaction networks via isochorismate pyruvate lyase (K04782), which consumes but does not produce the compound. In contrast, chorismate and prephenate are well-integrated metabolites produced by enzymes of the genomically complete shikimate pathway. Chorismate and prephenate then passed the “Ionizability” filter, which accounts for plausible chemical properties that make a compound suitable for measurement. Finally, the compounds passed the “Metabolic similarity” filter since they have similar metabolic roles as adjacent isomers in the shikimate pathway.

The next formula in the table, C₁₀H₁₁NO₂, matches 5-hydroxy-1H-indole-3-ethanol and indole-3-glycol. Both compounds passed the “Database isomer specificity” filter but were removed by the next “Annotation specificity” filter due to the broadness of the enzyme annotations that resulted in the compounds’ inclusion in the reaction networks. The compounds are involved in some of the numerous ModelSEED alcohol dehydrogenase reactions associated with EC 1.1.1.1, which annotates the enzymes S-(hydroxymethyl)glutathione dehydrogenase / alcohol dehydrogenase (K00121) and alcohol dehydrogenase, propanol-preferring (K13953). High uncertainty in the reaction specificity of the gene enzyme products reduces the likelihood that the particular compounds actually occur in the strains’ metabolomes.

Further down the table occurs C₁₂H₂₂O₁₁, which matches a number of disaccharides in the reaction networks. Before reaching the “Ionizability” filter, which would discard sugar compound matches due to the lack of an acidic proton for straightforward negative ionization, the compounds are filtered out by “Database isomer specificity.” A large number of possible compounds besides the matched disaccharides have the same formula in the reference databases (66 isomers in ModelSEED, 35 and 22 in the KEGG subsets), increasing the likelihood that other compounds produced by the strains actually represent the formula.

From this table, 53 molecular formulas had compound matches that passed the filters. These comprise SI Table 2f in our Füssel et al. publication.

The end

If you made it all the way here, you now know how we developed a computational workflow for “genomically guided metabolomics” and how we applied it to an intricate microbial co-culture experiment. In the end is the beginning, because this is only a first foray into a promising approach for plucking molecular needles out of the mass spectral haystack. Next we are exploring tandem mass spectrometry datasets to match structures instead of just molecular formulas and targeted metabolomes to further validate our approach. Best of luck, to you and us.