Redeconve — 空间转录组单细胞分辨率反卷积

概述

Redeconve 是用于空间转录组数据单细胞分辨率反卷积的 R 包。它通过整合 scRNA-seq 参考数据，估计每个空间 spot 中各细胞类型的组成，并提供可视化、共定位网络、基因过滤、熵困惑度计算等高级分析功能。

核心特性：
- 双层反卷积：单细胞级别（每个细胞）和细胞类型级别（每个细胞类型）
- 自动基因过滤：gene.filter() 自动筛选信息基因
- 共定位网络分析：coloc.network() + 子图提取
- 空间可视化：spatial.piechart()、spatial.cell.number()、spot.pie()
- 质量评估：残差 (resids)、熵与困惑度 (hNpp)
- Seurat 接口：extract.seurat.sc() 和 extract.seurat.st()
- 辅助分析：细胞采样 (cell.sampling)、空间基因表达 (spatial.gene)、DEG 点图 (deg.dotplot)

适用场景

当用户需要执行以下任务时使用此 Skill：

• 使用 Redeconve 对空间转录组数据进行反卷积
• 单细胞级别反卷积（每个细胞作为参考）
• 细胞类型级别反卷积（构建参考矩阵后）
• 共定位网络分析（细胞类型间空间相关性）
• 残差/熵困惑度等质量评估
• 通过 Seurat 接口处理数据

安装

# 从 GitHub 安装
devtools::install_github("ZxZhou4150/Redeconve")

# 加载依赖
library(Redeconve); library(Seurat); library(Matrix)
library(ggplot2); library(igraph); library(reshape2); library(jsonlite)

验证安装

library(Redeconve)
ver <- as.character(packageVersion("Redeconve"))
print(paste("Redeconve版本:", ver))

Redeconve 核心函数

反卷积

Seurat 接口

可视化

高级分析

关键参数

完整可运行脚本

输入数据要求

配置说明

所有脚本顶部都包含统一配置区域：

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR      <- "your/data/directory"        # 数据根目录
SPATIAL_SUBDIR     <- "binned_outputs/square_016um" # 空间数据子目录（相对 USER_DATA_DIR）
CELL_TYPE_COL      <- "cell_type"                   # scRNA-seq metadata 中细胞类型列名
SC_FILENAME        <- "scRNA.rds"                   # 单细胞 RDS 文件名
N_CELLS_TARGET     <- 10000                         # scRNA 采样目标（0=不采样）
N_SPOTS_TARGET     <- 10000                         # 空间采样目标（0=不采样）
OUTPUT_DIR         <- "Redeconve_results"           # 输出目录
# =====================================================

路径兼容性

• Windows：USER_DATA_DIR <- "D:/your/data/dir"
• WSL/Linux：USER_DATA_DIR <- "/mnt/d/your/data/dir" 或 "/home/user/data"
• macOS：USER_DATA_DIR <- "/Users/yourname/data"

> 注意：以下脚本均假设单细胞注释列名为 cell_type，脚本内部会自动回退到 Clusters、celltype。若你的数据列名不同，请修改 CELL_TYPE_COL。

脚本 1 — Redeconve 单细胞级别反卷积

#!/usr/bin/env Rscript
# 01_redeconve_single_cell_level.R
# Redeconve 单细胞级别反卷积分析
# 完全按照官方流程: https://github.com/ZxZhou4150/Redeconve
#
# 官方示例:
#   res <- deconvoluting(sc, st, genemode = "def", hpmode = "def", dopar = T, ncores = 8)

library(Redeconve); library(Seurat); library(Matrix)

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR  <- "your/data/directory"             # TODO: 修改为你的数据根目录
SC_FILENAME    <- "scRNA.rds"                       # 单细胞 RDS 文件名
SPATIAL_SUBDIR <- "binned_outputs/square_016um"    # 空间数据子目录
CELL_TYPE_COL  <- "cell_type"                       # scRNA-seq metadata 中细胞类型列名
N_CELLS_TARGET <- 1000                              # 单细胞采样目标（单细胞级别计算量大，建议 ≤1000）
N_SPOTS_TARGET <- 10000                             # 空间采样目标（0=不采样）
OUTPUT_DIR     <- "Redeconve_results"               # 输出目录
# =====================================================

print("Redeconve 单细胞级别反卷积分析")
print("官方流程: deconvoluting(sc, st, genemode='"'"'def'"'"', hpmode='"'"'def'"'"')")

# 路径组装
sc_data_path <- file.path(USER_DATA_DIR, SC_FILENAME)
spatial_dir  <- file.path(USER_DATA_DIR, SPATIAL_SUBDIR)
output_dir   <- OUTPUT_DIR
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)

# 步骤1: 加载单细胞参考数据
sc_seurat <- readRDS(sc_data_path)
print(paste("  ✓ 单细胞数据:", ncol(sc_seurat), "cells ×", nrow(sc_seurat), "genes"))

