774. 多组学数据标准化策略¶

不同组学的数据范围：
  RNA-seq: 0 ~ 1,000,000 (raw counts)
  蛋白组: 10 ~ 10,000,000 (强度值)
  代谢组: 0.001 ~ 100,000 (峰面积)
  甲基化: 0 ~ 1 (beta值)
  16S/宏基因组: 0 ~ 100,000 (reads数)

问题：
  直接合并 → 数值大的组学主导分析结果
  PCA/聚类/机器学习都会被"大数字"牵着走

解决：标准化 → 所有组学数据在同一尺度

标准化 vs 归一化：
  归一化(Normalization)：消除技术偏差（样本间可比）
    → 先做：CPM/TPM/RPKM等
  标准化(Standardization)：统一数值范围（特征间可比）
    → 后做：Z-score/MinMax等
  顺序：先归一化 → 再标准化

# ===== RNA-seq归一化方法对比 =====
import numpy as np  # 导入numpy
import pandas as pd  # 导入pandas

# 读取原始计数矩阵
counts = pd.read_csv("gene_counts.csv", index_col=0)  # 基因×样本

# ===== 方法一：CPM (Counts Per Million) =====
def cpm_normalize(counts):
    """CPM归一化：每百万读段的计数"""
    lib_sizes = counts.sum(axis=0)  # 每个样本的总reads数
    cpm = counts.div(lib_sizes, axis=1) * 1e6  # 除以总数×100万
    return cpm

cpm = cpm_normalize(counts)
print(f"CPM范围: {cpm.min().min():.2f} ~ {cpm.max().max():.2f}")

# ===== 方法二：TPM (Transcripts Per Million) =====
def tpm_normalize(counts, gene_lengths):
    """TPM归一化：考虑基因长度"""
    # Step 1: 除以基因长度(kb) → RPK
    rpk = counts.div(gene_lengths / 1000, axis=0)  # 每千碱基读段
    # Step 2: 除以RPK总和 × 100万 → TPM
    scaling = rpk.sum(axis=0) / 1e6  # 缩放因子
    tpm = rpk.div(scaling, axis=1)  # 除以缩放因子
    return tpm

# ===== 方法三：DESeq2的VST/rlog =====
# R代码：
# library(DESeq2)
# dds <- DESeqDataSetFromMatrix(countData=counts, colData=coldata, design=~condition)
# vsd <- vst(dds, blind=TRUE)    # 方差稳定变换（推荐，快）
# rld <- rlog(dds, blind=TRUE)   # 正则化log变换（慢但更准）
# normalized_matrix <- assay(vsd)

# ===== 方法四：edgeR的TMM归一化 =====
# R代码：
# library(edgeR)
# dge <- DGEList(counts=counts)
# dge <- calcNormFactors(dge, method="TMM")  # TMM归一化因子
# logCPM <- cpm(dge, log=TRUE, prior.count=1)  # log2 CPM

# ===== 组成数据归一化（宏基因组/16S）=====
# 宏基因组数据是组成数据(compositional data)
# 特点：各物种比例之和=100%，不独立
# → 不能直接用RNA-seq的方法

import numpy as np  # 导入numpy
import pandas as pd  # 导入pandas
from scipy import stats  # 导入scipy统计

# ===== CLR变换 (Centered Log-Ratio，推荐) =====
def clr_transform(abundance):
    """CLR变换：组成数据的标准变换"""
    # 处理零值（加伪计数）
    data = abundance.replace(0, 0.5)  # 零替换为0.5

    # 计算几何平均值
    log_data = np.log(data)  # 对数变换
    geo_mean = log_data.mean(axis=0)  # 每列的几何平均(log空间)

    # CLR = log(x) - 几何平均的log
    clr = log_data.subtract(geo_mean, axis=1)  # 中心化
    return clr

