RCTD-HDdata — RCTD 反卷积 Visium HD 多分辨率数据

概述

RCTD（Robust Cell Type Decomposition）是基于最大似然估计的空间转录组反卷积方法（来自 spacexr 包），将单细胞参考数据与空间表达数据整合，估计每个 spot 的细胞类型组成。本 Skill 提供对 Visium HD 多种分辨率（8µm、16µm、cellbin）的完整 RCTD 分析流程。

核心特性：

• 多分辨率支持：Visium HD 8µm / 16µm / cellbin（单细胞分辨率）
• 真实数据示例：scRNA-seq 参考 + Visium HD（示例：胰腺癌数据，可替换为任意组织）
• 完整可视化：分布图、箱线图、单细胞类型热图、综合图、气泡图、对比图
• Cellbin 真实坐标提取：从 cell_segmentations.geojson 解析细胞中心
• 统计摘要：自动生成各细胞类型平均比例与主导类型分布

---

适用场景

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

• 使用 RCTD 对 Visium HD 数据进行反卷积
• 在 8µm / 16µm / cellbin 多种分辨率间比较细胞类型组成
• 处理 cellbin 单细胞分辨率数据
• 从 cell_segmentations.geojson 提取真实空间坐标
• 生成 RCTD 高级可视化（热图、气泡图、对比图）

---

安装

环境要求

• R: >= 4.0.0
• 操作系统: Linux (推荐), Windows, macOS

安装 RCTD R 包

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("RCTD")

安装依赖包

install.packages(c("ggplot2", "dplyr", "Matrix"))

验证安装

library(RCTD)
print("RCTD 安装成功！")

常见安装问题

| 问题 | 解决方案 | |---|---| | Bioconductor 安装失败 | 更新 BiocManager: BiocManager::install(version = "3.18") | | 依赖包冲突 | 使用 clean R 环境，逐个安装依赖 |

---

RCTD 反卷积原理

RCTD 通过最大似然估计将空间 spot 的混合表达分解为各细胞类型的贡献，使用：

• Reference：单细胞表达矩阵 + 细胞类型标签
• SpatialRNA (Puck)：空间表达矩阵 + 坐标
• doublet_mode：
  • "full"：估计所有细胞类型比例（推荐）
  • "doublet"：仅检测双细胞

核心步骤：

1. 构建 Reference 对象（单细胞 + 细胞类型）
2. 构建 SpatialRNA 对象（空间表达 + 坐标）
3. create.RCTD() 创建 RCTD 对象
4. run.RCTD() 运行反卷积
5. normalize_weights() 归一化权重
6. 可视化与下游分析

---

Visium HD 三种分辨率对比

| 分辨率 | 用途 | 内存 | 速度 | |---|---|---|---| | 8µm | 高精度空间分布 | 高 | 慢 | | 16µm | 平衡精度与性能 | 中 | 中 | | Cellbin | 单细胞分辨率 | 极高（≥100GB） | 很慢 |

---

完整可运行脚本

输入数据要求

| 数据 | 格式 | 说明 | |---|---|---| | scRNA-seq 表达矩阵 | Seurat 对象 / dgCMatrix | 行=基因，列=细胞；meta.data 中需包含细胞类型注释列 | | 空间表达数据（16µm/8µm） | 10X h5 格式 | binned_outputs/square_016um/filtered_feature_bc_matrix.h5 | | 空间表达数据（Cellbin） | 10X h5 格式 | segmented_outputs/filtered_feature_cell_matrix.h5 | | Cellbin 真实坐标（可选） | geojson | segmented_outputs/cell_segmentations.geojson |

配置说明

所有脚本的前 15 行为可配置区域，包含以下参数（统一约定）：

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR   <- "your/data/directory"        # 数据根目录
CELL_TYPE_COL   <- "Clusters"                    # scRNA-seq metadata 中细胞类型列名
N_SPOTS_SAMPLE  <- 10000                         # 空间采样数（0 表示不采样）
MAX_CORES       <- 4                             # 并行核心数
OUTPUT_DIR      <- "RCTD_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"

注意：以下脚本均假设单细胞注释列名为 Clusters，若你的数据为 cell_type、celltype 等，请同步修改 CELL_TYPE_COL。

脚本 1 — RCTD 16µm 全流程

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# rctd_16um_workflow.R
# RCTD 空间反卷积分析 - 16µm 分辨率
# 使用单细胞参考数据对空间转录组数据进行细胞类型反卷积

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR   <- "your/data/directory"          # TODO: 修改为你的数据根目录
SPATIAL_SUBDIR  <- "binned_outputs/square_016um"  # 空间数据子目录（相对 USER_DATA_DIR）
CELL_TYPE_COL   <- "Clusters"                     # TODO: 修改为你的细胞类型注释列名
N_SPOTS_SAMPLE  <- 10000                          # 空间采样数（0=不采样）
MAX_CORES       <- 4                              # 并行核心数
OUTPUT_DIR      <- "RCTD_16um_results"            # 输出目录
# =====================================================

# 加载必要的包
suppressPackageStartupMessages({
  library(spacexr)
  library(Seurat)
  library(SeuratDisk)
  library(ggplot2)
  library(dplyr)
})

# 路径组装
sc_file        <- file.path(USER_DATA_DIR, "scRNA.rds")  # 单细胞 RDS
spatial_h5_dir <- file.path(USER_DATA_DIR, SPATIAL_SUBDIR)

cat("=== RCTD 16µm 空间反卷积分析 ===\n\n")

# 1. 读取单细胞参考数据
cat("1. 读取单细胞参考数据...\n")
cat("   文件路径:", sc_file, "\n")

sc_data <- readRDS(sc_file)
cat("   ✓ 单细胞数据读取成功\n")
cat("   - 细胞数量:", ncol(sc_data), "\n")
cat("   - 基因数量:", nrow(sc_data), "\n")

# 检查细胞类型注释
cat("\n2. 检查细胞类型注释...\n")
cat("   使用注释列:", CELL_TYPE_COL, "\n")

cell_type_table <- table(sc_data@meta.data[[CELL_TYPE_COL]])
cat("   细胞类型分布:\n")
for (ct in names(cell_type_table)) {
  cat(sprintf("     - %s: %d 个细胞\n", ct, cell_type_table[ct]))
}

# 3. 准备RCTD参考数据
cat("\n3. 准备RCTD参考数据...\n")
raw_counts <- GetAssayData(sc_data, layer = "counts")
cell_types <- as.factor(sc_data@meta.data[[CELL_TYPE_COL]])
names(cell_types) <- colnames(sc_data)

reference <- Reference(
  counts = raw_counts,
  cell_types = cell_types
)

cat("   ✓ Reference对象创建成功\n")
cat("   - 参考细胞类型数:", length(unique(cell_types)), "\n")

# 4. 读取空间转录组数据
cat("\n4. 读取空间转录组数据...\n")

spatial_data <- Load10X_Spatial(
  data.dir = spatial_h5_dir,
  filename = "filtered_feature_bc_matrix.h5"
)

cat("   ✓ 空间数据读取成功\n")
cat("   - Spot数量:", ncol(spatial_data), "\n")
cat("   - 基因数量:", nrow(spatial_data), "\n")

# 5. 准备RCTD空间数据 (Puck)
cat("\n5. 准备RCTD空间数据...\n")
spatial_counts <- GetAssayData(spatial_data, layer = "counts")
spatial_coords <- GetTissueCoordinates(spatial_data, scale = NULL)

if (ncol(spatial_coords) > 2) {
  spatial_coords <- spatial_coords[, c("x", "y")]
}

puck <- SpatialRNA(
  coords = spatial_coords,
  counts = spatial_counts
)

cat("   ✓ SpatialRNA (Puck) 对象创建成功\n")
cat("   - Spot数量:", ncol(puck@counts), "\n")

# 6. 创建RCTD对象
cat("\n6. 创建RCTD对象...\n")
my_rctd <- create.RCTD(
  spatialRNA = puck,
  reference = reference,
  max_cores = MAX_CORES
)

cat("   ✓ RCTD对象创建成功\n")

# 7. 运行RCTD反卷积
cat("\n7. 运行RCTD反卷积分析...\n")
cat("   这可能需要一些时间...\n")

my_rctd <- run.RCTD(
  my_rctd,
  doublet_mode = "full"
)

cat("   ✓ RCTD反卷积完成\n")

# 8. 提取反卷积结果
cat("\n8. 提取反卷积结果...\n")
results <- my_rctd@results
weights <- results$weights
norm_weights <- normalize_weights(weights)

weights_df <- as.data.frame(norm_weights)
weights_df$spot_id <- rownames(weights_df)

cat("   ✓ 细胞类型比例提取完成\n")
cat("   - Spot数量:", nrow(weights_df), "\n")
cat("   - 细胞类型数:", ncol(norm_weights), "\n")

