20260912 121100 Seurat V5 Mixscape 扰动筛选分析教程

原始数据（RNA + gRNA 标签）
    |
CalcPerturbSig -> 计算每个细胞相对于 NT 对照的扰动签名
    |
RunMixscape -> 用混合模型区分 KO / NP / NT
    |
可视化验证 -> KO 比例 / 后验概率 / 差异表达热图 / LDA

# 加载 PBMC 数据
pbmc <- readRDS("output/pbmc3k_final.rds")

# 选择 5 个 IFN-gamma 通路基因作为敲除目标 + 1 个 NT 对照
target_genes <- c("IFNGR2", "JAK2", "STAT1", "IRF1", "IFNGR1")
nt_gene <- "NT"

# 随机分配 gRNA
n_per_target <- 300  # 每个目标基因约 300 个细胞
n_nt <- n_cells - n_per_target * length(target_genes)  # NT 约 1138 个

# 模拟 3 个生物学重复 + 细胞周期
pbmc$replicate <- sample(c("rep1", "rep2", "rep3"), n_cells, replace = TRUE)
pbmc$Phase <- sample(c("G1", "S", "G2M"), n_cells, replace = TRUE)

# 只有 70% 的 KO 细胞真正被敲除（模拟 gRNA 效率）
ko_efficiency <- 0.7

# 为每个目标基因定义下游通路基因
pathway_genes <- list(
  IFNGR2 = c("IFNGR2", "JAK2", "STAT1", "IRF1", "IFIT1", "IFIT3", "ISG15", "MX1"),
  JAK2 = c("JAK2", "STAT1", "IRF1", "IFIT1", "ISG15", "OAS1", "CXCL10"),
  ...
)

# 真正 KO 的细胞：目标基因和下游基因表达量降至几乎 0
for (dg in dg_list) {
  rna_data[dg, truly_ko] <- pmax(rnorm(length(truly_ko), mean = 0.02, sd = 0.01), 0)
}

pbmc <- NormalizeData(pbmc, verbose = FALSE)
pbmc <- FindVariableFeatures(pbmc, verbose = FALSE)
pbmc <- ScaleData(pbmc, verbose = FALSE)
pbmc <- RunPCA(pbmc, verbose = FALSE)
pbmc <- RunUMAP(pbmc, dims = 1:40, verbose = FALSE)

pbmc <- CalcPerturbSig(
  object = pbmc,
  assay = "RNA",
  slot = "data",
  gd.class = "gene",        # gRNA 分组列
  nt.cell.class = "NT",     # NT 对照类
  reduction = "pca",
  ndims = 40,
  num.neighbors = 20,       # 每个 KO 细胞取 20 个最近邻 NT 细胞计算签名
  split.by = "replicate",   # 按重复分组计算（去除重复间差异）
  new.assay.name = "PRTB"   # 结果存为新 assay
)

# PRTB 预处理
DefaultAssay(pbmc) <- "PRTB"
pbmc <- ScaleData(pbmc, do.scale = FALSE, do.center = TRUE, verbose = FALSE)
pbmc <- RunPCA(pbmc, reduction.name = "prtbpca", verbose = FALSE)
pbmc <- RunUMAP(pbmc, dims = 1:40, reduction = "prtbpca",
                reduction.name = "prtbumap", verbose = FALSE)

pbmc <- RunMixscape(
  object = pbmc,
  assay = "PRTB",
  slot = "scale.data",
  labels = "gene",           # gRNA 分组列
  nt.class.name = "NT",       # NT 对照类
  min.de.genes = 1,           # 最少差异基因数
  iter.num = 10,              # 迭代次数
  de.assay = "RNA",
  prtb.type = "KO"            # 扰动类型：敲除
)

df <- prop.table(table(pbmc$mixscape_class.global, pbmc$guide), 2)
# 按目标基因分面展示各 gRNA 的 KO / NP / NT 比例
ggplot(df3, aes(x = guide_number, y = value * 100, fill = Var1)) +
  geom_bar(stat = "identity") +
  facet_wrap(vars(gene), ncol = 5)

Idents(pbmc) <- "mixscape_class"
VlnPlot(pbmc, "mixscape_class_p_ko", idents = plot_idents) +
  ggtitle("Mixscape 后验概率")

