空间转录组简介
一、什么是空间转录组？
空间转录组技术是一种可以在组织切片的空间背景下，高通量地测量基因表达的新型技术。它结合了单细胞转录组测序（scRNA-seq）和组织学图像信息，使我们不仅知道 哪些基因被哪些细胞表达，还知道这些细胞位于组织的什么位置。
二、为什么需要空间转录组？
单细胞转录组的局限：
- scRNA-seq 提供的是细胞内的转录信息，但丢失了细胞在组织中的空间信息。
- 对于像
肿瘤微环境、神经系统、胚胎发育等研究，空间位置是非常关键的。
空间转录的优势：
- 将转录数据映射回组织切片，可视化组织结构中基因表达的空间分布。
- 实现 空间维度 + 分子信息 的整合。
三、技术原理与流程
组织切片
↓
在玻片上固定
↓
玻片上预设有空间条形码的探针（如10X Visium）
↓
组织渗透、细胞裂解 → mRNA被探针捕获
↓
cDNA合成 + 建库测序
↓
结合条形码 → 获得每个空间点（spot）对应的转录数据
↓
图像 + 基因表达数据整合 → 空间转录图谱
四、空间转录组的数据结构
空间转录组常用的数据结构是：
- 每个点（spot）对应一组：
- 空间坐标 (x, y)
- 某些图像信息（染色图像）
- 一个转录组（通常是counts矩阵）
通常用 AnnData 对象（scanpy, squidpy）来组织数据，包括：
adata.obs # spot的meta信息（如位置、聚类、组织类型）
adata.var # 基因信息
adata.X # 表达矩阵（稀疏矩阵）
adata.obsm['spatial'] # 空间坐标 (x, y)
五、空间转录组的数据分析步骤
- 图像预处理: H&E图像/荧光图像，与空间点匹配
- 表达量标准化: 类似 scRNA-seq，如 log-normalize
- 空间聚类: 识别不同区域（如肿瘤核心 vs 边缘）
- 差异表达分析: 比较不同空间区域的基因表达差异
- 空间变异基因识别: 哪些基因表达呈现空间趋势
- 图像特征融合: 将图像特征与基因表达整合分析（如通过 squidpy）
- 注释与解读: 结合已有标注或 scRNA-seq 进行空间注释
六、常用分析工具
Scanpy: 基础分析框架，适用于空间转录组Squidpy空间数据专用分析（图像 + 空间邻近）Seurat V5: R 生态下空间转录组支持（SpatialFeaturePlot）stLearn: 深度学习 + 空间转录整合分析Giotto: 空间表达模式识别、可视化等
七、应用场景
- 肿瘤微环境分析: 找到免疫细胞、肿瘤细胞在组织中的空间分布
- 神经科学研究: 映射大脑区域基因表达
- 胚胎发育研究: 跟踪发育阶段基因表达空间变化
- 器官再生研究: 分析再生过程中不同区域活跃基因
- 疾病机制研究: 对比正常与病变组织的空间表达差异
八、单样本【点亮】最小工作流
假设，我们现在已经有10X Visium 的解析数据（outs），并希望点亮空间表达图，从 10x Visium 的 outs/ 读取数据 → “点亮”更鲜艳的空间图 → 单独点亮 AURKB（或任意基因）→（可选）多样本合并与统一绘图 → 批量保存 PNG 与 .h5ad。
保存为 st_visium_pointlight.py
# -*- coding: utf-8 -*-
"""
st_visium_pointlight.py
一体化：读取 10x Visium outs → 更鲜艳的空间点亮 → 单独点亮 AURKB → 多样本可选合并
使用示例：
单样本：
python st_visium_pointlight.py --single "D:\Visium\SampleA\outs" --outdir "D:\Visium\SampleA\plots"
多样本（根目录下每个子目录为样本，子目录内有 outs/）：
python st_visium_pointlight.py --multi-root "D:\Visium" --outdir "D:\Visium\plots" --merge
指定要点亮的基因（默认 AURKB）：
python st_visium_pointlight.py --single ".../outs" --gene AURKB --outdir ".../plots"
注意：
- 需要已安装：scanpy>=1.10, squidpy>=1.4, anndata, numpy, pandas, scipy, matplotlib
- 仅“点亮”与常规分析不需要 GPU；若后续用到 cell2location 等再配置 CUDA。
"""
import os
import sys
import glob
import argparse
from typing import List, Optional
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import anndata as ad
import scanpy as sc
import squidpy as sq
# ---------- 画图全局参数（更清晰） ----------
sc.set_figure_params(
dpi=140, dpi_save=300,
facecolor="white",
transparent=True
)
sc.settings.verbosity = 3 # 日志详细度
# ---------- 实用函数 ----------
def ensure_dir(path: str):
os.makedirs(path, exist_ok=True)
return path
def get_library_id_from_path(path: str) -> str:
"""
根据 outs 上级目录名作为 library_id，例如 D:/Visium/SRR22352787/outs → SRR22352787
"""
p = os.path.normpath(path)
parent = os.path.basename(os.path.dirname(p))
return parent if parent else "library"
def find_gene_name(adata: ad.AnnData, target: str = "AURKB") -> Optional[str]:
"""
在 adata 内尽力找到指定基因的可用变量名（兼容大小写、Aurkb、小鼠/人名、符号列等）
返回：匹配到的 var_name 或 None
"""
candidates = {target, target.upper(), target.lower(), "Aurkb", "aurkb", "AURKB"}
# 1) var_names 直接匹配（不区分大小写）
upper_map = {v.upper(): v for v in adata.var_names}
for c in list(candidates):
if c.upper() in upper_map:
return upper_map[c.upper()]
# 2) 常见符号列
for col in ["gene_symbols", "gene_symbol", "feature_name", "gene_name", "symbols", "Symbol", "SYMBOL"]:
if col in adata.var.columns:
col_upper = adata.var[col].astype(str).str.upper()
hit = np.where(col_upper.values == target.upper())[0]
if hit.size > 0:
return adata.var_names[hit[0]]
# 再试候选
for c in list(candidates):
hit2 = np.where(col_upper.values == c.upper())[0]
if hit2.size > 0:
return adata.var_names[hit2[0]]
return None
def vivid_palette(n: int, s: float = 0.9, v: float = 0.98, seed: int = 7) -> List[str]:
"""
生成 n 个高饱和度、对比强的 hex 颜色（HSV 均匀取色并打乱）
"""
import colorsys
rng = np.random.default_rng(seed)
hues = np.linspace(0, 1, n, endpoint=False)
rng.shuffle(hues)
hexes = []
for h in hues:
r, g, b = colorsys.hsv_to_rgb(h, s, v)
hexes.append('#%02x%02x%02x' % (int(r*255), int(g*255), int(b*255)))
return hexes
def plot_gene_vivid_single(
adata: ad.AnnData,
gene: str = "AURKB",
low: float = 2,
high: float = 98,
cmap: str = "turbo",
use_raw: bool = False,
spot_size: float = 1.35,
alpha_img: float = 0.4,
img_cmap: str = "gray",
img_key: str = "hires",
title_suffix: str = "vivid",
save: Optional[str] = None
):
"""
在空间图上单独点亮一个基因（更鲜艳）：
- 百分位裁剪 vmin/vmax 提升对比
- 灰底 + 降低底图透明度
- 放大 spot
"""
g = find_gene_name(adata, gene)
if g is None:
print(f"[WARN] 未在此对象中找到基因 {gene}，请检查基因名/物种/注释。")
return
# 获取表达向量
X = (adata.raw[:, g].X if (use_raw and adata.raw is not None) else adata[:, g].X)
if hasattr(X, "toarray"):
vec = X.toarray().ravel()
else:
vec = np.asarray(X).ravel()
# 计算百分位裁剪范围
if vec.size == 0:
vmin = vmax = None
else:
pos = vec[vec > 0]
vmin = float(np.percentile(pos, low)) if pos.size > 0 else float(np.percentile(vec, low))
vmax = float(np.percentile(vec, high))
if vmin >= vmax:
vmin = vmax = None
sc.pl.spatial(
adata,
img_key=img_key,
color=g,
cmap=cmap, # 高饱和 colormap：'turbo' / 'Spectral_r' / 'plasma' / 'inferno'
vmin=vmin, vmax=vmax,
alpha_img=alpha_img,
img_cmap=img_cmap, # 灰度底图
spot_size=spot_size,
frameon=False,
title=f"{g} — {title_suffix}",
save=save # 例如 "/AURKB_vivid.png"
)
def plot_clusters_vivid(
adata: ad.AnnData,
cluster_key: str = "leiden",
spot_size: float = 1.35,
alpha_img: float = 0.4,
img_cmap: str = "gray",
img_key: str = "hires",
save: Optional[str] = None,
seed: int = 7
):
"""
用高饱和随机调色板对类别（如 leiden）进行空间上色，并降低底图干扰
"""
if cluster_key not in adata.obs:
print(f"[WARN] {cluster_key} 不在 adata.obs 中，跳过聚类上色。")
return
clusters = list(map(str, adata.obs[cluster_key].astype(str).unique()))
try:
clusters_sorted = sorted(clusters, key=lambda x: int(x))
except Exception:
clusters_sorted = sorted(clusters)
adata.uns[f"{cluster_key}_colors"] = vivid_palette(len(clusters_sorted), s=0.95, v=0.98, seed=seed)
sc.pl.spatial(
adata,
img_key=img_key,
color=[cluster_key],
alpha_img=alpha_img,
img_cmap=img_cmap,
spot_size=spot_size,
frameon=False,
title=f"Spatial {cluster_key} (vivid)",
save=save # 如 "/clusters_vivid.png"
)
def basic_preprocess(adata: ad.AnnData, batch_key: Optional[str] = None) -> ad.AnnData:
"""
标准预处理：归一化→log1p→HVG（可带批次）→scale→PCA→邻居→Leiden→UMAP
返回副本（避免原地覆盖）
"""
a = adata.copy()
sc.pp.calculate_qc_metrics(a, inplace=True)
sc.pp.normalize_total(a, target_sum=1e4)
sc.pp.log1p(a)
a.raw = a
if batch_key is None:
sc.pp.highly_variable_genes(a, n_top_genes=3000, flavor="seurat_v3")
else:
sc.pp.highly_variable_genes(a, n_top_genes=3000, flavor="seurat_v3", batch_key=batch_key)
a = a[:, a.var['highly_variable']].copy()
sc.pp.scale(a, max_value=10)
sc.tl.pca(a, n_comps=50, svd_solver="arpack")
sc.pp.neighbors(a, n_neighbors=15, n_pcs=30)
sc.tl.leiden(a, resolution=0.6, key_added="leiden")
sc.tl.umap(a)
return a
def read_visium_outs(outs_path: str, library_id: Optional[str] = None) -> ad.AnnData:
"""
读取单个样本的 outs/ 目录并返回 AnnData；自动 var 去重
"""
if library_id is None:
library_id = get_library_id_from_path(outs_path)
a = sc.read_visium(
path=outs_path,
count_file="filtered_feature_bc_matrix.h5", # 若没有 .h5，则去掉该参数，read_visium 会读 mtx 目录
load_images=True,
library_id=library_id
)
a.var_names_make_unique()
a.obs["batch"] = library_id
return a
# ---------- 单样本流程 ----------
def run_single(outs_path: str, outdir: str, gene: str = "AURKB", img_key: str = "hires"):
ensure_dir(outdir)
sc.settings.figdir = outdir # 让 sc.pl.* 的 save= 生效到该目录
print(f"[INFO] 读取单样本：{outs_path}")
adata = read_visium_outs(outs_path)
# 1) 基础“点亮”：H&E + Spots（不着色，仅检查）
sc.pl.spatial(adata, img_key=img_key, color=None, alpha_img=1.0, spot_size=1.1,
title="H&E + Spots (check)", save="/he_spots.png")
# 2) 更鲜艳的基因点亮（AURKB 或自定义）
plot_gene_vivid_single(
adata, gene=gene, low=2, high=98, cmap="turbo", use_raw=False,
spot_size=1.4, alpha_img=0.35, img_cmap="gray", img_key=img_key,
title_suffix="vivid", save=f"/{gene}_vivid.png"
)
# 3) 预处理 + 聚类
a_proc = basic_preprocess(adata)
# 4) UMAP + 空间上色（鲜艳调色板）
sc.pl.umap(a_proc, color=["leiden"], legend_loc="on data", frameon=False, save="/umap_leiden.png")
plot_clusters_vivid(a_proc, cluster_key="leiden", spot_size=1.35, alpha_img=0.4,
img_cmap="gray", img_key=img_key, save="/clusters_vivid.png")
# 5) 空间变异基因（Moran's I） + “点亮”Top 基因
sq.gr.spatial_neighbors(a_proc, coord_type="grid")
sq.gr.spatial_autocorr(a_proc, mode="moran", n_perms=100)
if "moranI" in a_proc.uns:
top_genes = a_proc.uns["moranI"].sort_values("I", ascending=False).head(6).index.tolist()
if len(top_genes) > 0:
for g in top_genes:
plot_gene_vivid_single(a_proc, gene=g, low=2, high=98, cmap="turbo",
use_raw=False, spot_size=1.3, alpha_img=0.35,
img_cmap="gray", img_key=img_key,
title_suffix="Top Moran", save=f"/moran_{g}.png")
# 6) 保存对象
a_proc.write(os.path.join(outdir, "processed.h5ad"))
print(f"[DONE] 单样本完成，结果在：{outdir}")
# ---------- 多样本流程 ----------
def discover_samples(root: str) -> List[str]:
"""
在 root 下寻找形如 */outs 目录的样本路径列表
"""
candidates = []
for p in glob.glob(os.path.join(root, "*")):
outs = os.path.join(p, "outs")
if os.path.isdir(outs) and os.path.exists(outs):
candidates.append(outs)
return sorted(candidates)
def run_multi(root: str, outdir: str, gene: str = "AURKB", merge: bool = False, img_key: str = "hires"):
ensure_dir(outdir)
sc.settings.figdir = outdir
outs_list = discover_samples(root)
if not outs_list:
print(f"[ERROR] 在 {root} 下未发现任何包含 outs/ 的样本目录。")
sys.exit(1)
adatas = []
print(f"[INFO] 发现样本数：{len(outs_list)}")
for outs in outs_list:
sid = get_library_id_from_path(outs)
print(f" - 读取 {sid}")
a = read_visium_outs(outs, library_id=sid)
# 单片快速“点亮”该基因
plot_gene_vivid_single(
a, gene=gene, low=2, high=98, cmap="turbo",
use_raw=False, spot_size=1.35, alpha_img=0.35,
img_cmap="gray", img_key=img_key, title_suffix=f"{sid} vivid",
save=f"/{sid}_{gene}_vivid.png"
)
# 单片聚类并空间上色
a_proc = basic_preprocess(a)
plot_clusters_vivid(a_proc, cluster_key="leiden", spot_size=1.35, alpha_img=0.4,
img_cmap="gray", img_key=img_key, save=f"/{sid}_clusters_vivid.png")
# 保存每片处理结果
subdir = ensure_dir(os.path.join(outdir, sid))
a_proc.write(os.path.join(subdir, f"{sid}_processed.h5ad"))
adatas.append(a)
if merge and len(adatas) >= 2:
print("[INFO] 开始合并多样本……")
A = ad.concat(adatas, join="outer", label="batch",
keys=[a.obs["batch"][0] for a in adatas], fill_value=0)
A.var_names_make_unique()
# 合并对象预处理（带批次 HVG）
A_proc = basic_preprocess(A, batch_key="batch")
# 联合 UMAP
sc.pl.umap(A_proc, color=["batch", "leiden"], wspace=0.4, frameon=False, save="/merged_umap.png")
# 分样本分别空间上色（鲜艳）
for sid in sorted(A_proc.obs["batch"].unique()):
a_s = A_proc[A_proc.obs["batch"] == sid].copy()
plot_clusters_vivid(a_s, cluster_key="leiden", spot_size=1.35, alpha_img=0.4,
img_cmap="gray", img_key=img_key, save=f"/{sid}_clusters_vivid_merged.png")
# 每片再单独点亮目标基因（来自合并对象）
plot_gene_vivid_single(
a_s, gene=gene, low=2, high=98, cmap="turbo",
use_raw=False, spot_size=1.35, alpha_img=0.35,
img_cmap="gray", img_key=img_key, title_suffix=f"{sid} vivid (merged)",
save=f"/{sid}_{gene}_vivid_merged.png"
)
# 空间统计（多片会自动按片构图）
sq.gr.spatial_neighbors(A_proc, coord_type="grid")
sq.gr.spatial_autocorr(A_proc, mode="moran", n_perms=100)
# 保存合并对象
A_proc.write(os.path.join(outdir, "visium_merged.h5ad"))
print(f"[DONE] 多样本合并完成，结果在：{outdir}")
else:
if merge and len(adatas) < 2:
print("[WARN] 样本不足 2 个，跳过合并。")
print(f"[DONE] 多样本单片处理完成，结果在：{outdir}")
# ---------- CLI ----------
def parse_args():
ap = argparse.ArgumentParser(description="10x Visium 空间转录组更鲜艳点亮 & AURKB 专用脚本")
mode = ap.add_mutually_exclusive_group(required=True)
mode.add_argument("--single", type=str, help="单样本 outs 目录路径（包含 spatial/ 与矩阵）")
mode.add_argument("--multi-root", type=str, help="多样本根目录（其子目录内含 outs/）")
ap.add_argument("--outdir", type=str, required=True, help="输出目录（保存 PNG 与 .h5ad）")
ap.add_argument("--gene", type=str, default="AURKB", help="要单独点亮的基因名（默认 AURKB）")
ap.add_argument("--img-key", type=str, default="hires", help="使用的底图关键字（hires/lowres）")
ap.add_argument("--merge", action="store_true", help="多样本时是否执行合并分析")
return ap.parse_args()
def main():
args = parse_args()
if args.single:
run_single(outs_path=args.single, outdir=args.outdir, gene=args.gene, img_key=args.img_key)
else:
run_multi(root=args.multi_root, outdir=args.outdir, gene=args.gene,
merge=args.merge, img_key=args.img_key)
if __name__ == "__main__":
main()