Spotlight-linux-R — SPOTlight 空间反卷积

概述

SPOTlight 是基于 NMF（非负矩阵分解）和 Seurat 的空间转录组反卷积方法，结合 scRNA-seq 参考数据，估计每个空间 spot 的细胞类型组成。本 Skill 提供完整的 R 端 SPOTlight 分析流程，针对 Visium 和 Visium HD 数据进行了优化。

核心特性：
- NMF + Seurat：基于 marker 基因的 NMF 分解
- Visium / Visium HD：支持多种空间平台
- 自动降采样：大数据集自动降采样到 10000 spots
- 多种可视化：空间分布、饼图、热图、条形图、小提琴图
- 完整结果整合：将反卷积结果直接写入 Seurat 对象 metadata

适用场景

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

• 使用 SPOTlight 对 Visium/Visium HD 数据进行反卷积
• 处理大规模空间数据（自动降采样）
• 生成 SPOTlight 空间分布图
• 比较不同 cell type 的比例分布

安装

环境要求

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

安装 SPOTlight R 包

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

安装依赖包

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

验证安装

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

常见安装问题

SPOTlight 反卷积原理

SPOTlight 通过以下步骤反卷积空间数据：

1. marker 基因识别：使用 scRNA-seq 数据 FindAllMarkers 找到每种细胞类型的 marker 基因
2. NMF 分解：对 scRNA-seq 表达矩阵进行 NMF，得到每种细胞类型的特征基因权重
3. 基分解：对空间 spot 表达使用相同的基进行分解
4. 比例估计：通过基的系数估计每种细胞类型在每个 spot 中的贡献

核心步骤：
1. 预处理 scRNA-seq（NormalizeData → FindVariableFeatures → ScaleData → RunPCA）
2. 找 marker 基因（FindAllMarkers）
3. 运行 SPOTlight(x, y, groups, mgs, hvg, weight_id, group_id, gene_id)
4. 提取 $mat 矩阵作为反卷积结果

完整可运行脚本

输入数据要求

配置说明

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

# ========== 配置区域（请根据实际数据修改） ==========
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_PER_TYPE   <- 5000                          # 单细胞每类型最大采样数
N_SPOTS_SAMPLE     <- 10000                         # 空间采样数（0=不采样）
N_TOP_TYPES_PLOT   <- 6                             # 可视化显示前 N 种细胞类型
OUTPUT_DIR         <- "SPOTlight_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 — SPOTlight 反卷积（含自动降采样）

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# spotlight_deconvolution_downsampled.R
# SPOTlight Deconvolution Script with Downsampling
# Usage: Rscript spotlight_deconvolution_downsampled.R

library(Seurat); library(SPOTlight); library(ggplot2); library(dplyr); library(pheatmap)

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR    <- "your/data/directory"          # TODO: 修改为你的数据根目录
SC_FILENAME      <- "scRNA.rds"                    # TODO: 单细胞 RDS 文件名
SPATIAL_SUBDIR   <- "binned_outputs/square_016um"  # 空间数据子目录（相对 USER_DATA_DIR）
CELL_TYPE_COL    <- "cell_type"                    # TODO: 修改为你的细胞类型列名
N_CELLS_PER_TYPE <- 5000                           # 单细胞每类型最大采样数（避免不平衡）
N_SPOTS_SAMPLE   <- 10000                          # 空间采样数（0=不采样）
N_TOP_TYPES_PLOT <- 6                              # 可视化显示前 N 种细胞类型
OUTPUT_DIR       <- "SPOTlight_results"            # 输出目录
# =====================================================

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

# 自动检测细胞类型列名（按 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 参数。")
}

cat("============================================================`n")
cat("SPOTlight Deconvolution - With Downsampling`n")
cat("============================================================`n`n")

# 步骤1: 加载单细胞参考数据
sc_data <- readRDS(sc_data_path)
CELL_TYPE_COL <- resolve_cell_type_col(sc_data, CELL_TYPE_COL)
print(paste("  使用细胞类型列:", CELL_TYPE_COL))
sc_data$cell_type <- sc_data[[CELL_TYPE_COL, drop = TRUE]]
cat(sprintf("  Loaded: %d cells x %d genes`n", ncol(sc_data), nrow(sc_data)))