# 自动检测细胞类型列名（按 cell_type -> Clusters -> celltype 顺序回退）
resolve_cell_type_col <- function(seurat_obj, default_col) {
  candidates <- c(default_col, "cell_type", "Clusters", "celltype", "cell.type", "annot")
  for (c in candidates) {
    if (c %in% colnames(seurat_obj@meta.data)) return(c)
  }
  stop("无法找到细胞类型列。请修改 CELL_TYPE_COL 参数。")
}
CELL_TYPE_COL <- resolve_cell_type_col(sc_seurat, CELL_TYPE_COL)
print(paste("  使用细胞类型列:", CELL_TYPE_COL))

# 下采样到目标细胞数（单细胞级别计算量大，建议 ≤1000）
if (N_CELLS_TARGET > 0 && ncol(sc_seurat) > N_CELLS_TARGET) {
  set.seed(42)
  cells_to_keep <- sample(colnames(sc_seurat), N_CELLS_TARGET)
  sc_seurat <- sc_seurat[, cells_to_keep]
}
sc <- GetAssayData(sc_seurat, layer = "counts")

# 步骤2: 加载空间转录组数据
st_seurat <- Load10X_Spatial(data.dir = spatial_dir)
print(paste("  ✓ 空间数据:", ncol(st_seurat), "spots ×", nrow(st_seurat), "genes"))

if (ncol(st_seurat) >= N_SPOTS_TARGET) {
  set.seed(42)
  spots_to_keep <- sample(colnames(st_seurat), N_SPOTS_TARGET)
  st_seurat <- st_seurat[, spots_to_keep]
}
st <- GetAssayData(st_seurat, layer = "counts")

# 步骤3: 找到共同基因
common_genes <- intersect(rownames(sc), rownames(st))
if (length(common_genes) < 100) stop("共同基因太少")
sc <- sc[common_genes, ]
st <- st[common_genes, ]

# 步骤4: 运行 Redeconve 反卷积（单细胞级别）
res <- deconvoluting(
  ref = sc,
  st = st,
  genemode = "def",
  hpmode = "def",
  normalize = TRUE,
  dopar = FALSE,
  thre = 1e-10
)
print(paste("  ✓ 结果维度:", nrow(res), "cells ×", ncol(res), "spots"))

# 步骤5: 保存结果
res_prop <- to.proportion(res)

saveRDS(list(
  raw_results = res, proportions = res_prop,
  sc = sc, st = st
), file = file.path(output_dir, "01_redeconve_single_cell_results.rds"))

write.csv(as.data.frame(t(res_prop)),
          file = file.path(output_dir, "01_redeconve_single_cell_proportions.csv"))

# 统计摘要
cell_names <- rownames(res_prop)
mean_props <- rowMeans(res_prop)
stats_df <- data.frame(Cell = cell_names, MeanProportion = mean_props)
stats_df <- stats_df[order(-stats_df$MeanProportion), ]
write.csv(stats_df, file = file.path(output_dir, "01_redeconve_single_cell_statistics.csv"), row.names = FALSE)

print(paste("输出目录:", output_dir))

脚本 2 — Redeconve 细胞类型级别反卷积

#!/usr/bin/env Rscript
# 02_redeconve_cell_type_level.R
# Redeconve 细胞类型级别反卷积分析
# 官方流程: get.ref() -> deconvoluting()

library(Redeconve); library(Seurat); library(Matrix)

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR  <- "your/data/directory"
SC_FILENAME    <- "scRNA.rds"
SPATIAL_SUBDIR <- "binned_outputs/square_016um"
CELL_TYPE_COL  <- "cell_type"
N_CELLS_TARGET <- 10000
N_SPOTS_TARGET <- 10000
OUTPUT_DIR     <- "Redeconve_results"
# =====================================================

# 路径组装
sc_data_path <- file.path(USER_DATA_DIR, SC_FILENAME)
spatial_dir  <- file.path(USER_DATA_DIR, SPATIAL_SUBDIR)
output_dir   <- OUTPUT_DIR
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)

# 自动检测细胞类型列名
resolve_cell_type_col <- function(seurat_obj, default_col) {
  candidates <- c(default_col, "cell_type", "Clusters", "celltype", "cell.type", "annot")
  for (c in candidates) {
    if (c %in% colnames(seurat_obj@meta.data)) return(c)
  }
  stop("无法找到细胞类型列。请修改 CELL_TYPE_COL 参数。")
}

# 步骤1: 加载单细胞参考数据
sc_seurat <- readRDS(sc_data_path)
CELL_TYPE_COL <- resolve_cell_type_col(sc_seurat, CELL_TYPE_COL)
print(paste("  使用细胞类型列:", CELL_TYPE_COL))

if (N_CELLS_TARGET > 0 && ncol(sc_seurat) > N_CELLS_TARGET) {
  set.seed(42)
  cells_to_keep <- sample(colnames(sc_seurat), N_CELLS_TARGET)
  sc_seurat <- sc_seurat[, cells_to_keep]
}

