677 微生物组组成数据分析（CoDA）¶

# 示例：3个物种的相对丰度
真实绝对丰度: A=100, B=100, C=100 → 相对丰度: 33%, 33%, 33%
A翻倍后:      A=200, B=100, C=100 → 相对丰度: 50%, 25%, 25%
# B和C的绝对量没变，但相对丰度"下降"了！→ 假负相关

# 组成数据的对数比变换
library(compositions)  # 组成数据分析包
library(ALDEx2)        # 基于CLR的差异分析

# ============ 准备数据 ============
# 假设otu是物种丰度矩阵（行=物种，列=样本）
# 添加伪计数处理零值
pseudocount <- 0.5                       # 常用伪计数
otu_pseudo <- otu + pseudocount          # 添加伪计数

# ============ 1. CLR变换（最常用） ============
# CLR(x_i) = ln(x_i / g(x))
# g(x) = 几何均值
clr_transform <- function(x) {
  gm <- exp(mean(log(x)))               # 几何均值
  log(x / gm)                           # 中心对数比
}

# 对每个样本做CLR变换
otu_clr <- apply(otu_pseudo, 2, clr_transform)  # 按列（样本）变换
cat("CLR变换后维度:", dim(otu_clr), "\n")

# 使用compositions包
otu_comp <- acomp(t(otu_pseudo))         # 创建组成数据对象
otu_clr2 <- clr(otu_comp)               # CLR变换

# CLR的特点：
# - 维度不变（p个物种→p个CLR值）
# - 不需要选择参考物种
# - CLR值之和=0（每个样本内）
# - 可直接用于PCA、相关分析等

# ============ 2. ALR变换 ============
# ALR(x_i) = ln(x_i / x_ref)
# 需要选择一个参考物种
alr_transform <- function(x, ref_idx) {
  log(x[-ref_idx] / x[ref_idx])         # 相对于参考的对数比
}

# 选择最稳定的物种作为参考
# 方法：选择变异系数最小的物种
cv <- apply(otu_pseudo, 1, function(x) sd(x) / mean(x))  # 变异系数
ref_species <- which.min(cv)             # CV最小的物种
cat("参考物种:", rownames(otu)[ref_species], "\n")

# ALR变换
otu_alr <- apply(otu_pseudo, 2, alr_transform, ref_idx = ref_species)

# ALR的特点：
# - 维度减少1（p个物种→p-1个ALR值）
# - 需要选择参考物种（结果受选择影响）
# - 适合回归分析

# ============ 3. ILR变换（等距对数比） ============
# ILR使用正交基，统计上最正确
otu_ilr <- ilr(acomp(t(otu_pseudo)))     # ILR变换

# ILR的特点：
# - 维度减少1（p→p-1）
# - 不需要选择参考
# - 保持欧几里得距离（Aitchison距离）
# - 解释性差（坐标是正交组合）
# - 适合多元统计检验

# ============ 4. 使用ALDEx2做差异分析 ============
# ALDEx2基于CLR + 蒙特卡洛采样
# 安装: BiocManager::install("ALDEx2")
library(ALDEx2)

# ALDEx2分析
aldex_result <- aldex(
  reads = otu,                           # 原始计数（不加伪计数！）
  conditions = meta$Group,               # 分组
  mc.samples = 128,                      # MC采样次数
  test = "t",                            # t检验
  effect = TRUE,                         # 计算效应量
  include.sample.summary = FALSE,        # 不输出样本摘要
  denom = "all"                          # 使用所有物种作为参考
)

# ALDEx2结果解读
head(aldex_result[order(aldex_result$we.eBH), ])  # 按FDR排序
# 关键列：
# - rab.all: 所有样本的CLR中位数
# - diff.btw: 组间CLR差异
# - diff.win: 组内CLR差异
# - effect: 效应量（diff.btw / diff.win的中位数）
# - we.ep: Welch t检验p值
# - we.eBH: BH校正后的p值（FDR）

