630_基因组去污染Decontam¶

一句话概述： 实验室试剂和环境中的微生物DNA会严重污染低生物量样本（如皮肤、空气、血液）的微生物组数据，decontam是R/Bioconductor包通过频率法和流行度法识别污染物，micRoclean（2025年发布）则专门针对宏基因组组装后的MAG去污染。

核心知识点速查表¶

概念	白话解释
Contamination	污染——样本中混入了不是来自研究对象的微生物DNA
Kitome	试剂盒污染——DNA提取试剂自带的微生物DNA
Negative control	阴性对照——不加样本的空白提取（检测试剂污染）
Frequency method	频率法——污染物在高DNA浓度样本中占比低
Prevalence method	流行度法——污染物在阴性对照中更常见
decontam	R包，用统计方法识别和去除污染序列
micRoclean	2025年新工具，用于MAG去污染
DNA concentration	DNA浓度——样本总DNA量（频率法的关键输入）
Low biomass	低生物量——微生物含量很少的样本类型

一、安装¶

# === decontam安装（Bioconductor）===
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install("decontam")  # 安装decontam

# === 依赖包 ===
BiocManager::install("phyloseq")   # 微生物组数据管理
install.packages("ggplot2")        # 可视化
install.packages("gridExtra")      # 图形排列

library(decontam)   # 加载decontam
library(phyloseq)   # 加载phyloseq
library(ggplot2)    # 加载ggplot2

二、数据准备¶

2.1 构建phyloseq对象¶

# === 准备输入数据 ===
# decontam需要：
# 1. ASV/OTU丰度表
# 2. 样本元数据（含DNA浓度和/或阴性对照标记）
# 3. 分类注释（可选但推荐）

# 读取数据
otu <- read.csv("otu_table.csv", row.names = 1)  # OTU表
tax <- read.csv("taxonomy.csv", row.names = 1)    # 分类表
meta <- read.csv("sample_metadata.csv", row.names = 1)  # 元数据

# 元数据必须包含：
# DNA_concentration — DNA浓度（ng/µL），频率法需要
# is_neg — 是否是阴性对照（TRUE/FALSE），流行度法需要

# 构建phyloseq对象
ps <- phyloseq(
    otu_table(as.matrix(otu), taxa_are_rows = TRUE),  # OTU表
    tax_table(as.matrix(tax)),                          # 分类表
    sample_data(meta)                                   # 元数据
)

cat("样本数:", nsamples(ps), "\n")
cat("特征数:", ntaxa(ps), "\n")
cat("阴性对照数:", sum(sample_data(ps)$is_neg), "\n")

2.2 数据质量检查¶

# === 检查DNA浓度分布 ===
df <- data.frame(sample_data(ps))  # 提取元数据

# 画DNA浓度分布
ggplot(df, aes(x = DNA_concentration, fill = is_neg)) +
    geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
    scale_fill_manual(values = c("FALSE" = "steelblue", "TRUE" = "red"),
                      labels = c("Sample", "Negative Control")) +
    labs(x = "DNA Concentration (ng/µL)", 
         y = "Count",
         title = "DNA Concentration Distribution",
         fill = "Sample Type") +
    theme_minimal()

ggsave("dna_concentration_dist.pdf", width = 8, height = 5)

# === 检查阴性对照中的reads数 ===
neg_samples <- sample_names(ps)[sample_data(ps)$is_neg == TRUE]
neg_reads <- sample_sums(ps)[neg_samples]
cat("阴性对照reads数:\n")
print(neg_reads)

# 如果阴性对照reads太少（<100），流行度法可能不可靠
# 如果阴性对照reads太多（>10000），说明污染严重

三、频率法（Frequency Method）¶

3.1 基本用法¶

# === 频率法：利用DNA浓度识别污染物 ===
# 原理：真正的微生物在所有样本中的频率不随DNA浓度变化
#       污染物的频率随DNA浓度增加而降低（因为被稀释了）

contamdf_freq <- isContaminant(
    ps,                            # phyloseq对象
    method = "frequency",          # 使用频率法
    conc = "DNA_concentration",    # DNA浓度列名
    threshold = 0.1                # 阈值（p<0.1视为污染物）
)

