代谢组学归一化与质量控制

适用于统计分析前校正技术变异。推荐先读取 feature_table.csv 和 sample_info.csv，明确样本名、样本类型、实验组、批次和进样顺序。使用前核对 xcms、statTarget、sva、impute 等软件包的实际版本和 API。

处理顺序

1. 数据读取和预处理

# 读取数据
feature_table <- read.csv("feature_table.csv", row.names = 1)
sample_info <- read.csv("sample_info.csv", row.names = 1)

# 确保样本顺序一致
common_samples <- intersect(rownames(feature_table), rownames(sample_info))
feature_table <- feature_table[common_samples, ]
sample_info <- sample_info[common_samples, ]

# 分离QC和生物样本
qc_samples <- rownames(sample_info)[sample_info$sample_type == "QC"]
bio_samples <- rownames(sample_info)[sample_info$sample_type == "biological"]

2. 缺失值计算和特征过滤

# 计算缺失率
feature_missing_rate <- apply(is.na(feature_table), 2, mean)
sample_missing_rate <- apply(is.na(feature_table), 1, mean)

# 按实验组计算缺失率
calculate_group_missing <- function(feature, group_col) {
  tapply(feature, sample_info[[group_col]], function(x) {
    sum(is.na(x)) / length(x)
  })
}

# 保留在至少一个实验组中缺失率不超过20%的特征
keep_features <- sapply(colnames(feature_table), function(feat) {
  group_missing <- calculate_group_missing(feature_table[, feat], "experimental_group")
  any(group_missing <= 0.2, na.rm = TRUE)
})

filtered_table <- feature_table[, keep_features]
cat("保留特征数:", sum(keep_features), "/", ncol(feature_table), "\n")

3. QC-RSC信号漂移校正

library(stats)

# QC-RSC校正函数
qc_rsc_correction <- function(data, sample_info, qc_samples, span = 0.75) {
  corrected_data <- data
  injection_order <- sample_info$injection_order[match(rownames(data), rownames(sample_info))]
  
  for (feat in colnames(data)) {
    qc_values <- data[qc_samples, feat]
    qc_order <- injection_order[qc_samples]
    
    # LOESS拟合QC趋势
    loess_fit <- loess(qc_values ~ qc_order, span = span)
    
    # 预测所有样本的趋势
    all_predicted <- predict(loess_fit, injection_order)
    
    # 校正：除以趋势后乘以QC中位数
    qc_median <- median(qc_values, na.rm = TRUE)
    corrected_data[, feat] <- (data[, feat] / all_predicted) * qc_median
  }
  
  return(corrected_data)
}

# 应用校正（如果有足够的QC样本）
if (length(qc_samples) >= 6) {
  corrected_table <- qc_rsc_correction(filtered_table, sample_info, qc_samples)
} else {
  warning("QC样本不足，跳过QC-RSC校正")
  corrected_table <- filtered_table
}

4. 归一化（TIC或PQN）

# TIC归一化
tic_normalization <- function(data) {
  tic <- rowSums(data, na.rm = TRUE)
  normalized <- data / tic
  return(normalized)
}

# PQN归一化
pqn_normalization <- function(data, reference = NULL) {
  if (is.null(reference)) {
    # 使用中位数样本作为参考
    reference <- apply(data, 1, median, na.rm = TRUE)
  }
  
  # 计算每个样本与参考的比率
  quotients <- data / reference
  
  # 使用中位数比率进行归一化
  normalization_factors <- apply(quotients, 2, median, na.rm = TRUE)
  normalized <- sweep(data, 2, normalization_factors, "/")
  
  return(normalized)
}

# 选择归一化方法
normalized_table <- tic_normalization(corrected_table)
# 或者使用PQN：normalized_table <- pqn_normalization(corrected_table)

5. 缺失值插补

# 半最小值插补
half_min_impute <- function(data) {
  for (feat in colnames(data)) {
    min_value <- min(data[, feat], na.rm = TRUE)
    half_min <- min_value / 2
    data[is.na(data[, feat]), feat] <- half_min
  }
  return(data)
}

# KNN插补
knn_impute <- function(data, k = 5) {
  library(impute)
  # 转换为矩阵
  data_matrix <- as.matrix(data)
  # KNN插补
  imputed_data <- impute::impute.knn(data_matrix, k = k)
  return(imputed_data$data)
}

# 根据数据情况选择插补方法
if (sum(is.na(normalized_table)) / length(normalized_table) < 0.1) {
  imputed_table <- half_min_impute(normalized_table)
} else {
  imputed_table <- knn_impute(normalized_table)
}

6. 数据转换和缩放

# Log2转换
log_transform <- function(data) {
  log2_data <- log2(data + 1)
  return(log2_data)
}

