01 - Squidpy Visium Hello

目标：用 Squidpy 跑通一个 10x Visium H&E 示例数据，把「空间坐标 + 表达矩阵 + 组织图像」连起来。
这不是追求漂亮图，而是练三件事：
看懂 AnnData 里空间信息放在哪里。
构建空间邻接图，并做邻域富集分析。
把结果整理成后续 Agent tool 能返回的结构化报告。

第一次运行会下载 Squidpy 官方示例数据，体积约数百 MB。

总结

空间转录组是在每个空间 spot 上测基因表达。
先根据表达模式把 spot 聚成 cluster(简单说就是将基因表达特征相似的spot归位同一类)，再结合真实空间坐标看这些 cluster 在组织里怎么分布。
空间相邻富集看不同 cluster 是否异常靠近。
空间可变基因看某些基因的表达是否在空间上成片分布。
这些结果都是候选线索，不是机制验证。

自己理解：数据形状 spot × gene的意思是测了spot个区域，每个区域有gene种基因。之后对所有spot进行分类，分成若干个cluster。空间相邻富集是分析出每两个cluster之间是否比随机概率更容易相邻。空间可变基因分析是判断某些基因的表达是否与空间区域有关。

空间相邻富集对AI制药的作用：判断某些细胞群/组织区域是否在空间上形成了特殊关系。
空间可变基因对AI制药的作用：把“高表达基因”升级成“在关键空间区域高表达的候选靶点”。

1. 导入依赖

这里用 scanpy + squidpy：

scanpy 负责单细胞/表达矩阵的常规处理。
squidpy 负责空间坐标、空间邻接图、空间可视化和空间统计。

from pathlib import Path

# 【依赖分工】
# pathlib：只负责路径管理，避免在代码里到处手写字符串路径。
# matplotlib/pandas：辅助画图和表格整理，不承担空间组学算法。
# scanpy：处理 AnnData、表达矩阵、基础单细胞流程。
# squidpy：在 scanpy 的 AnnData 之上增加空间组学能力，例如空间坐标、邻接图、邻域富集。
import matplotlib.pyplot as plt
import pandas as pd
import scanpy as sc
import squidpy as sq

# 【工程习惯】
# notebook 很容易把图片、缓存、临时文件散落到当前目录。
# 先统一定义输出目录，后面封装成 Agent tool 时也能把所有 artifact 路径集中返回。
OUTPUT_DIR = Path("outputs")
OUTPUT_DIR.mkdir(exist_ok=True)

# 【可复现性】
# 固定图像大小和 dpi，不是为了美观，而是为了让每次运行产出的图尺寸稳定。
# 将来如果 Agent 把图路径写进报告，稳定尺寸能减少“同一段代码每次输出都不一样”的调试成本。
sc.settings.set_figure_params(figsize=(6, 5), dpi=120)

# 【环境观察 Observation】
# 先打印库版本。空间组学库 API 变化较快，版本号是排查问题的第一条证据。
print("scanpy", sc.__version__)
print("squidpy", sq.__version__)

2. 读取 Visium 示例数据

Squidpy 官方教程使用 visium_hne_adata()，这是一个预处理好的 10x Genomics Visium H&E 数据。

重点看三类对象：

位置	含义
`adata.X`	spot × gene 表达矩阵
`adata.obsm["spatial"]`	每个 spot 的组织切片坐标
`adata.obs["cluster"]`	官方预注释的空间/表达聚类标签

# 【输入数据】
# 这里用 Squidpy 官方自带的 10x Visium H&E 示例数据。
# 它适合入门，因为表达矩阵、空间坐标、组织图像和 cluster 标签都已经整理在 AnnData 里。
# 第一次运行可能会下载数据；这一步相当于未来 tool 的 data loader。
adata = sq.datasets.visium_hne_adata()

# 【基因名去重】
# AnnData 要求 var_names 最好唯一。重复基因名会让后续按基因取表达值时产生歧义。
# 这一步不是生物学分析，而是数据工程层面的“输入清洗”。
adata.var_names_make_unique()