sc <- GetAssayData(sc_seurat, layer = "counts")
annotations <- cbind(colnames(sc_seurat), as.character(sc_seurat[[CELL_TYPE_COL, drop = TRUE]]))
colnames(annotations) <- c("cell", "cell_type")
cell_types <- unique(sc_seurat[[CELL_TYPE_COL, drop = TRUE]])
print(paste("  ✓ 细胞类型:", length(cell_types), "种"))

# 步骤2: 加载空间转录组数据
st_seurat <- Load10X_Spatial(data.dir = spatial_dir)
if (N_SPOTS_TARGET > 0 && ncol(st_seurat) > N_SPOTS_TARGET) {
  set.seed(42)
  spots_to_keep <- sample(colnames(st_seurat), N_SPOTS_TARGET)
  st_seurat <- st_seurat[, spots_to_keep]
}
st <- GetAssayData(st_seurat, layer = "counts")

# 步骤3: 找到共同基因
common_genes <- intersect(rownames(sc), rownames(st))
sc <- sc[common_genes, ]
st <- st[common_genes, ]

# 步骤4: 构建细胞类型参考矩阵 (get.ref)
ref <- get.ref(sc = sc, annotations = annotations, dopar = FALSE)
print(paste("  ✓ 参考矩阵维度:", nrow(ref), "genes ×", ncol(ref), "cell types"))

# 步骤5: 运行 Redeconve 反卷积（细胞类型级别）
res.ct <- deconvoluting(
  ref = ref,
  st = st,
  genemode = "def",
  hpmode = "auto",
  normalize = TRUE,
  dopar = FALSE,
  thre = 1e-10
)
print(paste("  ✓ 结果维度:", nrow(res.ct), "cell types ×", ncol(res.ct), "spots"))

# 步骤6: 保存结果和生成可视化
res_prop <- to.proportion(res.ct)

saveRDS(list(
  raw_results = res.ct, proportions = res_prop,
  ref = ref, sc = sc, st = st, annotations = annotations
), file = file.path(output_dir, "02_redeconve_cell_type_results.rds"))

write.csv(as.data.frame(t(res_prop)),
          file = file.path(output_dir, "02_redeconve_cell_type_proportions.csv"))

# 统计摘要
celltype_names <- rownames(res_prop)
mean_props <- rowMeans(res_prop)
stats_df <- data.frame(CellType = celltype_names, MeanProportion = mean_props)
stats_df <- stats_df[order(-stats_df$MeanProportion), ]
write.csv(stats_df, file = file.path(output_dir, "02_redeconve_cell_type_statistics.csv"), row.names = FALSE)

# 空间分布图（前4种细胞类型）
coords <- GetTissueCoordinates(st_seurat)
common_spots <- intersect(rownames(coords), colnames(res_prop))
if (length(common_spots) > 0 && "x" %in% colnames(coords)) {
  coords <- coords[common_spots, ]
  x_coord <- coords$x; y_coord <- coords$y
  top_celltypes <- head(stats_df$CellType, 4)

  png(file.path(output_dir, "02_redeconve_celltype_spatial.png"), width = 1200, height = 1000, res = 150)
  par(mfrow = c(2, 2))
  for (ct in top_celltypes) {
    if (ct %in% rownames(res_prop)) {
      values <- res_prop[ct, common_spots]
      if (all(is.finite(values)) && length(unique(values)) > 1) {
        plot(x_coord, y_coord,
             col = colorRampPalette(c("white", "red"))(100)[cut(values, breaks = 100)],
             pch = 20, cex = 2,
             main = paste(ct, "\n(mean:", sprintf("%.3f", mean(values)), ")"),
             xlab = "X", ylab = "Y")
      }
    }
  }
  dev.off()
}

# 细胞类型比例条形图
png(file.path(output_dir, "02_redeconve_celltype_barplot.png"), width = 1000, height = 600, res = 150)
par(mar = c(10, 4, 4, 2))
barplot(stats_df$MeanProportion, names.arg = stats_df$CellType,
        las = 2, cex.names = 0.8,
        main = "Redeconve: Mean Cell Type Proportions",
        ylab = "Mean Proportion",
        col = rainbow(nrow(stats_df)))
dev.off()

print(paste("输出目录:", output_dir))

脚本 3 — Redeconve Seurat 接口

#!/usr/bin/env Rscript
# 03_redeconve_seurat_interface.R
# Redeconve Seurat 接口测试
# 官方流程: extract.seurat.sc() 和 extract.seurat.st()

library(Redeconve); library(Seurat); library(Matrix)

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR  <- "your/data/directory"
SC_FILENAME    <- "scRNA.rds"
SPATIAL_SUBDIR <- "binned_outputs/square_016um"
CELL_TYPE_COL  <- "cell_type"
N_CELLS_TARGET <- 10000
N_SPOTS_TARGET <- 10000
OUTPUT_DIR     <- "Redeconve_results"
# =====================================================

