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.