# 9. 保存结果
cat("\n9. 保存反卷积结果...\n")
if (!dir.exists(OUTPUT_DIR)) dir.create(OUTPUT_DIR, recursive = TRUE)

output_csv <- file.path(OUTPUT_DIR, "rctd_cell_type_proportions_16um.csv")
write.csv(weights_df, file = output_csv, row.names = FALSE)
cat("   ✓ 细胞类型比例已保存:", output_csv, "\n")

output_rdata <- file.path(OUTPUT_DIR, "rctd_results_16um.RData")
save(my_rctd, norm_weights, weights_df, file = output_rdata)
cat("   ✓ RCTD结果已保存:", output_rdata, "\n")

# 10. 可视化
cat("\n10. 生成可视化...\n")
dominant_cell_type <- colnames(norm_weights)[apply(norm_weights, 1, which.max)]
weights_df$dominant_cell_type <- dominant_cell_type

p1 <- ggplot(weights_df, aes(x = dominant_cell_type, fill = dominant_cell_type)) +
  geom_bar() +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "RCTD反卷积: 主要细胞类型分布",
       subtitle = "空间转录组数据 (16µm分辨率)",
       x = "主要细胞类型", y = "Spot数量") +
  guides(fill = "none")

ggsave(file.path(OUTPUT_DIR, "rctd_dominant_cell_types_16um.png"), p1, width = 10, height = 6, dpi = 150)

weights_long <- reshape2::melt(weights_df[, c("spot_id", colnames(norm_weights))],
                               id.vars = "spot_id",
                               variable.name = "cell_type",
                               value.name = "proportion")

p2 <- ggplot(weights_long, aes(x = cell_type, y = proportion, fill = cell_type)) +
  geom_boxplot() +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "RCTD反卷积: 细胞类型比例分布",
       x = "细胞类型", y = "比例") +
  guides(fill = "none")

ggsave(file.path(OUTPUT_DIR, "rctd_cell_type_proportions_boxplot_16um.png"), p2, width = 12, height = 6, dpi = 150)

# 11. 统计摘要
cat("\n11. 反卷积结果统计摘要\n")
mean_proportions <- colMeans(norm_weights)
cat("   各细胞类型平均比例:\n")
for (ct in names(mean_proportions)) {
  cat(sprintf("     - %s: %.3f\n", ct, mean_proportions[ct]))
}

dominant_stats <- table(dominant_cell_type)
cat("\n   各细胞类型作为主要类型的Spot数:\n")
for (ct in names(dominant_stats)) {
  pct <- dominant_stats[ct] / length(dominant_cell_type) * 100
  cat(sprintf("     - %s: %d spots (%.1f%%)\n", ct, dominant_stats[ct], pct))
}

# 12. 将结果添加到Seurat对象
cat("\n12. 将反卷积结果添加到Seurat对象...\n")
common_spots <- intersect(colnames(spatial_data), rownames(norm_weights))
spatial_data_subset <- spatial_data[, common_spots]

for (ct in colnames(norm_weights)) {
  spatial_data_subset@meta.data[[paste0("RCTD_", ct)]] <- norm_weights[common_spots, ct]
}
spatial_data_subset@meta.data$RCTD_dominant_type <- dominant_cell_type[common_spots]

saveRDS(spatial_data_subset, file = file.path(OUTPUT_DIR, "spatial_with_rctd_16um.rds"))

# 13. 生成空间可视化
cat("\n13. 生成空间可视化...\n")
top_cell_types <- names(sort(mean_proportions, decreasing = TRUE))[1:min(4, length(mean_proportions))]

p_list <- list()
for (i in seq_along(top_cell_types)) {
  ct <- top_cell_types[i]
  p <- SpatialFeaturePlot(spatial_data_subset,
    features = paste0("RCTD_", ct),
    pt.size.factor = 1.5, crop = FALSE) +
    scale_fill_gradientn(colors = c("grey90", "yellow", "red", "darkred"),
                         name = "Proportion") +
    ggtitle(paste0(ct, " (RCTD)"))
  p_list[[i]] <- p
}

if (length(p_list) >= 4) {
  combined_plot <- cowplot::plot_grid(
    p_list[[1]], p_list[[2]], p_list[[3]], p_list[[4]],
    ncol = 2, labels = c("A", "B", "C", "D")
  )
  ggsave(file.path(OUTPUT_DIR, "rctd_spatial_cell_types_16um.png"),
         combined_plot, width = 14, height = 12, dpi = 150)
}

p_dominant <- SpatialDimPlot(spatial_data_subset,
  group.by = "RCTD_dominant_type",
  pt.size.factor = 1.5, crop = FALSE, label = FALSE) +
  ggtitle("RCTD: 主要细胞类型空间分布")

ggsave(file.path(OUTPUT_DIR, "rctd_spatial_dominant_types_16um.png"),
       p_dominant, width = 10, height = 8, dpi = 150)

cat("\n=== RCTD反卷积分析完成 ===\n")

---

脚本 2 — RCTD 8µm 分辨率

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# rctd_deconvolution_8um.R
# RCTD 空间反卷积分析 - 8µm 分辨率

# 设置工作目录
setwd("OUTPUT_DIR")

suppressPackageStartupMessages({
  library(spacexr)
  library(Seurat)
  library(SeuratDisk)
  library(ggplot2)
  library(dplyr)
  library(reshape2)
  library(cowplot)
})

cat("=== RCTD 8µm 空间反卷积分析 ===\n\n")

# 1. 读取单细胞参考数据
sc_file <- "file.path(USER_DATA_DIR, "scRNA.rds")"
sc_data <- readRDS(sc_file)
cat("   ✓ 单细胞数据读取成功:", ncol(sc_data), "cells\n")

# 2. 采样单细胞数据以加快分析
cat("\n2. 采样单细胞数据 (每类型最多5000个细胞)...\n")
cell_type_col <- "Clusters"
cell_types <- sc_data@meta.data[[cell_type_col]]
set.seed(42)
selected_cells <- c()
for (ct in unique(cell_types)) {
  ct_cells <- colnames(sc_data)[cell_types == ct]
  if (length(ct_cells) > 5000) {
    ct_cells <- sample(ct_cells, 5000)
  }
  selected_cells <- c(selected_cells, ct_cells)
}
sc_data_subset <- sc_data[, selected_cells]

# 3. 准备RCTD参考数据
raw_counts <- GetAssayData(sc_data_subset, layer = "counts")
cell_types_factor <- as.factor(sc_data_subset@meta.data[[cell_type_col]])
names(cell_types_factor) <- colnames(sc_data_subset)

reference <- Reference(counts = raw_counts, cell_types = cell_types_factor)

# 4. 读取8um空间转录组数据
spatial_dir <- "USER_DATA_DIR"
spatial_data <- Load10X_Spatial(
  data.dir = file.path(spatial_dir, "binned_outputs/square_008um"),
  filename = "filtered_feature_bc_matrix.h5"
)

# 采样spots以加快分析
set.seed(123)
if (ncol(spatial_data) > 10000) {
  selected_spots <- sample(colnames(spatial_data), 10000)
  spatial_data <- spatial_data[, selected_spots]
}

# 5. 准备RCTD空间数据 (Puck)
spatial_counts <- GetAssayData(spatial_data, layer = "counts")
spatial_coords <- GetTissueCoordinates(spatial_data, scale = NULL)
if (ncol(spatial_coords) > 2) {
  spatial_coords <- spatial_coords[, c("x", "y")]
}

puck <- SpatialRNA(coords = spatial_coords, counts = spatial_counts)

# 6. 创建RCTD对象
my_rctd <- create.RCTD(spatialRNA = puck, reference = reference, max_cores = 4)

# 7. 运行RCTD反卷积
cat("\n8. 运行RCTD反卷积分析...\n")
my_rctd <- run.RCTD(my_rctd, doublet_mode = "full")

# 9. 提取反卷积结果
results <- my_rctd@results
weights <- results$weights
norm_weights <- normalize_weights(weights)
weights_df <- as.data.frame(norm_weights)
weights_df$spot_id <- rownames(weights_df)

# 10. 保存结果
output_csv <- "rctd_cell_type_proportions_8um.csv"
write.csv(weights_df, file = output_csv, row.names = FALSE)

save(my_rctd, norm_weights, weights_df, file = "rctd_results_8um.RData")

# 11. 基础可视化
dominant_cell_type <- colnames(norm_weights)[apply(norm_weights, 1, which.max)]
weights_df$dominant_cell_type <- dominant_cell_type

p1 <- ggplot(weights_df, aes(x = dominant_cell_type, fill = dominant_cell_type)) +
  geom_bar() + theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "RCTD反卷积: 主要细胞类型分布 (8um)",
       subtitle = "空间转录组数据 (8µm分辨率)",
       x = "主要细胞类型", y = "Spot数量") +
  guides(fill = "none")
ggsave("rctd_dominant_cell_types_8um.png", p1, width = 10, height = 6, dpi = 150)