sc_data_path <- file.path(USER_DATA_DIR, SC_FILENAME)
spatial_dir  <- file.path(USER_DATA_DIR, SPATIAL_SUBDIR)
output_dir   <- OUTPUT_DIR
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)

# 步骤1: 加载单细胞 Seurat 对象
sc_seurat <- readRDS(sc_data_path)

# 自动检测细胞类型列名
resolve_cell_type_col <- function(seurat_obj, default_col) {
  candidates <- c(default_col, "cell_type", "Clusters", "celltype", "cell.type", "annot")
  for (c in candidates) {
    if (c %in% colnames(seurat_obj@meta.data)) return(c)
  }
  stop("无法找到细胞类型列。请修改 CELL_TYPE_COL 参数。")
}
CELL_TYPE_COL <- resolve_cell_type_col(sc_seurat, CELL_TYPE_COL)
print(paste("  使用细胞类型列:", CELL_TYPE_COL))
Idents(sc_seurat) <- sc_seurat[[CELL_TYPE_COL, drop = TRUE]]

if (N_CELLS_TARGET > 0 && ncol(sc_seurat) >= N_CELLS_TARGET) {
  set.seed(42)
  cells_to_keep <- sample(colnames(sc_seurat), N_CELLS_TARGET)
  sc_seurat <- sc_seurat[, cells_to_keep]
}

# 步骤2: 使用 extract.seurat.sc 提取单细胞数据
sclist <- extract.seurat.sc(sc_seurat)
sc <- sclist$expr
annotations <- sclist$annotations
print(paste("  ✓ 表达矩阵:", nrow(sc), "genes ×", ncol(sc), "cells"))

# 步骤3: 加载空间 Seurat 对象
st_seurat <- Load10X_Spatial(data.dir = spatial_dir)
if (ncol(st_seurat) >= N_SPOTS_TARGET) {
  set.seed(42)
  spots_to_keep <- sample(colnames(st_seurat), N_SPOTS_TARGET)
  st_seurat <- st_seurat[, spots_to_keep]
}

# 步骤4: 使用 extract.seurat.st 提取空间数据
extract_success <- FALSE
stlist <- NULL
st <- NULL
coords <- NULL
tryCatch({
  stlist <- extract.seurat.st(st_seurat)
  st <- stlist$expr
  coords <- stlist$coords
  extract_success <- TRUE
}, error = function(e) {
  print(paste("  ⚠ extract.seurat.st失败:", e$message))
})

# 失败时手动提取
if (!extract_success) {
  st <- GetAssayData(st_seurat, layer = "counts")
  coords <- GetTissueCoordinates(st_seurat)
  stlist <- list(expr = st, coords = coords)
}

# 步骤5: 找到共同基因并保存
common_genes <- intersect(rownames(sc), rownames(st))
if (length(common_genes) < 100) stop("共同基因太少")
sc <- sc[common_genes, ]
st <- st[common_genes, ]

saveRDS(list(
  sc = sc, st = st, annotations = annotations,
  coords = coords, sclist = sclist, stlist = stlist
), file = file.path(output_dir, "03_redeconve_seurat_extracted.rds"))

# 显示细胞类型统计
cell_types <- table(annotations[, 2])
print("细胞类型统计:")
for (ct in names(cell_types)) {
  print(paste("  ", ct, ":", cell_types[ct], "cells"))
}

print("说明: extract.seurat.sc() 提取sc数据和注释；extract.seurat.st() 提取st数据和坐标")
print("后续可使用 01/02 脚本运行反卷积")

脚本 4 — Redeconve 可视化

#!/usr/bin/env Rscript
# 04_redeconve_visualization.R
# Redeconve 可视化功能测试
# 测试：cell.occur, spatial.cell.number, spatial.piechart, spot.pie

library(Redeconve); library(Seurat); library(Matrix); library(ggplot2)

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR  <- "your/data/directory"
SC_FILENAME    <- "scRNA.rds"
SPATIAL_SUBDIR <- "binned_outputs/square_016um"
CELL_TYPE_COL  <- "cell_type"
N_SPOTS_TARGET <- 10000
OUTPUT_DIR     <- "Redeconve_results"
# 依赖脚本 2 生成的 02_redeconve_cell_type_results.rds
# =====================================================

spatial_dir <- file.path(USER_DATA_DIR, SPATIAL_SUBDIR)
output_dir  <- OUTPUT_DIR
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)

# 步骤1: 加载反卷积结果
results <- readRDS(file.path(output_dir, "02_redeconve_cell_type_results.rds"))
nums <- results$raw_results
res_prop <- results$proportions

# 步骤2: 加载空间数据
st_seurat <- Load10X_Spatial(data.dir = spatial_dir)
if (ncol(st_seurat) >= N_SPOTS_TARGET) {
  set.seed(42)
  spots_to_keep <- sample(colnames(st_seurat), N_SPOTS_TARGET)
  st_seurat <- st_seurat[, spots_to_keep]
}
coords <- GetTissueCoordinates(st_seurat)

