195_微生物组绝对定量¶

# 组成性偏差示例
# 样本A：物种X=1000, 物种Y=1000, 物种Z=1000
#   相对丰度：X=33%, Y=33%, Z=33%
# 样本B：物种X=1000, 物种Y=1000, 物种Z=5000
#   相对丰度：X=14%, Y=14%, Z=71%
# 
# 相对丰度分析会认为X和Y在B中"下降"了，但实际数量没变！
# 绝对定量能正确反映：X和Y不变，Z增加了。

# Spike-in实验设计
# 1. 选择Spike-in标准品
#    - 合成DNA（如ZymoBIOMICS Spike-in Control）
#    - 已知浓度的非环境来源细菌（如Imtechella halotolerans）
#    - Sequins（合成的全长16S序列）

# 2. 在DNA提取前加入已知量的标准品
# 例如：加入10^6 copies的合成16S到1g土壤样本中

# 3. 正常进行16S扩增子或宏基因组测序

# 4. 数据分析：从测序结果中识别spike-in reads

import pandas as pd
import numpy as np

# 分析spike-in数据
def absolute_quantification(otu_table, spike_in_id, spike_in_copies):
    """
    将相对丰度转化为绝对定量

    Parameters:
    - otu_table: 样本×OTU的count矩阵
    - spike_in_id: spike-in OTU的ID
    - spike_in_copies: 加入的spike-in拷贝数
    """
    # 计算相对丰度
    total_reads = otu_table.sum(axis=1)
    spike_reads = otu_table[spike_in_id]
    spike_fraction = spike_reads / total_reads

    # 计算每个样本的总绝对量
    # spike_copies / spike_fraction = total_copies
    total_absolute = spike_in_copies / spike_fraction

    # 计算每个OTU的绝对丰度
    relative = otu_table.div(total_reads, axis=0)
    absolute = relative.multiply(total_absolute, axis=0)

    # 移除spike-in OTU
    absolute = absolute.drop(columns=spike_in_id)

    return absolute, total_absolute

# 使用示例
otu_table = pd.read_csv("otu_table_with_spikein.csv", index_col=0)
absolute, total = absolute_quantification(
    otu_table, 
    spike_in_id="SyntheticStandard", 
    spike_in_copies=1e6
)

print("各样本总绝对菌量:")
print(total)
print("\nTop OTU的绝对丰度:")
print(absolute.mean().sort_values(ascending=False).head(10))

# qPCR数据与相对丰度整合
library(ggplot2)

# qPCR结果：每个样本的16S拷贝数
qpcr_data <- data.frame(
  Sample = paste0("S", 1:10),
  copies_per_g = c(1e8, 2e8, 5e7, 3e8, 1e9, 2e7, 8e8, 4e8, 6e8, 1e8)
)

# 16S相对丰度数据
relative_abundance <- read.csv("relative_abundance.csv", row.names = 1)

# 计算绝对丰度
absolute_abundance <- sweep(relative_abundance, 1, qpcr_data$copies_per_g, "*")

# 对比相对vs绝对丰度
species_of_interest <- "Bacteroides_fragilis"

comparison <- data.frame(
  Sample = qpcr_data$Sample,
  Relative = relative_abundance[, species_of_interest],
  Absolute = absolute_abundance[, species_of_interest],
  Group = rep(c("Control", "Treatment"), each = 5)
)

# 相对丰度可能显示"下降"
p1 <- ggplot(comparison, aes(x = Group, y = Relative, fill = Group)) +
  geom_boxplot() + ggtitle("Relative Abundance") + theme_minimal()

# 绝对丰度显示"不变"
p2 <- ggplot(comparison, aes(x = Group, y = Absolute, fill = Group)) +
  geom_boxplot() + ggtitle("Absolute Abundance") + theme_minimal()

library(patchwork)
p1 | p2

# 流式细胞术数据整合
# 使用SYBR Green染色后流式计数

flow_data <- data.frame(
  Sample = paste0("S", 1:10),
  cells_per_ml = c(5e6, 8e6, 3e6, 1e7, 2e7, 4e6, 1.5e7, 9e6, 7e6, 6e6)
)

# 同样的计算方法
absolute_flow <- sweep(relative_abundance, 1, flow_data$cells_per_ml, "*")

# 比较不同定量方法的一致性
comparison_methods <- data.frame(
  qPCR_total = qpcr_data$copies_per_g,
  Flow_total = flow_data$cells_per_ml
)
cor(comparison_methods$qPCR_total, comparison_methods$Flow_total, method = "spearman")

# qPCR分析（使用通用16S引物）
# 引物：515F/806R 或 341F/785R
# 标准曲线：已知浓度的16S质粒梯度稀释

# ddPCR绝对定量（更准确，不需标准曲线）
# Bio-Rad QX200系统

# ZymoBIOMICS Spike-in标准品添加
# 在DNA提取前加入已知量到样本中

# 绝对定量计算
absolute <- sweep(relative, 1, total_counts, "*")

# 差异分析（使用绝对丰度）
library(Maaslin2)
Maaslin2(input_data = absolute_abundance, input_metadata = metadata,
         output = "absolute_diff", fixed_effects = "group",
         normalization = "NONE", transform = "LOG")