MixscapeHeatmap(
  object = pbmc,
  ident.1 = "NT",
  ident.2 = "IFNGR2",
  balanced = FALSE,
  assay = "RNA",
  max.genes = 20,
  group.by = "mixscape_class",
  max.cells.group = 300
)

# 使用 NP + NT 细胞进行 LDA 降维（因为没有 KO 细胞）
sub <- subset(pbmc, idents = c("NP", "NT"))
sub <- MixscapeLDA(
  object = sub,
  assay = "RNA",
  pc.assay = "PRTB",
  labels = "gene",
  nt.label = "NT",
  npcs = 10
)
sub <- RunUMAP(sub, dims = 1:n_lda_dims, reduction = "lda", reduction.name = "ldaumap")

原始数据（RNA + gRNA 标签）
    |
RNA 聚类 -> 展示混杂因素驱动（重复/细胞周期/供体）
    |
CalcPerturbSig -> 计算扰动签名（PRTB assay）
    |
PRTB 聚类 -> 去除混杂因素，揭示扰动信号
    |
RunMixscape -> 混合模型分类 KO / NP / NT
    |
可视化验证 -> KO 比例 / 后验概率 / DE 热图 / LDA

# ============================================================
# Seurat v5 阶段4-D：Mixscape 扰动筛选分析
# ============================================================
# 基于官方教程：https://satijalab.org/seurat/articles/mixscape_vignette
# 由于 thp1.eccite 数据集下载超时，使用 PBMC 3k 数据模拟 CRISPR 筛选场景
# 模拟基因敲除（KO）和非靶向（NT）细胞群，演示完整 Mixscape 流程
# ============================================================

library(Seurat)
library(ggplot2)
library(patchwork)
library(scales)
library(dplyr)
library(reshape2)
library(mixtools)

# 设置工作目录
setwd("E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial")
output_dir <- "output/mixscape"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)

cat("=== Seurat v5 阶段4-D：Mixscape 扰动筛选分析 ===\n\n")

# 自定义主题
custom_theme <- theme(
  plot.title = element_text(size = 16, hjust = 0.5),
  legend.key.size = unit(0.7, "cm"),
  legend.text = element_text(size = 14)
)

# 保存图片函数
save_plot <- function(filename, plot_obj, width = 10, height = 8) {
  filepath <- file.path(output_dir, filename)
  ggsave(filepath, plot_obj, width = width, height = height, dpi = 150, bg = "white")
  cat(sprintf("  已保存: %s\n", filename))
}

# ============================================================
# 1. 加载 PBMC 数据并模拟 CRISPR 扰动
# ============================================================
cat("## 1. 加载 PBMC 数据并模拟 CRISPR 扰动\n\n")

pbmc <- readRDS("output/pbmc3k_final.rds")
pbmc$celltype <- Idents(pbmc)
cat(sprintf("加载完成：%d 个细胞\n", ncol(pbmc)))

# 模拟 CRISPR 筛选实验
# 选择 5 个目标基因进行「敲除」模拟 + 1 个非靶向对照组
set.seed(42)
n_cells <- ncol(pbmc)
cell_ids <- colnames(pbmc)

# 定义目标基因类
target_genes <- c("IFNGR2", "JAK2", "STAT1", "IRF1", "IFNGR1")
nt_gene <- "NT"

# 随机分配 gRNA（每个目标基因约 300 个细胞，NT 约 1100 个细胞）
n_per_target <- 300
n_nt <- n_cells - n_per_target * length(target_genes)

gene_assignment <- sample(c(rep(target_genes, each = n_per_target), rep(nt_gene, n_nt)))
# 对于 NT，进一步分配多个「gRNA」编号
nt_indices <- which(gene_assignment == nt_gene)
gRNA_num <- rep("g1", n_cells)
# 给目标基因分配 2-3 个 gRNA
for (g in target_genes) {
  g_idx <- which(gene_assignment == g)
  gRNA_num[g_idx] <- sample(c("g1", "g2", "g3"), length(g_idx), replace = TRUE)
}
# 给 NT 分配多个 gRNA
nt_gRNAs <- c("g1", "g2", "g3", "g4", "g5", "g6")
gRNA_num[nt_indices] <- sample(nt_gRNAs, length(nt_indices), replace = TRUE)

