Dry lab · Software

Software

Any software tools we built to make our work easier or more reproducible.

In plain language

Every result on our Model page is reproducible from public data and the scripts below.

Reproducibility

Code base

Every result on this page is reproducible from public data and the scripts below. This section walks through the pipeline in the order we ran it, with the real code for each stage. The complete scripts, the intermediate tables behind each figure, and a README are provided as a downloadable bundle with this page. This is analysis code for reproducing our results, not a released software tool.

Get the code and data

Everything on this page, every script, every results table, every figure, and a README that maps each figure to its script and table, is in one download.

TexasiGEM_BestModel_Bundle.zip 6.8 MB, 66 files

Opens the file in our iGEM GitLab repository, where the Download button saves it. The repository is public, so no account is needed.

Inputs

Raw 16S reads: ENA accession PRJEB38494 (939 samples, 118 children). Reference database for metabolic modeling: AGORA2. Everything else is produced by the code here.

The analysis these commands produce, with the figures and the argument built on them, is on our Model page.

Stage 1

Stage 1: sequence processing (QIIME 2)

Import, denoise to ASVs with DADA2, and filter. This produces the feature table and representative sequences used everywhere downstream.

import, denoise and filterQIIME 2

qiime tools import --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path manifest.tsv --output-path demux.qza \
  --input-format PairedEndFastqManifestPhred33V2

qiime dada2 denoise-paired --i-demultiplexed-seqs demux.qza \
  --p-trunc-len-f 220 --p-trunc-len-r 220 --p-pooling-method pseudo \
  --o-table table.qza --o-representative-sequences rep-seqs.qza \
  --o-denoising-stats denoising-stats.qza

# drop samples below 2,000 reads (removes 3 of 942 -> 939 samples)
qiime feature-table filter-samples --i-table table.qza \
  --p-min-frequency 2000 --o-filtered-table table-filt.qza

qiime feature-table filter-features --i-table table-filt.qza \
  --p-min-frequency 10 --p-min-samples 2 --o-filtered-table table-final.qza

Stage 2

Stage 2: taxonomy and phylogeny

Assign names with a naive Bayes classifier against SILVA, then (for the metabolic step) re-classify against GTDB r207 so names match AGORA2. Build a tree for the diversity analyses.

classify and build treeQIIME 2

qiime feature-classifier classify-sklearn \
  --i-classifier silva-138-99-nb-classifier.qza \
  --i-reads rep-seqs.qza --o-classification taxonomy.qza

# re-classify against GTDB r207 for AGORA2 name compatibility (coverage 53% -> 95%)
qiime feature-classifier classify-sklearn \
  --i-classifier gtdb-r207-nb-classifier.qza \
  --i-reads rep-seqs.qza --o-classification taxonomy-gtdb.qza

qiime phylogeny align-to-tree-mafft-fasttree \
  --i-sequences rep-seqs.qza \
  --o-rooted-tree rooted-tree.qza --o-alignment aligned.qza

Stage 3

Stage 3: diversity and the community test

Core diversity metrics and the PERMANOVA that tests whether a grouping (arm, timepoint) separates the communities.

diversity + PERMANOVAQIIME 2

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza --i-table table-final.qza \
  --p-sampling-depth 15000 --m-metadata-file metadata.tsv \
  --output-dir core-metrics

qiime diversity beta-group-significance \
  --i-distance-matrix core-metrics/weighted_unifrac_distance_matrix.qza \
  --m-metadata-file metadata.tsv --m-metadata-column treatment_arm \
  --o-visualization permanova-arm.qzv --p-method permanova --p-pairwise

Stage 4

Stage 4: the mixed-effects differential abundance

This is the model of Equation (1): a per-child random intercept so the eight samples per child are not double-counted. We export the tables from QIIME, then fit the model in R with MaAsLin2 at both ASV and genus level. The core is one function call.

export tablesbash

qiime tools export --input-path table-final.qza --output-path exp-asv
biom convert -i exp-asv/feature-table.biom -o asv-table.tsv --to-tsv
qiime tools export --input-path table-genus-final.qza --output-path exp-genus
biom convert -i exp-genus/feature-table.biom -o genus-table.tsv --to-tsv

fit the per-child mixed modelR (MaAsLin2)

library(Maaslin2)

# MDCF-2 is the reference, so a negative treatment_arm coefficient = higher in MDCF-2
metadata$treatment_arm <- relevel(factor(metadata$treatment_arm), ref = "MDCF-2")