# 12. 添加到Seurat对象并保存
common_spots <- intersect(colnames(spatial_data), rownames(norm_weights))
spatial_data_subset <- spatial_data[, common_spots]
for (ct in colnames(norm_weights)) {
  spatial_data_subset@meta.data[[paste0("RCTD_", ct)]] <- norm_weights[common_spots, ct]
}
spatial_data_subset@meta.data$RCTD_dominant_type <- dominant_cell_type[common_spots]
saveRDS(spatial_data_subset, file = "spatial_pancreas_with_rctd_8um.rds")

# 13. 高级可视化（无HE背景的热图）
cat("\n14. 生成高级可视化（这次的4类图）...\n")
coords <- GetTissueCoordinates(spatial_data_subset, scale = NULL)
if (ncol(coords) > 2) coords <- coords[, c("x", "y")]
weights_subset <- weights_df[match(rownames(coords), weights_df$spot_id), ]
cell_type_cols <- setdiff(colnames(weights_df), c("spot_id", "dominant_cell_type"))

# 单个细胞类型热图（无HE背景）
for (ct in cell_type_cols) {
  plot_data <- data.frame(
    x = coords$x, y = coords$y,
    proportion = weights_subset[[ct]]
  )
  p <- ggplot(plot_data, aes(x = x, y = y, color = proportion)) +
    geom_point(size = 1.5, alpha = 0.8) +
    scale_color_gradientn(colors = c("grey90", "yellow", "orange", "red", "darkred"),
                         name = "Proportion",
                         limits = c(0, max(plot_data$proportion))) +
    coord_fixed() + theme_void() +
    theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
          legend.position = "right") +
    ggtitle(paste0(ct, " (RCTD 8um)"))
  output_file <- paste0("rctd_heatmap_", gsub(" ", "_", gsub("\\.", "_", ct)), "_no_bg_8um.png")
  ggsave(output_file, p, width = 8, height = 7, dpi = 150, bg = "white")
}

# 14.2 所有细胞类型综合图
dominant_types <- apply(weights_subset[, cell_type_cols], 1, function(x) cell_type_cols[which.max(x)])
combined_data <- data.frame(x = coords$x, y = coords$y, dominant_type = dominant_types)

cell_type_colors <- c(
  "ACINAR" = "#E41A1C", "B CELLS" = "#377EB8", "CYCLING DUCTAL" = "#4DAF4A",
  "CYCLING TNK" = "#984EA3", "CYCLING. MYELOID" = "#FF7F00", "DUCTAL" = "#FFFF33",
  "ENDOCRINE" = "#A65628", "ENDOTHELIAL" = "#F781BF", "FIBROBLASTS" = "#999999",
  "MAST" = "#66C2A5", "MYELOID" = "#FC8D62", "PERICYTES" = "#8DA0CB",
  "PLASMA" = "#E78AC3", "TNK" = "#A6D854"
)

used_colors <- cell_type_colors[unique(dominant_types)]
p_combined <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type)) +
  geom_point(size = 1.5, alpha = 0.8) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  coord_fixed() + theme_void() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "right") +
  ggtitle("RCTD: All Cell Types - Dominant (8um)")
ggsave("rctd_all_cell_types_combined_8um.png", p_combined, width = 10, height = 8, dpi = 150, bg = "white")

# 14.3 气泡图
combined_data$max_proportion <- apply(weights_subset[, cell_type_cols], 1, max)
p_bubble <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type, size = max_proportion)) +
  geom_point(alpha = 0.6) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  scale_size_continuous(name = "Proportion", range = c(0.5, 4)) +
  coord_fixed() + theme_void() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "right") +
  ggtitle("RCTD: All Cell Types - Size = Proportion (8um)")
ggsave("rctd_all_cell_types_bubble_8um.png", p_bubble, width = 10, height = 8, dpi = 150, bg = "white")

# 14.4 前4种细胞类型对比图
mean_proportions <- colMeans(norm_weights)
top_cell_types <- names(sort(mean_proportions, decreasing = TRUE))[1:min(4, length(mean_proportions))]
comparison_data <- NULL
for (ct in top_cell_types) {
  temp_data <- data.frame(x = coords$x, y = coords$y,
                         proportion = weights_subset[[ct]], cell_type = ct)
  comparison_data <- rbind(comparison_data, temp_data)
}
p_comparison <- ggplot(comparison_data, aes(x = x, y = y, color = proportion)) +
  geom_point(size = 1.2, alpha = 0.8) +
  scale_color_gradientn(colors = c("grey90", "yellow", "orange", "red", "darkred"),
                       name = "Proportion") +
  facet_wrap(~cell_type, ncol = 2) +
  coord_fixed() + theme_void() +
  theme(strip.text = element_text(size = 12, face = "bold"),
        legend.position = "right")
ggsave("rctd_top4_cell_types_comparison_8um.png", p_comparison, width = 12, height = 10, dpi = 150, bg = "white")

cat("\n=== RCTD反卷积分析完成 (8um分辨率) ===\n")

---

脚本 3 — RCTD Cellbin 单细胞分辨率

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# rctd_deconvolution_cellbin.R
# RCTD 空间反卷积分析 - Cellbin/单细胞分辨率

setwd("OUTPUT_DIR")

suppressPackageStartupMessages({
  library(spacexr); library(Seurat); library(SeuratDisk)
  library(ggplot2); library(dplyr); library(reshape2); library(cowplot); library(Matrix)
})

cat("=== RCTD Cellbin 空间反卷积分析 ===\n\n")

# 1. 读取单细胞参考数据
sc_file <- "file.path(USER_DATA_DIR, "scRNA.rds")"
sc_data <- readRDS(sc_file)
cat("   ✓ 单细胞数据读取成功:", ncol(sc_data), "cells\n")

# 2. 采样单细胞数据
set.seed(42)
cell_types <- sc_data@meta.data[["Clusters"]]
selected_cells <- c()
for (ct in unique(cell_types)) {
  ct_cells <- colnames(sc_data)[cell_types == ct]
  if (length(ct_cells) > 5000) ct_cells <- sample(ct_cells, 5000)
  selected_cells <- c(selected_cells, ct_cells)
}
sc_data_subset <- sc_data[, selected_cells]

# 3. 准备RCTD参考数据
raw_counts <- GetAssayData(sc_data_subset, layer = "counts")
cell_types_factor <- as.factor(sc_data_subset@meta.data[["Clusters"]])
names(cell_types_factor) <- colnames(sc_data_subset)
reference <- Reference(counts = raw_counts, cell_types = cell_types_factor)

# 4. 读取Cellbin空间转录组数据
spatial_dir <- "USER_DATA_DIR"
cellbin_h5_file <- file.path(spatial_dir, "segmented_outputs/filtered_feature_cell_matrix.h5")
cellbin_counts <- Read10X_h5(filename = cellbin_h5_file)
cat("   ✓ Cellbin表达矩阵读取成功:", ncol(cellbin_counts), "cells\n")

# 5. 采样细胞
set.seed(123)
if (ncol(cellbin_counts) > 10000) {
  selected_cells_spatial <- sample(colnames(cellbin_counts), 10000)
  cellbin_counts <- cellbin_counts[, selected_cells_spatial]
}

# 6. 创建模拟空间坐标（如无实际坐标文件）
n_cells <- ncol(cellbin_counts)
set.seed(456)
x_coords <- runif(n_cells, min = 0, max = 1000)
y_coords <- runif(n_cells, min = 0, max = 1000)
cellbin_coords <- data.frame(x = x_coords, y = y_coords)
rownames(cellbin_coords) <- colnames(cellbin_counts)

# 7. 准备RCTD空间数据 (Puck)
puck <- SpatialRNA(coords = cellbin_coords, counts = cellbin_counts)

# 8. 创建RCTD对象
my_rctd <- create.RCTD(spatialRNA = puck, reference = reference, max_cores = 4)

# 9. 运行RCTD反卷积
cat("\n9. 运行RCTD反卷积分析...\n")
my_rctd <- run.RCTD(my_rctd, doublet_mode = "full")

# 10. 提取反卷积结果
results <- my_rctd@results
weights <- results$weights
norm_weights <- normalize_weights(weights)
weights_df <- as.data.frame(norm_weights)
weights_df$cell_id <- rownames(weights_df)

# 11. 保存结果
write.csv(weights_df, file = "rctd_cell_type_proportions_cellbin.csv", row.names = FALSE)
save(my_rctd, norm_weights, weights_df, file = "rctd_results_cellbin.RData")

# 12. 基础可视化（与8um/16um相同的逻辑）
dominant_cell_type <- colnames(norm_weights)[apply(norm_weights, 1, which.max)]
weights_df$dominant_cell_type <- dominant_cell_type

