空间转录组Visium分析¶

# 安装Space Ranger（需要从10x Genomics官网下载）
# 下载地址：https://www.10xgenomics.com/support/software/space-ranger

# 运行spaceranger count（核心命令）
spaceranger count \                      # 主命令：计数分析
  --id=sample_01 \                       # 样本ID，自己命名
  --transcriptome=/ref/refdata-gex-GRCh38 \  # 参考基因组路径
  --fastqs=/data/fastqs/ \               # 测序FASTQ文件目录
  --sample=Sample1 \                     # FASTQ文件中的样本名
  --image=/data/tissue_image.tif \       # 组织H&E染色图像
  --slide=V19L01-041 \                   # 载玻片序列号
  --area=A1 \                            # 捕获区域编号
  --localcores=16 \                      # 使用CPU核心数
  --localmem=64                          # 最大内存（GB）

# 安装必要的R包
# install.packages("Seurat")           # 安装Seurat（如果没有）
library(Seurat)                         # 加载Seurat包

# 加载Visium数据
# Load10X_Spatial()读取spaceranger输出，自动关联图像
spatial_data <- Load10X_Spatial(
  data.dir = "spaceranger_output/outs/",  # Space Ranger输出目录
  filename = "filtered_feature_bc_matrix.h5",  # 表达矩阵文件
  assay = "Spatial",                     # 检测类型命名为Spatial
  slice = "tissue_slice"                 # 切片名称
)

# 查看数据基本信息
spatial_data                             # 显示基因数、spot数等信息

# 质控：计算线粒体基因比例（高比例说明细胞可能受损）
spatial_data[["percent.mt"]] <- PercentageFeatureSet(
  spatial_data,                          # 输入Seurat对象
  pattern = "^MT-"                       # 线粒体基因以MT-开头
)

# 空间可视化质控指标
SpatialFeaturePlot(
  spatial_data,                          # Seurat对象
  features = "nCount_Spatial"            # 可视化每个spot的UMI总数
)

# 过滤低质量spot
spatial_data <- subset(
  spatial_data,                          # 原始数据
  subset = nFeature_Spatial > 200 &      # 至少检测到200个基因
    nFeature_Spatial < 6000 &            # 最多6000个基因（去除doublet）
    percent.mt < 20                      # 线粒体基因比例不超过20%
)

# SCTransform归一化（Seurat推荐方法）
# 白话：校正技术噪音，让生物学差异凸显出来
spatial_data <- SCTransform(
  spatial_data,                          # 输入Seurat对象
  assay = "Spatial",                     # 指定Spatial assay
  verbose = FALSE                        # 不显示冗余信息
)

# PCA降维（把高维数据压缩到少数几个主成分）
spatial_data <- RunPCA(
  spatial_data,                          # 输入数据
  assay = "SCT",                         # 使用SCTransform后的数据
  verbose = FALSE                        # 静默模式
)

# 构建近邻图 + 聚类
spatial_data <- FindNeighbors(
  spatial_data,                          # 输入数据
  reduction = "pca",                     # 基于PCA结果
  dims = 1:30                            # 使用前30个主成分
)
spatial_data <- FindClusters(
  spatial_data,                          # 输入数据
  verbose = FALSE,                       # 静默模式
  resolution = 0.8                       # 分辨率，越大cluster越多
)

# UMAP降维（用于可视化）
spatial_data <- RunUMAP(
  spatial_data,                          # 输入数据
  reduction = "pca",                     # 基于PCA
  dims = 1:30                            # 使用前30个主成分
)

# 在组织切片上可视化聚类结果
SpatialDimPlot(
  spatial_data,                          # Seurat对象
  label = TRUE,                          # 显示cluster编号
  label.size = 3                         # 标签字体大小
)

# 使用Moran's I检测空间可变基因
# Moran's I衡量空间自相关性：值越高说明空间分布越不随机
spatial_data <- FindSpatiallyVariableFeatures(
  spatial_data,                          # 输入数据
  assay = "SCT",                         # 使用SCTransform后的数据
  features = VariableFeatures(spatial_data),  # 候选基因列表
  selection.method = "moransi"           # 使用Moran's I方法
)

# 查看top空间可变基因
top_svg <- head(
  SpatiallyVariableFeatures(spatial_data, selection.method = "moransi"),
  n = 6                                  # 前6个SVG
)
print(top_svg)                            # 打印结果

# 可视化空间可变基因的表达分布
SpatialFeaturePlot(
  spatial_data,                          # Seurat对象
  features = top_svg[1:4],               # 显示前4个SVG
  ncol = 2,                              # 每行2个图
  alpha = c(0.1, 1)                      # 透明度范围
)

# 需要一个单细胞RNA-seq参考数据集
# 假设已有配套的scRNA-seq Seurat对象：sc_ref
# 1. 准备参考数据
library(spacexr)                         # 加载spacexr包（包含RCTD）

# 提取参考数据的计数矩阵和细胞类型信息
counts_ref <- GetAssayData(sc_ref, slot = "counts")  # 参考表达矩阵
cell_types <- sc_ref$celltype            # 每个细胞的类型标签
names(cell_types) <- colnames(sc_ref)    # 确保有细胞barcode名字

# 创建Reference对象
reference <- Reference(
  counts = counts_ref,                   # 参考表达矩阵
  cell_types = cell_types                # 细胞类型标签
)

# 2. 准备空间数据
counts_st <- GetAssayData(spatial_data, assay = "Spatial", slot = "counts")
coords <- GetTissueCoordinates(spatial_data)  # 获取空间坐标
query <- SpatialRNA(
  coords = coords,                       # 空间坐标
  counts = counts_st                     # spot表达矩阵
)

# 3. 运行RCTD
rctd <- create.RCTD(query, reference, max_cores = 8)  # 创建RCTD对象
rctd <- run.RCTD(rctd, doublet_mode = "doublet")      # 运行反卷积

# 4. 将结果加入Seurat对象
spatial_data <- AddMetaData(
  spatial_data,                          # 原始Seurat对象
  metadata = rctd@results$weights        # RCTD估计的细胞类型比例
)

# ========== 快速命令速查 ==========

# Space Ranger流程
spaceranger count --id=ID --transcriptome=REF --fastqs=DIR --image=IMG --slide=SN --area=A1

# R分析核心流程（6步）
data <- Load10X_Spatial("outs/")         # 1. 加载
data <- SCTransform(data, assay="Spatial")  # 2. 归一化
data <- RunPCA(data)                     # 3. PCA
data <- FindNeighbors(data, dims=1:30)   # 4. 建图
data <- FindClusters(data, resolution=0.8)  # 5. 聚类
SpatialDimPlot(data, label=TRUE)         # 6. 可视化

# 空间可变基因
FindSpatiallyVariableFeatures(data, selection.method="moransi")

# 关键可视化函数
SpatialFeaturePlot(data, features="GENE")  # 单基因空间表达图
SpatialDimPlot(data)                     # 聚类空间分布图