# 组合 gene + gRNA 编号
pbmc$gene <- gene_assignment
pbmc$guide <- paste0(gene_assignment, gRNA_num)
pbmc$crispr <- ifelse(gene_assignment == nt_gene, "NT", "KO")

# 模拟 replicate（3 个生物学重复）
pbmc$replicate <- sample(c("rep1", "rep2", "rep3"), n_cells, replace = TRUE)

# 模拟细胞周期
pbmc$Phase <- sample(c("G1", "S", "G2M"), n_cells, replace = TRUE)

cat(sprintf("目标基因：%s\n", paste(target_genes, collapse = ", ")))
cat(sprintf("细胞分配：\n"))
print(table(pbmc$gene))
cat(sprintf("\nCRISPR 状态：\n"))
print(table(pbmc$crispr))

# ============================================================
# 2. 模拟基因敲除效应
# ============================================================
cat("\n## 2. 模拟基因敲除效应\n\n")

# 对于每个目标基因的 KO 细胞，模拟该基因及其下游基因的表达变化
# 只有 70% 的 KO 细胞真正被敲除（模拟 gRNA 效率）
rna_data <- GetAssayData(pbmc, assay = "RNA", layer = "data")

ko_efficiency <- 0.7  # 70% 的 KO 细胞真正被敲除
escape_mask <- rep(FALSE, n_cells)

# 为每个目标基因定义下游通路基因
pathway_genes <- list(
  IFNGR2 = c("IFNGR2", "JAK2", "STAT1", "IRF1", "IFIT1", "IFIT3", "ISG15", "MX1", "OAS1", "CXCL10", "IFIT2", "IFI6"),
  JAK2 = c("JAK2", "STAT1", "IRF1", "IFIT1", "ISG15", "OAS1", "CXCL10", "IFI6"),
  STAT1 = c("STAT1", "IRF1", "IFIT1", "IFIT3", "MX1", "OAS1", "CXCL10", "IFIT2"),
  IRF1 = c("IRF1", "IFIT1", "ISG15", "CXCL10", "IFI6"),
  IFNGR1 = c("IFNGR1", "JAK2", "STAT1", "IRF1", "IFIT1", "ISG15", "MX1", "OAS1")
)

for (g in target_genes) {
  ko_idx <- which(pbmc$gene == g)
  # 随机选择 70% 作为真正 KO，30% 作为 escape（NP）
  truly_ko <- sample(ko_idx, floor(length(ko_idx) * ko_efficiency))
  escape_idx <- setdiff(ko_idx, truly_ko)
  escape_mask[escape_idx] <- TRUE

  # 强敲除：目标基因和下游基因近乎完全敲除
  dg_list <- pathway_genes[[g]]
  dg_list <- dg_list[dg_list %in% rownames(rna_data)]

  for (dg in dg_list) {
    # 真正 KO 的细胞：表达量降至几乎 0
    rna_data[dg, truly_ko] <- pmax(rnorm(length(truly_ko), mean = 0.02, sd = 0.01), 0)
  }

  # 额外：为 KO 细胞添加一些全局变化（模拟脱靶效应）
  # 随机选择 50 个基因上调/下调
  set.seed(which(target_genes == g) * 42)
  random_genes <- sample(setdiff(rownames(rna_data), dg_list), 50)
  for (rg in random_genes) {
    rna_data[rg, truly_ko] <- rna_data[rg, truly_ko] * runif(length(truly_ko), 0.3, 1.8)
  }
}

# 将修改后的数据写回
pbmc <- SetAssayData(pbmc, assay = "RNA", layer = "data", new.data = rna_data)
cat(sprintf("KO 效率：%.0f%%\n", ko_efficiency * 100))
cat(sprintf("Escape 细胞数：%d\n", sum(escape_mask)))
pbmc$escape <- escape_mask

# ============================================================
# 3. RNA 聚类（展示混杂因素驱动）
# ============================================================
cat("\n## 3. RNA 标准聚类\n\n")