# 查看结果
table(contamdf_freq$contaminant)  # TRUE=污染物, FALSE=非污染物
cat("识别的污染物数:", sum(contamdf_freq$contaminant), "\n")

# 阈值选择：
# 0.1 — 默认值，较保守
# 0.5 — 更积极去除污染（可能误伤真实微生物）
# 0.05 — 非常保守（可能漏掉污染物）

3.2 可视化频率模式¶

# === 可视化：污染物vs真实微生物的频率模式 ===
# 选几个疑似污染物看频率-浓度关系

# 找到所有被标记为污染物的ASV
contam_asvs <- rownames(contamdf_freq[contamdf_freq$contaminant == TRUE, ])

# 画频率-浓度图（前6个污染物）
plot_frequency(
    ps,                            # phyloseq对象
    taxa = head(contam_asvs, 6),   # 展示前6个污染物
    conc = "DNA_concentration"     # DNA浓度列
) + 
    labs(title = "Contaminant Frequency vs DNA Concentration")

ggsave("contaminant_frequency.pdf", width = 12, height = 8)

# 解读：
# 污染物：频率随DNA浓度增加而下降（负相关）
# 真实微生物：频率不随DNA浓度变化（无相关）

四、流行度法（Prevalence Method）¶

4.1 基本用法¶

# === 流行度法：利用阴性对照识别污染物 ===
# 原理：污染物在阴性对照中出现的频率>在真实样本中出现的频率

contamdf_prev <- isContaminant(
    ps,                            # phyloseq对象
    method = "prevalence",         # 使用流行度法
    neg = "is_neg",                # 阴性对照标记列名
    threshold = 0.1                # 阈值
)

table(contamdf_prev$contaminant)  # 查看结果
cat("识别的污染物数:", sum(contamdf_prev$contaminant), "\n")

4.2 可视化流行度模式¶

# === 可视化流行度比较 ===
# 计算每个ASV在阴性对照和真实样本中的流行度
ps_neg <- prune_samples(sample_data(ps)$is_neg == TRUE, ps)   # 阴性对照
ps_pos <- prune_samples(sample_data(ps)$is_neg == FALSE, ps)  # 真实样本

# 流行度数据框
prev_df <- data.frame(
    prev_neg = apply(otu_table(ps_neg), 1, function(x) sum(x > 0)),
    prev_pos = apply(otu_table(ps_pos), 1, function(x) sum(x > 0)),
    contaminant = contamdf_prev$contaminant
)
prev_df$prev_neg_frac <- prev_df$prev_neg / nsamples(ps_neg)
prev_df$prev_pos_frac <- prev_df$prev_pos / nsamples(ps_pos)

# 画散点图
ggplot(prev_df, aes(x = prev_neg_frac, y = prev_pos_frac, 
                     color = contaminant)) +
    geom_point(alpha = 0.6, size = 2) +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    scale_color_manual(values = c("FALSE" = "grey60", "TRUE" = "red"),
                       labels = c("Non-contaminant", "Contaminant")) +
    labs(x = "Prevalence in Negative Controls",
         y = "Prevalence in True Samples",
         title = "Prevalence Comparison") +
    theme_minimal()

ggsave("prevalence_comparison.pdf", width = 8, height = 7)

五、联合方法¶

5.1 两种方法联合使用¶

# === 联合频率法和流行度法 ===
contamdf_combined <- isContaminant(
    ps,
    method = "combined",           # 联合两种方法
    conc = "DNA_concentration",    # DNA浓度
    neg = "is_neg",                # 阴性对照标记
    threshold = 0.1                # 阈值
)

table(contamdf_combined$contaminant)  # 查看结果
cat("联合方法识别的污染物数:", sum(contamdf_combined$contaminant), "\n")

# 比较三种方法的结果
cat("\n=== 方法比较 ===\n")
cat("频率法:", sum(contamdf_freq$contaminant), "个\n")
cat("流行度法:", sum(contamdf_prev$contaminant), "个\n")
cat("联合法:", sum(contamdf_combined$contaminant), "个\n")

5.2 去除污染物并保存¶