p1 <- ggplot(weights_df, aes(x = dominant_cell_type, fill = dominant_cell_type)) +
  geom_bar() + theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "RCTD反卷积: 主要细胞类型分布 (Cellbin)",
       subtitle = "空间转录组数据 (Cellbin/单细胞分辨率)",
       x = "主要细胞类型", y = "细胞数量") +
  guides(fill = "none")
ggsave("rctd_dominant_cell_types_cellbin.png", p1, width = 10, height = 6, dpi = 150)
ggsave("rctd_dominant_cell_types_cellbin.pdf", p1, width = 10, height = 6)

# 13. 高级可视化
weights_subset <- weights_df[match(rownames(cellbin_coords), weights_df$cell_id), ]
cell_type_cols <- setdiff(colnames(weights_df), c("cell_id", "dominant_cell_type"))

# 单个细胞类型热图
for (ct in cell_type_cols) {
  plot_data <- data.frame(x = cellbin_coords$x, y = cellbin_coords$y,
                         proportion = weights_subset[[ct]])
  p <- ggplot(plot_data, aes(x = x, y = y, color = proportion)) +
    geom_point(size = 1.2, alpha = 0.8) +
    scale_color_gradientn(colors = c("grey90", "yellow", "orange", "red", "darkred"),
                         name = "Proportion",
                         limits = c(0, max(plot_data$proportion))) +
    coord_fixed() + theme_void() +
    theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
          legend.position = "right") +
    ggtitle(paste0(ct, " (RCTD Cellbin)"))
  output_file_png <- paste0("rctd_heatmap_", gsub(" ", "_", gsub("\\.", "_", ct)), "_no_bg_cellbin.png")
  output_file_pdf <- paste0("rctd_heatmap_", gsub(" ", "_", gsub("\\.", "_", ct)), "_no_bg_cellbin.pdf")
  ggsave(output_file_png, p, width = 8, height = 7, dpi = 150, bg = "white")
  ggsave(output_file_pdf, p, width = 8, height = 7)
}

# 综合图 / 气泡图 / 对比图（与8um逻辑相同）
cell_type_colors <- c(
  "ACINAR" = "#E41A1C", "B CELLS" = "#377EB8", "CYCLING DUCTAL" = "#4DAF4A",
  "CYCLING TNK" = "#984EA3", "CYCLING. MYELOID" = "#FF7F00", "DUCTAL" = "#FFFF33",
  "ENDOCRINE" = "#A65628", "ENDOTHELIAL" = "#F781BF", "FIBROBLASTS" = "#999999",
  "MAST" = "#66C2A5", "MYELOID" = "#FC8D62", "PERICYTES" = "#8DA0CB",
  "PLASMA" = "#E78AC3", "TNK" = "#A6D854"
)

dominant_types <- apply(weights_subset[, cell_type_cols], 1, function(x) cell_type_cols[which.max(x)])
combined_data <- data.frame(x = cellbin_coords$x, y = cellbin_coords$y, dominant_type = dominant_types)
used_colors <- cell_type_colors[unique(dominant_types)]

p_combined <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type)) +
  geom_point(size = 1.2, alpha = 0.8) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  coord_fixed() + theme_void() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "right") +
  ggtitle("RCTD: All Cell Types - Dominant (Cellbin)")
ggsave("rctd_all_cell_types_combined_cellbin.png", p_combined, width = 10, height = 8, dpi = 150, bg = "white")
ggsave("rctd_all_cell_types_combined_cellbin.pdf", p_combined, width = 10, height = 8)

combined_data$max_proportion <- apply(weights_subset[, cell_type_cols], 1, max)
p_bubble <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type, size = max_proportion)) +
  geom_point(alpha = 0.6) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  scale_size_continuous(name = "Proportion", range = c(0.5, 3)) +
  coord_fixed() + theme_void() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "right") +
  ggtitle("RCTD: All Cell Types - Size = Proportion (Cellbin)")
ggsave("rctd_all_cell_types_bubble_cellbin.png", p_bubble, width = 10, height = 8, dpi = 150, bg = "white")
ggsave("rctd_all_cell_types_bubble_cellbin.pdf", p_bubble, width = 10, height = 8)

cat("\n=== RCTD反卷积分析完成 (Cellbin/单细胞分辨率) ===\n")

---

脚本 4 — RCTD Visium HD 16µm 最终优化版

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# rctd_visium_hd_16um.R
# RCTD Visium HD 16um 空间反卷积分析 - 最终优化版

# 配置参数
SC_DATA_PATH <- "file.path(USER_DATA_DIR, "scRNA.rds")"
SPATIAL_DATA_DIR <- "USER_DATA_DIR/binned_outputs/square_016um"
OUTPUT_DIR <- "RCTD_16um_results"
N_SPOTS_SAMPLE <- 10000
MAX_CORES <- 4
CELL_TYPE_COL <- "Clusters"

setwd(dirname(parent.frame(2)$ofile))
if (!dir.exists(OUTPUT_DIR)) dir.create(OUTPUT_DIR, recursive = TRUE)

suppressPackageStartupMessages({
  library(spacexr); library(Seurat); library(SeuratDisk)
  library(ggplot2); library(dplyr); library(reshape2)
})

# 步骤1: 读取单细胞参考数据
if (!file.exists(SC_DATA_PATH)) stop("错误: 找不到单细胞数据文件: ", SC_DATA_PATH)
sc_data <- readRDS(SC_DATA_PATH)
cat(sprintf("   - 细胞数量: %d, 基因数量: %d\n", ncol(sc_data), nrow(sc_data)))

# 步骤2: 检查细胞类型注释
if (!(CELL_TYPE_COL %in% colnames(sc_data@meta.data))) stop("错误: 找不到细胞类型注释列: ", CELL_TYPE_COL)
cell_type_table <- table(sc_data@meta.data[[CELL_TYPE_COL]])

# 步骤3: 准备RCTD参考数据
set.seed(42)
cell_types_all <- unique(sc_data@meta.data[[CELL_TYPE_COL]])
cells_to_keep <- c()
for (ct in cell_types_all) {
  ct_cells <- WhichCells(sc_data, expression = !!sym(CELL_TYPE_COL) == ct)
  if (length(ct_cells) > 5000) ct_cells <- sample(ct_cells, 5000)
  cells_to_keep <- c(cells_to_keep, ct_cells)
}
sc_data <- subset(sc_data, cells = cells_to_keep)

raw_counts <- GetAssayData(sc_data, layer = "counts")
cell_types <- as.factor(sc_data@meta.data[[CELL_TYPE_COL]])
names(cell_types) <- colnames(sc_data)
reference <- Reference(counts = raw_counts, cell_types = cell_types)

# 步骤4: 读取空间转录组数据
spatial_data <- Load10X_Spatial(data.dir = SPATIAL_DATA_DIR, filename = "filtered_feature_bc_matrix.h5")
cat(sprintf("   - Spot数量: %d\n", ncol(spatial_data)))

# 步骤5: 数据预处理 - 采样
if (N_SPOTS_SAMPLE > 0 && ncol(spatial_data) > N_SPOTS_SAMPLE) {
  set.seed(42)
  spots_to_keep <- sample(colnames(spatial_data), N_SPOTS_SAMPLE)
  spatial_data <- subset(spatial_data, cells = spots_to_keep)
}

# 步骤6: 准备RCTD空间数据
spatial_counts <- GetAssayData(spatial_data, layer = "counts")
spatial_coords <- GetTissueCoordinates(spatial_data, scale = NULL)
if (ncol(spatial_coords) > 2) spatial_coords <- spatial_coords[, c("x", "y")]
puck <- SpatialRNA(coords = spatial_coords, counts = spatial_counts)

# 步骤7: 创建并运行RCTD
cat("\n【步骤7/9】创建RCTD对象并运行反卷积...\n")
cat("   这可能需要15-30分钟，请耐心等待...\n")
my_rctd <- create.RCTD(spatialRNA = puck, reference = reference, max_cores = MAX_CORES)
my_rctd <- run.RCTD(my_rctd, doublet_mode = "full")

# 步骤8: 提取和保存结果
results <- my_rctd@results
weights <- results$weights
norm_weights <- normalize_weights(weights)
weights_df <- as.data.frame(norm_weights)
weights_df$spot_id <- rownames(weights_df)

write.csv(weights_df, file = file.path(OUTPUT_DIR, "rctd_cell_type_proportions_16um.csv"), row.names = FALSE)
save(my_rctd, norm_weights, weights_df, file = file.path(OUTPUT_DIR, "rctd_results_16um.RData"))

# 步骤9: 生成可视化
dominant_cell_type <- colnames(norm_weights)[apply(norm_weights, 1, which.max)]
weights_df$dominant_cell_type <- dominant_cell_type

p1 <- ggplot(weights_df, aes(x = dominant_cell_type, fill = dominant_cell_type)) +
  geom_bar() + theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "RCTD 16um: Dominant Cell Types",
       subtitle = "Visium HD 16um",
       x = "Dominant Cell Type", y = "Number of Spots") +
  guides(fill = "none")