# ===== CSS归一化 (Cumulative Sum Scaling) =====
# metagenomeSeq包的方法
# R代码：
# library(metagenomeSeq)
# mgseq <- newMRexperiment(counts)
# mgseq <- cumNorm(mgseq, p=cumNormStat(mgseq))
# css_normalized <- MRcounts(mgseq, norm=TRUE, log=TRUE)

# ===== 稀释(Rarefaction)=====
def rarefy(counts, depth=None):
    """稀释到统一深度"""
    if depth is None:
        depth = counts.sum(axis=0).min()  # 最小样本深度

    rarefied = pd.DataFrame(index=counts.index, columns=counts.columns)
    for sample in counts.columns:  # 遍历每个样本
        pool = []  # 读段池
        for taxon, count in counts[sample].items():  # 展开计数
            pool.extend([taxon] * int(count))

        # 随机抽取depth个读段
        np.random.seed(42)  # 设置随机种子
        if len(pool) >= depth:
            sampled = np.random.choice(pool, size=depth, replace=False)
            for taxon in counts.index:
                rarefied.loc[taxon, sample] = np.sum(sampled == taxon)

    return rarefied.astype(float)

print("归一化方法选择:")
print("  CLR → 统计分析(推荐，保留组成信息)")
print("  CSS → 差异分析(metagenomeSeq)")
print("  稀释 → 多样性分析(alpha/beta多样性)")
print("  TSS → 简单相对丰度(总和缩放)")

# ===== 蛋白组/代谢组归一化 =====

# ===== 中位数归一化 =====
def median_normalize(data):
    """中位数归一化：使所有样本中位数相同"""
    medians = data.median(axis=0)  # 每个样本的中位数
    global_median = medians.median()  # 全局中位数
    normalized = data.div(medians, axis=1) * global_median  # 缩放
    return normalized

# ===== 分位数归一化（跨样本分布一致）=====
def quantile_normalize(data):
    """分位数归一化：使所有样本分布相同"""
    rank_mean = data.stack().groupby(
        data.rank(method="first").stack().astype(int)
    ).mean()  # 按排名取均值

    normalized = data.rank(method="min").stack().astype(int).map(rank_mean).unstack()
    return normalized

# ===== Log2变换（必须在归一化后）=====
def log2_transform(data, pseudocount=1):
    """Log2变换：压缩动态范围"""
    return np.log2(data + pseudocount)  # 加伪计数后取log2

# ===== 跨组学标准化方法 =====
from sklearn.preprocessing import StandardScaler, MinMaxScaler  # 导入缩放器
from sklearn.preprocessing import RobustScaler  # 导入稳健缩放器

# 假设已经完成了各组学的内部归一化
rna_data = pd.read_csv("rna_normalized.csv", index_col=0)   # RNA已归一化
protein_data = pd.read_csv("protein_normalized.csv", index_col=0)  # 蛋白已归一化
metab_data = pd.read_csv("metab_normalized.csv", index_col=0)  # 代谢已归一化

# ===== 方法一：Z-score标准化（推荐）=====
def zscore_standardize(data):
    """Z-score: (x - mean) / std，均值0标准差1"""
    scaler = StandardScaler()  # 创建标准化器
    scaled = pd.DataFrame(
        scaler.fit_transform(data),  # 标准化
        index=data.index, columns=data.columns
    )
    return scaled

rna_z = zscore_standardize(rna_data)  # RNA Z-score
protein_z = zscore_standardize(protein_data)  # 蛋白 Z-score
metab_z = zscore_standardize(metab_data)  # 代谢 Z-score

# ===== 方法二：Min-Max缩放（0-1范围）=====
def minmax_scale(data):
    """Min-Max: (x - min) / (max - min)，范围[0,1]"""
    scaler = MinMaxScaler()  # 创建缩放器
    scaled = pd.DataFrame(
        scaler.fit_transform(data),
        index=data.index, columns=data.columns
    )
    return scaled