# 步骤2: 降采样单细胞数据
cat("`nStep 2: Downsampling single-cell data...`n")
set.seed(42)
cell_types <- unique(sc_data$cell_type)
cells_to_keep <- c()
for (ct in cell_types) {
  ct_cells <- WhichCells(sc_data, expression = cell_type == ct)
  if (N_CELLS_PER_TYPE > 0 && length(ct_cells) > N_CELLS_PER_TYPE) {
    ct_cells <- sample(ct_cells, N_CELLS_PER_TYPE)
  }
  cells_to_keep <- c(cells_to_keep, ct_cells)
}
sc_data <- subset(sc_data, cells = cells_to_keep)
cat(sprintf("  Downsampled to %d cells`n", ncol(sc_data)))

# 步骤3: 预处理单细胞数据
cat("`nStep 3: Preprocessing single-cell data...`n")
sc_data <- subset(sc_data, nFeature_RNA > 200 & nFeature_RNA < 7000 & percent.mt < 20)
sc_data <- NormalizeData(sc_data)
sc_data <- FindVariableFeatures(sc_data, selection.method = "vst", nfeatures = 2000)
sc_data <- ScaleData(sc_data)
sc_data <- RunPCA(sc_data, features = VariableFeatures(object = sc_data))

# 步骤4: 加载空间数据
cat("`nStep 4: Loading spatial data...`n")
spatial_data <- Load10X_Spatial(
  data.dir = spatial_base_path,
  filename = "filtered_feature_bc_matrix.h5"
)
cat(sprintf("  Loaded: %d spots x %d genes`n", ncol(spatial_data), nrow(spatial_data)))

# 步骤5: 预处理空间数据
cat("`nStep 5: Preprocessing spatial data...`n")
spatial_data <- subset(spatial_data, nFeature_Spatial > 10)
spatial_data <- NormalizeData(spatial_data, assay = "Spatial")

