seurat 读入10X spaceranger 细胞分割的结果并显示
seurat 读入10X spaceranger 细胞分割的结果并显示
下面是代码：
#!/usr/bin/env Rscript
############################################################################
# 读取 10x Visium HD 细胞分割结果并展示
# 依赖: Seurat(>=5.x, 含 Read10X_Segmentations), ggplot2
############################################################################
# ---------------------- spaceranger 输出 outs/ 目录说明 ----------------------
# outs/
# ├── web_summary.html QC 汇总网页(浏览器打开即可看)
# ├── metrics_summary.csv 关键指标表(细胞/spot 数、测序质量、比对率等)
# ├── molecule_info.h5 每个分子的位置与归属(分子级信息)
# ├── barcode_mappings.parquet 2µm bin 与细胞分割 barcode 的对应关系
# ├── feature_slice.h5 / probe_set.csv 探针组(feature)定义 / 探针列表
# ├── cloupe_*.cloupe Loupe Browser 可视化文件
# │
# ├── spatial/ 【顶层图像】全分辨率组织图与对齐信息(此层无坐标表)
# │ ├── tissue_hires_image.png 高分辨率组织图
# │ ├── tissue_lowres_image.png 低分辨率组织图(出图常用)
# │ ├── cytassist_image.tiff CytAssist 原始成像
# │ ├── aligned_fiducials.jpg 基准点对齐图(人工核对用)
# │ ├── aligned_tissue_image.jpg 组织图与载玻片对齐结果
# │ ├── detected_tissue_image.jpg 自动识别出的组织区域
# │ └── final_alignment.json 图像对齐参数
# │
# ├── binned_outputs/ 【方阵(bin)结果】每档 bin 一个子目录
# │ ├── square_002um/ 2µm(最细, 最大, 一般不用)
# │ ├── square_008um/ 8µm
# │ ├── square_016um/ 16µm
# │ └── square_050um/ 50µm
# │ ├── filtered_feature_bc_matrix/ 过滤后表达矩阵(仅组织上的 bin), 10x 三件套:
# │ │ ├── barcodes.tsv.gz bin 名称
# │ │ ├── features.tsv.gz 基因 ID / 名称 / 类型
# │ │ └── matrix.mtx.gz 表达计数(稀疏矩阵)
# │ ├── filtered_feature_bc_matrix.h5 同上, 合并的 HDF5 格式
# │ ├── raw_feature_bc_matrix[.h5] 未过滤(全部 bin, 含背景)
# │ ├── analysis/ Space Ranger 自带聚类/PCA/UMAP/差异分析
# │ └── spatial/
# │ ├── tissue_positions.parquet ★每个 bin 的行列号 + 像素坐标(坐标表)
# │ ├── scalefactors_json.json 图像缩放系数 + spot_diameter_fullres
# │ └── tissue_hires/lowres_image.png
# │
# └── segmented_outputs/ 【细胞分割结果】← 本代码读取这里
# ├── filtered_feature_cell_matrix/ 基于分割细胞的表达矩阵(10x 三件套)
# │ ├── barcodes.tsv.gz 细胞名 cellid_000000001-1
# │ ├── features.tsv.gz
# │ └── matrix.mtx.gz
# ├── filtered_feature_cell_matrix.h5 同上, HDF5 格式
# ├── raw_feature_cell_matrix[.h5] 未过滤的细胞矩阵
# ├── cell_segmentations.geojson ★细胞分割: 每个细胞的 polygon 顶点 + 质心 + id
# ├── nucleus_segmentations.geojson ★细胞核分割(结构同上)
# ├── graphclust_annotated_*.geojson 带 graphclust 聚类注释的分割文件
# ├── cloupe.cloupe Loupe Browser 文件
# ├── analysis/ 分割数据的聚类/PCA/UMAP/差异分析
# └── spatial/ 分割数据用的组织图与 scalefactors
# ├── scalefactors_json.json 缩放系数(注意: 不含 spot_diameter_fullres)
# └── tissue_hires/lowres_image.png 组织图(与顶层相同)
#
# 关键点: 细胞分割没有 tissue_positions, 每个细胞的坐标只能从 *_segmentations.geojson 取。
# ---------------------------------------------------------------------------
library(Seurat)
library(SeuratObject)
library(ggplot2)
## ---- 1. 指定 spaceranger 的输出目录(按需修改) ----
outs_dir <- "Visium_HD_Mouse_Brain_Fixed_Frozen/outs" # spaceranger 的 outs 目录
seg_dir <- file.path(outs_dir, "segmented_outputs") # 分割结果目录
image_dir <- file.path(seg_dir, "spatial") # 组织图 + scalefactors_json.json
## ---- 2. 读入细胞计数矩阵(filtered_feature_cell_matrix) ----
counts <- Read10X(file.path(seg_dir, "filtered_feature_cell_matrix"), gene.column = 2)
obj <- CreateSeuratObject(counts = counts, assay = "Spatial")
message("读入细胞数: ", ncol(obj))
## ---- 3. 读入细胞分割(多边形边界 + 质心) ----
# Read10X_Segmentations 固定读取 <data.dir>/segmented_outputs/<type>_segmentations.geojson
# 返回含 segmentations(多边形) 与 centroids(质心) 两种边界的 VisiumV2 图像对象
v2 <- Read10X_Segmentations(
image.dir = image_dir, # 组织图与 scalefactors 目录
data.dir = outs_dir, # 其下含 segmented_outputs/
image.name = "tissue_lowres_image.png", # 也可用 tissue_hires_image.png
assay = "Spatial",
slice = "cellbin", # 图像名(自定义, 别用特殊字符)
segmentation.type = "cell", # "cell" 或 "nucleus"
compact = TRUE # 多边形只存 sf.data, 省内存
)
## ---- 4. 矩阵与分割取交集后, 把分割图像挂到 Seurat 对象上 ----
common <- intersect(Cells(obj), Cells(v2))
obj <- subset(obj, cells = common)
v2 <- v2[common]
obj[["cellbin"]] <- v2
## ---- 5. 展示分割结果 ----
# (a) 直接画分割多边形轮廓(颜色为 ident)
p1 <- SpatialPlot(obj, image.alpha = 1, image.scale = "lowres",
crop = FALSE, plot_segmentations = TRUE) + theme(aspect.ratio = 1)
# (b) 按指标给每个细胞的多边形上色(更常用), 这里用每个细胞的总 counts
p2 <- SpatialFeaturePlot(obj, features = "nCount_Spatial", image.scale = "lowres",
crop = FALSE, plot_segmentations = TRUE, stroke = NA)
#ggsave("cell_segmentation.png", p2, width = 6, height = 6, dpi = 300)
print(p2)
python 版本：
#!/usr/bin/env python3
############################################################################
# 读取 10x Visium HD 细胞分割结果(geojson)并展示分割多边形
# 依赖: numpy, matplotlib
############################################################################
# ---------------------- spaceranger 输出 outs/ 目录说明 ----------------------
# outs/
# ├── web_summary.html QC 汇总网页(浏览器打开即可看)
# ├── metrics_summary.csv 关键指标表(细胞/spot 数、测序质量、比对率等)
# ├── molecule_info.h5 每个分子的位置与归属(分子级信息)
# ├── barcode_mappings.parquet 2µm bin 与细胞分割 barcode 的对应关系
# ├── feature_slice.h5 / probe_set.csv 探针组(feature)定义 / 探针列表
# ├── cloupe_*.cloupe Loupe Browser 可视化文件
# │
# ├── spatial/ 【顶层图像】全分辨率组织图与对齐信息(此层无坐标表)
# │ ├── tissue_hires_image.png 高分辨率组织图
# │ ├── tissue_lowres_image.png 低分辨率组织图(出图常用)
# │ ├── cytassist_image.tiff CytAssist 原始成像
# │ ├── aligned_fiducials.jpg 基准点对齐图(人工核对用)
# │ ├── aligned_tissue_image.jpg 组织图与载玻片对齐结果
# │ ├── detected_tissue_image.jpg 自动识别出的组织区域
# │ └── final_alignment.json 图像对齐参数
# │
# ├── binned_outputs/ 【方阵(bin)结果】每档 bin 一个子目录
# │ ├── square_002um/ 2µm(最细, 最大, 一般不用)
# │ ├── square_008um/ 8µm
# │ ├── square_016um/ 16µm
# │ └── square_050um/ 50µm
# │ ├── filtered_feature_bc_matrix/ 过滤后表达矩阵(仅组织上的 bin), 10x 三件套:
# │ │ ├── barcodes.tsv.gz bin 名称
# │ │ ├── features.tsv.gz 基因 ID / 名称 / 类型
# │ │ └── matrix.mtx.gz 表达计数(稀疏矩阵)
# │ ├── filtered_feature_bc_matrix.h5 同上, 合并的 HDF5 格式
# │ ├── raw_feature_bc_matrix[.h5] 未过滤(全部 bin, 含背景)
# │ ├── analysis/ Space Ranger 自带聚类/PCA/UMAP/差异分析
# │ └── spatial/
# │ ├── tissue_positions.parquet ★每个 bin 的行列号 + 像素坐标(坐标表)
# │ ├── scalefactors_json.json 图像缩放系数 + spot_diameter_fullres
# │ └── tissue_hires/lowres_image.png
# │
# └── segmented_outputs/ 【细胞分割结果】← 本代码读取这里
# ├── filtered_feature_cell_matrix/ 基于分割细胞的表达矩阵(10x 三件套)
# │ ├── barcodes.tsv.gz 细胞名 cellid_000000001-1
# │ ├── features.tsv.gz
# │ └── matrix.mtx.gz
# ├── filtered_feature_cell_matrix.h5 同上, HDF5 格式
# ├── raw_feature_cell_matrix[.h5] 未过滤的细胞矩阵
# ├── cell_segmentations.geojson ★细胞分割: 每个细胞的 polygon 顶点 + 质心 + id
# ├── nucleus_segmentations.geojson ★细胞核分割(结构同上)
# ├── graphclust_annotated_*.geojson 带 graphclust 聚类注释的分割文件
# ├── cloupe.cloupe Loupe Browser 文件
# ├── analysis/ 分割数据的聚类/PCA/UMAP/差异分析
# └── spatial/ 分割数据用的组织图与 scalefactors
# ├── scalefactors_json.json 缩放系数(注意: 不含 spot_diameter_fullres)
# └── tissue_hires/lowres_image.png 组织图(与顶层相同)
#
# 关键点: 细胞分割没有 tissue_positions, 每个细胞的坐标只能从 *_segmentations.geojson 取。
# ---------------------------------------------------------------------------
import json
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection
from matplotlib.image import imread
## ---- 1. 指定 spaceranger 的输出目录(按需修改) ----
outs_dir = Path("path/to/spaceranger/outs")
seg_dir = outs_dir / "segmented_outputs"
image_dir = seg_dir / "spatial"
geojson = seg_dir / "cell_segmentations.geojson" # 或 nucleus_segmentations.geojson
## ---- 2. 读取 geojson: 每个细胞的多边形顶点 ----
with open(geojson) as f:
gj = json.load(f)
polys = [] # 每个细胞一个 (N,2) 顶点数组
for feat in gj["features"]:
coords = feat["geometry"]["coordinates"] # Polygon: [[[x, y], ...]]
ring = coords[0] # 取外环顶点
if ring:
polys.append(np.asarray(ring, dtype=float))
print("细胞数:", len(polys))
## ---- 3. 读取组织图与缩放系数(fullres 像素 -> 图像显示像素) ----
img = imread(str(image_dir / "tissue_lowres_image.png"))
if img.ndim == 2: # 灰度图转 RGB, 便于叠加
img = np.stack([img] * 3, axis=-1)
sf = json.loads((image_dir / "scalefactors_json.json").read_text())
scale = sf["tissue_lowres_scalef"] # lowres 图相对 fullres 的缩放比
H, W = img.shape[:2]
polys_disp = [p * scale for p in polys] # 多边形换算到 lowres 坐标
## ---- 4. 绘制: 组织图 + 分割多边形(全景 + 局部放大) ----
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
def draw(ax, xlim, ylim, title):
ax.imshow(img) # 背景组织图
pc = PolyCollection(polys_disp, # 叠加分割多边形
facecolor=(0.2, 0.6, 1.0, 0.5), # 半透明填充
edgecolors="black", linewidths=0.05) # 黑色细胞边界
ax.add_collection(pc)
ax.set_xlim(*xlim); ax.set_ylim(*ylim) # y 轴反向(原点在左上)
ax.set_aspect("equal"); ax.axis("off"); ax.set_title(title)
draw(axes[0], (0, W), (H, 0), "cell segmentation (full)")
cx, cy, half = W * 0.42, H * 0.52, 180 # 放大区域(可调)
draw(axes[1], (cx - half, cx + half), (cy + half, cy - half), "zoom")
plt.tight_layout()
plt.savefig("cell_segmentation.png", dpi=300)
plt.show()
- 发表于 1天前
- 阅读 ( 29 )
- 分类：转录组