# ===== 方法三：稳健标准化（对异常值不敏感）=====
def robust_standardize(data):
    """Robust: (x - median) / IQR"""
    scaler = RobustScaler()  # 基于中位数和IQR
    scaled = pd.DataFrame(
        scaler.fit_transform(data),
        index=data.index, columns=data.columns
    )
    return scaled

# ===== 方法四：Pareto缩放（代谢组常用）=====
def pareto_scale(data):
    """Pareto: (x - mean) / sqrt(std)"""
    mean = data.mean(axis=0)  # 均值
    std = data.std(axis=0)  # 标准差
    pareto = (data - mean) / np.sqrt(std)  # Pareto缩放
    return pareto

# ===== 多组学数据整合 =====

# 合并标准化后的多组学数据
combined = pd.concat([rna_z, protein_z, metab_z], axis=1)  # 列拼接
print(f"整合矩阵: {combined.shape}")  # 样本×(RNA基因+蛋白+代谢物)

# 检查标准化效果
import matplotlib.pyplot as plt  # 导入matplotlib

fig, axes = plt.subplots(1, 3, figsize=(15, 4))  # 3个子图

# 标准化前分布
for i, (name, data) in enumerate([
    ("RNA-seq", rna_data), ("Protein", protein_data), ("Metabolome", metab_data)
]):
    axes[i].hist(data.values.flatten(), bins=50, alpha=0.7)  # 直方图
    axes[i].set_title(f"{name} (标准化前)")
    axes[i].set_xlabel("Value")

plt.tight_layout()
plt.savefig("before_standardization.png", dpi=300)

# 标准化后分布
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for i, (name, data) in enumerate([
    ("RNA-seq", rna_z), ("Protein", protein_z), ("Metabolome", metab_z)
]):
    axes[i].hist(data.values.flatten(), bins=50, alpha=0.7, color="orange")
    axes[i].set_title(f"{name} (Z-score后)")
    axes[i].set_xlabel("Z-score")

plt.tight_layout()
plt.savefig("after_standardization.png", dpi=300)

选择标准化方法的决策树：

你的数据是什么类型？
├── RNA-seq (count data)
│   ├── 差异表达分析 → DESeq2 VST或TMM+logCPM
│   ├── 聚类/PCA → VST或rlog
│   └── 机器学习 → VST → Z-score
│
├── 16S/宏基因组 (组成数据)
│   ├── 统计分析 → CLR变换
│   ├── 多样性分析 → 稀释(Rarefaction)
│   └── 差异分析 → CSS或ALDEx2 CLR
│
├── 蛋白组/代谢组 (强度值)
│   ├── 先做 → 中位数归一化或分位数归一化
│   └── 再做 → Log2 → Z-score或Pareto
│
├── 甲基化 (beta值, 0-1)
│   ├── 差异分析 → M值 = log2(beta/(1-beta))
│   └── 可视化 → beta值直接用
│
└── 多组学整合
    ├── 各组学先内部归一化
    └── 统一Z-score或Robust标准化

# ===== 标准化速查 =====

# RNA-seq
# 差异分析 → DESeq2(raw counts)
# PCA/聚类 → VST: vst(dds)
# 跨样本比较 → TPM

# 宏基因组/16S
# 统计分析 → CLR: log(x) - mean(log(x))
# 多样性 → Rarefaction
# 差异分析 → ALDEx2 CLR

# 蛋白组/代谢组
# 归一化 → 中位数归一化 → Log2
# 标准化 → Z-score或Pareto

# 跨组学整合
# 步骤：内部归一化 → Z-score → 合并
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
data_z = scaler.fit_transform(data)

# 方法对比
# Z-score: (x-mean)/std → 均值0, 标准差1
# Min-Max: (x-min)/(max-min) → 范围[0,1]
# Robust: (x-median)/IQR → 对异常值稳健
# Pareto: (x-mean)/sqrt(std) → 代谢组常用
# CLR: log(x/geometric_mean) → 组成数据必用