# 步骤6: 找 marker 基因
cat("`nStep 6: Finding marker genes...`n")
Idents(sc_data) <- sc_data$cell_type
markers <- FindAllMarkers(sc_data, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top_markers <- markers %>%
  group_by(cluster) %>%
  top_n(n = 50, wt = avg_log2FC) %>%
  ungroup()
cat(sprintf("  Found %d markers`n", nrow(markers)))

# 步骤7: 准备数据
cat("`nStep 7: Preparing data...`n")
sc_counts <- GetAssayData(sc_data, layer = "counts")
spatial_counts <- GetAssayData(spatial_data, layer = "counts", assay = "Spatial")
sc_cell_types <- sc_data$cell_type
names(sc_cell_types) <- colnames(sc_data)

common_genes <- intersect(rownames(sc_counts), rownames(spatial_counts))
common_genes <- intersect(common_genes, top_markers$gene)
cat(sprintf("  Common genes: %d`n", length(common_genes)))

# 步骤8: 运行 SPOTlight
cat("`nStep 8: Running SPOTlight...`n")
spotlight_results <- SPOTlight(
  x = sc_counts,
  y = spatial_counts,
  groups = sc_cell_types,
  mgs = top_markers,
  hvg = common_genes,
  weight_id = "avg_log2FC",
  group_id = "cluster",
  gene_id = "gene"
)

cat("`n✓ SPOTlight completed!`n")
decon_mat <- spotlight_results$mat
spatial_data@meta.data <- cbind(spatial_data@meta.data, decon_mat)

# 步骤9: 可视化
cat("`nStep 9: Creating visualizations...`n")
cell_types <- colnames(decon_mat)
decon_props <- decon_mat / rowSums(decon_mat)

# 空间分布图（前 6 种细胞类型）
n_types <- min(N_TOP_TYPES_PLOT, length(cell_types))
for (i in 1:n_types) {
  ct <- cell_types[i]
  p <- SpatialFeaturePlot(spatial_data, features = ct, pt.size.factor = 1.5) +
    ggtitle(paste("Cell Type:", ct))
  ggsave(file.path(output_dir, sprintf("spotlight_%s.png", gsub(" ", "_", ct))),
         p, width = 8, height = 8, dpi = 150)
}

# 热图
if (nrow(decon_props) > 500) {
  sample_idx <- sample(1:nrow(decon_props), 500)
  decon_props_plot <- decon_props[sample_idx, ]
} else {
  decon_props_plot <- decon_props
}
png(file.path(output_dir, "spotlight_heatmap.png"), width = 1200, height = 800, res = 150)
pheatmap(decon_props_plot, cluster_rows = TRUE, cluster_cols = TRUE,
         show_rownames = FALSE, main = "Cell Type Proportions",
         color = colorRampPalette(c("white", "yellow", "red", "darkred"))(100))
dev.off()

# 步骤10: 保存结果
cat("`nStep 10: Saving results...`n")
write.csv(decon_mat, file.path(output_dir, "spotlight_deconvolution.csv"), row.names = TRUE)
write.csv(decon_props, file.path(output_dir, "spotlight_proportions.csv"), row.names = TRUE)
saveRDS(spatial_data, file.path(output_dir, "spatial_with_spotlight.rds"))

cat("`n✓ Analysis completed!`n")
cat(sprintf("Output directory: %s`n", output_dir))

脚本 2 — SPOTlight 反卷积（任意路径版本）

> 与脚本 1 功能相同，但去掉了内部降采样，假设单细胞输入已经过预处理/降采样。路径完全通过 USER_DATA_DIR 配置。

# Author: LKP <kunpeng.liao@abiosciences.com>
# Date:   2026-08-18
#
# spotlight_deconvolution.R
# SPOTlight Deconvolution Script (path-agnostic)
# Usage: Rscript spotlight_deconvolution.R

library(Seurat); library(SPOTlight); library(ggplot2); library(dplyr); library(pheatmap)

# ========== 配置区域（请根据实际数据修改） ==========
USER_DATA_DIR    <- "your/data/directory"          # Windows: "D:/data"; WSL: "/mnt/d/data"; Linux: "/home/user/data"
SC_FILENAME      <- "scRNA.rds"                    # TODO: 单细胞 RDS 文件名（已预处理的）
SPATIAL_SUBDIR   <- "binned_outputs/square_016um"  # 空间数据子目录
CELL_TYPE_COL    <- "cell_type"                    # TODO: 修改为你的细胞类型列名
N_SPOTS_SAMPLE   <- 10000                          # 空间采样数（0=不采样）
N_TOP_TYPES_PLOT <- 6                              # 可视化显示前 N 种细胞类型
OUTPUT_DIR       <- "SPOTlight_results"            # 输出目录
# =====================================================

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

# 自动检测细胞类型列名（按 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 参数。")
}

cat("============================================================`n")
cat("SPOTlight Deconvolution`n")
cat("============================================================`n`n")

# 步骤1: 加载并预处理单细胞参考数据
sc_data <- readRDS(sc_data_path)
CELL_TYPE_COL <- resolve_cell_type_col(sc_data, CELL_TYPE_COL)
print(paste("  使用细胞类型列:", CELL_TYPE_COL))
sc_data$cell_type <- sc_data[[CELL_TYPE_COL, drop = TRUE]]
sc_data <- subset(sc_data, nFeature_RNA > 200 & percent.mt < 20)
sc_data <- NormalizeData(sc_data)
sc_data <- FindVariableFeatures(sc_data, selection.method = "vst", nfeatures = 2000)
sc_data <- ScaleData(sc_data)
sc_data <- RunPCA(sc_data, features = VariableFeatures(object = sc_data))
cat(sprintf("  SC data: %d cells`n", ncol(sc_data)))

