[批次效应评估] 泛癌多样本来源的数据

针对不同样本来源的泛癌数据，pan-cancer怎么做

发现各组有各组的风格，比如张泽民组更偏向R seurat流程，pan-cancer文章普遍用的是harmony；没有明确说去批次是按样本还是数据集；评估批次效应的方法也没有提到

而王凌华组23年naruew medicine 这一篇就说了他们比较了rPCA和harmony方法，使用轮廓系数评估

总之做法不同，有理即可

评估批次矫正效果：Silhouette Score

王凌华组这一篇Pan-cancer T cell atlas links a cellular stress response state to immunotherapy resistance，用Silhouette Score（轮廓系数）评估harmony和rPCA两种方法的批次矫正效果 （最后他们选了Seurat rPCA）

核心思想是：如果批次效应很强，同一个批次的细胞会抱团（即使细胞类型不同），那么Silhouette Score会很高

如果矫正得好，相同批次的细胞应该分散开，Silhouette Score会降低

对下面公式的解释

对于每个细胞i：

a(i) = 这个细胞和同一批次其他细胞的平均距离

• 批次矫正好时：a(i) 大（同批次分散） ；同批次细胞聚集不好，如 T细胞、B细胞、髓系细胞都因为”来自批次1”而挤在一起

b(i) = 这个细胞到最近的另一个批次的平均距离

• 批次矫正好时：b(i) 小（批次间混合好）

Silhouette Score = (b - a) / max(a, b)

• 范围：-1 到 +1
• **接近+1**：细胞和自己批次很像，和其他批次很不像 → 批次效应强
• **接近0或负数**：细胞和自己批次不那么像，批次混合好 → 批次矫正好

library(Seurat)
library(cluster)
library(ggplot2)

# =============================================================================
# 计算Silhouette Score评估批次效应
# =============================================================================

# 假设你已经有整合好的Seurat对象
# seu_integrated: 批次矫正后的对象
# seu_raw: 批次矫正前的对象（用PCA降维）

# 步骤1: 提取降维后的坐标
# 批次矫正前（PCA空间）
pca_coords <- Embeddings(seu_raw, reduction = "pca")[, 1:30]

# 批次矫正后（Harmony或rPCA空间）
harmony_coords <- Embeddings(seu_integrated, reduction = "harmony")[, 1:30]
# 或者
# rpca_coords <- Embeddings(seu_integrated, reduction = "integrated.dr")[, 1:30]

# 步骤2: 计算Silhouette Score
# 批次信息
batch_labels <- seu_raw$dataset  # 或 seu_raw$sample_id

# 计算距离矩阵（使用欧氏距离）
dist_before <- dist(pca_coords)
dist_after <- dist(harmony_coords)

# 计算Silhouette Score
sil_before <- silhouette(as.numeric(factor(batch_labels)), dist_before)
sil_after <- silhouette(as.numeric(factor(batch_labels)), dist_after)

# 步骤3: 可视化对比
# 提取分数
sil_scores_before <- data.frame(
  score = sil_before[, "sil_width"],
  batch = batch_labels,
  stage = "Before correction"
)

sil_scores_after <- data.frame(
  score = sil_after[, "sil_width"],
  batch = batch_labels,
  stage = "After correction"
)

sil_scores <- rbind(sil_scores_before, sil_scores_after)

# 箱线图对比
ggplot(sil_scores, aes(x = stage, y = score, fill = stage)) +
  geom_boxplot() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Batch Effect Evaluation",
    subtitle = "Lower score = Better batch mixing",
    y = "Silhouette Score",
    x = ""
  ) +
  theme_bw() +
  theme(legend.position = "none")

# 统计摘要
cat("=== Silhouette Score 统计 ===\n")
cat("矫正前平均分数:", mean(sil_scores_before$score), "\n")
cat("矫正后平均分数:", mean(sil_scores_after$score), "\n")
cat("分数降低:", mean(sil_scores_before$score) - mean(sil_scores_after$score), "\n")

# =============================================================================
# 完整流程：对比Harmony和rPCA
# =============================================================================

# 准备数据
seu_combined <- NormalizeData(seu_combined)
seu_combined <- FindVariableFeatures(seu_combined, nfeatures = 2000)
seu_combined <- ScaleData(seu_combined)
seu_combined <- RunPCA(seu_combined, npcs = 30)

# 方法1: Harmony矫正
library(harmony)
seu_harmony <- RunHarmony(
  seu_combined,
  group.by.vars = "dataset",
  reduction = "pca",
  dims.use = 1:30
)