# 步骤3: 运行可视化函数

# 1. cell.occur - 细胞出现频率图
png(file.path(output_dir, "04_cell_occur.png"), width = 1000, height = 600, res = 150)
p <- cell.occur(nums)
print(p)
dev.off()

# 2. spatial.cell.number - 特定细胞类型的空间分布
top_celltypes <- rownames(res_prop)[order(rowMeans(res_prop), decreasing = TRUE)[1:4]]
common_spots <- intersect(rownames(coords), colnames(nums))
if (length(common_spots) > 0) {
  coords_subset <- coords[common_spots, ]
  nums_subset <- nums[, common_spots]
  coords_subset$x <- as.numeric(coords_subset$x)
  coords_subset$y <- as.numeric(coords_subset$y)

  png(file.path(output_dir, "04_spatial_cell_number.png"), width = 1200, height = 1000, res = 150)
  par(mfrow = c(2, 2))
  for (ct in top_celltypes) {
    if (ct %in% rownames(nums_subset)) {
      values <- as.numeric(nums_subset[ct, ])
      colors <- colorRampPalette(c("white", "red"))(100)[cut(values, breaks = 100, include.lowest = TRUE)]
      plot(coords_subset$x, coords_subset$y, col = colors, pch = 20, cex = 1.5,
           main = paste(ct, "\n(mean:", sprintf("%.3f", mean(values)), ")"),
           xlab = "X", ylab = "Y")
    }
  }
  dev.off()
}

# 3. spatial.piechart - 空间饼图
subset_spots <- sample(colnames(res_prop), min(100, ncol(res_prop)))
nums_subset <- nums[, subset_spots]
coords_subset <- coords[subset_spots, ]
p <- spatial.piechart(nums = nums_subset, coords = coords_subset, title = "Cell Type Composition")
png(file.path(output_dir, "04_spatial_piechart.png"), width = 1200, height = 1000, res = 150)
print(p)
dev.off()

# 4. spot.pie - 单个 spot 的饼图
spot_sums <- colSums(nums)
top_spot <- names(sort(spot_sums, decreasing = TRUE))[1]
p <- spot.pie(num = nums[, top_spot], title = paste("Spot:", top_spot))
png(file.path(output_dir, "04_spot_pie.png"), width = 800, height = 600, res = 150)
print(p)
dev.off()

# 5. show.cellsofinterest
if (exists("occurred.cells")) {
  p <- show.cellsofinterest(nums = nums, occurred.cells = occurred.cells, outliernum = 10)
  png(file.path(output_dir, "04_cells_of_interest.png"), width = 1000, height = 600, res = 150)
  print(p)
  dev.off()
}

# 汇总报告
summary_text <- c(
  "Redeconve 可视化功能测试报告",
  "================================", "",
  "测试的可视化函数:",
  "1. cell.occur() - 细胞出现频率图",
  "2. spatial.cell.number() - 特定细胞类型的空间分布",
  "3. spatial.piechart() - 空间饼图",
  "4. spot.pie() - 单个 spot 的饼图",
  "5. show.cellsofinterest() - 显示感兴趣的细胞",
  "", paste("测试时间:", Sys.time()), ""
)
writeLines(summary_text, file.path(output_dir, "04_visualization_summary.txt"))

print(paste("输出目录:", output_dir))

脚本 5 — Redeconve 共定位分析

#!/usr/bin/env Rscript
# 05_redeconve_colocalization.R
# Redeconve 共定位分析测试
# 测试：coloc.network, extract.subgraph

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR  <- "your/data/directory"
SC_FILENAME    <- "scRNA.rds"
SPATIAL_SUBDIR <- "binned_outputs/square_016um"
CELL_TYPE_COL  <- "cell_type"
N_SPOTS_TARGET <- 10000
OUTPUT_DIR     <- "Redeconve_results"
# 依赖脚本 2 生成的 02_redeconve_cell_type_results.rds
# =====================================================

library(Redeconve); library(Seurat); library(Matrix); library(igraph); library(ggplot2)

sc_data_path <- file.path(USER_DATA_DIR, SC_FILENAME)
spatial_dir  <- file.path(USER_DATA_DIR, SPATIAL_SUBDIR)
output_dir   <- OUTPUT_DIR
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)

# 步骤1: 加载反卷积结果
results <- readRDS(file.path(output_dir, "02_redeconve_cell_type_results.rds"))
nums <- results$raw_results
res_prop <- results$proportions

# 步骤2: 计算细胞类型相关性
prop_t <- t(as.matrix(res_prop))
corr_matrix <- cor(prop_t, method = "pearson")
corr_matrix[is.na(corr_matrix)] <- 0

# 步骤3: 构建共定位网络
cell_types <- rownames(res_prop)
n_types <- length(cell_types)
threshold <- 0.3