# Pareto缩放
pareto_scaling <- function(data) {
  scaled <- scale(data, center = TRUE, scale = apply(data, 2, sd))
  scaled <- sqrt(scaled)
  return(scaled)
}

# Z-score自动缩放
autoscaling <- function(data) {
  scaled <- scale(data, center = TRUE, scale = apply(data, 2, sd))
  return(scaled)
}

# 应用转换和缩放
log_table <- log_transform(imputed_table)
scaled_table <- pareto_scaling(log_table)  # 或使用 autoscaling(log_table)

7. 批次效应校正

library(sva)

# ComBat批次校正
combat_correction <- function(data, sample_info, batch_col = "batch", biological_factors = NULL) {
  batch <- sample_info[[batch_col]]
  
  # 创建设计矩阵，保留生物学分组信息
  if (!is.null(biological_factors)) {
    design <- model.matrix(~ 0 + factor(sample_info[[biological_factors]]))
    colnames(design) <- levels(factor(sample_info[[biological_factors]]))
  } else {
    design <- NULL
  }
  
  # 应用ComBat
  corrected_data <- ComBat(dat = as.matrix(data), 
                          batch = batch, 
                          mod = design,
                          par.prior = TRUE,
                          prior.plots = FALSE)
  
  return(corrected_data)
}

# 检查批次和生物学分组是否混杂
if ("batch" %in% colnames(sample_info)) {
  # 检查混杂情况
  contingency_table <- table(sample_info$batch, sample_info$experimental_group)
  print(contingency_table)
  
  # 如果不完全混杂，应用ComBat
  final_table <- combat_correction(scaled_table, sample_info, 
                                   batch_col = "batch", 
                                   biological_factors = "experimental_group")
} else {
  final_table <- scaled_table
}

QC 判定

QC指标计算

# 计算QC样本的RSD
calculate_qc_rsd <- function(data, qc_samples) {
  qc_data <- data[qc_samples, ]
  rsd_values <- apply(qc_data, 2, function(x) {
    sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE) * 100
  })
  return(rsd_values)
}

# 计算QC指标
qc_rsd_before <- calculate_qc_rsd(corrected_table, qc_samples)
qc_rsd_after <- calculate_qc_rsd(final_table, qc_samples)

# 统计QC质量
qc_stats <- data.frame(
  metric = c("样本数", "QC数", "特征数", "总体缺失率(%)", 
            "QC RSD中位数", "RSD<30%特征数"),
  before_correction = c(nrow(final_table), length(qc_samples), ncol(final_table),
                       round(mean(is.na(corrected_table)) * 100, 2),
                       round(median(qc_rsd_before), 2), sum(qc_rsd_before < 30)),
  after_correction = c(nrow(final_table), length(qc_samples), ncol(final_table),
                      round(mean(is.na(final_table)) * 100, 2),
                      round(median(qc_rsd_after), 2), sum(qc_rsd_after < 30))
)

print(qc_stats)

QC可视化

# PCA比较
pca_comparison <- function(data_before, data_after, sample_info) {
  pca_before <- prcomp(t(data_before), scale. = TRUE)
  pca_after <- prcomp(t(data_after), scale. = TRUE)
  
  par(mfrow = c(1, 2))
  plot(pca_before$x[,1], pca_before$x[,2], 
       col = as.factor(sample_info$sample_type),
       main = "校正前PCA", xlab = "PC1", ylab = "PC2")
  
  plot(pca_after$x[,1], pca_after$x[,2], 
       col = as.factor(sample_info$sample_type),
       main = "校正后PCA", xlab = "PC1", ylab = "PC2")
  
  par(mfrow = c(1, 1))
}

# RSD分布比较
rsd_distribution <- function(qc_rsd_before, qc_rsd_after) {
  par(mfrow = c(1, 2))
  hist(qc_rsd_before, main = "校正前QC RSD分布", 
       xlab = "RSD (%)", col = "lightblue", breaks = 30)
  abline(v = 30, col = "red", lty = 2)
  
  hist(qc_rsd_after, main = "校正后QC RSD分布", 
       xlab = "RSD (%)", col = "lightgreen", breaks = 30)
  abline(v = 30, col = "red", lty = 2)
  
  par(mfrow = c(1, 1))
}

# 执行可视化
pca_comparison(corrected_table, final_table, sample_info)
rsd_distribution(qc_rsd_before, qc_rsd_after)

保存处理结果

# 保存最终处理的数据
write.csv(final_table, "normalized_feature_table.csv")

# 保存QC统计信息
write.csv(qc_stats, "qc_summary.csv")