# 【分组字段】
# cluster 是示例数据里已有的区域/簇标签。
# 后面的邻域富集不是直接比较每个 spot，而是比较“不同 cluster 是否更常相邻”。
cluster_key = "cluster"

# 【类型转换】
# category 类型会告诉 scanpy/squidpy：这是离散分组标签，不是连续数值。
# 如果不转成 category，某些绘图或统计函数可能无法按“类别”正确处理。
adata.obs[cluster_key] = adata.obs[cluster_key].astype("category")

# 【先观察输入，不急着分析】
# Agent tool 设计里，这些打印信息未来应变成结构化返回值：n_spots、n_genes、obs/obsm keys。
# 原因是：模型不能只看图片，必须知道自己分析的数据规模和字段是否符合预期。
print("数据形状 spot × gene:", adata.shape)
print("obsm keys:", list(adata.obsm.keys()))
print("obs columns 示例:", list(adata.obs.columns[:10]))
print("cluster 数量:", adata.obs[cluster_key].nunique())

3. 空间可视化

这一步不是 UMAP。它直接使用真实组织坐标，把每个 spot 的 cluster 画回切片位置。

要记住：

UMAP = 表达相似性的二维压缩
spatial plot = 真实组织坐标

# 【空间图，不是 UMAP】
# spatial_scatter 会读取 adata.obsm["spatial"] 里的真实组织坐标。
# 如果 adata.uns["spatial"] 里有组织切片图像，Squidpy 也可以把 spot 叠到图片背景上。
# 这张图回答的问题是：“不同 cluster 在组织切片的哪里？”
# 它不回答：“哪些 spot 的表达最相似？”后者才是 UMAP/PCA 这类降维图的问题。
sq.pl.spatial_scatter(adata, color=cluster_key, title="Visium H&E clusters")

# 【artifact 记录】
# 现在图片显示在 notebook 里；未来封装成 Agent tool 时，最好把图片保存到文件并返回路径。
# 这样 LLM 负责解释，Python tool 负责稳定产出证据，职责更清楚。
figure_paths = {
    "spatial_cluster": "notebook_inline_spatial_cluster",
}

4. 构建空间邻接图

空间分析的关键不是“点在图上好看”，而是先定义谁和谁是邻居。

Squidpy 1.8 以后推荐使用更明确的邻接构建函数。Visium 是规则网格，所以这里用 spatial_neighbors_grid()。

后续分析会用到：

adata.obsp["spatial_connectivities"]
adata.obsp["spatial_distances"]

# 【构建空间邻接图】
# 空间分析的核心不是“ spot 有坐标”这么简单，而是要定义谁和谁是邻居。
# Visium spot 排列近似六边形网格，所以这里使用 spatial_neighbors_grid。
# 这一步会在 adata.obsp 中写入空间连接矩阵和距离矩阵，供后续统计函数复用。

# n_neighs=6：Visium 六边形网格中，一个 spot 理论上最多有 6 个一圈邻居。
# n_rings=1：只看第一圈邻居；如果调成 2，就会把更远的第二圈 spot 也纳入空间邻域。
# key_added="spatial"：让输出字段命名为 spatial_connectivities / spatial_distances，便于后续函数默认读取。
sq.gr.spatial_neighbors_grid(
    adata,
    spatial_key="spatial",
    n_neighs=6,
    n_rings=1,
    key_added="spatial",
)

# 【邻接矩阵解释】
# adata.obsp["spatial_connectivities"] 是 spot × spot 的稀疏矩阵。
# 某个位置非 0，表示两个 spot 被定义为空间邻居。
# getnnz(axis=1) 统计每个 spot 有多少邻居；平均值可作为邻接图是否异常的粗略 QC 指标。
connectivities = adata.obsp["spatial_connectivities"]
avg_neighbors = connectivities.getnnz(axis=1).mean()

# 【Observation】
# 对 Agent 来说，这些是工具调用后的观察结果。
# 如果平均邻居数极低或极高，后续邻域富集结论就不可靠，必须进入 quality_flags。
print("obsp keys:", list(adata.obsp.keys()))
print("平均空间邻居数:", round(float(avg_neighbors), 2))