ggsave(file.path(OUTPUT_DIR, "rctd_dominant_cell_types_16um.png"), p1, width = 10, height = 6, dpi = 150)
ggsave(file.path(OUTPUT_DIR, "rctd_dominant_cell_types_16um.pdf"), p1, width = 10, height = 6)

weights_long <- melt(weights_df[, c("spot_id", colnames(norm_weights))],
                     id.vars = "spot_id", variable.name = "cell_type", value.name = "proportion")
p2 <- ggplot(weights_long, aes(x = cell_type, y = proportion, fill = cell_type)) +
  geom_boxplot() + theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "RCTD 16um: Cell Type Proportions",
       x = "Cell Type", y = "Proportion") +
  guides(fill = "none")
ggsave(file.path(OUTPUT_DIR, "rctd_cell_type_proportions_boxplot_16um.png"), p2, width = 12, height = 6, dpi = 150)
ggsave(file.path(OUTPUT_DIR, "rctd_cell_type_proportions_boxplot_16um.pdf"), p2, width = 12, height = 6)

# 保存Seurat对象
common_spots <- intersect(colnames(spatial_data), rownames(norm_weights))
spatial_data_subset <- spatial_data[, common_spots]
for (ct in colnames(norm_weights)) {
  spatial_data_subset@meta.data[[paste0("RCTD_", ct)]] <- norm_weights[common_spots, ct]
}
spatial_data_subset@meta.data$RCTD_dominant_type <- dominant_cell_type[common_spots]
saveRDS(spatial_data_subset, file = file.path(OUTPUT_DIR, "spatial_pancreas_with_rctd_16um.rds"))

# 统计摘要
mean_proportions <- colMeans(norm_weights)
cat("\n各细胞类型平均比例:\n")
for (ct in names(mean_proportions)) {
  cat(sprintf("  - %s: %.4f (%.2f%%)\n", ct, mean_proportions[ct], mean_proportions[ct]*100))
}

stats_file <- file.path(OUTPUT_DIR, "rctd_statistics_16um.txt")
sink(stats_file)
cat("=== RCTD 16µm反卷积结果统计 ===\n\n")
cat(sprintf("分析日期: %s\n", format(Sys.time(), "%Y-%m-%d %H:%M:%S")))
cat(sprintf("样本: Visium HD 16µm\n"))
cat(sprintf("Spots: %d\n", nrow(norm_weights)))
cat(sprintf("细胞类型: %d\n\n", ncol(norm_weights)))
cat("各细胞类型平均比例:\n")
for (ct in names(mean_proportions)) {
  cat(sprintf("  - %s: %.4f (%.2f%%)\n", ct, mean_proportions[ct], mean_proportions[ct]*100))
}
sink()

cat("\n=== RCTD 16um 分析完成! ===\n")

---

脚本 5 — RCTD Visium HD Cellbin（GeoJSON 坐标提取版）

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# rctd_visium_hd_cellbin.R
# RCTD Visium HD Cellbin 空间反卷积分析 - 最终优化版
# 从 cell_segmentations.geojson 提取真实坐标

# 配置参数
SC_DATA_PATH <- "file.path(USER_DATA_DIR, "scRNA.rds")"
SPATIAL_DATA_DIR <- "USER_DATA_DIR/segmented_outputs"
OUTPUT_DIR <- "RCTD_cellbin_results"
N_CELLS_SAMPLE <- 10000
MAX_CORES <- 4
CELL_TYPE_COL <- "Clusters"

setwd(dirname(parent.frame(2)$ofile))
if (!dir.exists(OUTPUT_DIR)) dir.create(OUTPUT_DIR, recursive = TRUE)

suppressPackageStartupMessages({
  library(spacexr); library(Seurat); library(SeuratDisk)
  library(ggplot2); library(dplyr); library(reshape2); library(jsonlite)
})

# 步骤1-3: 与脚本4相同（读取单细胞、准备Reference）
sc_data <- readRDS(SC_DATA_PATH)
cat(sprintf("   - 细胞数量: %d\n", ncol(sc_data)))

set.seed(42)
cells_to_keep <- c()
for (ct in unique(sc_data@meta.data[[CELL_TYPE_COL]])) {
  ct_cells <- WhichCells(sc_data, expression = !!sym(CELL_TYPE_COL) == ct)
  if (length(ct_cells) > 5000) ct_cells <- sample(ct_cells, 5000)
  cells_to_keep <- c(cells_to_keep, ct_cells)
}
sc_data <- subset(sc_data, cells = cells_to_keep)

raw_counts <- GetAssayData(sc_data, layer = "counts")
cell_types <- as.factor(sc_data@meta.data[[CELL_TYPE_COL]])
names(cell_types) <- colnames(sc_data)
reference <- Reference(counts = raw_counts, cell_types = cell_types)

# 步骤4: 读取空间转录组数据 (Cellbin)
spatial_data <- Load10X_Spatial(data.dir = SPATIAL_DATA_DIR, filename = "filtered_feature_cell_matrix.h5")
cat(sprintf("   - Cell数量: %d\n", ncol(spatial_data)))

# 步骤5: 从geojson提取坐标
cat("\n【步骤5/10】从cell_segmentations.geojson提取坐标...\n")
geojson_path <- file.path(SPATIAL_DATA_DIR, "cell_segmentations.geojson")
if (!file.exists(geojson_path)) stop("错误: 找不到geojson文件: ", geojson_path)

geojson_data <- fromJSON(geojson_path)
cell_coords <- data.frame(cell_id = integer(), x = numeric(), y = numeric(), stringsAsFactors = FALSE)
for (feature in geojson_data$features) {
  cell_id <- feature$properties$cell_id
  coords <- feature$geometry$coordinates[[1]]
  x_coords <- sapply(coords, function(c) c[1])
  y_coords <- sapply(coords, function(c) c[2])
  cell_coords <- rbind(cell_coords, data.frame(
    cell_id = cell_id, x = mean(x_coords), y = mean(y_coords)
  ))
}
cat(sprintf("   ✓ 从geojson提取: %d cells\n", nrow(cell_coords)))

# 步骤6: 匹配cell ID
extract_cell_number <- function(cell_id_str) {
  match <- regmatches(cell_id_str, regexpr("cellid_(\\d+)-", cell_id_str))
  if (length(match) > 0) {
    num_str <- gsub("cellid_|-", "", match)
    return(as.integer(num_str))
  }
  return(NA)
}

spatial_data$cell_num <- sapply(colnames(spatial_data), extract_cell_number)
valid_cells <- !is.na(spatial_data$cell_num)
spatial_data <- subset(spatial_data, cells = colnames(spatial_data)[valid_cells])

cell_coords <- cell_coords %>% distinct(cell_id, .keep_all = TRUE)
rownames(cell_coords) <- cell_coords$cell_id

matched_coords <- data.frame(
  x = numeric(nrow(spatial_data@meta.data)),
  y = numeric(nrow(spatial_data@meta.data))
)
for (i in 1:nrow(spatial_data@meta.data)) {
  cell_num <- spatial_data$cell_num[i]
  if (cell_num %in% rownames(cell_coords)) {
    matched_coords[i, ] <- cell_coords[as.character(cell_num), c("x", "y")]
  } else {
    matched_coords[i, ] <- c(NA, NA)
  }
}

valid_mask <- !is.na(matched_coords$x)
spatial_data <- subset(spatial_data, cells = colnames(spatial_data)[valid_mask])
matched_coords <- matched_coords[valid_mask, ]
cat(sprintf("   ✓ 坐标匹配成功: %d cells\n", ncol(spatial_data)))

# 步骤7-8: 采样 + 准备Puck
if (N_CELLS_SAMPLE > 0 && ncol(spatial_data) > N_CELLS_SAMPLE) {
  set.seed(42)
  cells_to_keep <- sample(colnames(spatial_data), N_CELLS_SAMPLE)
  spatial_data <- subset(spatial_data, cells = cells_to_keep)
  matched_coords <- matched_coords[cells_to_keep, ]
}

spatial_counts <- GetAssayData(spatial_data, layer = "counts")
puck <- SpatialRNA(coords = matched_coords, counts = spatial_counts)

# 步骤9: 创建并运行RCTD
my_rctd <- create.RCTD(spatialRNA = puck, reference = reference, max_cores = MAX_CORES)
my_rctd <- run.RCTD(my_rctd, doublet_mode = "full")

# 步骤10: 提取和保存结果
results <- my_rctd@results
weights <- results$weights
norm_weights <- normalize_weights(weights)
weights_df <- as.data.frame(norm_weights)
weights_df$cell_id <- rownames(weights_df)

write.csv(weights_df, file = file.path(OUTPUT_DIR, "rctd_cell_type_proportions_cellbin.csv"), row.names = FALSE)
save(my_rctd, norm_weights, weights_df, matched_coords, file = file.path(OUTPUT_DIR, "rctd_results_cellbin.RData"))