g <- coloc.network(
  corr = corr_matrix,
  thresh = threshold,
  cell.type = TRUE,
  annotations = cell_types,
  ntypes = n_types
)
print(paste("    ✓ 网络节点数:", vcount(g), ", 边数:", ecount(g)))

# 保存网络信息
network_info <- data.frame(
  Node = V(g)$name,
  CellType = V(g)$annotations,
  Color = V(g)$color
)
write.csv(network_info, file.path(output_dir, "05_network_nodes.csv"), row.names = FALSE)

if (ecount(g) > 0) {
  edges <- as.data.frame(get.edgelist(g))
  colnames(edges) <- c("From", "To")
  edges$Weight <- E(g)$weight
  write.csv(edges, file.path(output_dir, "05_network_edges.csv"), row.names = FALSE)
}

# 绘制共定位网络图
if (ecount(g) > 0) {
  png(file.path(output_dir, "05_coloc_network.png"), width = 1200, height = 1000, res = 150)
  layout <- layout_in_circle(g)
  plot(g, layout = layout,
       vertex.size = 30, vertex.label.cex = 0.8, vertex.label.color = "black",
       vertex.color = V(g)$color, edge.width = E(g)$weight * 5,
       edge.label = round(E(g)$weight, 2), edge.label.cex = 0.7, edge.label.color = "red",
       main = "Cell Type Co-localization Network",
       sub = paste("Threshold:", threshold))
  dev.off()
}

# 提取子网络（hub节点）
if (ecount(g) > 0) {
  degrees <- degree(g)
  hub_node <- names(which.max(degrees))

  neighbors_nodes <- neighbors(g, hub_node)
  subg_nodes <- c(hub_node, names(neighbors_nodes))
  subg <- induced_subgraph(g, subg_nodes)

  png(file.path(output_dir, "05_subgraph.png"), width = 800, height = 600, res = 150)
  subg_layout <- layout_in_circle(subg)
  vertex_colors <- ifelse(V(subg)$name == hub_node, "red", V(subg)$color)
  vertex_sizes <- ifelse(V(subg)$name == hub_node, 40, 30)
  plot(subg, layout = subg_layout,
       vertex.size = vertex_sizes, vertex.label.cex = 0.8,
       vertex.color = vertex_colors, edge.width = E(subg)$weight * 5,
       main = paste("Subgraph centered at", hub_node))
  dev.off()
}

# 相关性热图
library(reshape2)
corr_df <- melt(corr_matrix)
colnames(corr_df) <- c("CellType1", "CellType2", "Correlation")
p <- ggplot(corr_df, aes(x = CellType1, y = CellType2, fill = Correlation)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white",
                       midpoint = 0, limit = c(-1, 1)) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
        axis.text.y = element_text(size = 8), axis.title = element_blank()) +
  labs(title = "Cell Type Correlation Heatmap")

png(file.path(output_dir, "05_correlation_heatmap.png"), width = 1000, height = 900, res = 150)
print(p)
dev.off()

# 网络统计
network_stats <- data.frame(
  Metric = c("Number of Nodes", "Number of Edges", "Density", "Average Degree", "Clustering Coefficient"),
  Value = c(
    vcount(g), ecount(g),
    round(edge_density(g), 4),
    round(mean(degree(g)), 2),
    if (vcount(g) > 2) round(transitivity(g), 4) else NA
  )
)
write.csv(network_stats, file.path(output_dir, "05_network_statistics.csv"), row.names = FALSE)

# 汇总报告
summary_text <- c(
  "Redeconve 共定位分析测试报告",
  "================================", "",
  "测试的共定位函数:",
  "1. coloc.network() - 构建共定位网络",
  "2. extract.subgraph() - 提取子网络",
  "", "分析参数:",
  paste("  - 相关系数阈值:", threshold),
  paste("  - 细胞类型数:", nrow(nums)),
  "",
  paste("  - 节点数:", vcount(g)),
  paste("  - 边数:", ecount(g)),
  paste("  - 网络密度:", round(edge_density(g), 4)),
  "", paste("测试时间:", Sys.time()), ""
)
writeLines(summary_text, file.path(output_dir, "05_colocalization_summary.txt"))

print(paste("输出目录:", output_dir))

脚本 6 — Redeconve 高级分析

#!/usr/bin/env Rscript
# 06_redeconve_advanced_analysis.R
# Redeconve 高级分析测试
# 测试：gene.filter, resids, hNpp, cell.sampling, extract.seurat.sc/st, spatial.gene, deg.dotplot

library(Redeconve); library(Seurat); library(Matrix); library(ggplot2)

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR  <- "your/data/directory"
SC_FILENAME    <- "scRNA.rds"
SPATIAL_SUBDIR <- "binned_outputs/square_016um"
CELL_TYPE_COL  <- "cell_type"
N_CELLS_TARGET <- 10000
N_SPOTS_TARGET <- 10000
OUTPUT_DIR     <- "Redeconve_results"
# =====================================================

