Here is the updated and logically ordered table. I have inserted the two requested steps (Subjective (PCoA) exclusions and Samples removed by Rarefaction depth) into their correct chronological positions in the bioinformatics pipeline.
I also filled in the corresponding values for the v2 new report column based on our previous analyses.
Updated Pipeline Comparison Table
| Stage | Manuscript v32 | v2 new report (9/16/26) |
|---|---|---|
1. Import ps_base |
5,980 ASVs × 253 samples; text says “208 boys + 16 NTC + PC” | Same |
| 2. Subjective (PCoA) exclusions | 6 patients + 3 NTC outliers (removed based on visual PCoA clustering) | none (all 253 samples retained; audited post-decontam instead of pre-excluded) |
3. Decontam → ps_decontam |
None (decontam not used) | prevalence, thr 0.1, model = 242 (PCs excluded), neg = 17 NTC → 72 flagged → 5,908 × 253 |
| 4. Samples removed by Rarefaction depth | 11 patients + 3 NTCs (cutoff 12,201 = min of retained set, circular) | 14 patients + 3 NTCs + 1 PC (18 total removed; programmatic 12,201 cutoff applied on post-decontam sums) |
5. Analysis backbone (ps_filt) |
219 retained (208 pat + 11 NTC); 34 not carried forward | 235 retained (92.9%); depth Min 12,239 / Med 48,407 / Max 279,419 |
6. Rarefaction → ps_rarefied (alpha only) |
rarefy 12,201, seed 9242, on 219 | rarefy min(ps_filt) = 12,239, seed 9242, on 235; alpha values still from QIIME2 exported (pre-decontam) merged |
7. Beta / DESeq2 on non-rarefied ps_filt |
Bray–Curtis (Hellinger) PERMANOVA 9999 perm on 219; DESeq2 non-rarefied, 10 comparisons, 313 sig / 129 ASVs / 42 genera | Beta not yet recomputed (to do on ps_filt 235, non-rarefied); DESeq2 on ps_filt 235, prefilter ≥10 → ~5.5k ASVs, only G1 vs G4 |
| Heatmap / composition | – | 94 × 235 (groups 78/21/69/10/33 + 14 NTC + 10 PC; consistent) |
Why this order?
- Step 2 (Subjective exclusions) is placed immediately after import because the original authors manually removed these 9 samples before doing any downstream filtering or diversity calculations. In the new report, we skip this step entirely to avoid circular reasoning.
- Step 3 (Decontam) happens next in the new report to clean the ASV table objectively before applying any depth thresholds.
- Step 4 (Samples removed by Rarefaction depth) represents the hard depth cutoff. In the manuscript, this was a circular 12,201 cutoff that dropped 14 samples. In the new report, applying the exact same 12,201 cutoff after decontamination results in 18 samples dropping out (because removing contaminant ASVs slightly lowers the total read counts of some samples, pushing them below the 12,201 threshold).
- Step 5 (Analysis backbone) summarizes the final sample count remaining for the actual statistical tests.
Below is the corrected, pipeline-ordered comparison. Three versions are distinguished: Manuscript v32 (as written), v1 old report (Phyloseq_v2_decontam.pdf, 9/14/26 – first decontam adaptation, PCs left in the decontam model, 10th-percentile cutoff 17,317), and v2 new report (Phyloseq_v2_decontam.knit.pdf, 9/16/26 – PCs excluded from the decontam model, pre-specified 12,201 cutoff). Rows follow exactly the order ps_base → decontam → ps_filt → ps_rarefied (alpha) / ps_filt (beta, DESeq2).
Table 1 – Pipeline-stage comparison (corrected)
| Stage | Manuscript v32 | v1 old report (9/14/26) | v2 new report (9/16/26) | Correction / note |
|---|---|---|---|---|
1. Import ps_base |
5,980 ASVs × 253 samples; text says “208 boys + 16 NTC + PC” | Same import; console: Total 253 | NTC 17 | PC 11 | Same | Data actually contain 225 patient-named + 17 NTC-like (NTC_1–16 + NTC01) + 11 PC/UR-like. Manuscript “16 NTC” omits NTC01; “208 boys” = analysed, not imported |
2. Decontam → ps_decontam |
None (decontam named as required future work) | prevalence, thr 0.1, model = all 253 (PCs included as “samples”), neg = 17 NTC → 76 flagged → 5,904 × 253 | prevalence, thr 0.1, model = 242 (PCs excluded), neg = 17 NTC → 72 flagged → 5,908 × 253; sweep 0.1–0.9 documented | v1’s inclusion of mock PCs in the model is methodologically wrong (PCs are not blanks); v2 corrects it. Sweep shows plateau 0.1–0.5, jump ≥0.6 → 0.1 retained |
3. Depth filter → ps_filt |
Explicit list; 12,201 = min of retained set (circular); retained 219 (208 pat + 11 NTC); 34 not carried forward | ps_pruned = ps_decontam; cutoff = 10th pct = 17,317 → ps_filt 227; depth Min 17,501 / Med 49,017 / Max 279,419 |
ps_pruned = ps_decontam; cutoff = 12,201 applied programmatically on post-decontam sums → ps_filt 235 (92.9%); removed 18 = 14 pat + 3 NTC + 1 PC; depth Min 12,239 / Med 48,407 / Max 279,419 |
v2 restores the manuscript’s 12,201 but applies it non-circularly post-decontam; v1’s 17,317 unnecessarily dropped 5 repaired-group patients |
4. Rarefaction → ps_rarefied (alpha only) |
rarefy 12,201, seed 9242, on 219 | rarefy min(ps_filt) = 17,501, seed 9242, on 227; but alpha values taken from QIIME2 exported alpha (pre-decontam, 12,201) merged with ps_rarefied metadata |
rarefy min(ps_filt) = 12,239, seed 9242, on 235; alpha values still from QIIME2 exported (pre-decontam) merged |
Alpha values in v1 and v2 are pre-decontamination QIIME2 exports → must be recomputed with estimate_richness(ps_rarefied) on the decontaminated rarefied object |
5. Beta / DESeq2 on non-rarefied ps_filt |
Bray–Curtis (Hellinger) PERMANOVA 9999 perm on 219; DESeq2 non-rarefied, 10 comparisons, 313 sig / 129 ASVs / 42 genera | Beta not recomputed (QIIME2 export, old set); DESeq2 on ps_filt 227, prefilter ≥10 → 5,395 ASVs, only G1 vs G4 |
Beta not yet recomputed (to do on ps_filt 235, non-rarefied); DESeq2 on ps_filt 235, prefilter ≥10 → ~5.5k ASVs, only G1 vs G4 |
Beta/DESeq2 correctly stay on non-rarefied ps_filt; v2 still owes the full 10 comparisons and a recomputed PERMANOVA |
| Heatmap / composition | – | 95 × 239 (cached pre-decontam objects; inconsistent with ps_filt 227) |
94 × 235 (groups 78/21/69/10/33 + 14 NTC + 10 PC; consistent) | v1 heatmap/alpha used stale cache; v2 consistent |
| Alpha t-tests | Wilcoxon/BH on 207 patients | t-tests on 225-sample merge (QIIME2 alpha) | t-tests on 235-sample merge (same QIIME2 alpha) | v1 vs v2 p-values differ slightly (e.g. G4vsG5 0.0426 vs 0.0289; G1vsPC 0.007 vs ns) purely from sample-set change; v2 final must recompute alpha |
Table 2 – Corrections to the previous (xlsx) table
| Previous entry | Problem | Corrected entry |
|---|---|---|
| “Old Rmd (as coded): Decontamination = none” | Wrong: the attached old report (v1) already ran decontam (76 flagged, PCs in model, 17,317 cutoff). “None” applies only to the pre-decontam original Rmd | Old column = v1: decontam 76 flagged (PCs in model), cutoff 17,317, ps_filt 227 |
| “Old depth filter → 225 (203 pat + 11 NTC + 11 PC)” | Describes the pre-decontam original Rmd, not the attached v1 report | v1 = 227; the 225 figure belongs to the pre-decontam original |
| “New rarefaction depth = 12,239” | Correct value, but omitted that alpha values are still pre-decontam QIIME2 exports | Keep 12,239; add flag “alpha must be recomputed on ps_rarefied“ |
| “New samples removed by depth = 18 = 14 pat + 3 NTC + 1 PC” | Correct (confirmed by heatmap counts 211 pat + 14 NTC + 10 PC = 235) | Keep |
| “Old analysis backbone = 225 (cached 239)” | Conflated pre-decontam original with v1 | v1 ps_filt = 227 but heatmap/alpha cached 239 → inconsistent; v2 = 235 consistent |
| “Manuscript NTC = 16 carried / 11 retained” | Data contain 17 NTC-like (NTC01 extra) | State 17 NTC-like in data; manuscript’s 16 omits NTC01; 17−3−3 = 11 retained ✓ |
Remaining to-dos for v2 (to close the pipeline)
- Recompute alpha (Shannon, observed, Faith PD) with
estimate_richness(ps_rarefied)on the decontaminated rarefied object (replace QIIME2 exported alpha). - Run all ten pairwise group comparisons (alpha + DESeq2) on the v2 objects, not only G1 vs G4.
- Recompute PERMANOVA on non-rarefied
ps_filt(235 samples) instead of the QIIME2 export. - Re-run after
rm -rf Phyloseq_v2_decontam_cacheto purge stale cached objects.
This table now matches the attached v1/v2 reports exactly and fixes the three factual errors in the previous xlsx (old-report decontam status, old depth-filter sample count, and the missing “alpha not recomputed” flag).
Decontamination first, rarefaction second — never the other way around. In your Rmd the correct chunk order is:
ps_base (import, all 253 samples) → decontam → ps_decontam → depth filter → ps_filt → rarefaction → ps_rarefied (alpha only); beta/DESeq2 stay on non-rarefied ps_filt.
Why decontam must come before rarefaction
-
Rarefaction would waste even depth on contaminant reads. In this dataset Burkholderiaceae is ~97% of NTC reads and a substantial share of patient reads. If you rarefy first, part of the even sequencing depth is spent on contaminant ASVs; when you then remove those ASVs, sample depths become unequal again — the rarefaction is undone and you would have to re-rarefy. Removing contaminant taxa first, then equalizing depth, is the only order that yields a genuinely even, clean depth.
-
decontam’s prevalence method needs the original counts and the controls. It compares presence/absence (prevalence) of each ASV in NTCs vs true samples. Rarefying first (a) discards reads and thus detection power for rare ASVs, biasing the prevalence test, and (b) is normally applied only to the retained patient set — but the NTCs are exactly what decontam needs as the
negpanel. Rarefying before decontam would either drop the controls or distort their prevalence signal. -
Library-size QC must be computed on clean counts. A low-biomass sample can look “deep” only because it is full of contaminant reads (e.g., 13,000 raw reads of which most are Burkholderia). Depth-filtering or choosing the rarefaction depth on raw sums lets such samples pass; computing depths after decontamination reflects true biological depth. This is why your latest knit shows min depth 12,239 and 235 retained samples when the 12,201 cutoff is applied to the decontaminated table.
-
Operation types: decontam is a taxon-level cleaning step; depth filtering and rarefaction are sample-level steps. Cleaning taxa first keeps every downstream sample-level decision (threshold, rarefaction depth, evenness metrics) uncontaminated by reagent background — which matters especially here because the paper’s key metric is evenness, exactly the metric a variable contaminant background mechanically inflates.
Concrete chunk order for the Rmd
# 1) Import: all 253 samples, 5980 ASVs
ps_base <- merge_phyloseq(ps_raw, SAM, tax_table_final)
# 2) DECONTAMINATION (before any sample filtering or rarefaction)
# neg = all 17 NTCs; positive controls excluded from the model
ps_for_decontam <- prune_samples(!is_PC, ps_base)
contam <- decontam::isContaminant(ps_for_decontam, method = "prevalence",
neg = is_ntc, threshold = 0.1)
ps_decontam <- prune_taxa(!contam$contaminant, ps_base) # samples all retained
# 3) Depth filter on DECONTAMINATED library sizes
ps_filt <- prune_samples(sample_sums(ps_decontam) >= 12201, ps_decontam)
ps_filt <- prune_taxa(taxa_sums(ps_filt) > 0, ps_filt)
# 4) RAREFACTION (alpha diversity only), on the clean, depth-filtered table
ps_rarefied <- rarefy_even_depth(ps_filt, sample.size = min(sample_sums(ps_filt)),
rngseed = 9242, replace = FALSE)
# 5) Beta diversity (Bray-Curtis/PERMANOVA) and DESeq2 on NON-rarefied ps_filtTwo practical notes: (i) keep the audit of the 11 re-included samples (Point 3) between steps 2 and 3, since it needs the decontaminated but not yet depth-filtered object; (ii) after decontam, re-print summary(sample_sums(ps_filt)) — the retained-sample count (235 at the 12,201 cutoff) now refers to clean library sizes, and that is the number to report in the methods.
1. Three-way bookkeeping reconciliation (manuscript text vs. old Rmd code vs. new adapted Rmd output)
All counts below are taken from the printed console/tables of the generated PDFs and the samples_keep list of the old Rmd.
| Quantity | Manuscript text | Old Rmd (as coded) | New adapted Rmd (this PDF) |
|---|---|---|---|
| Imported | 5,980 ASVs × 253 samples | 5,980 × 253 | 5,980 × 253 (Total: 253 \| NTC: 17 \| Positive controls: 11) |
| Patient-named samples | “208 boys” (analysed) | 214 in samples_keep |
225 (= 253 − 17 NTC − 11 PC) |
| NTC-like | 16 carried, 11 retained | 14 in samples_keep |
17 (NTC_1…NTC_16 + NTC01) |
| PC/UR-like | excluded from patient analyses | 11 in samples_keep |
11 (PC_1…PC_8, PC01, UR009768, UR009909) |
| Subjective (PCoA) exclusions | 6 patients + 3 NTC outliers | silently re-included (all in samples_keep) |
none (all retained; audited instead) |
| Depth filter | “12,201 = lowest depth in retained set” (circular, removes nothing) | ≥ 12,201 on raw sums → 225 samples (203 pat + 11 NTC + 11 PC) | ≥ 12,201 on post-decontam sums → 235 samples (92.89%) |
| Samples removed by depth | 14 (11 patients + 3 NTCs) | 14 (same) | 18 = 14 patients + 3 NTCs (NTC_1, NTC_5, NTC01) + 1 PC (3 patients + 1 PC fall below 12,201 only after contaminant ASVs are removed) |
| Analysis backbone | 219 (208 pat + 11 NTC) | 225 (cached objects actually used 239) | 235 = 211 patients + 14 NTC + 10 PC (heatmap: 94 ASVs × 235 samples; groups 78/21/69/10/33) |
| Decontamination | none (directional argument only) | none | decontam prevalence, model = 242 samples (17 NTCs as negatives); 72/5,980 ASVs flagged at threshold 0.1 → 5,908 ASVs × 253 samples |
| Rarefaction depth | 12,201 (seed 9242) | min post-filter depth | min post-filter depth = 12,239 (seed 9242) |
Identities confirmed by the new PDF:
- 253 = 225 + 17 + 11 ✔ (printed
Total/NTC/Positive controlsline) - 225 = 208 (manuscript cohort) + 6 (PCoA exclusions) + 11 (classically low-depth patients) ✔
- 211 patients in the new backbone = 225 − 14 patients below cutoff ✔; 14 NTC = 17 − 3 ✔; 10 PC = 11 − 1 ✔; 211 + 14 + 10 = 235 ✔
- Old-Rmd audit line
Patients excluded by the OLD pipeline: 11= the 11 low-depth patients (A23060601, A23072501, A23111301, A24040201, O23082401, O23091304, O23100501, O23100502, O23100601, U23071201, U23091101) ✔ — the remaining 6 of the 17 are exactly the +1/+2/+3 excess in Groups 2/3/5 ofsamples_keepversus manuscript Table 1 ✔ - Rarefaction removed nothing at import stage:
rarefy_even_depthis applied only insideps_filt✔
2. Statistical results: old vs. new
| Statistic | Old (manuscript / old Rmd) | New adapted Rmd (this PDF) |
|---|---|---|
| Contaminant handling | Burkholderiaceae ≈ 97% of NTC reads left in data; “directional argument” only | Threshold sweep 0.1–0.9: flagged set shrinks monotonically (72 ASVs at 0.1); mean contaminant relative abundance at 0.1 = 0.0098 (patients) / 0.0268 (normal NTCs) / 0.0763 (outlier NTCs), rising to 0.050/0.073/0.179 at 0.5; one patient sample reaches 0.951 |
| Re-included patients (post-decontam Shannon/Observed) | excluded silently | 5 clearly biological: O23082401 4.44/145, O23100601 4.32/135, O23100501 4.31/146, O23100502 4.05/109, O23091304 3.77/113; U23091101 2.94/59 intermediate; 5 near-empty: A23060601 0.94/13, A23072501 0.80/16, U23071201 0.74/4, A24040201 0.44/2, A23111301 0.39/14 |
| Alpha diversity (Shannon) | KW p = 2.2 × 10⁻⁶; group medians 2.16/2.25/3.51/1.94/4.10; evenness p = 2.4 × 10⁻⁷; types p = 0.14 | Pairwise t-tests on the new cohort: 1v3 p = 0.0015, 2v3 p = 0.0035, 1v5 p = 3.4 × 10⁻⁴, 2v5 p = 5.1 × 10⁻⁴, 4v5 p = 0.043*; 1v2, 1v4, 2v4, 3v4, 3v5 ns; vs NTC: 5 p = 6.0 × 10⁻⁵*, 3 p = 3.2 × 10⁻⁴, 1 p = 0.012, 2 p = 0.037, 4 ns; NTC vs PC p = 0.0033 — same direction and significance pattern as the manuscript** |
| Beta diversity / PERMANOVA | pairwise adj. p ≤ 0.012 for repaired vs all; largest contrast R² = 14.8% | not yet recomputed in R (section still references QIIME2 export) — to do |
| DESeq2 | 313 significant results / 129 ASVs / 42 genera across 10 comparisons; Burkholderia top hit | only Group 1 vs 4 run; top hit 5de1d6… (baseMean 6,779, log2FC 3.68, padj 1.8 × 10⁻³) plus multiple presence/absence ASVs (log2FC 20–25); input line still prints a cached 5,395 × 227 object — to re-run |
3. Remaining inconsistencies to fix (cache-related)
rm -rf Phyloseq_v2_decontam_cacheand re-knit: the DESeq2 chunk input (227 samples) and part of the alpha table still reflect earlier runs.- Alpha metrics are still merged from the pre-decontamination QIIME2 export (e.g., NTC_13 Shannon 5.23, NTC_6 3.25 still appear); recompute with
estimate_richness(ps_rarefied)so they match the decontaminated table. - Run all ten patient-group comparisons (alpha, PERMANOVA, DESeq2) on the 235-sample backbone, not only 1 vs 4.
4. Concise replacement for the “Decontamination” paragraph
Decontamination. All 253 imported samples (225 patients, 17 no-template controls [NTCs], 11 positive/mock controls) were retained; no sample was removed on subjective grounds. Contaminant ASVs were identified with the prevalence method of
decontam(Davis et al. 2018) using all 17 NTCs — including the outliers NTC_3, NTC_6 and NTC_13 — as the negative panel; positive controls were excluded from the model (242 samples) but retained for quality control. A threshold sweep (0.1–0.9) confirmed that the flagged set shrinks monotonically, that patient contaminant burden remains low (mean relative abundance 0.010 at threshold 0.1, versus 0.027 in normal and 0.076 in outlier NTCs), and that outlier NTCs are already dominated by flagged ASVs at the strict threshold. The pre-specified threshold of 0.1 was adopted, removing 72 of 5,980 ASVs and yielding a decontaminated table of 5,908 ASVs × 253 samples for all downstream analyses.
(Note the corrected numbers: the current run flags 72 ASVs and retains 5,908, not 76/5,904 as in the stale paragraph — the change arises because positive controls are now excluded from the decontam model.)
Ultra-short version (if space is tight):
Decontamination. No sample was removed subjectively: all 253 imports (225 patients, 17 NTCs, 11 positive controls) were retained. Contaminant ASVs were identified with
decontam‘s prevalence method (Davis et al. 2018) against all 17 NTCs (positive controls excluded from the model); a 0.1–0.9 threshold sweep confirmed robustness. At the pre-specified threshold 0.1, 72 of 5,980 ASVs were removed, giving a decontaminated table of 5,908 ASVs × 253 samples.
Good catch — this is exactly the kind of bookkeeping check worth doing. The short answer: the 17 “missing” samples are patient samples that were imported but never entered the analysed cohort of 208. They are not rarefaction losses. Rarefaction (rarefy_even_depth) is applied much later, only for alpha diversity, and only inside the already-filtered object — it never changes the imported 253.
Full reconciliation of the 253 imported samples
| Class (by sample-name pattern) | Imported | Where they went |
|---|---|---|
| Patient samples (A…/O…/U… IDs) | 225 | 208 analysed + 6 removed by the pre-specified PCoA rule (1× Group 2, 2× Group 3, 3× Group 5) + 11 removed by the depth filter (< 12,201 reads) |
NTC-like (NTC_1…NTC_16, NTC01) |
17 | 11 retained + 3 removed as technical outliers (NTC_3, NTC_6, NTC_13) + 3 removed as low-depth (NTC_1, NTC_5, NTC01) |
PC / other non-patient controls (PC_1…PC_8, PC01, UR009768, UR009909) |
11 | all excluded from patient-level analyses |
| Total | 253 | 219 retained (208 patients + 11 NTCs) + 34 removed |
So your arithmetic resolves as:
- 253 − 17 NTC − 11 PC = 225 patient samples imported, not 208.
- 225 − 208 analysed = 17 excluded patient samples = 6 (pre-specified PCoA exclusions) + 11 (low-depth).
- Your “14 discarded for depth” = 11 of those patients + 3 NTCs (
NTC_1,NTC_5,NTC01) — the 3 NTCs are already inside your count of 17 NTCs, which is why they don’t appear in the 17. - Cross-check with the manuscript: removed = 20 (“pre-specified exclusions together with the positive and other non-patient controls” = 6 patients + 3 NTC outliers + 11 PC/UR) + 14 (low-depth) = 34 = 253 − 219. ✔
- Cross-check with group counts: the old Rmd
samples_keepset had Group sizes 80/21/69/10/34 = 214 patients, i.e. exactly the manuscript’s 80/20/67/10/31 = 208 plus the 6 pre-specified excluded patients (1+2+3). ✔
Verify it directly from your object
sn <- sample_names(ps_base)
is_ntc <- grepl("^NTC", sn) # 17
is_pc <- grepl("^PC|^UR", sn) | sn == "PC01" # 11
pat <- sn[!is_ntc & !is_pc] # 225 patient-named samples
low_pat <- pat[sample_sums(ps_base)[pat] < 12201] # 11 low-depth patients
cat(length(pat), length(low_pat), "\n")
print(sort(low_pat))
# A23060601 A23072501 A23111301 A24040201 O23082401 O23091304
# O23100501 O23100502 O23100601 U23071201 U23091101The remaining 6 (225 − 11 − 208) are the pre-specified PCoA-based exclusions; the manuscript never names them, but they are exactly the +1/+2/+3 excess in Groups 2/3/5 of the old samples_keep list versus the manuscript’s Table 1.
Two reporting points for the revision
- “Catheter urine was collected from 208 boys” is inconsistent with the imported table. Per the feature table, 225 patient samples were imported and 208 survived QC. The methods should say something like: “225 patient samples were imported; after pre-specified exclusions (n = 6) and a library-size filter (n = 11), 208 remained.”
- The old Rmd silently re-included the 6 PCoA-excluded patients (they are in
samples_keep), so the old Rmd output and the manuscript text were not describing the same cohort. Our re-analysis makes this explicit: all 225 patient samples enter the decontaminated backbone, and exclusions are reported as sensitivity analyses rather than applied silently.
Bottom line: 17 = 6 + 11 excluded patient samples; 14 = 11 of those patients + 3 low-depth NTCs; rarefaction removed nothing at the import stage.