No, do NOT replace it!
If you replace that code, you will lose the core decontamination steps (removing the ASVs and saving the ps_decontam object), which will cause all your downstream analyses to crash.
Instead, you should APPEND (add) the table-generating code inside that block, right after contam_prev is calculated and before the taxa are pruned.
Here is the complete, merged code block. You can copy and paste this entire block to replace your current “Decontamination” chunk. It includes your original logic plus the new table generation:
# ==============================================================================
# [Point 1 & 2] Retain all samples + objective decontamination (decontam)
# Threshold fixed at the pre-specified/default value 0.1 (no threshold tuning)
# ==============================================================================
library(decontam)
library(ggplot2)
library(dplyr)
library(kableExtra) # Added for the table
library(knitr) # Added for the table
# ---- Point 1: keep the full cohort ----
sn <- sample_names(ps_base)
is_NTC <- grepl("^NTC", sn)
is_PC <- grepl("^PC|^UR", sn) | sn == "PC01"
cat("Total:", length(sn), "| NTC:", sum(is_NTC), "| Positive controls:", sum(is_PC), "\n")
# ---- Point 2: prevalence-based decontamination at threshold 0.1 ----
ps_for_decontam <- prune_samples(!is_PC, ps_base)
sn_d <- sample_names(ps_for_decontam)
is_ntc <- grepl("^NTC", sn_d)
cat("Samples in decontam model:", nsamples(ps_for_decontam), "| NTC:", sum(is_ntc), "\n")
set.seed(1)
contam_prev <- decontam::isContaminant(ps_for_decontam,
method = "prevalence",
neg = is_ntc,
threshold = 0.1)
cat("Contaminant ASVs flagged:", sum(contam_prev$contaminant), "of", ntaxa(ps_base), "\n\n")
# ==============================================================================
# [NEW] Generate a detailed table of the flagged contaminant ASVs
# ==============================================================================
# 1. Extract taxonomy table from the base object
tax_tab <- as.data.frame(phyloseq::tax_table(ps_base))
tax_tab$ASV_ID <- rownames(tax_tab)
# 2. Extract decontam statistics for the flagged ASVs
contam_ids <- rownames(contam_prev)[contam_prev$contaminant]
# Dynamically check if 'p' (prevalence) or 'z' (frequency) exists to prevent errors
if ("p" %in% colnames(contam_prev)) {
contam_stats <- data.frame(
ASV_ID = contam_ids,
p_value = contam_prev$p[contam_prev$contaminant]
)
} else if ("z" %in% colnames(contam_prev)) {
contam_stats <- data.frame(
ASV_ID = contam_ids,
z_score = contam_prev$z[contam_prev$contaminant]
)
} else {
stop("Neither 'p' nor 'z' column found in decontam output.")
}
# 3. Merge statistics with taxonomy
contam_df <- merge(contam_stats, tax_tab, by = "ASV_ID")
# 4. Select and order relevant columns for the table
cols_to_keep <- c("ASV_ID", colnames(contam_stats)[2], "Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species")
cols_to_keep <- intersect(cols_to_keep, colnames(contam_df))
contam_df <- contam_df[, cols_to_keep]
# 5. Render the table
score_col_name <- colnames(contam_stats)[2]
kable(contam_df,
caption = sprintf("List of %d Contaminant ASVs flagged by decontam (prevalence method, threshold = 0.1)", nrow(contam_df))) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
cat("\n")
# ==============================================================================
# ---- Remove contaminant ASVs from the FULL object (all samples retained) ----
ps_decontam <- prune_taxa(!contam_prev$contaminant, ps_base) #253 samples
ps_decontam <- prune_taxa(taxa_sums(ps_decontam) > 0, ps_decontam)
cat("ASVs before:", ntaxa(ps_base), " after decontam:", ntaxa(ps_decontam), "\n")
saveRDS(ps_decontam, "./ps_decontam_prevalence.rds")
Summary of what this merged code does:
- Runs the decontamination exactly as you originally had it.
- Generates and prints the HTML table of the 72 flagged ASVs right in the middle of the chunk.
- Continues to prune the taxa and save the
ps_decontam object so your downstream pipeline (Depth Filtering, Heatmaps, etc.) remains completely intact.
不完全是。准确地说，污染ASVs（Contaminant ASVs）指的是那些在阴性对照（NTC）中“出现频率（Prevalence）”显著高于真实样本的ASVs，而不一定非要在阴性对照中“占主导地位（Dominate / 丰度最高）”。
这里需要区分微生物组分析中的两个核心概念：出现频率（Prevalence） 和 相对丰度/主导地位（Abundance/Dominance）。
结合你代码中使用的 decontam 包的 method = "prevalence"（基于频率的方法），具体解释如下：
1. 算法判定的核心是“频率”，而不是“主导地位”
- 出现频率（Prevalence）：指的是一个ASV在多少个样本中出现（即reads数 > 0）。
- 主导地位（Dominance）：指的是一个ASV在某个样本中的相对丰度或reads数是否最高。
decontam 的判定逻辑是：如果一个ASV在大多数阴性对照（比如你的17个NTC中有15个）中都出现了，但在真实患者样本中很少出现（比如200个样本中只有10个有），它就会被判定为污染物。
它不需要在NTC中“占主导”（即不需要是NTC里reads数最多的物种），只要它“普遍存在于NTC中”即可。
2. 为什么不用“占主导”来定义污染物？
- 低丰度污染物（Low-abundance contaminants）：很多来自试剂或实验室环境的污染物DNA含量很低，测序后产生的reads数并不多（在NTC中不占主导），但它们几乎出现在每一个NTC中。如果用“占主导”来筛选，就会漏掉这些低丰度但普遍存在的“隐形”污染物。
- 偶然污染 vs 系统性污染：有时某个真实样本中的优势菌可能偶然污染了某一个NTC，导致它在那一个NTC中占主导，但它不会在大多数NTC中出现。基于频率的方法可以排除这种偶然事件，只抓出系统性污染。
3. 实际数据中的现象
虽然算法是基于“频率”判定的，但在你的实际数据（以及大多数16S测序数据）中，很多典型的污染物（例如 Acinetobacter 不动杆菌、Pseudomonas 假单胞菌、Sphingomonas 鞘氨醇单胞菌等）确实会在NTC中占据极高的比例，甚至在NTC中占主导。
- 这是因为真实样本的DNA浓度通常远高于NTC（NTC中几乎全是污染物DNA），所以在NTC中，这些污染物ASVs的相对丰度会被放大，从而显得“占主导”。
总结
污染ASVs的准确定义是：在阴性对照中出现的概率/频率，显著高于在真实样本中出现的概率的ASVs。
“在阴性对照中占主导”只是很多高浓度污染物的表现现象，而不是 decontam 算法判定污染的根本依据。
是的，您的理解非常准确。
污染 ASVs（Contaminant ASVs） 确实是指那些在阴性对照（Negative Controls, NTCs）中高频出现，甚至在群落中占据主导地位（dominate） 的微生物序列。
为了让您更清晰地理解，可以从以下三个方面来详细说明：
1. 为什么它们会在阴性对照中“占据主导”？
在微生物组测序中，阴性对照（NTC）里面是没有真实生物样本的。因此，NTC 中测出来的所有 DNA 几乎都来源于实验试剂（如提取试剂盒、PCR酶）、实验室环境或操作过程中的背景污染（学术界常称为 “Kitome”）。因为 NTC 本身的生物量极低，这些污染 DNA 就会在 NTC 的测序结果中占据绝对的主导地位（dominate）。
2. decontam 是如何识别它们的？
在您提供的 Rmd 代码中，使用的是 decontam 包的 流行率/出现频率方法（prevalence method）：
contam_prev <- decontam::isContaminant(ps_for_decontam, method = "prevalence", neg = is_ntc, threshold = 0.1)
- 判断逻辑：它并不是单纯看“谁在 NTC 里丰度最高”，而是看 “出现频率（Prevalence）”。
- 如果一个 ASV 在 NTC 样本中出现的频率（比如在 17 个 NTC 中有 15 个都测到了它），显著高于它在真实病人样本中出现的频率，
decontam 就会在统计学上将其判定为“污染 ASV”。
- 当然，由于它们在 NTC 中丰度极高，它们通常也符合“在 NTC 中占据主导”这一直观特征。
3. 剔除它们意味着什么？
一旦某个 ASV 被 decontam 标记为污染（Contaminant），代码中的 prune_taxa(!contam_prev$contaminant, ps_base) 会将这些 ASV 从所有样本（包括真实的病人样本和 NTC）中彻底删除。
- 目的：防止这些背景噪音干扰后续的分析（比如 Alpha 多样性、Beta 多样性、差异丰度分析等），确保您看到的微生物差异是真实的生物学差异，而不是试剂盒批次差异。
总结：
您的理解是对的。污染 ASVs 就是那些主要来源于实验背景、在阴性对照中大量存在（占据主导）且高频出现的序列，必须通过 decontam 等算法将它们从数据集中清洗掉。
Here are a few ways to improve and polish that sentence, depending on where you plan to use it in your report:
Option 1: Formal & Scientific (Best for the main text or Methods section)
“Samples that passed the pre-decontamination depth filter (≥ 12,201 reads) but were subsequently excluded because the removal of contaminant ASVs reduced their final library size below the threshold (e.g., O23092004, O24010402, O23121305).”
Option 2: Concise & Clear (Best for figure legends or table footnotes)
“Samples excluded solely due to post-decontamination read loss. These samples met the initial depth threshold (≥ 12,201 reads) but fell below it after contaminant ASVs were removed (e.g., O23092004, O24010402, O23121305).”
Option 3: Direct Polish (Closest to your original structure, but with fixed punctuation and flow)
“Samples that passed the initial depth check but were dropped solely because decontamination removed enough reads to push their library size below the threshold (Original ≥ 12,201; Post-decontam < 12,201), such as O23092004, O24010402, and O23121305."
Option 4: Optimized for your RMarkdown Table Footnote
If you are putting this in the add_footnote() function from kableExtra, this version is highly readable and fits well:
“Yellow: Excluded solely due to post-decontamination read loss. These samples met the initial depth threshold (≥ 12,201 reads) but fell below it after contaminant ASVs were removed (e.g., O23092004, O24010402, O23121305).”
Key Improvements Made:
- Removed the awkward comma before the parenthesis: Changed
...threshold), (e.g.,... to ...threshold (e.g.,... or ...threshold, such as....
- Clarified the mechanism: Replaced “dropped only because” with “excluded solely due to” or “subsequently excluded because”, which sounds more professional.
- Improved the shorthand: Changed
(Original ≥ 12,201, Post < 12,201) to (Original ≥ 12,201; Post-decontam < 12,201) to make it explicitly clear what “Post” refers to.
Based on your Phyloseq_reanalysis.Rmd, you are actually only filtering once (post-decontamination). The object ps_pruned is assigned directly from ps_decontam (ps_pruned <- ps_decontam), and the depth filter is applied to ps_pruned to create ps_filt. There is no pre-decontamination depth filtering in your current pipeline.
To clearly visualize why the three samples (O23092004, O24010402, O23121305) are dropped, we can generate a table comparing the read counts before (ps_base) and after (ps_decontam) decontamination.
Here is the optimized code to replace your “Preprocessing statistics for each sample” chunk. It uses the phyloseq objects directly (which is more accurate than reading the external TSV) and color-codes the table to distinguish between the two types of exclusions:
library(kableExtra)
library(dplyr)
min_depth <- 12201
# 1. Extract sample sums (total reads) before and after decontamination
sums_before <- sample_sums(ps_base)
sums_after <- sample_sums(ps_decontam)
# 2. Combine into a data frame
depth_df <- data.frame(
SampleID = names(sums_before),
Original_NonChimeric = as.numeric(sums_before),
PostDecontam_NonChimeric = as.numeric(sums_after)
)
# 3. Determine filtering status
depth_df$Filter_Status <- ifelse(
depth_df$Original_NonChimeric < min_depth,
"Excluded (Pre-decontam low depth)",
ifelse(
depth_df$PostDecontam_NonChimeric < min_depth,
"Excluded (Post-decontam low depth)",
"Passed"
)
)
# 4. Identify rows for coloring
rows_pre_excluded <- which(depth_df$Filter_Status == "Excluded (Pre-decontam low depth)")
rows_post_excluded <- which(depth_df$Filter_Status == "Excluded (Post-decontam low depth)")
# 5. Print summary log
cat(sprintf("[INFO] Total samples: %d\n", nrow(depth_df)))
cat(sprintf("[INFO] Excluded due to pre-decontam low depth (< %d): %d\n", min_depth, length(rows_pre_excluded)))
cat(sprintf("[INFO] Excluded ONLY due to decontamination (Original >= %d, Post < %d): %d\n", min_depth, min_depth, length(rows_post_excluded)))
cat(sprintf("[INFO] Passed: %d\n\n", nrow(depth_df) - length(rows_pre_excluded) - length(rows_post_excluded)))
if (length(rows_post_excluded) > 0) {
cat("[NOTE] Samples excluded ONLY after decontamination:\n")
cat(paste(depth_df$SampleID[rows_post_excluded], collapse = ", "), "\n\n")
}
# 6. Render table with color coding
kable(depth_df, caption = sprintf("Sample depth before and after decontamination (Threshold: %d reads)", min_depth)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
row_spec(rows_pre_excluded, background = "#ffcccc", color = "black", bold = TRUE) %>% # Red for pre-decontam
row_spec(rows_post_excluded, background = "#ffeb99", color = "black", bold = TRUE) %>% # Yellow for post-decontam
add_footnote(
sprintf("Red: Failed depth filter before decontamination. Yellow: Passed before, but failed after decontamination (e.g., O23092004, O24010402, O23121305). Threshold: %d", min_depth),
notation = "none"
)
Key Features of this Code:
- Direct Phyloseq Extraction: Uses
sample_sums(ps_base) and sample_sums(ps_decontam) to get the exact read counts used by the pipeline, avoiding any parsing issues with external TSV files.
- Dual Color Coding:
- Red (
#ffcccc): Samples that would have been excluded regardless of decontamination (Original < 12,201).
- Yellow (
#ffeb99): Samples that passed the initial depth check but were dropped only because decontamination removed enough reads to push them below the threshold (Original ≥ 12,201, Post < 12,201). This will highlight your three specific samples.
- Clear Logging: The console output explicitly lists the three samples that were lost purely due to the decontamination step.
The “+3 patients” (11 → 14) are the samples whose library size falls below 12,201 only after the 72 contaminant ASVs are removed — because in the new pipeline the cutoff is applied to clean (post-decontam) sums, whereas in the old pipeline the cutoff was applied to raw, contaminant-inflated sums (and was circular, i.e. set equal to the minimum of the retained set, so it removed nothing by construction).
The 3 additional patients (and the 1 PC)
| Sample |
Group |
Raw depth |
Post-decontam depth |
Why it drops out only in the new pipeline |
| O23092004 |
1 |
12,230 |
< 12,201 |
Marginally above the cutoff pre-decontam; losing a few dozen contaminant reads pushes it below |
| O24010402 |
1 |
12,363 |
< 12,201 |
Same marginal-depth mechanism |
| O23121305 |
5 |
75,260 |
collapses far below 12,201 |
Library was almost pure contaminant (Shannon 0.35, 11 ASVs, i.e. ~97% Burkholderia); after removing contaminant ASVs its clean library collapses |
| UR009768 (the +1 PC) |
PC |
21,031 |
collapses far below 12,201 |
Nearly pure contaminant (Shannon 0.50, 5 ASVs); collapses after decontam |
This is exactly consistent with the heatmap group counts: old samples_keep-based backbone 80/21/69/10/34 (+14 NTC +11 PC) vs new 78/21/69/10/33 (+14 NTC +10 PC) — i.e. −2 in Group 1 (O23092004, O24010402), −1 in Group 5 (O23121305), −1 PC (UR009768).
The 11 patients + 3 NTCs removed in both pipelines are the raw-low-depth set: patients A23060601, A23072501, A23111301, A24040201, O23082401, O23091304, O23100501, O23100502, O23100601, U23071201, U23091101, plus NTC_1, NTC_5, NTC01.
Interpretation
- The old pipeline’s “11 patients + 3 NTCs” were removed before/without decontamination, on raw sums, with a circular threshold.
- The new pipeline’s “14 patients + 3 NTCs + 1 PC” is the methodologically correct order (decontam → depth filter → rarefaction): the threshold now acts on clean library sizes, so samples whose apparent depth was created by contaminant reads are correctly dropped.
- The two Group-1 dropouts are pure threshold-margin cases; the Group-5 sample and the PC are contaminant-dominated libraries whose “depth” was an artifact — dropping them is a QC improvement, not a loss of biological signal.
Verify directly in R
removed <- setdiff(sample_names(ps_pruned), sample_names(ps_filt))
print(sort(removed)) # should print the 18 IDs
sort(sample_sums(ps_pruned)[removed]) # their post-decontam depths, all < 12,201
sort(sample_sums(ps_decontam)[removed]) # same, post-decontam sums
sort(sample_sums(ps_base)[removed]) # their raw sums (shows the 3 marginal/contaminant-dominated cases)
If the printed list matches the 18 IDs above, the table row is confirmed; if any ID differs, substitute the printed IDs into the table (the mechanism — post-decontam sum < 12,201 — is unchanged).
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_cache to 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 neg panel. 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_filt
Two 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 controls line)
- 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 of samples_keep versus manuscript Table 1 ✔
- Rarefaction removed nothing at import stage:
rarefy_even_depth is applied only inside ps_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_cache and 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_keep set 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 U23091101
The 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.