673 微生物组数据标准化¶

# 微生物组数据标准化方法合集
library(phyloseq)       # 微生物组数据管理
library(metagenomeSeq)  # CSS标准化
library(edgeR)          # TMM标准化
library(DESeq2)         # DESeq2标准化
library(vegan)          # rarefaction

# ============ 准备数据 ============
# 假设有一个phyloseq对象ps
# ps <- readRDS("phyloseq_object.rds")

# 查看原始数据
cat("样本数:", nsamples(ps), "\n")       # 样本数
cat("物种数:", ntaxa(ps), "\n")          # 物种数
cat("测序深度范围:", range(sample_sums(ps)), "\n")  # 深度范围

# ============ 1. TSS (Total Sum Scaling) ============
# 最简单的方法——除以总reads数得到相对丰度
tss_normalize <- function(ps) {
  otu <- otu_table(ps)                   # 提取OTU表
  if (!taxa_are_rows(ps)) otu <- t(otu)  # 确保行是物种
  tss <- apply(otu, 2, function(x) x / sum(x))  # 每列除以列总和
  otu_table(ps) <- otu_table(tss, taxa_are_rows = TRUE)  # 替换OTU表
  return(ps)                             # 返回标准化后的对象
}
ps_tss <- tss_normalize(ps)              # TSS标准化

# ============ 2. Rarefaction (抽平) ============
# 经典方法——抽样到最小测序深度
min_depth <- min(sample_sums(ps))        # 最小测序深度
cat("抽平深度:", min_depth, "\n")         # 打印深度
ps_rare <- rarefy_even_depth(            # 抽平函数
  ps,
  sample.size = min_depth,               # 抽平到最小深度
  rngseed = 42,                          # 随机种子（可重复）
  replace = FALSE,                       # 不放回抽样
  trimOTUs = TRUE                        # 移除零丰度OTU
)
cat("抽平后样本数:", nsamples(ps_rare), "\n")  # 可能移除了低深度样本
cat("抽平后物种数:", ntaxa(ps_rare), "\n")     # 抽平后物种数

# ============ 3. CSS (Cumulative Sum Scaling) ============
# metagenomeSeq包的方法——基于累积分布的缩放
css_normalize <- function(ps) {
  otu <- as(otu_table(ps), "matrix")     # 转为矩阵
  if (!taxa_are_rows(ps)) otu <- t(otu)  # 确保行是物种

  # 创建MRexperiment对象
  mr <- newMRexperiment(otu)             # metagenomeSeq对象
  p <- cumNormStatFast(mr)               # 计算CSS百分位参数
  mr <- cumNorm(mr, p = p)              # CSS标准化
  css_counts <- MRcounts(mr, norm = TRUE, log = TRUE)  # 提取标准化计数

  otu_table(ps) <- otu_table(css_counts, taxa_are_rows = TRUE)
  return(ps)
}
ps_css <- css_normalize(ps)              # CSS标准化

# ============ 4. TMM (Trimmed Mean of M-values) ============
# edgeR的方法——计算样本间的缩放因子
tmm_normalize <- function(ps) {
  otu <- as(otu_table(ps), "matrix")     # 转为矩阵
  if (!taxa_are_rows(ps)) otu <- t(otu)  # 确保行是物种

  dge <- DGEList(counts = otu)           # 创建DGEList对象
  dge <- calcNormFactors(dge, method = "TMM")  # 计算TMM因子
  tmm_counts <- cpm(dge, log = TRUE)    # 转为log-CPM

  otu_table(ps) <- otu_table(tmm_counts, taxa_are_rows = TRUE)
  return(ps)
}
ps_tmm <- tmm_normalize(ps)             # TMM标准化

# ============ 5. DESeq2归一化 ============
# 基于几何均值比的中位数因子
deseq2_normalize <- function(ps) {
  # 需要非零行
  ps_nz <- prune_taxa(taxa_sums(ps) > 0, ps)  # 移除零丰度物种
  dds <- phyloseq_to_deseq2(ps_nz, ~ 1)       # 转为DESeq2对象
  dds <- estimateSizeFactors(dds, type = "poscounts")  # 估计大小因子
  normalized_counts <- counts(dds, normalized = TRUE)  # 标准化计数

  otu_table(ps_nz) <- otu_table(normalized_counts, taxa_are_rows = TRUE)
  return(ps_nz)
}
ps_deseq <- deseq2_normalize(ps)         # DESeq2标准化

# ============ 6. CLR (Centered Log-Ratio) ============
# 组成数据分析推荐方法
clr_normalize <- function(ps, pseudocount = 0.5) {
  otu <- as(otu_table(ps), "matrix")     # 转为矩阵
  if (!taxa_are_rows(ps)) otu <- t(otu)  # 确保行是物种

  # 添加伪计数处理零值
  otu_pseudo <- otu + pseudocount        # 加伪计数避免log(0)

  # 计算CLR
  gm <- apply(otu_pseudo, 2, function(x) exp(mean(log(x))))  # 几何均值
  clr_values <- log(otu_pseudo / gm)     # CLR变换

  otu_table(ps) <- otu_table(clr_values, taxa_are_rows = TRUE)
  return(ps)
}
ps_clr <- clr_normalize(ps)             # CLR标准化

# Python实现标准化方法选择和比较
import numpy as np   # 数值计算
import pandas as pd  # 数据处理
from scipy import stats  # 统计检验
import matplotlib.pyplot as plt  # 绑图

# 1. TSS标准化
def tss_normalize(count_matrix):
    """TSS标准化：除以列总和得到相对丰度"""
    col_sums = count_matrix.sum(axis=0)  # 每个样本的总reads
    return count_matrix / col_sums       # 除以总reads