5. 邻域富集分析

邻域富集回答的是：

某两个 cluster 是否比随机情况更常挨在一起？

注意：这只是空间共现/邻近的统计结果，不等于真实生物学互作。

# 【邻域富集分析】
# nhood_enrichment 的问题不是“两个 cluster 是否真的发生细胞通讯”，而是：
# 在已经定义好的空间邻接图上，cluster A 和 cluster B 是否比随机置换更常相邻。
# 因此它输出的是空间共现统计线索，不是机制证明。

# n_perms：随机置换次数。次数越大，统计更稳，但运行更慢。
# 学习阶段用 100 是为了快速看懂流程；正式报告建议提高，并记录 seed 方便复现。
# seed：固定随机种子，避免每次运行 zscore 有轻微变化。
sq.gr.nhood_enrichment(
    adata,
    cluster_key=cluster_key,
    n_perms=100,
    seed=42,
    show_progress_bar=False,
)

# 【可视化】
# 热图颜色通常表示 z-score：越高说明两个 cluster 比随机情况更常相邻。
# 但“更常相邻”仍然只是空间统计，不要写成“发生了相互作用”。
sq.pl.nhood_enrichment(adata, cluster_key=cluster_key, title="Neighborhood enrichment")

# 【结果落点】
# Squidpy 会把结果写进 adata.uns，而不是直接返回一个表。
# 这体现了 AnnData 的工作方式：分析步骤不断往 adata 里追加结果，adata 就是整个 pipeline 的状态容器。
nhood_result = adata.uns[f"{cluster_key}_nhood_enrichment"]
print("邻域富集结果字段:", list(nhood_result.keys()))
print("zscore 矩阵形状:", nhood_result["zscore"].shape)

6. 空间可变基因：Moran's I

空间可变基因回答的是：

这个基因的表达是否呈现空间聚集/空间模式？

如果基因 A 的高表达 spot 在空间上成片聚集，而不是随机散落，

说明基因 A 可能和某个空间区域 / 组织结构 / 生物学 domain 有关。它仍然是统计线索，不是靶点验证。

邻域富集说的是"谁和谁挨着"（cluster 之间），空间可变基因说的是"一个基因在哪里高"（基因×空间）。

# 【空间可变基因】
# Moran's I 衡量一个基因表达是否呈现空间自相关：
# 如果高表达 spot 倾向于挨在一起，Moran's I 往往更高。
# 这一步回答的是“表达是否有空间模式”，不是“这个基因是否是有效药物靶点”。

# 示例数据是小鼠脑组织，这里挑几个常见脑组织相关基因做入门观察。
# 将来封装成 tool 时，这个列表应来自用户入参 target_genes。
candidate_genes = ["Mbp", "Hpca", "Ttr", "Snap25", "Pcp4"]

# 【输入护栏】
# 真实数据里用户给的基因名可能不存在，比如大小写不一致、人鼠同源基因混淆、面板数据未检测该基因。
# 这里先过滤存在的基因，避免整个流程因为一个缺失基因直接失败。  
genes = [gene for gene in candidate_genes if gene in adata.var_names]

# 【异常处理】
# 如果一个目标基因都不存在，就明确报错。
# 这比静默返回空结果更适合 Agent tool，因为 LLM 需要知道失败原因，才能给用户解释或换参数。
if not genes:
    raise ValueError("候选基因都不在 adata.var_names 中，请换一组 target_genes。")

# 【统计计算】
# mode="moran" 使用 Moran's I。
# n_perms=None 表示这里先不做置换检验，学习阶段只看空间自相关分数本身。
# 如果要写正式报告，应加入置换检验或多重检验校正信息。
sq.gr.spatial_autocorr(
    adata,
    genes=genes,
    mode="moran",
    n_perms=None,
    show_progress_bar=False,
)

# 【结果排序】
# adata.uns["moranI"] 是 Squidpy 写回的结果表。
# 按 I 从高到低排序，方便找空间模式最明显的候选基因。
moran_table = adata.uns["moranI"].sort_values("I", ascending=False)
display(moran_table)