# 基础可视化
dominant_cell_type <- colnames(norm_weights)[apply(norm_weights, 1, which.max)]
weights_df$dominant_cell_type <- dominant_cell_type

p1 <- ggplot(weights_df, aes(x = dominant_cell_type, fill = dominant_cell_type)) +
  geom_bar() + theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "RCTD Cellbin: Dominant Cell Types",
       subtitle = "Visium HD Cellbin",
       x = "Dominant Cell Type", y = "Number of Cells") +
  guides(fill = "none")
ggsave(file.path(OUTPUT_DIR, "rctd_dominant_cell_types_cellbin.png"), p1, width = 10, height = 6, dpi = 150)
ggsave(file.path(OUTPUT_DIR, "rctd_dominant_cell_types_cellbin.pdf"), p1, width = 10, height = 6)

cat("\n=== RCTD Cellbin 分析完成! ===\n")

---

脚本 6 — RCTD 高级可视化（单细胞类型热图 / 综合图 / 气泡图 / 对比图）

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# rctd_visualization_advanced.R
# RCTD 高级可视化（基于已有Seurat对象）

setwd("OUTPUT_DIR")

suppressPackageStartupMessages({
  library(Seurat); library(ggplot2); library(dplyr); library(reshape2)
})

cat("=== RCTD 高级可视化 ===\n\n")

# 1. 读取数据
spatial_data <- readRDS("spatial_pancreas_with_rctd_subset.rds")
weights_df <- read.csv("rctd_cell_type_proportions_subset.csv")

# 2. 获取细胞类型列名
cell_type_cols <- setdiff(colnames(weights_df), "spot_id")

# 3. 单个细胞类型热图（无HE背景）
coords <- GetTissueCoordinates(spatial_data, scale = NULL)
if (ncol(coords) > 2) coords <- coords[, c("x", "y")]
common_spots <- intersect(rownames(coords), weights_df$spot_id)
coords <- coords[common_spots, ]
weights_subset <- weights_df[match(common_spots, weights_df$spot_id), ]

for (ct in cell_type_cols) {
  plot_data <- data.frame(x = coords$x, y = coords$y, proportion = weights_subset[[ct]])
  p <- ggplot(plot_data, aes(x = x, y = y, color = proportion)) +
    geom_point(size = 1.5, alpha = 0.8) +
    scale_color_gradientn(colors = c("grey90", "yellow", "orange", "red", "darkred"),
                         name = "Proportion",
                         limits = c(0, max(plot_data$proportion))) +
    coord_fixed() + theme_void() +
    theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
          legend.position = "right") +
    ggtitle(paste0(ct, " (RCTD)"))
  output_file <- paste0("rctd_heatmap_", gsub(" ", "_", ct), "_no_bg.png")
  ggsave(output_file, p, width = 8, height = 7, dpi = 150, bg = "white")
}

# 4. 所有细胞类型综合图
dominant_types <- apply(weights_subset[, cell_type_cols], 1, function(x) cell_type_cols[which.max(x)])
combined_data <- data.frame(x = coords$x, y = coords$y, dominant_type = dominant_types)

cell_type_colors <- c(
  "ACINAR" = "#E41A1C", "B CELLS" = "#377EB8", "CYCLING DUCTAL" = "#4DAF4A",
  "CYCLING TNK" = "#984EA3", "CYCLING. MYELOID" = "#FF7F00", "DUCTAL" = "#FFFF33",
  "ENDOCRINE" = "#A65628", "ENDOTHELIAL" = "#F781BF", "FIBROBLASTS" = "#999999",
  "MAST" = "#66C2A5", "MYELOID" = "#FC8D62", "PERICYTES" = "#8DA0CB",
  "PLASMA" = "#E78AC3", "TNK" = "#A6D854"
)
used_colors <- cell_type_colors[unique(dominant_types)]

p_combined <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type)) +
  geom_point(size = 1.5, alpha = 0.8) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  coord_fixed() + theme_void() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "right") +
  ggtitle("RCTD: All Cell Types (Dominant)")
ggsave("rctd_all_cell_types_combined.png", p_combined, width = 10, height = 8, dpi = 150, bg = "white")

# 5. 气泡图
combined_data$max_proportion <- apply(weights_subset[, cell_type_cols], 1, max)
p_bubble <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type, size = max_proportion)) +
  geom_point(alpha = 0.6) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  scale_size_continuous(name = "Proportion", range = c(0.5, 4)) +
  coord_fixed() + theme_void() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "right") +
  ggtitle("RCTD: All Cell Types (Size = Proportion)")
ggsave("rctd_all_cell_types_bubble.png", p_bubble, width = 10, height = 8, dpi = 150, bg = "white")

# 6. 前4种细胞类型对比图
mean_props <- colMeans(weights_subset[, cell_type_cols])
top4_types <- names(sort(mean_props, decreasing = TRUE))[1:4]
comparison_data <- NULL
for (ct in top4_types) {
  temp_data <- data.frame(x = coords$x, y = coords$y,
                         proportion = weights_subset[[ct]], cell_type = ct)
  comparison_data <- rbind(comparison_data, temp_data)
}
p_comparison <- ggplot(comparison_data, aes(x = x, y = y, color = proportion)) +
  geom_point(size = 1.2, alpha = 0.8) +
  scale_color_gradientn(colors = c("grey90", "yellow", "orange", "red", "darkred"),
                       name = "Proportion") +
  facet_wrap(~cell_type, ncol = 2) +
  coord_fixed() + theme_void() +
  theme(strip.text = element_text(size = 12, face = "bold"),
        legend.position = "right")
ggsave("rctd_top4_cell_types_comparison.png", p_comparison, width = 12, height = 10, dpi = 150, bg = "white")

cat("\n=== 高级可视化完成 ===\n")

---

脚本 7 — 从 Cellbin geojson 提取真实坐标

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# read_cellbin_coordinates.R
# 读取 cellbin 真实坐标并重新生成可视化

library(jsonlite); library(ggplot2); library(dplyr)
setwd("OUTPUT_DIR")

cat("=== 读取Cellbin真实坐标并重新生成可视化 ===\n\n")

# 1. 读取geojson文件
geojson_file <- "USER_DATA_DIR/segmented_outputs/cell_segmentations.geojson"
geo_data <- fromJSON(geojson_file)
cat("   - Features数量:", length(geo_data$features), "\n")

# 2. 提取细胞坐标
n_features <- length(geo_data$features)
cell_ids <- integer(n_features)
x_coords <- numeric(n_features)
y_coords <- numeric(n_features)

for (i in 1:n_features) {
  cell_ids[i] <- geo_data$features[[i]]$properties$cell_id
  coords <- geo_data$features[[i]]$geometry$coordinates[[1]]
  x_vals <- sapply(coords, function(c) c[[1]])
  y_vals <- sapply(coords, function(c) c[[2]])
  x_coords[i] <- mean(x_vals)
  y_coords[i] <- mean(y_vals)
  if (i %% 10000 == 0) cat(sprintf("   已处理 %d/%d 个细胞...\n", i, n_features))
}

cell_coords <- data.frame(cell_id = cell_ids, x = x_coords, y = y_coords, stringsAsFactors = FALSE)
write.csv(cell_coords, file = "cellbin_coordinates_real.csv", row.names = FALSE)

# 3. 读取RCTD结果并合并
weights_df <- read.csv("rctd_cell_type_proportions_cellbin.csv")
weights_df$cell_num <- as.integer(sub(".*-", "", weights_df$cell_id))
merged_data <- merge(cell_coords, weights_df, by.x = "cell_id", by.y = "cell_num", all.y = TRUE)
merged_data <- merged_data[!is.na(merged_data$x), ]

# 4. 重新生成可视化（使用真实坐标）
cell_type_cols <- setdiff(colnames(weights_df), c("cell_id", "cell_num"))

# 单个细胞类型热图
for (ct in cell_type_cols) {
  plot_data <- data.frame(x = merged_data$x, y = merged_data$y, proportion = merged_data[[ct]])
  p <- ggplot(plot_data, aes(x = x, y = y, color = proportion)) +
    geom_point(size = 0.8, alpha = 0.8) +
    scale_color_gradientn(colors = c("grey90", "yellow", "orange", "red", "darkred"),
                         name = "Proportion",
                         limits = c(0, max(plot_data$proportion, na.rm = TRUE))) +
    coord_fixed() + theme_void() +
    theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
          legend.position = "right") +
    ggtitle(paste0(ct, " (RCTD Cellbin Real Coords)"))
  output_file_png <- paste0("rctd_heatmap_", gsub(" ", "_", gsub("\\.", "_", ct)), "_real_coords_cellbin.png")
  output_file_pdf <- paste0("rctd_heatmap_", gsub(" ", "_", gsub("\\.", "_", ct)), "_real_coords_cellbin.pdf")
  ggsave(output_file_png, p, width = 10, height = 8, dpi = 150, bg = "white")
  ggsave(output_file_pdf, p, width = 10, height = 8)
}