# 2. Rarefaction（随机抽样）
def rarefy(count_matrix, depth=None, seed=42):
    """Rarefaction：抽样到指定深度"""
    np.random.seed(seed)                 # 设随机种子
    if depth is None:
        depth = int(count_matrix.sum(axis=0).min())  # 最小深度
    rarefied = pd.DataFrame(            # 创建空矩阵
        0, index=count_matrix.index, columns=count_matrix.columns
    )
    for sample in count_matrix.columns:  # 遍历每个样本
        counts = count_matrix[sample].values  # 该样本的计数
        total = int(counts.sum())        # 总reads数
        if total < depth:                # 跳过深度不够的样本
            continue
        # 构造reads池
        pool = np.repeat(np.arange(len(counts)), counts.astype(int))
        chosen = np.random.choice(pool, size=depth, replace=False)  # 随机抽样
        unique, cnt = np.unique(chosen, return_counts=True)  # 统计
        rarefied.iloc[unique, rarefied.columns.get_loc(sample)] = cnt
    return rarefied

# 3. CSS标准化（Python简化版）
def css_normalize(count_matrix, quantile=0.5):
    """CSS标准化：基于累积和缩放"""
    normalized = count_matrix.copy()     # 复制数据
    for sample in count_matrix.columns:  # 遍历样本
        counts = count_matrix[sample].sort_values()  # 排序计数
        nonzero = counts[counts > 0]     # 非零值
        cumsum = nonzero.cumsum()         # 累积和
        threshold = nonzero.quantile(quantile)  # 分位数阈值
        scaling_factor = counts[counts <= threshold].sum()  # 缩放因子
        if scaling_factor > 0:
            normalized[sample] = count_matrix[sample] / scaling_factor * 1e6
    return normalized

# 4. CLR变换
def clr_transform(count_matrix, pseudocount=0.5):
    """CLR变换：中心对数比"""
    data = count_matrix + pseudocount    # 加伪计数
    log_data = np.log(data)              # 取对数
    gm = log_data.mean(axis=0)           # 几何均值（对数空间的算术均值）
    return log_data - gm                 # 减去几何均值

# 5. 比较不同标准化方法的效果
def compare_normalizations(count_matrix, metadata, group_col="Group"):
    """比较不同标准化方法对PCoA结果的影响"""
    from sklearn.decomposition import PCA  # PCA降维
    from scipy.spatial.distance import pdist, squareform  # 距离计算

    methods = {
        "Raw": count_matrix,                    # 原始数据
        "TSS": tss_normalize(count_matrix),     # TSS
        "CLR": clr_transform(count_matrix),     # CLR
        "CSS": css_normalize(count_matrix),     # CSS
    }

    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    for idx, (name, data) in enumerate(methods.items()):
        ax = axes[idx // 2, idx % 2]     # 子图位置
        # PCA降维
        pca = PCA(n_components=2)        # 2维PCA
        coords = pca.fit_transform(data.T)  # 转置后降维（样本在行）

        # 按组着色
        groups = metadata[group_col]     # 分组信息
        for group in groups.unique():    # 遍历每个组
            mask = groups == group       # 该组的样本
            ax.scatter(coords[mask, 0], coords[mask, 1],
                      label=group, alpha=0.7, s=60)  # 散点图
        ax.set_xlabel(f"PC1 ({pca.explained_variance_ratio_[0]:.1%})")
        ax.set_ylabel(f"PC2 ({pca.explained_variance_ratio_[1]:.1%})")
        ax.set_title(f"{name} 标准化")    # 标题
        ax.legend()                       # 图例

    plt.suptitle("不同标准化方法的PCA比较", fontsize=14)
    plt.tight_layout()
    plt.savefig("normalization_comparison.png", dpi=150)
    print("比较图已保存: normalization_comparison.png")

下游分析是什么？
│
├─ α多样性 → Rarefaction（标准做法）或 不标准化（用原始计数）
│
├─ β多样性/排序
│   ├─ Bray-Curtis → TSS（相对丰度）
│   ├─ UniFrac → Rarefaction 或 不标准化
│   └─ Aitchison → CLR变换（推荐）
│
├─ 差异丰度分析
│   ├─ DESeq2 → DESeq2内部标准化（不需额外处理）
│   ├─ edgeR → TMM标准化
│   ├─ metagenomeSeq → CSS标准化
│   ├─ ANCOM-BC → 内部bias校正
│   └─ ALDEx2 → CLR + MC采样
│
├─ 相关性/网络分析 → CLR变换（避免组成偏差）
│
└─ 机器学习 → TSS或CLR（视模型而定）

# 标准化方法速查
TSS: 简单直观，适合初步探索
Rarefaction: 经典方法，适合多样性分析
CSS: 适合差异分析（metagenomeSeq）
TMM: 适合差异分析（edgeR）
DESeq2: 适合差异分析（自动处理）
CLR: 组成数据分析推荐（相关性/网络）
GMPR: 零值多的稀疏数据

# 什么时候不该用Rarefaction？
- 差异丰度分析（浪费统计功效）
- 样本量少时（每条read都珍贵）
- 机器学习特征提取（需要最大信息量）

# 伪计数策略
CLR变换: 0.5 或 min(非零值)/2
DESeq2: 自动处理（poscounts方法）
log变换: 1（log(x+1)）

# 2024-2025新方法
TaxaNorm: 考虑批次效应的归一化
ANCOM-BC2: 自带bias校正的差异分析
wrench: 基于分组信息的归一化因子估计