# 筛选差异物种
sig_species <- aldex_result[aldex_result$we.eBH < 0.05 &
                             abs(aldex_result$effect) > 1, ]
cat("\n显著差异物种（FDR<0.05, |effect|>1）:\n")
print(sig_species[, c("diff.btw", "effect", "we.eBH")])

# ALDEx2效应量图
aldex.plot(aldex_result, type = "MA",    # MA图
           test = "welch",               # Welch检验
           cutoff.pval = 0.05,           # p值阈值
           cutoff.effect = 1)            # 效应量阈值

# 基于CLR的Aitchison PCA
library(compositions)  # 组成数据包
library(ggplot2)       # 绑图

# 1. Aitchison距离
# = CLR变换后的欧几里得距离
# = 组成数据的"真实距离"
aitchison_dist <- dist(otu_clr2)         # CLR后的欧几里得距离

# 等价于：
# library(robCompositions)
# aDist(acomp(t(otu_pseudo)))

# 2. Aitchison PCA（=对CLR数据做PCA）
pca_result <- prcomp(otu_clr2, scale. = FALSE)  # PCA
summary(pca_result)                      # 方差解释比例

# 可视化
pca_df <- data.frame(
  PC1 = pca_result$x[, 1],              # 第1主成分
  PC2 = pca_result$x[, 2],              # 第2主成分
  Group = meta$Group                     # 分组
)
var_explained <- summary(pca_result)$importance[2, 1:2] * 100

ggplot(pca_df, aes(x = PC1, y = PC2, color = Group)) +
  geom_point(size = 3, alpha = 0.7) +   # 散点
  stat_ellipse(level = 0.95) +           # 95%置信椭圆
  labs(x = paste0("PC1 (", round(var_explained[1], 1), "%)"),
       y = paste0("PC2 (", round(var_explained[2], 1), "%)"),
       title = "Aitchison PCA（CLR变换后）") +
  theme_minimal()                        # 简洁主题
ggsave("aitchison_pca.png", width = 8, height = 6, dpi = 150)

# 3. PERMANOVA with Aitchison距离
library(vegan)
adonis2(aitchison_dist ~ Group + Age, data = meta, permutations = 999)

# Python组成数据分析
import numpy as np   # 数值计算
import pandas as pd  # 数据处理
from scipy.spatial.distance import pdist, squareform  # 距离计算
import matplotlib.pyplot as plt  # 绑图

# 1. CLR变换
def clr_transform(counts, pseudocount=0.5):
    """中心对数比变换"""
    data = counts + pseudocount          # 加伪计数
    log_data = np.log(data)              # 取对数
    gm = log_data.mean(axis=0)           # 几何均值（列方向）
    clr = log_data - gm                  # 减去几何均值
    return clr

# 2. ALR变换
def alr_transform(counts, ref_idx=-1, pseudocount=0.5):
    """加性对数比变换"""
    data = counts + pseudocount          # 加伪计数
    ref = data.iloc[ref_idx]             # 参考物种
    alr = np.log(data.drop(data.index[ref_idx]) / ref)  # 对数比
    return alr

# 3. ILR变换
def ilr_transform(counts, pseudocount=0.5):
    """等距对数比变换"""
    data = (counts + pseudocount).values  # 加伪计数
    n_components = data.shape[0]         # 物种数
    # 构造Helmert矩阵作为正交基
    basis = np.zeros((n_components - 1, n_components))
    for i in range(n_components - 1):
        basis[i, :i+1] = 1.0 / (i + 1)
        basis[i, i+1] = -(i + 1.0) / (i + 1)
        basis[i] *= np.sqrt((i + 1.0) / (i + 2.0))
    # ILR = Helmert矩阵 × log(数据)
    ilr = basis @ np.log(data)           # 矩阵乘法
    return pd.DataFrame(ilr, columns=counts.columns)