# 保存处理记录
processing_log <- list(
  original_features = ncol(feature_table),
  filtered_features = sum(keep_features),
  samples = nrow(feature_table),
  qc_samples = length(qc_samples),
  normalization_method = "TIC",  # 或"PQN"
  imputation_method = "half_min",  # 或"KNN"
  scaling_method = "pareto",  # 或"autoscaling"
  batch_correction = "batch" %in% colnames(sample_info)
)

saveRDS(processing_log, "processing_log.rds")

原有处理步骤说明

分离 QC 和生物学样本，计算每个特征与每个样本的缺失率。
先过滤高缺失特征。默认可保留至少一个实验组内缺失率不超过 20% 的特征；须在报告中说明阈值。
对有足够 QC 注射的数据，优先基于进样顺序使用 LOESS/QC-RSC 校正信号漂移：以 QC 趋势预测全部样本，除以趋势后乘以 QC 中位数。
按研究设计选择 TIC/总和归一化或 PQN。TIC 适用于总信号差异；PQN 对少数大幅变化的特征更稳健。
对保留的缺失值按缺失机制使用 KNN 或半最小值插补；避免在过滤前无区别插补。
对强度做 log2(x + 1) 转换；根据目的选择 Pareto 缩放或 z-score 自动缩放。
当存在处理批次时，在对数转换数据上使用 sva::ComBat()，并在设计矩阵中保留关注的生物学分组。

QC 判定要点

RSD计算标准：sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE) * 100；常用目标为 RSD < 30%，但应按平台和研究预先定义。
校正验证：比较校正前后的 QC RSD 分布、PCA、进样顺序趋势和缺失率，确认技术校正没有消除真实分组差异。
报告指标：样本数、QC 数、特征数、总体缺失率、QC RSD 中位数和 RSD < 30% 的特征数。

具体计算和可视化代码参见上方"QC指标计算"和"QC可视化"部分。

注意事项

不要在无 QC 或 QC 过少时声称完成了 QC-RSC。ComBat 的批次不能与生物学分组完全混杂；若混杂，应报告该限制而非强行校正。所有变换、过滤与插补均需在结果中可复现。

依赖包安装

使用前请确保安装了以下R包：

# 核心包
install.packages(c("impute", "sva"))

# 如果使用MetaboAnalystR
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("MetaboAnalystR")

# 可选包（用于高级功能）
install.packages(c("ggplot2", "pheatmap", "FactoMineR", "factoextra"))

完整工作流程示例

# 完整流程：从原始数据到归一化数据
# 假设已有 feature_table.csv 和 sample_info.csv

library(impute)
library(sva)

# 1. 数据加载和预处理
feature_table <- read.csv("feature_table.csv", row.names = 1)
sample_info <- read.csv("sample_info.csv", row.names = 1)

# 2. 质量控制流程
qc_samples <- rownames(sample_info)[sample_info$sample_type == "QC"]

# 3. 特征过滤
keep_features <- sapply(colnames(feature_table), function(feat) {
  group_missing <- tapply(feature_table[, feat], 
                         sample_info$experimental_group,
                         function(x) sum(is.na(x)) / length(x))
  any(group_missing <= 0.2)
})

filtered_data <- feature_table[, keep_features]

# 4. QC校正（如适用）
if (length(qc_samples) >= 6) {
  corrected_data <- qc_rsc_correction(filtered_data, sample_info, qc_samples)
} else {
  corrected_data <- filtered_data
}

# 5. 归一化
normalized_data <- tic_normalization(corrected_data)

# 6. 缺失值插补
imputed_data <- knn_impute(normalized_data)

# 7. 数据转换和缩放
log_data <- log_transform(imputed_data)
scaled_data <- autoscaling(log_data)

# 8. 批次校正（如适用）
if ("batch" %in% colnames(sample_info)) {
  final_data <- combat_correction(scaled_data, sample_info, 
                                batch_col = "batch",
                                biological_factors = "experimental_group")
} else {
  final_data <- scaled_data
}

# 9. 保存结果
write.csv(final_data, "final_normalized_data.csv")

常见问题和解决方案

1. QC样本数量不足

问题：QC样本少于6个，无法进行QC-RSC校正
解决：跳过QC-RSC步骤，在报告中说明局限性

2. 批次和生物学分组混杂

问题：ComBat无法区分批次效应和生物学差异
解决：检查混杂矩阵，如完全混杂则不进行批次校正

3. 缺失值过多

问题：缺失率超过50%的样本或特征
解决：提高过滤阈值，删除低质量样本或特征

4. 异常值影响

问题：极端值影响归一化效果
解决：考虑使用鲁棒统计方法（如中位数代替均值）

东方天意代谢组学归一化与质量控制技能