# 综合图 / 气泡图 / 对比图（与前面脚本相同的逻辑，使用真实坐标）
dominant_types <- apply(merged_data[, cell_type_cols], 1, function(x) cell_type_cols[which.max(x)])
combined_data <- data.frame(x = merged_data$x, y = merged_data$y, dominant_type = dominant_types)

cell_type_colors <- c(
  "ACINAR" = "#E41A1C", "B CELLS" = "#377EB8", "CYCLING DUCTAL" = "#4DAF4A",
  "CYCLING TNK" = "#984EA3", "CYCLING. MYELOID" = "#FF7F00", "DUCTAL" = "#FFFF33",
  "ENDOCRINE" = "#A65628", "ENDOTHELIAL" = "#F781BF", "FIBROBLASTS" = "#999999",
  "MAST" = "#66C2A5", "MYELOID" = "#FC8D62", "PERICYTES" = "#8DA0CB",
  "PLASMA" = "#E78AC3", "TNK" = "#A6D854"
)
used_colors <- cell_type_colors[unique(dominant_types)]

p_combined <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type)) +
  geom_point(size = 0.8, alpha = 0.8) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  coord_fixed() + theme_void() +
  ggtitle("RCTD: All Cell Types - Dominant (Cellbin Real Coords)")
ggsave("rctd_all_cell_types_combined_real_coords_cellbin.png", p_combined, width = 12, height = 10, dpi = 150, bg = "white")
ggsave("rctd_all_cell_types_combined_real_coords_cellbin.pdf", p_combined, width = 12, height = 10)

combined_data$max_proportion <- apply(merged_data[, cell_type_cols], 1, max)
p_bubble <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type, size = max_proportion)) +
  geom_point(alpha = 0.6) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  scale_size_continuous(name = "Proportion", range = c(0.3, 2.5)) +
  coord_fixed() + theme_void() +
  ggtitle("RCTD: All Cell Types - Size = Proportion (Cellbin Real Coords)")
ggsave("rctd_all_cell_types_bubble_real_coords_cellbin.png", p_bubble, width = 12, height = 10, dpi = 150, bg = "white")
ggsave("rctd_all_cell_types_bubble_real_coords_cellbin.pdf", p_bubble, width = 12, height = 10)

cat("\n=== 使用真实坐标的可视化完成 ===\n")

---

脚本 8 — Cellbin 完整可视化（综合图/气泡图/对比图）

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# rctd_cellbin_visualization_remaining.R
# RCTD Cellbin 剩余可视化生成（综合图/气泡图/对比图）

setwd("OUTPUT_DIR")

suppressPackageStartupMessages({
  library(ggplot2); library(dplyr); library(reshape2)
})

cat("=== RCTD Cellbin 剩余可视化生成 ===\n\n")

# 1. 读取数据
weights_df <- read.csv("rctd_cell_type_proportions_cellbin.csv")
cell_type_cols <- setdiff(colnames(weights_df), "cell_id")

# 2. 计算主要细胞类型
dominant_cell_type <- apply(weights_df[, cell_type_cols], 1, function(x) cell_type_cols[which.max(x)])
weights_df$dominant_cell_type <- dominant_cell_type

# 3. 创建模拟空间坐标（用真实坐标请见 read_cellbin_coordinates.R）
set.seed(456)
n_cells <- nrow(weights_df)
x_coords <- runif(n_cells, min = 0, max = 1000)
y_coords <- runif(n_cells, min = 0, max = 1000)

# 4. 所有细胞类型综合图
combined_data <- data.frame(x = x_coords, y = y_coords, dominant_type = dominant_cell_type)
cell_type_colors <- c(
  "ACINAR" = "#E41A1C", "B CELLS" = "#377EB8", "CYCLING DUCTAL" = "#4DAF4A",
  "CYCLING TNK" = "#984EA3", "CYCLING. MYELOID" = "#FF7F00", "DUCTAL" = "#FFFF33",
  "ENDOCRINE" = "#A65628", "ENDOTHELIAL" = "#F781BF", "FIBROBLASTS" = "#999999",
  "MAST" = "#66C2A5", "MYELOID" = "#FC8D62", "PERICYTES" = "#8DA0CB",
  "PLASMA" = "#E78AC3", "TNK" = "#A6D854"
)
used_colors <- cell_type_colors[unique(dominant_cell_type)]

p_combined <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type)) +
  geom_point(size = 1.2, alpha = 0.8) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  coord_fixed() + theme_void() +
  ggtitle("RCTD: All Cell Types - Dominant (Cellbin)")
ggsave("rctd_all_cell_types_combined_cellbin.png", p_combined, width = 10, height = 8, dpi = 150, bg = "white")
ggsave("rctd_all_cell_types_combined_cellbin.pdf", p_combined, width = 10, height = 8)

# 5. 气泡图
combined_data$max_proportion <- apply(weights_df[, cell_type_cols], 1, max)
p_bubble <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type, size = max_proportion)) +
  geom_point(alpha = 0.6) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  scale_size_continuous(name = "Proportion", range = c(0.5, 3)) +
  coord_fixed() + theme_void() +
  ggtitle("RCTD: All Cell Types - Size = Proportion (Cellbin)")
ggsave("rctd_all_cell_types_bubble_cellbin.png", p_bubble, width = 10, height = 8, dpi = 150, bg = "white")
ggsave("rctd_all_cell_types_bubble_cellbin.pdf", p_bubble, width = 10, height = 8)

# 6. 前4种细胞类型对比图
mean_proportions <- colMeans(weights_df[, cell_type_cols])
top4_types <- names(sort(mean_proportions, decreasing = TRUE))[1:min(4, length(mean_proportions))]
comparison_data <- NULL
for (ct in top4_types) {
  temp_data <- data.frame(x = x_coords, y = y_coords,
                         proportion = weights_df[[ct]], cell_type = ct)
  comparison_data <- rbind(comparison_data, temp_data)
}
p_comparison <- ggplot(comparison_data, aes(x = x, y = y, color = proportion)) +
  geom_point(size = 1.0, alpha = 0.8) +
  scale_color_gradientn(colors = c("grey90", "yellow", "orange", "red", "darkred"),
                       name = "Proportion") +
  facet_wrap(~cell_type, ncol = 2) +
  coord_fixed() + theme_void() +
  theme(strip.text = element_text(size = 12, face = "bold"))
ggsave("rctd_top4_cell_types_comparison_cellbin.png", p_comparison, width = 12, height = 10, dpi = 150, bg = "white")
ggsave("rctd_top4_cell_types_comparison_cellbin.pdf", p_comparison, width = 12, height = 10)

# 7. 统计摘要
cat("\n各细胞类型平均比例:\n")
for (ct in names(mean_proportions)) {
  cat(sprintf("     - %s: %.3f\n", ct, mean_proportions[ct]))
}

cat("\n=== 剩余可视化生成完成 ===\n")

---

脚本 9 — Cellbin 修正坐标可视化

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# generate_cellbin_plots_fixed.R
# 使用修正后的真实坐标生成 Cellbin 可视化图

library(ggplot2); library(dplyr)
setwd("OUTPUT_DIR")

cat("=== 使用修正后的真实坐标生成Cellbin可视化图 ===\n\n")

# 1. 读取修正后的数据
merged_data <- read.csv("cellbin_merged_real_coords.csv")
cat("   ✓ 数据读取成功:", nrow(merged_data), "个细胞\n")