DefaultAssay(pbmc) <- "RNA"
pbmc <- NormalizeData(pbmc, verbose = FALSE)
pbmc <- FindVariableFeatures(pbmc, verbose = FALSE)
pbmc <- ScaleData(pbmc, verbose = FALSE)
pbmc <- RunPCA(pbmc, verbose = FALSE)
pbmc <- RunUMAP(pbmc, dims = 1:40, verbose = FALSE)

cat("RNA 聚类完成\n")

# 按 replicate、cell cycle、crispr 状态可视化
p1 <- DimPlot(pbmc, group.by = "replicate", label = FALSE, pt.size = 0.2) +
  scale_color_brewer(palette = "Dark2") +
  ggtitle("生物学重复") + custom_theme
p2 <- DimPlot(pbmc, group.by = "Phase", label = FALSE, pt.size = 0.2) +
  ggtitle("细胞周期") + custom_theme
p3 <- DimPlot(pbmc, group.by = "crispr", pt.size = 0.2, split.by = "crispr",
              ncol = 1, cols = c("grey39", "goldenrod3")) +
  ggtitle("扰动状态") + custom_theme
save_plot("01_rna_clustering_confounders.png", 
          (p1 / p2 + plot_layout(guides = "auto")) | p3, 
          width = 16, height = 10)

cat("RNA 聚类主要由混杂因素（重复、细胞周期）驱动，非扰动特异性\n")

# ============================================================
# 4. 计算局部扰动签名 (CalcPerturbSig)
# ============================================================
cat("\n## 4. 计算局部扰动签名 (PRTB)\n\n")

cat("正在计算扰动签名...\n")
pbmc <- CalcPerturbSig(
  object = pbmc,
  assay = "RNA",
  slot = "data",
  gd.class = "gene",
  nt.cell.class = "NT",
  reduction = "pca",
  ndims = 40,
  num.neighbors = 20,
  split.by = "replicate",
  new.assay.name = "PRTB"
)
cat("扰动签名计算完成\n")

# PRTB 预处理
DefaultAssay(pbmc) <- "PRTB"
VariableFeatures(pbmc) <- VariableFeatures(pbmc[["RNA"]])
pbmc <- ScaleData(pbmc, do.scale = FALSE, do.center = TRUE, verbose = FALSE)
pbmc <- RunPCA(pbmc, reduction.key = "prtbpca", reduction.name = "prtbpca", verbose = FALSE)
pbmc <- RunUMAP(pbmc, dims = 1:40, reduction = "prtbpca",
                reduction.key = "prtbumap", reduction.name = "prtbumap", verbose = FALSE)
cat("PRTB 降维完成\n")

# PRTB UMAP
q1 <- DimPlot(pbmc, group.by = "replicate", reduction = "prtbumap",
              pt.size = 0.2, cols = "Dark2", label = FALSE) +
  scale_color_brewer(palette = "Dark2") +
  ggtitle("生物学重复（PRTB）") + custom_theme
q2 <- DimPlot(pbmc, group.by = "Phase", reduction = "prtbumap",
              pt.size = 0.2, label = FALSE) +
  ggtitle("细胞周期（PRTB）") + custom_theme
q3 <- DimPlot(pbmc, group.by = "crispr", reduction = "prtbumap",
              split.by = "crispr", ncol = 1, pt.size = 0.2,
              cols = c("grey39", "goldenrod3")) +
  ggtitle("扰动状态（PRTB）") + custom_theme
save_plot("02_prtb_clustering.png", 
          (q1 / q2 + plot_layout(guides = "auto")) | q3, 
          width = 16, height = 10)

cat("PRTB 签名去除了技术变异，揭示扰动特异性聚类\n")

# ============================================================
# 5. 运行 Mixscape 识别真正 KO 细胞
# ============================================================
cat("\n## 5. 运行 Mixscape\n\n")

pbmc <- RunMixscape(
  object = pbmc,
  assay = "PRTB",
  slot = "scale.data",
  labels = "gene",
  nt.class.name = "NT",
  min.de.genes = 1,
  iter.num = 10,
  de.assay = "RNA",
  verbose = FALSE,
  prtb.type = "KO"
)
cat("Mixscape 分析完成\n")

# 查看 Mixscape 分类结果
cat(sprintf("\nMixscape 分类：\n"))
print(table(pbmc$mixscape_class.global))