# 步骤2: 找 marker 基因
Idents(sc_data) <- sc_data$cell_type
markers <- FindAllMarkers(sc_data, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top_markers <- markers %>%
  group_by(cluster) %>%
  top_n(n = 50, wt = avg_log2FC) %>%
  ungroup()

# 步骤3: 加载空间数据并降采样
spatial_data <- Load10X_Spatial(data.dir = spatial_base_path, filename = "filtered_feature_bc_matrix.h5")
set.seed(42)
if (N_SPOTS_SAMPLE > 0 && ncol(spatial_data) > N_SPOTS_SAMPLE) {
  spots_to_keep <- sample(colnames(spatial_data), N_SPOTS_SAMPLE)
  spatial_data <- subset(spatial_data, cells = spots_to_keep)
}
spatial_data <- subset(spatial_data, nFeature_Spatial > 10)
spatial_data <- NormalizeData(spatial_data, assay = "Spatial")
cat(sprintf("  Spatial data: %d spots`n", ncol(spatial_data)))

# 步骤4: 准备数据
sc_counts <- GetAssayData(sc_data, layer = "counts")
spatial_counts <- GetAssayData(spatial_data, layer = "counts", assay = "Spatial")
sc_cell_types <- sc_data$cell_type
names(sc_cell_types) <- colnames(sc_data)

common_genes <- intersect(rownames(sc_counts), rownames(spatial_counts))
common_genes <- intersect(common_genes, top_markers$gene)
cat(sprintf("  Common genes: %d`n", length(common_genes)))

# 步骤5: 运行 SPOTlight
spotlight_results <- SPOTlight(
  x = sc_counts, y = spatial_counts,
  groups = sc_cell_types, mgs = top_markers, hvg = common_genes,
  weight_id = "avg_log2FC", group_id = "cluster", gene_id = "gene"
)

decon_mat <- spotlight_results$mat
cat(sprintf("Deconvolution matrix: %d spots x %d cell types`n",
            nrow(decon_mat), ncol(decon_mat)))
spatial_data@meta.data <- cbind(spatial_data@meta.data, decon_mat)

# 步骤6: 可视化
cell_types <- colnames(decon_mat)

# 1. 空间分布图（前 6 种细胞类型）
n_types <- min(N_TOP_TYPES_PLOT, length(cell_types))
for (i in 1:n_types) {
  ct <- cell_types[i]
  safe_name <- gsub("[ .]", "_", ct)
  p <- SpatialFeaturePlot(spatial_data, features = ct,
                          pt.size.factor = 1.5, alpha = c(0.1, 1)) +
    ggtitle(paste("Cell Type:", ct))
  ggsave(file.path(output_dir, sprintf("spotlight_%s.png", safe_name)),
         p, width = 8, height = 8, dpi = 150)
}

# 2. 综合图（主导细胞类型）
dominant_cell_type <- colnames(decon_mat)[apply(decon_mat, 1, which.max)]
spatial_data$dominant_cell_type <- dominant_cell_type
p_combined <- SpatialDimPlot(spatial_data, group.by = "dominant_cell_type",
                             pt.size.factor = 1.5, label = TRUE) +
  ggtitle("Dominant Cell Types")
ggsave(file.path(output_dir, "spotlight_combined.png"),
       p_combined, width = 12, height = 10, dpi = 150)

# 3. 比例热图
decon_props <- decon_mat / rowSums(decon_mat)
if (nrow(decon_props) > 500) {
  set.seed(42)
  sample_idx <- sample(1:nrow(decon_props), 500)
  decon_props_plot <- decon_props[sample_idx, ]
} else {
  decon_props_plot <- decon_props
}
png(file.path(output_dir, "spotlight_heatmap.png"), width = 1200, height = 800, res = 150)
pheatmap(decon_props_plot, cluster_rows = TRUE, cluster_cols = TRUE,
         show_rownames = FALSE, main = "Cell Type Proportions per Spot",
         color = colorRampPalette(c("white", "yellow", "red", "darkred"))(100))
dev.off()

# 4. 平均比例条形图
avg_props <- colMeans(decon_props)
prop_df <- data.frame(CellType = names(avg_props), Proportion = avg_props)
p <- ggplot(prop_df, aes(x = reorder(CellType, -Proportion), y = Proportion, fill = CellType)) +
  geom_bar(stat = "identity") + theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Average Cell Type Proportions", x = "Cell Type", y = "Average Proportion") +
  guides(fill = "none")
ggsave(file.path(output_dir, "spotlight_barplot.png"), p, width = 10, height = 6, dpi = 150)

# 5. 小提琴图
decon_long <- as.data.frame(decon_props)
decon_long$spot <- rownames(decon_long)
decon_long <- tidyr::pivot_longer(decon_long, cols = -spot,
                                  names_to = "CellType", values_to = "Proportion")
if (nrow(decon_long) > 50000) {
  set.seed(42)
  sample_idx <- sample(1:nrow(decon_long), 50000)
  decon_long_plot <- decon_long[sample_idx, ]
} else {
  decon_long_plot <- decon_long
}
p <- ggplot(decon_long_plot, aes(x = CellType, y = Proportion, fill = CellType)) +
  geom_violin() + geom_boxplot(width = 0.1, alpha = 0.5) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Distribution of Cell Type Proportions",
       x = "Cell Type", y = "Proportion") +
  guides(fill = "none")