# 2. 生成可视化图
cell_type_cols <- c("ACINAR", "B.CELLS
", "CYCLING.DUCTAL", "CYCLING.TNK", "CYCLING..MYELOID", "DUCTAL", "ENDOCRINE", "ENDOTHELIAL", "FIBROBLASTS", "MAST", "MYELOID", "PERICYTES", "PLASMA", "TNK")

# 2.1 单个细胞类型热图
cat("2.1 生成单个细胞类型热图...`n")
for (ct in cell_type_cols) {
  plot_data <- data.frame(x = merged_data$x, y = merged_data$y, proportion = merged_data[[ct]])
  p <- ggplot(plot_data, aes(x = x, y = y, color = proportion)) +
    geom_point(size = 0.6, alpha = 0.8) +
    scale_color_gradientn(colors = c("grey90", "yellow", "orange", "red", "darkred"),
                         name = "Proportion", limits = c(0, max(plot_data$proportion, na.rm = TRUE))) +
    coord_fixed() + theme_void() +
    theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"), legend.position = "right") +
    ggtitle(paste0(ct, " (RCTD Cellbin Fixed)"))
  output_file_png <- paste0("rctd_heatmap_", gsub("\\.", "_", ct), "_fixed_cellbin.png")
  output_file_pdf <- paste0("rctd_heatmap_", gsub("\\.", "_", ct), "_fixed_cellbin.pdf")
  ggsave(output_file_png, p, width = 12, height = 10, dpi = 150, bg = "white")
  ggsave(output_file_pdf, p, width = 12, height = 10)
}

# 2.2 所有细胞类型综合图
dominant_types <- apply(merged_data[, cell_type_cols], 1, function(x) cell_type_cols[which.max(x)])
combined_data <- data.frame(x = merged_data$x, y = merged_data$y, dominant_type = dominant_types)
cell_type_colors <- c("ACINAR" = "#E41A1C", "B.CELLS" = "#377EB8", "CYCLING.DUCTAL" = "#4DAF4A",
                     "CYCLING.TNK" = "#984EA3", "CYCLING..MYELOID" = "#FF7F00", "DUCTAL" = "#FFFF33",
                     "ENDOCRINE" = "#A65628", "ENDOTHELIAL" = "#F781BF", "FIBROBLASTS" = "#999999",
                     "MAST" = "#66C2A5", "MYELOID" = "#FC8D62", "PERICYTES" = "#8DA0CB",
                     "PLASMA" = "#E78AC3", "TNK" = "#A6D854")
used_colors <- cell_type_colors[unique(dominant_types)]

p_combined <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type)) +
  geom_point(size = 0.6, alpha = 0.8) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  coord_fixed() + theme_void() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "right") +
  ggtitle("RCTD: All Cell Types - Dominant (Cellbin Fixed)")
ggsave("rctd_all_cell_types_combined_fixed_cellbin.png", p_combined, width = 14, height = 12, dpi = 150, bg = "white")
ggsave("rctd_all_cell_types_combined_fixed_cellbin.pdf", p_combined, width = 14, height = 12)

# 2.3 气泡图
combined_data$max_proportion <- apply(merged_data[, cell_type_cols], 1, max)
p_bubble <- ggplot(combined_data, aes(x = x, y = y, color = dominant_type, size = max_proportion)) +
  geom_point(alpha = 0.6) +
  scale_color_manual(values = used_colors, name = "Cell Type") +
  scale_size_continuous(name = "Proportion", range = c(0.3, 2.5)) +
  coord_fixed() + theme_void() +
  ggtitle("RCTD: All Cell Types - Size = Proportion (Cellbin Fixed)")
ggsave("rctd_all_cell_types_bubble_fixed_cellbin.png", p_bubble, width = 14, height = 12, dpi = 150, bg = "white")
ggsave("rctd_all_cell_types_bubble_fixed_cellbin.pdf", p_bubble, width = 14, height = 12)

# 2.4 前4种细胞类型对比图
mean_proportions <- colMeans(merged_data[, cell_type_cols], na.rm = TRUE)
top4_types <- names(sort(mean_proportions, decreasing = TRUE))[1:min(4, length(mean_proportions))]
comparison_data <- NULL
for (ct in top4_types) {
  temp_data <- data.frame(x = merged_data$x, y = merged_data$y,
                         proportion = merged_data[[ct]], cell_type = ct)
  comparison_data <- rbind(comparison_data, temp_data)
}
p_comparison <- ggplot(comparison_data, aes(x = x, y = y, color = proportion)) +
  geom_point(size = 0.5, alpha = 0.8) +
  scale_color_gradientn(colors = c("grey90", "yellow", "orange", "red", "darkred"),
                       name = "Proportion") +
  facet_wrap(~cell_type, ncol = 2) +
  coord_fixed() + theme_void() +
  theme(strip.text = element_text(size = 12, face = "bold"))
ggsave("rctd_top4_cell_types_comparison_fixed_cellbin.png", p_comparison, width = 16, height = 14, dpi = 150, bg = "white")
ggsave("rctd_top4_cell_types_comparison_fixed_cellbin.pdf", p_comparison, width = 16, height = 14)

cat("`n=== 使用修正后真实坐标的可视化完成 ===`n")

---

RCTD 核心 API

主函数（来自 spacexr 包）

| 函数 | 说明 | |---|---| | Reference(counts, cell_types) | 构建 RCTD 单细胞参考对象 | | SpatialRNA(coords, counts) | 构建空间表达对象（Puck） | | create.RCTD(spatialRNA, reference, max_cores) | 创建 RCTD 对象 | | run.RCTD(rctd, doublet_mode) | 运行反卷积（"full" 或 "doublet"） | | normalize_weights(weights) | 归一化权重矩阵 |

RCTD 对象结果槽

• @results$weights：原始权重矩阵
• @results$cell_type_names：细胞类型列表
• @results$weights_doublet：双细胞权重（如 doublet_mode = "doublet"）

关键参数

| 参数 | 默认值 | 说明 | |---|---|---| | doublet_mode | "doublet" | "full" 估计所有细胞类型比例；"doublet" 仅检测双细胞 | | max_cores | 1 | 并行核心数 | | CELL_TYPE_COL | - | scRNA-seq metadata 中的细胞类型列名 |

---

输出文件清单

数据文件

• rctd_cell_type_proportions*.csv：归一化细胞类型比例（spots × cell types）
• rctd_results*.RData：完整 RCTD 对象（含 weights）
• spatial_pancreas_with_rctd*.rds：带 RCTD 注释的 Seurat 对象
• cellbin_coordinates_real.csv：Cellbin 真实坐标
• cellbin_merged_real_coords.csv：Cellbin 合并数据

图形文件

• rctd_dominant_cell_types*.png/pdf：主要细胞类型分布
• rctd_cell_type_proportions_boxplot*.png/pdf：箱线图
• rctd_spatial_cell_types*.png：前 4 种细胞类型空间分布组合图
• rctd_spatial_dominant_types*.png：主导细胞类型空间分布
• rctd_heatmap_*_no_bg*.png/pdf：单细胞类型热图
• rctd_all_cell_types_combined*.png/pdf：综合图
• rctd_all_cell_types_bubble*.png/pdf：气泡图
• rctd_top4_cell_types_comparison*.png/pdf：前 4 种对比图

报告

• rctd_statistics_*.txt：统计摘要

---

故障排除

数据格式错误

错误：输入数据格式不正确

解决：

1. 检查 scRNA-seq 数据格式（Seurat 对象或 count 矩阵）
2. 检查 spatial 数据格式
3. 确保基因名一致

内存不足

错误：MemoryError / cannot allocate vector

解决：

1. 减少数据量（采样 spots）
2. 使用分块处理
3. 增加系统内存
4. 8µm/Cellbin 推荐 ≥64GB RAM

计算速度慢

解决：

1. 减少迭代次数 / 采样 spots
2. 使用并行计算（max_cores）
3. 减少单细胞参考规模（每类型 ≤5000 cells）

结果不准确

解决：

1. 检查参考数据质量（细胞类型注释）
2. 调整参数（doublet_mode、max_cores）
3. 使用更多细胞类型 / 优化采样

Cellbin 坐标为空

原因：cell ID 格式不匹配

解决：

• 使用 extract_cell_number() 正则提取数字部分
• 确保 cell_segmentations.geojson 存在
• 检查 RCTD 结果的 cell_id 格式

---

性能优化

大数据集处理

• 分块处理
• 减少高变基因数量
• 使用降采样（每类型 5000 cells + 10000 spots）

参数调优

| 参数 | 调整建议 | |---|---| | N_SPOTS_SAMPLE | 根据内存调整（0 = 不采样） | | MAX_CORES | 设为 CPU 核心数 | | doublet_mode | "full" 更准确但更慢 |

---

推荐工作流

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

按使用场景推荐：

| 场景 | 推荐脚本 | |---|---| | 16µm 标准全流程 | 脚本 1 | | 16µm 最终优化版（含采样） | 脚本 4 | | 8µm 高精度 | 脚本 2 | | Cellbin 单细胞分辨率（模拟坐标） | 脚本 3 | | Cellbin 单细胞分辨率（GeoJSON 真实坐标） | 脚本 5 | | Cellbin 真实坐标再生成可视化 | 脚本 7 | | 高级可视化（基于已有 Seurat 对象） | 脚本 6 | | Cellbin 综合图/气泡图/对比图 | 脚本 8 | | Cellbin 修正坐标可视化 | 脚本 9 |

典型使用顺序：

1. 先运行脚本 1（16µm）或脚本 2（8µm）得到基础结果
2. 若需 cellbin 级别，运行脚本 5（提供 geojson 时）或脚本 3（无 geojson 时）
3. 用脚本 7 重新生成基于真实坐标的可视化
4. 用脚本 6 / 8 / 9 生成综合图、气泡图、对比图

---

资源链接

• spacexr 源码：github.com/dmcable/spacexr
• RCTD Bioconductor：bioconductor.org/packages/RCTD
• RCTD 论文：Cable DM, et al. Robust decomposition of cell type mixtures in spatial transcriptomics. Nature Biotechnology, 2022.

---

许可证

GPL-3 License