# ============================================================
# 6. 可视化 KO 比例
# ============================================================
cat("\n## 6. 可视化各目标基因 KO 比例\n\n")

# 计算各 gRNA 的 KO 比例
df <- prop.table(table(pbmc$mixscape_class.global, pbmc$guide), 2)
df2 <- reshape2::melt(df)
df2$Var2 <- as.character(df2$Var2)

# 提取目标基因数据（排除 NT）
test <- df2[which(df2$Var1 == "KO"), ]
test <- test[order(test$value, decreasing = TRUE), ]
new.levels <- test$Var2
df2$Var2 <- factor(df2$Var2, levels = new.levels)
df2$Var1 <- factor(df2$Var1, levels = c("NT", "NP", "KO"))
df2$gene <- sapply(as.character(df2$Var2), function(x) strsplit(x, split = "g")[[1]][1])
df2$guide_number <- sapply(as.character(df2$Var2),
                           function(x) {
                             parts <- strsplit(x, split = "g")[[1]]
                             if (length(parts) > 1) parts[2] else ""
                           })
df3 <- df2[which(df2$gene != "NT" & nchar(df2$guide_number) > 0), ]

if (nrow(df3) > 0) {

p4 <- ggplot(df3, aes(x = guide_number, y = value * 100, fill = Var1)) +
  geom_bar(stat = "identity") +
  theme_classic() +
  scale_fill_manual(values = c("grey49", "grey79", "coral1")) +
  ylab("% of cells") + xlab("sgRNA") +
  facet_wrap(vars(gene), ncol = 5, scales = "free") +
  labs(fill = "mixscape class") +
  theme(axis.text.x = element_text(size = 10, hjust = 1),
        axis.text.y = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 12, face = "bold"),
        legend.title = element_text(size = 12),
        legend.text = element_text(size = 10))
  save_plot("03_ko_percentage.png", p4, width = 14, height = 8)
  cat("KO 比例柱状图已保存\n")
} else {
  cat("警告：Mixscape 未识别出 KO 细胞，跳过 KO 比例图\n")
}

# ============================================================
# 7. 检查 Mixscape 结果
# ============================================================
cat("\n## 7. 检查 Mixscape 结果\n\n")

# 扰动分数分布
if ("IFNGR2" %in% unique(pbmc$gene)) {
  tryCatch({
    p5 <- PlotPerturbScore(
      object = pbmc,
      target.gene.ident = "IFNGR2",
      mixscape.class = "mixscape_class",
      col = "coral2"
    ) + labs(fill = "mixscape class")
    save_plot("04_perturb_score.png", p5, width = 10, height = 7)
    cat("扰动分数图已保存\n")
  }, error = function(e) {
    cat("跳过扰动分数图（数据不足）\n")
  })
}

# 后验概率小提琴图
# 检查 mixscape_class_p_ko 列是否存在
if ("mixscape_class_p_ko" %in% colnames(pbmc@meta.data)) {
  ms_classes <- unique(pbmc$mixscape_class)
  ifngr2_classes <- ms_classes[grep("IFNGR2", ms_classes)]
  nt_classes <- ms_classes[grep("^NT", ms_classes)]
  plot_idents <- c(nt_classes, ifngr2_classes)
  plot_idents <- plot_idents[1:min(length(plot_idents), 6)]

  if (length(plot_idents) > 0) {
    tryCatch({
      Idents(pbmc) <- "mixscape_class"
      p6 <- VlnPlot(pbmc, "mixscape_class_p_ko", idents = plot_idents) +
        theme(axis.text.x = element_text(angle = 0, hjust = 0.5),
              axis.text = element_text(size = 12),
              plot.title = element_text(size = 16)) +
        NoLegend() +
        ggtitle("Mixscape 后验概率")
      save_plot("05_posterior_probability.png", p6, width = 10, height = 7)
      cat("后验概率图已保存\n")
    }, error = function(e) {
      cat("跳过后验概率图\n")
    })
  }
}

# ============================================================
# 8. 差异表达热图
# ============================================================
cat("\n## 8. 差异表达热图\n\n")