# 4. Aitchison距离
def aitchison_distance(counts, pseudocount=0.5):
    """计算Aitchison距离矩阵"""
    clr = clr_transform(counts, pseudocount)  # CLR变换
    dist_matrix = squareform(pdist(clr.T, metric='euclidean'))  # 欧几里得
    return pd.DataFrame(dist_matrix,
                        index=counts.columns,
                        columns=counts.columns)

# 5. 组成偏差检测
def detect_compositional_bias(counts, method="correlation"):
    """检测组成偏差——用CLR前后的相关性比较"""
    # 直接计算相关性（受组成偏差影响）
    raw_corr = np.corrcoef(counts.values)  # 原始相关矩阵
    # CLR后计算相关性（消除组成偏差）
    clr = clr_transform(counts)
    clr_corr = np.corrcoef(clr.values)   # CLR相关矩阵

    # 比较
    n = len(counts)
    raw_neg = (raw_corr[np.triu_indices(n, k=1)] < 0).mean()
    clr_neg = (clr_corr[np.triu_indices(n, k=1)] < 0).mean()

    print(f"原始数据负相关比例: {raw_neg:.1%}")
    print(f"CLR后负相关比例: {clr_neg:.1%}")
    if raw_neg > clr_neg + 0.1:
        print("→ 检测到组成偏差！原始数据的负相关部分是假的")
    else:
        print("→ 组成偏差不严重")

    return raw_corr, clr_corr

# 6. 可视化对比
def plot_composition_effect(counts, output="composition_bias.png"):
    """可视化组成偏差效应"""
    clr = clr_transform(counts)

    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    # 原始相对丰度的相关矩阵
    tss = counts / counts.sum(axis=0)    # TSS标准化
    raw_corr = tss.T.corr()              # 相关矩阵
    im1 = axes[0].imshow(raw_corr, cmap='RdBu_r', vmin=-1, vmax=1)
    axes[0].set_title("TSS相对丰度相关矩阵\n(含组成偏差)")

    # CLR变换后的相关矩阵
    clr_corr = clr.T.corr()              # 相关矩阵
    im2 = axes[1].imshow(clr_corr, cmap='RdBu_r', vmin=-1, vmax=1)
    axes[1].set_title("CLR变换后相关矩阵\n(消除组成偏差)")

    fig.colorbar(im2, ax=axes, label="相关系数")
    plt.tight_layout()
    plt.savefig(output, dpi=150)
    print(f"对比图已保存: {output}")

# CoDA变换选择决策树
你要做什么？
│
├─ 差异丰度分析 → ALDEx2 (CLR + MC) 或 ANCOM-BC
│
├─ PCA/排序分析 → CLR变换 → Aitchison PCA
│
├─ 距离矩阵/PERMANOVA → Aitchison距离（CLR+欧几里得）
│
├─ 相关性/网络分析 → CLR变换 → 标准相关分析
│   (不要对TSS数据直接算Pearson相关！)
│
├─ 回归分析 → ALR变换（减少一个维度）
│
└─ 多元统计检验 → ILR变换（统计上最正确）

# 对数比变换对比
CLR: 维度不变(p→p), 不需参考, 最常用
ALR: 维度减1(p→p-1), 需选参考, 适合回归
ILR: 维度减1(p→p-1), 不需参考, 统计最正确

# 伪计数策略
经典: 0.5（最简单）
Bayesian: 每个成分加1/物种数
CZM: Compositional Zero Multiplier
ALDEx2: 蒙特卡洛采样（推荐）

# Aitchison距离 = CLR后的欧几里得距离
# 是组成数据的"正确"距离度量

# 常见错误
× 对TSS数据直接用Pearson相关 → 假负相关
× 对TSS数据做PCA → 不在正确空间
× 忽略零值直接取对数 → Inf值
✓ CLR + PCA = Aitchison PCA
✓ CLR + 欧几里得 = Aitchison距离
✓ ALDEx2处理差异分析