Species-level analysis asks “is Lactobacillus rhamnosus
present?”
Strain-level analysis asks “which specific strain of L. rhamnosus —
GG, LC705, or GR-1?”
This distinction is clinically critical for probiotic work because two strains of the same species can have opposite effects. L. rhamnosus GG prevents antibiotic-associated diarrhea; L. rhamnosus LC705 does not. You cannot see this difference at species level.
| Method | Resolution | Can do strain-level? |
|---|---|---|
| 16S rRNA amplicon | Genus → Species | ✗ No — not enough sequence diversity |
| Shotgun metagenomics (WGS) | Species → Strain | ✓ Yes — reads whole genome |
| Long-read (PacBio/Oxford Nanopore) | Strain → Sub-strain | ✓ Yes — best for novel strains |
16S sequencing, which is what most public datasets use, physically cannot resolve strains because the 16S gene is too conserved across strains of the same species. You need whole-genome shotgun sequencing as the baseline requirement for everything below.
These compare sequencing reads against a database of known strain genomes using k-mer fingerprints.
StrainGE (Strain Genome Explorer)
- Best for: tracking known strains across longitudinal samples
- How it works: StrainGST scores each reference strain by comparing the
k-mer profile of the reference to the sample k-mers, reporting the
reference closest to the highest abundant strain in the sample. StrainGR
then characterizes SNVs (single nucleotide variants) between the
detected strain and its closest reference.
- Input: FASTQ reads + reference genome database
- Output: strain identity + nucleotide-level difference map
- GitHub: https://github.com/broadinstitute/StrainGE
StrainScan
- Best for: single-species strain identification
- Limitation: one species at a time — not multi-species
Sylph
- Best for: fast species and strain-level profiling
- Speed: ~30 seconds per sample vs MetaPhlAn’s 20 minutes
These assemble raw reads into longer contiguous sequences (contigs) first, then classify at strain level from the assembled genome.
inStrain
- Best for: measuring within-population diversity and strain
co-existence
- How it works: maps reads to a reference, then calls SNVs and measures
population-level allele frequencies — detects when multiple strains
co-exist in one sample
- Input: BAM alignment file + reference genome
- Output: strain-level SNV profiles, population diversity metrics
DESMAN (Disentangling Strains from
Metagenomes)
- Best for: de novo strain inference without a reference
- Works even without a known reference strain — reconstructs strain
haplotypes from co-varying SNVs
STRONG
- Best for: strain resolution from assembled metagenomes with
uncertainty quantification
PanTax (published January 2026 — most current)
- Best for: multi-species strain-level classification from both short
and long reads
- Key advantage over others: handles multiple species simultaneously,
not one at a time
- How it works: builds pangenome graphs (a compressed representation of
all known genomic variation within a species) and maps reads against
those graphs to determine which strain path best explains the observed
reads
- Handles both Illumina short reads and Oxford Nanopore/PacBio long
reads
- Published: Genome Research, January 2026
Centrifuge / Centrifuger
- Currently the only tools that handle both NGS and TGS reads and
perform strain-level taxonomic profiling for multiple species at a
time
- Fast enough for clinical turnaround
Raw Data Prerequisite:
Shotgun WGS (Illumina, minimum 10M reads per sample)
OR long-read (PacBio/Nanopore, minimum 1M reads per sample)
Reference Database (choose based on organism of interest):
UHGG (Unified Human Gastrointestinal Genome catalog)
→ 204,938 reference genomes, gut-specific
NCBI RefSeq bacterial genomes
→ comprehensive but not gut-specific
Custom strain database (if proprietary strains)
Read QC tools:
fastp → quality trimming, adapter removal
FastQC → quality reports
Alignment tools (needed by assembly-based methods):
Bowtie2 → short-read alignment to reference genomes
Minimap2 → long-read alignment
Assembly tools (for assembly-based strain inference):
MEGAHIT → fast de novo metagenomic assembly
metaSPAdes → high-quality assembly, slower
Core strain-level tools (pick based on use case):
StrainGE → known strain tracking across samples
inStrain → within-sample strain diversity
PanTax → multi-species, short+long read, newest
Centrifuger → fast multi-species strain profiling
Downstream strain abundance table:
→ format identical to species OTU table
→ rows = samples, columns = strain IDs
→ values = relative abundance (0–1)
→ plug directly into your existing pipeline steps 3–7Once you have a strain-level abundance table from any of the above tools, the rest of your pipeline stays identical — just with more columns:
Species-level: 1,656 features (genus.species)
Strain-level: 5,000–20,000 features (genus.species.strainID)
Step 2 feature filtering applies identically
Step 3 CLR transform applies identically
Step 4 ComBat batch correction applies identically
Step 5 missingness flags — more structural zeros expected
(strains are rarer than species — prevalence threshold
should drop from 10% to 2–5%)
Steps 6–7 unchangedThe main difference is that some strains are extremely rare — present in <1% of samples. The prevalence filter in step 2 becomes more important, not less.
This section explains every mathematical concept used in the pipeline, in plain language with formulas. Written for a math/CS undergraduate who knows linear algebra and basic statistics but has no biology background.
Microbiome data is compositional: the numbers for all species in one sample always add up to 1 (or 100%). This is called the constant sum constraint.
Imagine you have three species: A, B, C with abundances:
Sample 1: A=0.33, B=0.33, C=0.33 (all equal)
Sample 2: A=0.50, B=0.25, C=0.25 (A doubled)In sample 2, only A changed biologically. But the numbers for B and C also changed — not because B and C changed, but because they had to, since everything must sum to 1. This creates spurious (fake) correlations between species that have nothing to do with each other biologically.
Standard Pearson correlation between two compositional variables is always biased:
\[\text{Cov}(x_i, x_j) = -\frac{1}{p^2} \sum_{k \neq i,j} \text{Var}(x_k) \neq 0 \text{ even when truly independent}\]
In plain terms: even if species A and species B have no biological relationship, their correlation will be non-zero just because they’re both measured as fractions of a whole.
CLR moves data from the simplex space (where everything sums to 1) into real Euclidean space (where standard statistics work normally).
For a sample with species abundances x₁, x₂, …, xₚ:
\[\text{clr}(x_i) = \log(x_i) - \frac{1}{p}\sum_{j=1}^{p} \log(x_j)\]
The second term is the geometric mean of all abundances in log space.
After CLR, the row means are always exactly zero:
\[\sum_{i=1}^{p} \text{clr}(x_i) = 0\]
This is the mathematical proof it’s working correctly — our step 3 code checks this.
We cannot take log(0). So before the CLR we add a tiny number (1e-6) to all values:
X_plus = X + 1e-6
clr = log(X_plus) - mean(log(X_plus), axis=1)Removes systematic differences between studies (batch effects) while preserving disease signal.
ComBat models each feature (taxon) value as:
\[Y_{ijg} = \alpha_g + X\beta_g + \gamma_{ig} + \delta_{ig}\epsilon_{ijg}\]
Where: - \(Y_{ijg}\) = observed CLR value for sample j, in batch i, for taxon g - \(\alpha_g\) = overall mean for taxon g - \(X\beta_g\) = biological signal to preserve (disease label) - \(\gamma_{ig}\) = additive batch effect to remove (study-level mean shift) - \(\delta_{ig}\) = multiplicative batch effect to remove (study-level variance scaling) - \(\epsilon_{ijg}\) = residual noise
Think of it as: observed = biology + lab artifacts + noise. ComBat estimates the lab artifact terms (γ and δ) using an empirical Bayes approach (borrows information across all taxa to get stable estimates even when n is small), then subtracts them out, leaving only biology + noise.
Instead of estimating γᵢ separately for each taxon (which is noisy), it assumes all the γᵢ values across taxa come from a common distribution, and uses that distribution to regularise each individual estimate. This is why it works even when n < 20 samples per batch.
XGBoost builds many simple decision trees one after another. Each tree tries to correct the mistakes of all the trees before it. The final prediction is a weighted sum of all trees.
Step 1 — Start with a constant prediction
For multi-class, initialise all predictions to equal probability: 1/K
for K classes.
Step 2 — Compute residuals (what we got wrong)
For each sample, compute the difference between what we predicted and
what was actually true. These are called
pseudo-residuals or gradients:
\[g_i = \frac{\partial L(y_i, \hat{y}_i)}{\partial \hat{y}_i}\]
where L is the loss function (log-loss for classification).
Step 3 — Fit a tree to the residuals
Build a decision tree that predicts the residuals, not the original
labels. The tree splits are chosen to maximise the
gain:
\[\text{Gain} = \frac{1}{2}\left[\frac{G_L^2}{H_L + \lambda} + \frac{G_R^2}{H_R + \lambda} - \frac{(G_L+G_R)^2}{H_L+H_R+\lambda}\right] - \gamma\]
Where: - G = sum of gradients in a leaf - H = sum of second-order
gradients (hessians) in a leaf - λ = L2 regularisation (our
reg_lambda=1.0) - γ = minimum gain required to make a split
(our min_child_weight)
Step 4 — Add the tree to the ensemble (with shrinkage)
\[\hat{y}_i^{(t)} = \hat{y}_i^{(t-1)} + \eta \cdot f_t(x_i)\]
Where η is the learning rate (our learning_rate=0.1). A
smaller learning rate means each tree contributes less — the model
learns more slowly but generalises better.
Step 5 — Repeat for t = 1 to n_estimators
reg_alpha = 0.1 → L1: pushes some feature weights to exactly 0
(like LASSO — automatic feature selection)
reg_lambda = 1.0 → L2: keeps all weights small
(like Ridge — prevents any one feature dominating)
min_child_weight=5 → minimum sum of sample weights in a leaf
(prevents overfitting on small subgroups)
subsample = 0.8 → use 80% of samples per tree (like bagging)
colsample_bytree=0.6 → use 60% of features per tree
(reduces feature correlation between trees)Our dataset has:
IBS: 305 samples (27.7%)
Healthy_Control: 353 samples (32.1%)
Crohn Disease: 265 samples (24.1%)
Ulcerative Colitis: 177 samples (16.1%) ← smallest classA naive model could achieve 32% accuracy by always predicting “Healthy” — without learning anything. This is the class imbalance problem.
\[w_k = \frac{N}{K \cdot n_k}\]
Where: - N = total samples (880 training) - K = number of classes (4) - nₖ = samples in class k
For Ulcerative Colitis (the smallest class): \[w_{UC} = \frac{880}{4 \times 142} = 1.549\]
For Healthy_Control (the largest class): \[w_{HC} = \frac{880}{4 \times 282} = 0.780\]
UC samples count 1.549/0.780 = 2× as much as Healthy samples during training. This forces the model to pay equal attention to all classes regardless of their frequency.
We need an honest estimate of how the model will perform on new patients it has never seen. If we evaluate on the same data we trained on, the model looks better than it really is (overfitting).
A simple random split might, by chance, put all UC samples in training and none in test. Stratified splitting guarantees each class appears in both sets at the same proportion:
Before split: UC = 16.1% of all 1100 samples
After split: UC = 16.1% of train AND 16.1% of testWe went further and stratified on disease × dataset jointly, because IBS only comes from one study — we needed MARS samples in both train and test.
Instead of one train/test split, we split the training data into K=5 equal parts (folds):
Fold 1: [val] [train] [train] [train] [train]
Fold 2: [train] [val] [train] [train] [train]
Fold 3: [train] [train] [val] [train] [train]
Fold 4: [train] [train] [train] [val] [train]
Fold 5: [train] [train] [train] [train] [val]Train 5 models, each evaluated on a different validation fold. Average the 5 accuracy scores. This gives a more reliable performance estimate than one split and uses all the data for training.
Critical rule: the held-out test set (220 samples) is never used during cross-validation. It is touched exactly once — at the very end, after all model decisions are finalised.
If 90% of patients are healthy and 10% have Crohn’s, a model that always predicts “healthy” gets 90% accuracy — without learning anything useful.
For each class, compute:
\[\text{Precision}_k = \frac{\text{True Positives}_k}{\text{True Positives}_k + \text{False Positives}_k}\]
\[\text{Recall}_k = \frac{\text{True Positives}_k}{\text{True Positives}_k + \text{False Negatives}_k}\]
\[F1_k = 2 \times \frac{\text{Precision}_k \times \text{Recall}_k}{\text{Precision}_k + \text{Recall}_k}\]
Macro F1 = simple average of F1 across all classes. Equal weight to each class regardless of size. Our result: 0.864.
A single number summarising a multi-class confusion matrix. Ranges from -1 (always wrong) to +1 (perfect) with 0 = no better than random:
\[\text{MCC} = \frac{\sum_{k}\sum_{l}\sum_{m} C_{kk}C_{lm} - C_{kl}C_{mk}}{\sqrt{\sum_k(\sum_l C_{kl})(\sum_{k'\neq k}\sum_l C_{k'l})} \times \sqrt{\sum_k(\sum_l C_{lk})(\sum_{k'\neq k}\sum_l C_{lk'})}}\]
In plain terms: MCC measures the correlation between predicted labels and true labels across all classes simultaneously. Our result: 0.822 — substantially better than random (0.0).
For each class, plot the True Positive Rate vs False Positive Rate at every possible probability threshold. AUC = area under this curve.
SHAP (SHapley Additive exPlanations) assigns each feature a contribution score to each prediction, based on game theory.
The SHAP value φᵢ for feature i is:
\[\phi_i = \sum_{S \subseteq F \setminus \{i\}} \frac{|S|!(|F|-|S|-1)!}{|F|!} \left[f(S \cup \{i\}) - f(S)\right]\]
Where: - F = set of all features - S = a subset of features not including i - f(S) = model prediction using only features in S
Think of each feature as a player in a game. The game’s prize is the model’s prediction. SHAP calculates each player’s fair share of the prize by averaging their contribution across all possible orderings of players joining the game.
A positive SHAP value for “Escherichia coli” in the Crohn’s Disease class means: “having high E. coli abundance pushed this patient’s Crohn’s risk score upward.”
For each feature, take the average of the absolute SHAP values across all test samples:
\[\text{Importance}_i = \frac{1}{n}\sum_{j=1}^{n} |\phi_i^{(j)}|\]
The higher this number, the more that feature consistently influenced predictions — in either direction.
The VAE trains by minimising two terms simultaneously:
\[\mathcal{L} = \underbrace{\mathbb{E}_{q(z|x)}[\log p(x|z)]}_{\text{Reconstruction loss}} - \underbrace{\beta \cdot D_{KL}(q(z|x) \| p(z))}_{\text{KL divergence}}\]
\[\mathcal{L}_{\text{recon}} = \frac{1}{n}\sum_{i=1}^{n}(x_i - \hat{x}_i)^2\]
Measures how well the decoder rebuilt the original CLR values. Minimising this forces the network to preserve information through the bottleneck.
\[D_{KL}(q \| p) = -\frac{1}{2}\sum_{j=1}^{d}\left(1 + \log\sigma_j^2 - \mu_j^2 - \sigma_j^2\right)\]
Where μ and σ are the encoder’s output for each latent dimension. This term pulls the latent space toward a standard normal distribution N(0,1). Without it, the encoder would memorise training data rather than learning a smooth, generalisable representation.
We cannot backpropagate through a random sampling operation. So instead of sampling Z ~ N(μ, σ²) directly, we write:
\[Z = \mu + \sigma \cdot \epsilon, \quad \epsilon \sim \mathcal{N}(0, I)\]
The randomness is in ε (which has no parameters), so gradients flow through μ and σ normally.
| Concept | Formula/Method | Pipeline step | Undergrad owns |
|---|---|---|---|
| Compositionality | Simplex constraint | Background | ✓ Explain in report |
| CLR transform | log(xᵢ) − mean(log(x)) | Step 3 | ✓ Implement + verify |
| Pseudocount | x + 1e-6 before log | Step 3 | ✓ Implement |
| ComBat | Empirical Bayes batch model | Step 4 | ✓ Run + validate PCA |
| Missingness flags | Binary indicator per study | Step 5 | ✓ Build mask logic |
| Stratified split | Joint disease×dataset key | Step 6 | ✓ Implement + verify |
| Class weights | wₖ = N/(K·nₖ) | Step 7 | ✓ Compute + apply |
| XGBoost gain | G²/(H+λ) split criterion | Step 7 | ✓ Tune hyperparameters |
| F1, MCC, AUC | Precision/recall curves | Step 7 | ✓ Compute + report |
| SHAP | Shapley value attribution | Step 7 | ✓ Visualise + interpret |
| VAE loss | Recon + β·KL | Step 3b | ✓ Implement + train |
| Reparameterisation | Z = μ + σ·ε | Step 3b | ✓ Implement |