381_单细胞数据伪批量分析¶

原始单细胞数据（细胞×基因）
    ↓ 按样本 + 细胞类型分组合并
伪批量矩阵（每个"样本-细胞类型"组合是一行）
    ↓ DESeq2 / edgeR 差异分析
差异基因结果（基于样本级别统计）

library(Seurat)      # 单细胞分析
library(DESeq2)      # 差异分析（伪批量最佳搭配）
library(dplyr)
library(ggplot2)

# 假设 seurat_obj 已有:
# - cell_type: 细胞类型注释
# - sample_id: 每个细胞来自哪个供体/样本
# - condition: 疾病状态 (Disease/Control)

# ============================================================
# 步骤1：按"样本 + 细胞类型"分组，对reads计数求和
# ============================================================
# 创建"样本-细胞类型"组合的标识
seurat_obj$sample_celltype <- paste(
  seurat_obj$sample_id,    # 样本ID（如 P001, P002...）
  seurat_obj$cell_type,    # 细胞类型（如 T_cell, B_cell...）
  sep = "_"
)

# 聚合reads：对每个"样本-细胞类型"组合的所有细胞求和
# 使用 AggregateExpression（Seurat v5专用函数）
pseudo_bulk <- AggregateExpression(
  seurat_obj,
  group.by = c("sample_id", "cell_type"),   # 按样本+细胞类型分组
  assays = "RNA",
  slot = "counts",                            # 使用raw counts（DESeq2要求）
  return.seurat = FALSE                       # 返回矩阵而非Seurat对象
)
pb_counts <- pseudo_bulk$RNA   # 取RNA层

# 查看结果维度
cat("伪批量矩阵维度：", dim(pb_counts), "\n")
# 格式：行=基因，列="sample_id_cell_type"

# ============================================================
# 步骤2：只分析某一细胞类型（以T细胞为例）
# ============================================================
target_celltype <- "T_cell"
tcell_cols <- grepl(paste0("_", target_celltype, "$"), colnames(pb_counts))
tcell_pb <- pb_counts[, tcell_cols]    # 只取T细胞列

# 构建样本元数据
tcell_meta <- data.frame(
  col_name = colnames(tcell_pb)
) %>%
  mutate(
    sample_id = sub(paste0("_", target_celltype), "", col_name),  # 提取样本ID
  ) %>%
  left_join(
    # 从Seurat对象提取样本级元数据（每个样本一行）
    seurat_obj@meta.data %>%
      group_by(sample_id, condition) %>%
      summarise(.groups = "drop"),
    by = "sample_id"
  )
rownames(tcell_meta) <- tcell_meta$col_name

# 检查：过滤掉T细胞数量太少的伪批量样本
cell_counts <- table(seurat_obj$sample_id[seurat_obj$cell_type == target_celltype])
keep_samples <- names(cell_counts)[cell_counts >= 20]   # 至少20个T细胞
tcell_pb <- tcell_pb[, rownames(tcell_meta)[tcell_meta$sample_id %in% keep_samples]]
tcell_meta <- tcell_meta[colnames(tcell_pb), ]

cat("保留样本数:", ncol(tcell_pb), "\n")

# ============================================================
# 步骤3：DESeq2 差异分析
# ============================================================
dds <- DESeqDataSetFromMatrix(
  countData = tcell_pb,           # 整数count矩阵（基因×伪批量样本）
  colData   = tcell_meta,         # 样本元数据
  design    = ~ condition         # 统计模型：按疾病状态比较
)

# 过滤低表达基因（加速分析，减少多重检验负担）
keep_genes <- rowSums(counts(dds) >= 10) >= 3   # 至少3个样本中有>=10个reads
dds <- dds[keep_genes, ]

# 运行 DESeq2
dds <- DESeq(dds)

# 提取结果（Disease vs Control）
res <- results(
  dds,
  contrast = c("condition", "Disease", "Control"),
  alpha = 0.05       # FDR 阈值
)
res <- as.data.frame(res)
res <- res[order(res$pvalue), ]    # 按 p 值排序

# 显著差异基因
sig_genes <- res[!is.na(res$padj) & res$padj < 0.05 & abs(res$log2FoldChange) > 1, ]
cat("显著差异基因数:", nrow(sig_genes), "\n")

# ============================================================
# 步骤4：可视化
# ============================================================
# 火山图
res$status <- "NS"
res$status[res$padj < 0.05 & res$log2FoldChange > 1] <- "Up"
res$status[res$padj < 0.05 & res$log2FoldChange < -1] <- "Down"

ggplot(res[!is.na(res$padj), ], aes(x = log2FoldChange, y = -log10(padj),
                                      color = status)) +
  geom_point(alpha = 0.6, size = 1) +
  scale_color_manual(values = c("Up" = "red", "Down" = "blue", "NS" = "gray")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  labs(title = paste("T细胞伪批量分析：Disease vs Control"),
       x = "Log2 Fold Change", y = "-Log10 adj.p") +
  theme_classic()