# 【把统计结果投回真实空间】
# 只看表格不知道表达区域在哪里，所以把 Moran's I 最高的基因画回组织坐标。
# 这一步是空间组学的关键：统计线索必须回到 tissue context 里解释。
top_gene = str(moran_table.index[0])
sq.pl.spatial_scatter(adata, color=top_gene, title=f"Spatial expression: {top_gene}")

7. 生成结构化质量报告

这一步开始向 Agent tool 靠拢：不要只把图留在 notebook 里，而要返回机器和人都能读的结构化结果。

# 【质量护栏 quality_flags】
# Agent 不能像湿实验一样立刻验证一个空间组学结论是否真实。
# 所以 tool 必须把“不确定性”和“不能过度解释的地方”显式返回给 LLM。
# 这些 flags 不是报错，而是防止模型把统计线索说成生物学定论。
quality_flags = []

# 护栏 1：Visium 是 spot 级，不允许直接当单细胞解释。
# 这是平台层面的限制，和代码是否跑通无关。
quality_flags.append("Visium 是 spot 级数据；一个 spot 可能混合多个细胞，不能直接解释成单细胞结论。")

# 护栏 2：空间邻居数异常时，提示邻接图可能有问题。
# 平均邻居数过低：空间图太稀，邻域统计不稳定。
# 平均邻居数过高：邻域定义太宽，可能把本不相邻的区域混在一起。
if avg_neighbors < 3:
    quality_flags.append("平均邻居数偏低，空间邻接图可能过稀。")
if avg_neighbors > 12:
    quality_flags.append("平均邻居数偏高，空间邻接图可能过密。")

# 护栏 3：邻域富集只能说明空间共现，不能说明真实互作。
# 这是把“统计结果”和“生物学机制”分开的关键，后续写报告必须保留。
quality_flags.append("邻域富集只能说明空间共现/邻近，不能证明细胞间真实互作。")

# 护栏 4：Moran's I 只能说明空间模式，不能直接推出药物靶点有效。
# 对 AI 制药来说，它最多提供候选靶点/候选区域线索，不能替代实验验证。
quality_flags.append("空间可变基因是候选线索，不等于靶点验证。")

# 【结构化输出】
# 这里开始按 tool schema 的思路组织结果：
# - 基础元数据：platform、n_spots、n_genes
# - 分析参数：cluster_key、spatial_key
# - 核心结果：avg_spatial_neighbors、moran_top_gene、moran_table
# - 可解释性材料：quality_flags、figure_paths
# 这样 LLM 不需要从 notebook 输出里猜信息，可以直接消费这个 dict 生成报告。
quality_report = {
    "platform": "10x Visium H&E",
    "n_spots": int(adata.n_obs),
    "n_genes": int(adata.n_vars),
    "cluster_key": cluster_key,
    "n_clusters": int(adata.obs[cluster_key].nunique()),
    "spatial_key": "spatial",
    "avg_spatial_neighbors": round(float(avg_neighbors), 2),
    "moran_top_gene": top_gene,
    "moran_table": moran_table.reset_index().rename(columns={"index": "gene"}).to_dict(orient="records"),
    "quality_flags": quality_flags,
    "figure_paths": figure_paths,
}

# 【最终 Observation】
# notebook 里直接显示 dict；封装成函数时，这一行应改成 return quality_report。
quality_report

8. 这一节要带走什么

你现在已经跑通了空转最小闭环：

读取 Visium 数据 -> 空间可视化 -> 空间邻接图 -> 邻域富集 -> 空间可变基因 -> 质量报告

这套流程后面可以封装成一个 tool：

analyze_spatial_transcriptomics(adata_path, platform, cluster_key, target_genes) -> dict

重点不是 Agent 会不会调用 Squidpy，而是 tool 的输出必须带护栏：

Visium 不能说成单细胞结果。
邻域富集不能说成真实互作。
Moran's I 不能说成靶点验证。
notebook 里的图要配合结构化结果一起返回。