# 方法2: rPCA矫正（Seurat整合）
seu_list <- SplitObject(seu_combined, split.by = "dataset")
seu_list <- lapply(seu_list, function(x) {
  x <- FindVariableFeatures(x, nfeatures = 2000)
})
features <- SelectIntegrationFeatures(seu_list, nfeatures = 2000)
seu_list <- lapply(seu_list, function(x) {
  x <- ScaleData(x, features = features)
  x <- RunPCA(x, features = features)
})

anchors <- FindIntegrationAnchors(seu_list, dims = 1:30, reduction = "rpca")
seu_rpca <- IntegrateData(anchors, dims = 1:30)

# 计算Silhouette Score
batch_labels <- seu_combined$dataset

# 原始PCA
pca_coords <- Embeddings(seu_combined, "pca")[, 1:30]
dist_pca <- dist(pca_coords)
sil_pca <- silhouette(as.numeric(factor(batch_labels)), dist_pca)

# Harmony
harmony_coords <- Embeddings(seu_harmony, "harmony")[, 1:30]
dist_harmony <- dist(harmony_coords)
sil_harmony <- silhouette(as.numeric(factor(batch_labels)), dist_harmony)

# rPCA
rpca_coords <- Embeddings(seu_rpca, "integrated")[, 1:30]
dist_rpca <- dist(rpca_coords)
sil_rpca <- silhouette(as.numeric(factor(batch_labels)), dist_rpca)

# 汇总结果
results <- data.frame(
  Method = c("No correction", "Harmony", "rPCA"),
  Mean_Silhouette = c(
    mean(sil_pca[, "sil_width"]),
    mean(sil_harmony[, "sil_width"]),
    mean(sil_rpca[, "sil_width"])
  ),
  Median_Silhouette = c(
    median(sil_pca[, "sil_width"]),
    median(sil_harmony[, "sil_width"]),
    median(sil_rpca[, "sil_width"])
  )
)

print(results)

# 可视化对比
sil_combined <- data.frame(
  score = c(
    sil_pca[, "sil_width"],
    sil_harmony[, "sil_width"],
    sil_rpca[, "sil_width"]
  ),
  method = rep(
    c("No correction", "Harmony", "rPCA"),
    each = nrow(sil_pca)
  )
)

ggplot(sil_combined, aes(x = method, y = score, fill = method)) +
  geom_violin(alpha = 0.7) +
  geom_boxplot(width = 0.2, fill = "white", outlier.shape = NA) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Comparison of Batch Correction Methods",
    subtitle = "Lower silhouette score = Better batch mixing",
    y = "Silhouette Score",
    x = "Method"
  ) +
  theme_bw() +
  theme(legend.position = "none")

# =============================================================================
# 下采样分析（像文章中那样）
# =============================================================================

# 对主要细胞类型进行20次下采样分析
calculate_silhouette_downsample <- function(coords, batch_labels, downsample_ratio = 0.2, n_iter = 20) {
  all_scores <- list()
  
  for (i in 1:n_iter) {
    # 下采样20%的细胞
    sample_idx <- sample(1:nrow(coords), size = floor(nrow(coords) * downsample_ratio))
    coords_sub <- coords[sample_idx, ]
    batch_sub <- batch_labels[sample_idx]
    
    # 计算Silhouette Score
    dist_sub <- dist(coords_sub)
    sil_sub <- silhouette(as.numeric(factor(batch_sub)), dist_sub)
    all_scores[[i]] <- sil_sub[, "sil_width"]
  }
  
  # 汇总所有迭代的分数
  mean(unlist(all_scores))
}

# 对CD4+ T细胞和CD8+ T细胞分别计算
cd4_cells <- subset(seu_harmony, subset = celltype == "CD4+ T")
cd8_cells <- subset(seu_harmony, subset = celltype == "CD8+ T")

cd4_harmony_coords <- Embeddings(cd4_cells, "harmony")[, 1:30]
cd8_harmony_coords <- Embeddings(cd8_cells, "harmony")[, 1:30]

cd4_sil <- calculate_silhouette_downsample(
  cd4_harmony_coords, 
  cd4_cells$dataset,
  downsample_ratio = 0.2,
  n_iter = 20
)

cd8_sil <- calculate_silhouette_downsample(
  cd8_harmony_coords, 
  cd8_cells$dataset,
  downsample_ratio = 0.2,
  n_iter = 20
)

cat("CD4+ T cells Silhouette Score:", cd4_sil, "\n")
cat("CD8+ T cells Silhouette Score:", cd8_sil, "\n")