https://huttenhower.sph.harvard.edu/biobakery_workflows/
Whole metagenome shotgun sequencing data can be processed through read-level quality control (KneadData), taxonomic profiling (MetaPhlAn), functional profiling (HUMAnN), and strain profiling (StrainPhlAn) to generate a report with publication-ready figures with two workflow commands.
-
Prepare the toy datasets
jhuang@WS-2290C:/mnt/md1/DATA/Data_Tam_Metagenomics_2026_Soil$ find . -name "*.fastq.gz"
#./biobakery_input/Soil_Loc4_2.fastq.gz
#./biobakery_input/Soil_Loc4_1.fastq.gz
#./biobakery_input/Soil_Loc1_1.fastq.gz
#./biobakery_input/Soil_Loc1_2.fastq.gz
-
Create Pseudo-replicates for Testing Pipelines by creating subsampled replicates
# For Soil_Loc1 - create two pseudo-replicates by random subsampling
# Install seqtk if not already installed
#conda install -c bioconda seqtk
# Create replicate 1 (50% of reads)
seqtk sample -s100 Soil_Loc1_1.fastq.gz 0.5 > Soil_Loc1_rep1_1.fastq
seqtk sample -s100 Soil_Loc1_2.fastq.gz 0.5 > Soil_Loc1_rep1_2.fastq
# Create replicate 2 (different 50% using different seed)
seqtk sample -s200 Soil_Loc1_1.fastq.gz 0.5 > Soil_Loc1_rep2_1.fastq
seqtk sample -s200 Soil_Loc1_2.fastq.gz 0.5 > Soil_Loc1_rep2_2.fastq
# Compress them
gzip Soil_Loc1_rep*_*.fastq
# Repeat for Soil_Loc4
seqtk sample -s100 Soil_Loc4_1.fastq.gz 0.5 > Soil_Loc4_rep1_1.fastq
seqtk sample -s100 Soil_Loc4_2.fastq.gz 0.5 > Soil_Loc4_rep1_2.fastq
seqtk sample -s200 Soil_Loc4_1.fastq.gz 0.5 > Soil_Loc4_rep2_1.fastq
seqtk sample -s200 Soil_Loc4_2.fastq.gz 0.5 > Soil_Loc4_rep2_2.fastq
gzip Soil_Loc4_rep*_*.fastq
-
Run docker
docker run -it \
-v /mnt/nvme4n1p1/biobakery_db:/biobakery_databases \
-v /mnt/md1/DATA/Data_Tam_Metagenomics_2026_Soil/biobakery_input_rep:/data \
biobakery/workflows:fixed \
/bin/bash
export BIOBAKERY_WORKFLOWS_DATABASES=/biobakery_databases
# ---- Configure databases (read-level quality control (1_KneadData), taxonomic profiling (2_MetaPhlAn), functional profiling (3_HUMAnN), and strain profiling (4_StrainPhlAn)) ----
# By default in the environment: 1_KneadData_databases 路径: /biobakery_databases/kneaddata_db_human_genome
# Check 2_MetaPhlAn_databases if correct
python3 -c "import metaphlan, os; print(os.path.join(os.path.dirname(metaphlan.__file__), 'metaphlan_databases'))"
#/usr/local/lib/python3.6/dist-packages/metaphlan/metaphlan_databases
ls -lh $(python3 -c "import metaphlan, os; print(os.path.join(os.path.dirname(metaphlan.__file__), 'metaphlan_databases'))")
# Check 3_HUMAnN
humann_config --print
#If not configured, using the following commands configuring them.
humann_config --update database_folders nucleotide /biobakery_databases/humann/chocophlan
humann_config --update database_folders protein /biobakery_databases/humann/uniref
humann_config --update database_folders utility_mapping /biobakery_databases/humann/utility_mapping
# By default in the environment: 4_StrainPhlAn_databases 路径: strainphlan_db_reference(empty) and strainphlan_db_markers (1.4G)
# If new running, optimally clean up the partial results from the failed run.
rm -rf /data/results/
rm -rf /data/results/*fastqc.zip _fastqc #IMPORTANT, so that no fastqc-related files existing under /data/results/
rm -rf /data/results/humann
# Run the workflow
biobakery_workflows wmgx \
--input /data \
--output /data/results \
--threads 64 \
--pair-identifier "_1"
biobakery_workflows wmgx_vis \
--input /data/results \
--output /data/results_vis \
--project-name wastewater_2026
-
Analyze metagenomics data from biobakery output using R (WITH REPLICATES) Project: Soil Metagenomics 2026 — Loc1 vs Loc4 comparison (Pseudo-replicates)
(r_env) Rscript analyze_biobakery_output.R
analyze_biobakery_output.R
Key Updates Made:
- Metadata Parsing: Updated to automatically detect the new pseudo-replicate naming convention (e.g.,
Soil_Loc1_rep1, Soil_Loc4_rep2) and extract both Location and Replicate information.
- Dynamic Excel Exports: Removed hardcoded sample names (
Soil_Loc1, Soil_Loc4). The script now dynamically calculates the mean abundance across replicates for each location to compute Diff and Log2FC.
- MaAsLin2 Setup: Removed all the leftover code from the hospital wastewater experiment (Treatment/TimePoint). Set up clean MaAsLin2 models to test the
Location effect for both species and pathways.
- Beta Diversity PCoA: Added a Principal Coordinates Analysis (PCoA) plot, which is the standard and most informative way to visualize beta diversity when you have replicates.
- Pathway File Path: Updated the HUMAnN pathway file path to point to the new
biobakery_input_rep directory.
⚠️ Important Statistical Note on Pseudo-Replicates
Pseudo-replicates created by subsampling FASTQ files are technical replicates, not biological replicates. They help you understand the technical variance of your sequencing pipeline and allow statistical models to run, but they do not capture true biological variance between different soil cores. Therefore, while MaAsLin2 and PERMANOVA will run and likely yield highly significant p-values, you should interpret these results as “technically distinct” rather than “biologically significant” until you sequence true biological replicates.
Yes, exactly! The qval is the FDR-adjusted p-value (specifically using the Benjamini-Hochberg procedure, since you set correction = "BH").
A qval < 0.05 means that after accounting for the fact that you are testing thousands of species/pathways simultaneously (multiple testing correction), this result has a False Discovery Rate of less than 5%. In other words, it is statistically significant.
Here is a breakdown of how to read your MaAsLin2 output table, along with an interpretation of your specific results.
📖 Glossary of Your Output Columns
| Column |
Meaning |
feature |
The species or pathway being tested (e.g., Bradyrhizobium.diazoefficiens). |
metadata |
The variable you are testing (in this case, Location). |
value |
The group being compared to the reference. Since your reference was Loc1, Loc4 means “Loc4 compared to Loc1”. |
coef |
The Coefficient (Effect Size). A negative value means the feature is less abundant in Loc4 than Loc1. A positive value means it is more abundant in Loc4 than Loc1. |
stderr |
The standard error of the coefficient. |
pval |
The raw, unadjusted p-value. |
qval |
The FDR-adjusted p-value. This is the most important metric for significance. |
N |
Total number of samples in the model (4 in your case: 2 reps per location). |
N.not.zero |
How many samples actually contained this feature. If N=4 and N.not.zero=2, it means the species was only detected in 2 of the 4 samples (e.g., present in Loc1, but completely absent in Loc4). |
🔬 Interpreting Your Top Hits
Because you set transform = "NONE", the coef represents the raw difference in mean relative abundance between Loc4 and Loc1.
1. Bradyrhizobium.diazoefficiens
coef: -0.0475qval: 0.021 (Significant!)N.not.zero: 2- Interpretation: This nitrogen-fixing bacterium is significantly depleted in Location 4 compared to Location 1. It was likely only detected in your Loc1 replicates (hence
N.not.zero = 2).
2. Candidatus Nitrosocosmicus oleophilus
coef: +0.7463qval: 0.021 (Significant!)N.not.zero: 4- Interpretation: This ammonia-oxidizing archaeon is significantly enriched in Location 4 compared to Location 1. It was detected across all 4 of your pseudo-replicates, but its abundance is substantially higher in Loc4.
3. Dyella marensis
coef: -0.2885qval: 0.053 (Marginally non-significant)- Interpretation: It is less abundant in Loc4, but because the
qval is just over 0.05, it does not strictly pass the 5% FDR threshold. It is a “trend” but you should not claim it as definitively different.
⚠️ A Crucial Reminder on “Pseudo-Replicates” for Publishing
When you write your methods or results section, you must be transparent about how these replicates were generated.
Because your replicates were created by computationally splitting the same FASTQ file (using seqtk) rather than sequencing independent soil cores, your statistical models (MaAsLin2) are calculating “technical variance,” not “biological variance.”
- What this means: The
qval proves that the differences between Loc1 and Loc4 are much larger than the sequencing noise/technical variance of the machine. It proves the differences are technically real and robustly detectable by your pipeline.
- What it doesn’t mean: It does not prove that every square foot of soil in Loc1 differs from Loc4, because you only sampled one physical soil core per location.
How to phrase this in a publication:
“To evaluate the technical robustness of differential abundance between the two sampling sites, in silico pseudo-replicates were generated via random subsampling of sequencing reads. Differential abundance testing was performed using MaAsLin2. While these models successfully control for technical and sequencing variance (yielding significant FDR-adjusted q-values < 0.05), the lack of independent biological replicates means these results reflect localized, site-specific differences rather than broad population-level biological variance."