Maaslin2(
  input_data      = feat,                              # samples x features
  input_metadata  = meta,                              # arm, timepoint, participant_id
  output          = out_dir,
  fixed_effects   = c("treatment_arm", "timepoint"),   # the beta1, beta2 terms
  random_effects  = c("participant_id"),              # the u_i term: one per child
  reference       = c("treatment_arm,MDCF-2"),
  normalization   = "TSS", transform = "LOG",          # compositional handling
  min_prevalence  = 0.1, max_significance = 0.05)       # BH-FDR q < 0.05

Compositional differential abundance was cross-checked with ANCOM-BC2 (Equation 1 of Section 2.3) through the QIIME q2-composition plugin, giving the same qualitative result.

Stage 5

Stage 5: community metabolic modeling (MICOM)

One master script builds the community models, runs every metabolic analysis, and writes every Results II figure and table. The genome-scale optimization of Equation (3) is handled inside MICOM; our script sets it up, extracts the fluxes, fits the same mixed model on predicted metabolism, and traces the cross-feeding. The key loading step rebuilds the growth result directly from the three tables inside the QIIME artifact, which is version-proof:

load models + run the metabolic analysesPython (MICOM)

from micom.workflows import build, grow
from micom.workflows.results import GrowthResults
from micom import interaction as it

# build and grow the community models (tradeoff 0.3; 0.5 stalled the solver)
manifest = build(abundances, model_db="agora201_gtdb207_genus_1.qza",
                 out_folder="models", cutoff=1e-4, threads=8)
growth = grow(manifest, "models", medium="western_diet_gut_agora.qza",
              tradeoff=0.3, threads=8)

# rebuild GrowthResults from the 3 CSVs inside growth.qza (version-proof)
res = GrowthResults(growth_rates, exchange_fluxes, annotations)

# cross-feeding: directed hand-offs among the growth-associated taxa
ints = it.interactions(res, taxa=wlz_taxa, progress=False)   # provided / received
mes  = it.MES(res)                                           # which metabolites are traded

Community production is the export flux in the shared compartment; per-arm and per-timepoint tests reuse the mixed model of Equation (1) with Benjamini-Hochberg control (Equation 2). The full script (model build, SCFA test, 48-metabolite scan, MES, cross-feeding network, roles, activity, and every figure) is micom_master_make_figures.py in the bundle.

Stage 6

Stage 6: the minimal cross-feeding model

The acetate-to-butyrate competition model of Equations (4)-(5), solved to its steady state. Running it reproduces the sign of the empirical butyrate-abundance relationship.

steady state of the minimal modelPython

def steady(P, F, C, qA=1.0, Va=1.0, Ka=1.5, Vc=1.3, Kc=1.5,
           Yb=0.6, kB=0.9, Kb=1.0, D=0.25, iters=500, dt=0.02):
    """P acetate producers, F butyrate makers, C non-butyrogenic consumers."""
    A, B = 1.0, 1.0
    for _ in range(iters):
        dA = qA*P - Va*F*A/(Ka+A) - Vc*C*A/(Kc+A) - D*A    # Eq. (4)
        dB = Yb*Va*F*A/(Ka+A) - kB*C*B/(Kb+B) - D*B        # Eq. (5)
        A = max(A + dt*dA, 1e-6)
        B = max(B + dt*dB, 1e-6)
    return A, B    # as the consumer-dominated community grows, B* falls

The bundle

Bundle contents

The bundle (TexasiGEM_BestModel_Bundle.zip) is organized so a reader can run the pipeline end to end:

What each directory in the released bundle contains
Path	Contents
01_qiime/	All QIIME 2 commands as runnable shell scripts (Stages 1-3), with the exact parameters above
02_diffabund/	1_prepare_tables.sh (export) and 2_run_maaslin2.R (the per-child mixed model, ASV and genus)
03_micom/	micom_master_make_figures.py (full metabolic pipeline and all figures), plus the standalone extraction and figure scripts
04_model/	crossfeeding_model.py (the minimal acetate-to-butyrate model)
tables/	Intermediate outputs behind every figure: per-taxon exchange fluxes, the 48-metabolite scan with FDR, cross-feeding hand-offs and producer/consumer roles, growth rates
README.md	Run order, software versions, and which script and table produce each figure

Software: QIIME 2 (amplicon 2026.1), R with MaAsLin2, Python with MICOM on COBRApy, statsmodels and matplotlib. Each figure maps to one script and one released table, listed in the README.