sc_data_path <- file.path(USER_DATA_DIR, SC_FILENAME)
spatial_dir  <- file.path(USER_DATA_DIR, SPATIAL_SUBDIR)
output_dir   <- OUTPUT_DIR
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)

# 步骤1: 加载数据
sc_seurat <- readRDS(sc_data_path)

# 自动检测细胞类型列名
resolve_cell_type_col <- function(seurat_obj, default_col) {
  candidates <- c(default_col, "cell_type", "Clusters", "celltype", "cell.type", "annot")
  for (c in candidates) {
    if (c %in% colnames(seurat_obj@meta.data)) return(c)
  }
  stop("无法找到细胞类型列。请修改 CELL_TYPE_COL 参数。")
}
CELL_TYPE_COL <- resolve_cell_type_col(sc_seurat, CELL_TYPE_COL)
print(paste("  使用细胞类型列:", CELL_TYPE_COL))

if (N_CELLS_TARGET > 0 && ncol(sc_seurat) >= N_CELLS_TARGET) {
  set.seed(42)
  cells_to_keep <- sample(colnames(sc_seurat), N_CELLS_TARGET)
  sc_seurat <- sc_seurat[, cells_to_keep]
}
sc <- GetAssayData(sc_seurat, layer = "counts")

st_seurat <- Load10X_Spatial(data.dir = spatial_dir)
if (N_SPOTS_TARGET > 0 && ncol(st_seurat) >= N_SPOTS_TARGET) {
  set.seed(42)
  spots_to_keep <- sample(colnames(st_seurat), N_SPOTS_TARGET)
  st_seurat <- st_seurat[, spots_to_keep]
}
st <- GetAssayData(st_seurat, layer = "counts")

common_genes <- intersect(rownames(sc), rownames(st))
sc <- sc[common_genes, ]
st <- st[common_genes, ]

# 步骤2: 高级分析功能测试
annotations <- cbind(colnames(sc_seurat), as.character(sc_seurat[[CELL_TYPE_COL, drop = TRUE]]))
colnames(annotations) <- c("cell", "cell_type")
ref <- get.ref(sc = sc, annotations = annotations, dopar = FALSE)

# 1. gene.filter
filtered_genes <- gene.filter(ref = ref, st = st, var_thresh = 0.025, exp_thresh = 0.003)
print(paste("    ✓ 过滤后基因数:", length(filtered_genes)))
write.csv(data.frame(Gene = filtered_genes),
          file.path(output_dir, "06_filtered_genes.csv"), row.names = FALSE)

# 加载反卷积结果
results <- readRDS(file.path(output_dir, "02_redeconve_cell_type_results.rds"))
nums <- results$raw_results
res_prop <- results$proportions

# 2. resids - 残差
spot_resids <- resids(nums = nums, ref = ref, st = st[, colnames(nums)],
                     mode = "spot", inside = "out")
write.csv(data.frame(Spot = names(spot_resids), Residual = spot_resids),
          file.path(output_dir, "06_spot_residuals.csv"), row.names = FALSE)

# 残差分布
resids_df <- data.frame(Spot = names(spot_resids), Residual = spot_resids)
p <- ggplot(resids_df, aes(x = Residual)) +
  geom_histogram(bins = 30, fill = "steelblue", alpha = 0.7) +
  theme_minimal() +
  labs(title = "Distribution of Spot-wise Residuals", x = "Residual", y = "Count")
png(file.path(output_dir, "06_residuals_distribution.png"), width = 800, height = 600, res = 150)
print(p); dev.off()

# 3. hNpp - 熵和困惑度
entropy_perplexity <- hNpp(nums = nums, thre = 1e-10)
hNpp_df <- data.frame(
  Spot = colnames(nums),
  Entropy = entropy_perplexity[1, ],
  Perplexity = entropy_perplexity[2, ]
)
write.csv(hNpp_df, file.path(output_dir, "06_entropy_perplexity.csv"), row.names = FALSE)

# 4. cell.sampling
sampled_cells <- cell.sampling(ncells = ncol(sc), annotations = annotations,
                              size = min(1000, ncol(sc)), prot = TRUE)
write.csv(sampled_cells, file.path(output_dir, "06_sampled_cells.csv"), row.names = FALSE)

# 5. Seurat接口演示
sc_list <- extract.seurat.sc(sc_seurat)
st_list <- extract.seurat.st(st_seurat)

# 6. spatial.gene
top_genes <- rownames(sc)[order(rowMeans(sc), decreasing = TRUE)[1:4]]
coords <- GetTissueCoordinates(st_seurat)
coords$x <- as.numeric(coords$x)
coords$y <- as.numeric(coords$y)