# === 从数据中去除污染物 ===
# 选择使用哪种方法的结果（推荐联合法或根据数据情况选择）
contaminants <- contamdf_combined$contaminant  # 污染物标记

# 去除污染物
ps_clean <- prune_taxa(!contaminants, ps)  # 保留非污染物

# 同时去除阴性对照样本
ps_clean <- prune_samples(
    sample_data(ps_clean)$is_neg == FALSE,  # 只保留真实样本
    ps_clean
)

cat("清理前: ", ntaxa(ps), "个特征,", nsamples(ps), "个样本\n")
cat("清理后: ", ntaxa(ps_clean), "个特征,", nsamples(ps_clean), "个样本\n")

# 保存清理后的数据
saveRDS(ps_clean, "phyloseq_decontaminated.rds")

# === 查看被去除的污染物是什么 ===
contam_taxa <- tax_table(ps)[contaminants, ]  # 污染物的分类信息
cat("\n=== 被去除的污染物 ===\n")
print(as.data.frame(contam_taxa))

六、常见污染微生物参考¶

# === 已知的常见试剂污染物 ===
# 参考：Salter et al., BMC Biology 2014

known_contaminants <- c(
    # 常见试剂污染属
    "Bradyrhizobium",      # DNA提取试剂盒常见
    "Ralstonia",           # 水和试剂中常见
    "Methylobacterium",    # 实验室环境常见
    "Sphingomonas",        # 水和试剂中常见
    "Pseudomonas",         # 无处不在的环境菌（部分是污染）
    "Acinetobacter",       # 皮肤和环境
    "Stenotrophomonas",    # 水系统
    "Herbaspirillum",      # DNA提取试剂盒
    "Burkholderia",        # 试剂和环境
    "Propionibacterium",   # 皮肤（现更名Cutibacterium）
    "Corynebacterium"      # 皮肤
)

# 检查数据中是否存在已知污染物
if (!is.null(tax_table(ps))) {
    genus_col <- "Genus"  # 属级分类列名
    taxa_df <- as.data.frame(tax_table(ps))
    found_contaminants <- taxa_df[taxa_df[[genus_col]] %in% known_contaminants, ]
    cat("数据中发现的已知污染属:\n")
    print(unique(found_contaminants[[genus_col]]))
}

七、micRoclean MAG去污染¶

# === micRoclean：宏基因组MAG去污染（2025年新工具）===
pip install microclean  # 安装micRoclean

# micRoclean检测MAG（Metagenome-Assembled Genome）中的污染contig
# 适用于宏基因组binning后的质量控制

# 基本用法
microclean \
    --input mag_bins/ \            # MAG bin目录
    --output cleaned_bins/ \       # 清理后的输出
    --threads 16 \                 # 16线程
    --min-completeness 50 \        # 最低完整性50%
    --max-contamination 5          # 最高污染度5%

# micRoclean使用机器学习方法：
# 基于tetranucleotide频率、GC含量、覆盖度等特征
# 识别每个contig是否属于该MAG
# 去除与MAG主体不一致的污染contig

八、最佳实践流程¶

# === 完整的去污染最佳实践 ===

# 1. 实验设计阶段
#    - 每批次至少包含2-3个阴性对照（空白提取）
#    - 记录每个样本的DNA浓度（Qubit测量）
#    - 阳性对照（mock community）验证流程正确性

# 2. 数据分析阶段
# Step 1: 数据导入和初步质控
ps <- readRDS("raw_phyloseq.rds")  # 读取原始数据

# Step 2: 检查阴性对照
neg_reads <- sample_sums(ps)[sample_data(ps)$is_neg == TRUE]
cat("阴性对照reads:\n")
print(neg_reads)  # 应远低于真实样本

# Step 3: 选择去污染方法
# 有DNA浓度 → 频率法或联合法
# 只有阴性对照 → 流行度法
# 两者都有 → 联合法（推荐）

# Step 4: 运行decontam
contamdf <- isContaminant(ps, method = "combined",
                          conc = "DNA_concentration",
                          neg = "is_neg",
                          threshold = 0.1)