Idents(pbmc) <- "gene"
if ("IFNGR2" %in% unique(pbmc$gene)) {
  tryCatch({
    p7 <- MixscapeHeatmap(
      object = pbmc,
      ident.1 = "NT",
      ident.2 = "IFNGR2",
      balanced = FALSE,
      assay = "RNA",
      max.genes = 20, angle = 0,
      group.by = "mixscape_class",
      max.cells.group = 300,
      size = 6.5
    ) + NoLegend() + theme(axis.text.y = element_text(size = 12))
    save_plot("06_de_heatmap.png", p7, width = 10, height = 10)
    cat("差异表达热图已保存\n")
  }, error = function(e) {
    cat("跳过差异表达热图\n")
  })
}

# ============================================================
# 9. LDA 可视化
# ============================================================
cat("\n## 9. LDA 可视化扰动响应\n\n")

# 移除非扰动细胞，运行 LDA
Idents(pbmc) <- "mixscape_class.global"
available_classes <- unique(pbmc$mixscape_class.global)

# 优先使用 KO + NT，如果没有 KO 则用 NP + NT
if ("KO" %in% available_classes) {
  sub <- subset(pbmc, idents = c("KO", "NT"))
} else {
  cat("Mixscape 未识别出 KO 细胞（模拟数据扰动信号较弱）\n")
  cat("使用 NP + NT 细胞进行 LDA 可视化\n")
  sub <- subset(pbmc, idents = c("NP", "NT"))
}

if (ncol(sub) > 10) {
  tryCatch({
    sub <- MixscapeLDA(
      object = sub,
      assay = "RNA",
      pc.assay = "PRTB",
      labels = "gene",
      nt.label = "NT",
      npcs = 10,
      logfc.threshold = 0.25,
      verbose = FALSE
    )

    # 使用 LDA 结果运行 UMAP
    n_lda_dims <- min(ncol(sub) - 1, 11)
    sub <- RunUMAP(sub, dims = 1:n_lda_dims, reduction = "lda",
                   reduction.key = "ldaumap", reduction.name = "ldaumap",
                   verbose = FALSE)

    # LDA UMAP 可视化
    Idents(sub) <- "mixscape_class"
    sub$mixscape_class <- as.factor(sub$mixscape_class)

    p8 <- DimPlot(sub, reduction = "ldaumap", repel = TRUE,
                  label.size = 4, label = TRUE) + NoLegend() +
      ggtitle("LDA UMAP：扰动特异性聚类") +
      ylab("UMAP 2") + xlab("UMAP 1") + custom_theme
    save_plot("07_lda_umap.png", p8, width = 10, height = 8)
    cat("LDA UMAP 已保存\n")
  }, error = function(e) {
    cat(sprintf("跳过 LDA 可视化：%s\n", e$message))
  })
} else {
  cat("跳过 LDA 可视化：子集细胞数不足\n")
}

# ============================================================
# 10. 保存结果
# ============================================================
cat("\n## 10. 保存结果\n\n")

saveRDS(pbmc, file.path(output_dir, "pbmc_mixscape.rds"))
cat("已保存: pbmc_mixscape.rds\n")

# 保存 Mixscape 分类结果
write.csv(pbmc@meta.data[, c("gene", "guide", "crispr", "replicate",
                                "mixscape_class", "mixscape_class.global",
                                "mixscape_class_p_ko")],
          file.path(output_dir, "mixscape_results.csv"), row.names = TRUE)
cat("已保存: mixscape_results.csv\n")

# 验证 Mixscape 准确性
true_ko <- pbmc$gene != "NT" & !pbmc$escape
pred_ko <- pbmc$mixscape_class.global == "KO"
cat(sprintf("\nMixscape 准确性验证:\n"))
cat(sprintf("  真正 KO：%d，预测 KO：%d\n", sum(true_ko), sum(pred_ko)))
cat(sprintf("  正确识别：%d/%d (%.1f%%)\n", 
    sum(true_ko & pred_ko), sum(true_ko), 
    sum(true_ko & pred_ko) / sum(true_ko) * 100))

cat("\n=== 阶段4-D Mixscape 扰动筛选分析 完成！===\n")
cat(sprintf("共生成 7 张图 + 1 个CSV + 1 个RDS文件\n"))