ggsave(file.path(output_dir, "spotlight_violin.png"), p, width = 12, height = 6, dpi = 150)

# 步骤7: 保存结果
write.csv(decon_mat, file.path(output_dir, "spotlight_deconvolution.csv"), row.names = TRUE)
write.csv(decon_props, file.path(output_dir, "spotlight_proportions.csv"), row.names = TRUE)
saveRDS(spatial_data, file.path(output_dir, "spatial_with_spotlight.rds"))

cat("`n✓ All results saved!`n")
cat("`n============================================================`n")
cat("SPOTlight analysis completed successfully!`n")
cat("============================================================`n")
cat("`nSummary:`n")
cat(sprintf("  - Reference cells: %d`n", ncol(sc_counts)))
cat(sprintf("  - Spatial spots: %d`n", ncol(spatial_counts)))
cat(sprintf("  - Cell types: %d`n", length(cell_types)))
cat(sprintf("  - Output directory: %s`n", output_dir))

SPOTlight 核心 API

主函数

SPOTlight(
  x,                  # scRNA-seq 表达矩阵（dgCMatrix）
  y,                  # 空间表达矩阵（dgCMatrix）
  groups,             # scRNA-seq 细胞类型向量
  mgs,                # marker 基因 data.frame（包含 cluster/gene/avg_log2FC）
  hvg,                # 高变基因名向量
  weight_id = "avg_log2FC",
  group_id = "cluster",
  gene_id = "gene",
  ...
)

返回值

• mat：spots × cell_types 的反卷积比例矩阵
• deconMots：每个 spot 的分解基
• res：拟合结果

关键参数

输出文件清单

数据文件

• spotlight_deconvolution.csv：反卷积矩阵（spots × cell types）
• spotlight_proportions.csv：归一化比例矩阵
• spatial_with_spotlight.rds：带 SPOTlight 注释的 Seurat 对象

图形文件

• spotlight_*.png：前 6 种细胞类型的空间分布图
• spotlight_combined.png：主导细胞类型综合空间图
• spotlight_heatmap.png：细胞类型比例热图
• spotlight_barplot.png：平均比例条形图
• spotlight_violin.png：细胞类型比例小提琴图

故障排除

数据格式错误

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

解决：
1. 检查 scRNA-seq 数据格式（Seurat 对象）
2. 检查 spatial 数据格式
3. 确保基因名一致
4. 确认 marker 基因 data.frame 包含正确的列名

内存不足

错误：MemoryError / cannot allocate vector

解决：
1. 减少数据量（采样 spots）
2. 使用降采样（自动到 10000 spots）
3. 增加系统内存

计算速度慢

解决：
1. 使用降采样（每类型 ≤ 5000 cells）
2. 减少 marker 基因数（top_n = 50）
3. 减少 HVG 数（nfeatures = 2000）

结果不准确

解决：
1. 检查参考数据质量（细胞类型注释）
2. 增加 marker 基因数（top_n = 100）
3. 调整 min.pct 和 logfc.threshold 参数

性能优化

大数据集处理

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

参数调优

推荐工作流

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

按使用场景推荐：

两个脚本的差异：
- 脚本 1 显式对 scRNA-seq 做降采样（每类型 ≤ N_CELLS_PER_TYPE），适合大数据集
- 脚本 2 依赖输入文件已经过预处理/降采样，运行更快

典型使用顺序：

1. 先用脚本 1（标准化全流程）
2. 若数据已预处理，改用脚本 2（更轻量）

资源链接

• SPOTlight 源码：github.com/MarcElosua/SPOTlight
• Bioconductor：bioconductor.org/packages/SPOTlight
• SPOTlight 论文：Elosua-Bayes A, et al. SPOTlight: seeded NMF regression to deconvolute spatial transcriptomics spots with single-cell transcriptomes. Nucleic Acids Research, 2021.

许可证

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