# Step 5: 检查结果合理性
# 查看被标记的是否包含已知污染物
# 确认没有误伤重要的生物学信号

# Step 6: 去除污染物
ps_clean <- prune_taxa(!contamdf$contaminant, ps)
ps_clean <- prune_samples(sample_data(ps_clean)$is_neg == FALSE, ps_clean)

# Step 7: 在方法部分报告去污染步骤
cat("\n=== 报告内容 ===\n")
cat("去污染方法: decontam (combined method)\n")
cat("阈值: 0.1\n")
cat("去除污染ASV数:", sum(contamdf$contaminant), "\n")
cat("保留ASV数:", ntaxa(ps_clean), "\n")

九、常见报错与解决¶

报错信息	原因	解决方案
`No concentration data`	元数据中没有DNA浓度列	添加DNA浓度数据或改用流行度法
`No negative controls`	没有标记阴性对照	在元数据中添加is_neg列
`All taxa are contaminants`	阈值太宽松	降低阈值（如0.05）
`No contaminants found`	阈值太严格或确实没污染	提高阈值（如0.5）或检查数据
`Negative control has 0 reads`	阴性对照没有reads	可能测序深度不够或确实没污染
`phyloseq object error`	数据格式不匹配	检查OTU表、分类表、元数据的行列对应

十、面试高频题¶

Q1：微生物组数据为什么需要去污染？¶

答： DNA提取试剂、PCR试剂、实验室环境中都含有微量微生物DNA。对于高生物量样本（如粪便），这些污染被淹没在大量目标DNA中影响很小。但对于低生物量样本（如血液、CSF、皮肤、空气），污染DNA可能占据大部分甚至全部检测到的微生物信号。经典例子：Salter等人2014年发现仅DNA提取试剂盒就贡献了数十个假阳性菌属。不去污染会导致：(1) 假阳性——报告不存在的微生物；(2) 丰度偏差——真实微生物的相对丰度被稀释；(3) 研究结论错误——在"无菌"部位发现的微生物可能全是污染。

Q2：decontam的频率法和流行度法分别适用于什么场景？¶

答：频率法需要DNA浓度数据——原理是污染物来自固定量的试剂，所以在高浓度样本中被稀释、频率低，在低浓度样本中被"放大"、频率高。适合：有精确DNA浓度测量（Qubit）的实验。流行度法需要阴性对照——原理是污染物在阴性对照中比在真实样本中出现得更频繁。适合：有足够阴性对照（每批次≥2个）的实验。如果两者都有，用联合法最可靠。注意：频率法对DNA浓度变异小的样本效果差；流行度法对阴性对照数量太少时不可靠。

Q3：如何判断decontam的阈值设多少合适？¶

答：默认阈值0.1是个合理的起点，但应该根据数据特征调整：(1) 画频率-浓度图看趋势是否合理——被标记的ASV是否确实呈现负相关趋势；(2) 检查被标记的是否包含已知污染物（如Bradyrhizobium、Ralstonia）——如果已知污染物没被标记，阈值可能太严格；(3) 检查是否误伤了生物学上重要的微生物——如果目标微生物被标记，阈值可能太宽松；(4) 用不同阈值（0.05、0.1、0.2、0.5）做敏感性分析，看核心结论是否稳定。

Q4：除了decontam，还有哪些去污染策略？¶

答： (1) 实验层面——最有效的方法是从源头减少污染：使用经过DNA去除处理的试剂、在超净工作台操作、设置多个阴性对照；(2) 简单过滤——去除在阴性对照中出现的所有ASV（粗暴但有效）；(3) 比例过滤——去除在阴性对照中丰度>在样本中丰度X%的ASV；(4) SourceTracker——用贝叶斯方法估计每个样本中污染源的贡献比例；(5) micRoclean（2025年新工具）——专门针对MAG的contig级去污染，用机器学习方法。最佳实践是多种方法结合：实验控制+统计去污染+结果敏感性分析。

参考资料：decontam: Davis et al., Microbiome 2018 | Salter et al., BMC Biology 2014 | Karstens et al., mSystems 2019 | micRoclean: 2025 | 实验室污染综述: Eisenhofer et al., Trends Microbiol 2019