png(file.path(output_dir, "06_spatial_gene_expression.png"), width = 1200, height = 1000, res = 150)
par(mfrow = c(2, 2))
for (gene in top_genes) {
  if (gene %in% rownames(st)) {
    expr_values <- as.numeric(st[gene, rownames(coords)])
    colors <- colorRampPalette(c("white", "blue"))(100)[cut(expr_values, breaks = 100, include.lowest = TRUE)]
    plot(coords$x, coords$y, col = colors, pch = 20, cex = 1.5,
         main = paste(gene, "\n(mean:", sprintf("%.2f", mean(expr_values)), ")"),
         xlab = "X", ylab = "Y")
  }
}
dev.off()

# 7. deg.dotplot
cell_types <- unique(annotations[, 2])
marker_genes <- c()
for (ct in cell_types) {
  ct_cells <- annotations[annotations[, 2] == ct, 1]
  if (length(ct_cells) > 0) {
    ct_mean <- rowMeans(sc[, ct_cells, drop = FALSE])
    top_genes <- names(sort(ct_mean, decreasing = TRUE))[1:3]
    marker_genes <- c(marker_genes, top_genes)
  }
}
marker_genes <- unique(marker_genes)

if (length(marker_genes) > 0) {
  p <- deg.dotplot(sc = sc, annotations = annotations, genelist = marker_genes)
  png(file.path(output_dir, "06_deg_dotplot.png"), width = 1200, height = 800, res = 150)
  print(p); dev.off()
}

# 汇总报告
summary_text <- c(
  "Redeconve 高级分析测试报告",
  "================================", "",
  "测试的高级分析函数:",
  "1. gene.filter() - 基因过滤",
  "2. resids() - 计算残差",
  "3. hNpp() - 计算熵和困惑度",
  "4. cell.sampling() - 细胞采样",
  "5. extract.seurat.sc/st() - Seurat接口",
  "6. spatial.gene() - 基因空间表达",
  "7. deg.dotplot() - 差异表达基因点图",
  "", paste("测试时间:", Sys.time()), ""
)
writeLines(summary_text, file.path(output_dir, "06_advanced_analysis_summary.txt"))

print(paste("输出目录:", output_dir))

推荐工作流

> 脚本已全部内嵌到本 SKILL.md 中，请按需复制对应代码块到本地 .R 文件运行。每个脚本顶部都有配置区域，按需修改即可。

按使用场景推荐：

典型使用顺序：

1. 脚本 2（细胞类型反卷积）→ 生成 02_redeconve_cell_type_results.rds
2. 脚本 4 / 5 读取该结果文件做可视化和共定位分析
3. 脚本 6 做独立的高级分析（残差、熵等）

输出文件清单

数据文件

• 01_redeconve_single_cell_results.rds：单细胞级别结果
• 01_redeconve_single_cell_proportions.csv：单细胞比例矩阵
• 02_redeconve_cell_type_results.rds：细胞类型级别结果
• 02_redeconve_cell_type_proportions.csv：细胞类型比例矩阵
• 03_redeconve_seurat_extracted.rds：Seurat 提取的数据
• 05_network_nodes.csv / 05_network_edges.csv：网络节点/边
• 06_filtered_genes.csv：过滤后基因列表
• 06_spot_residuals.csv：Spot 残差
• 06_entropy_perplexity.csv：熵和困惑度
• 06_sampled_cells.csv：采样细胞

图形文件

• 02_redeconve_celltype_spatial.png：前 4 种细胞类型空间分布
• 02_redeconve_celltype_barplot.png：细胞类型比例条形图
• 04_cell_occur.png：细胞出现频率图
• 04_spatial_cell_number.png：特定细胞类型空间分布
• 04_spatial_piechart.png：空间饼图
• 04_spot_pie.png：单 spot 饼图
• 05_coloc_network.png：共定位网络
• 05_subgraph.png：子网络
• 05_correlation_heatmap.png：相关性热图
• 06_residuals_distribution.png：残差分布
• 06_spatial_gene_expression.png：基因空间表达
• 06_deg_dotplot.png：DEG 点图

报告

• 01/02_*_statistics.csv：统计摘要
• 04_visualization_summary.txt：可视化汇总
• 05_colocalization_summary.txt：共定位汇总
• 06_advanced_analysis_summary.txt：高级分析汇总

故障排除

反卷积失败

解决：
1. 检查单细胞/空间表达矩阵的格式（dgCMatrix 或 matrix）
2. 确保共同基因 ≥ 100
3. 检查基因名是否一致

extract.seurat.sc/st 失败

解决：
- 使用手动提取作为替代（脚本 3 中已提供 fallback）

内存不足

解决：
- 单细胞级别限制 ≤ 1000 cells
- 空间采样到 10000 spots
- dopar = FALSE 关闭并行

共定位网络为空

解决：
- 降低阈值（默认 0.3，可改为 0.2）
- 检查细胞类型比例（去掉过低比例的类型）

资源链接

• Redeconve 源码：github.com/ZxZhou4150/Redeconve
• 官方文档：参考 GitHub README

许可证

遵循原作者的开源许可证（详见 GitHub 仓库）。