空间转录组数据拿到手之后最常被问到的一个问题就是某个基因的表达为什么会在空间上呈现梯度某个细胞群落为什么只在特定区域出现这些问题本质上都指向同一个分析目标——细胞的“空间关系”。我以前在课程里第一次接触 mistyR 的时候以为它只是又一个相关性分析工具后来真正用 Python 版本mistyPy跑完一轮数据才发现它解决的是普通相关系数完全回答不了的“多尺度拆解”问题。这篇文章就是把课程复习时整理的思路结合我自己跑数据的经验重新梳理了一遍从原理到代码到坑点都过一遍给打算做空间上下文分析的同学一个能直接上手的参考。1. 为什么要分析细胞的“空间关系”先搞清楚问题的边界1.1 拿到空间坐标之后数据分析不应该还停留在“找差异基因”做过单细胞转录组的人应该都有体会单细胞数据最大的遗憾是细胞在组织里的位置信息丢了。空间转录组技术比如 Visium、Slide-seq、Xenium 这类把“表达量”和“空间位置”对应起来数据维度一下子丰富了很多。但很多拿到空间数据的同学第一反应仍然是老一套找 marker、做聚类、跑差异表达。这些当然要做但空间数据最独特的信息——细胞在组织里的邻居是谁、处于什么微环境、周围信号的强弱如何——往往被忽略了。我在实际项目里就遇到过这种情况某个基因在整体表达量上没有明显差异但把它按空间位置画出来你会发现它只在肿瘤边缘的一个窄带里高表达而周围一圈细胞的配体基因恰好也在那里活跃。这种“位置决定的表达状态”如果只看全局表达均值根本看不出来。要回答这类问题就得把每个 spot或每个细胞的“空间背景”作为解释变量纳入模型看看它对目标基因的表达有没有贡献、贡献有多大。1.2 普通相关性分析为什么在这里不够用很多人第一个想法是直接算两个基因在空间上的相关系数不就行了吗比如 A 基因高表达的地方 B 基因也高表达就说明它们有空间共定位。这个思路在简单场景下能用但问题也不小。第一相关系数只描述“两者同步变化”无法区分因果关系和混杂因素。B 基因的高表达可能不是因为 A 基因而是因为第三个基因 C 在这个区域整体活性高。第二相关性分析没有“空间尺度”的概念。两个基因可能在相邻 spot 之间相互作用也可能在跨几个 spot 的较大范围内同步变化这两种关系的生物学含义完全不同。第三一个 spot 的基因表达从来不是被单个因子决定的而是自身调控、邻近细胞通讯、远处组织信号共同作用的结果。普通相关性分析给不出这种“多因子贡献分解”的答案。mistyR以及它的 Python 版本 mistyPy正是针对这个问题设计的它把目标基因的表达变异拆成几个空间尺度的贡献然后告诉你“自身基线解释了多少、邻域解释了多少、更大范围的空间背景解释了多少”。1.3 空间上下文spatial context到底是什么misty 的核心概念是“空间上下文”。简单理解每个 spot 的基因表达不仅由它自己的分子状态决定还受它周围一定范围内的其他 spot 影响。这个“周围”不是一个距离而是分层级的。我当时学这个的时候觉得特别像理解“人际关系”一个人现在的状态有一部分是自己长期以来的习惯自身基线有一部分是身边亲密朋友的影响物理相邻的 spot还有一部分是更大的社区氛围、远程协作或整个环境趋势的影响更大范围内的信号。这样分层之后你才能说清楚“到底是身边人把你带偏了还是整个大环境的问题”。2. mistyR 与 mistyPy 的多视图建模思想2.1 三种视图的直觉理解与数学表达misty 把空间上下文划分为三种视图view。intraview 代表 spot 自身的表达谱也就是作为基线的那部分juxtaview 代表物理上相邻的 spot 对该 spot 的影响通常定义为直接邻居例如距离最近的 K 个 spotparaview 代表更大范围内 spot 的贡献通过计算每个 spot 与周围一定距离内所有其他 spot 的平均表达来定义。用数学一点的表述就是对于目标基因 g 在 spot i 的表达 x_g,i模型写作x_g,i beta_intra * x_intra,i beta_juxta * x_juxta,i beta_para * x_para,i epsilon_i这里的 x_intra 是同一 spot 内其他基因的表达向量x_juxta 是相邻 spot 的表达汇总x_para 是远处 spot 的表达汇总。模型通过正则化线性回归可理解为加了惩罚项的线性拟合估计这些系数然后计算每个视图对方差的解释比例。这里有一个容易误解的点misty 并不是预测“某个基因在所有 spot 中的表达值”而是尝试用一个 spot 的“空间上下文”特征去解释该 spot 中目标基因表达的变异。所以输出的不是一张简单的相关性热图而是一套“方差分解表”。2.2 R2 分解、gain 与基因重要性三个最核心的输出第一次跑完 mistyPy我看着输出结果里的 R2、gain、importances 这几个字段有点懵后来才算彻底理清楚。R2 分解是最直观的它告诉你目标基因的表达变异中每个视图分别贡献了多少比例。比如 A 基因的 intraview 解释率是 20%juxtaview 是 35%paraview 是 10%。这意味着 A 基因的表达变异主要受邻近 spot 影响而不是自身基线决定的。gain 描述的是“增加某个视图后模型预测能力的提升幅度”。它衡量某个视图的增量价值。如果一个视图的 gain 很小说明去掉这个视图对预测几乎没影响。importance 则对应具体基因的重要性。它综合了回归系数和表达变异反映了“哪些基因的空间分布模式对目标基因最有解释力”。这比单纯看相关系数要更合理因为它是在控制其他因素之后计算的。2.3 Python 版 mistyPy 和 R 版 mistyR 怎么选mistyR 最早是 R 包发表在 2022 年前后社区用过的人相对多教程和 issue 也多。mistyPy 是 Python 实现核心算法思想一致包名在 PyPI 上叫 mistypyGitHub 仓库是 saezlab/mistyPy。对于习惯 Python 分析流程的人来说mistyPy 可以直接嵌入 scanpy 生态处理 AnnData 对象非常方便这是它最大的吸引力。但也要说实话mistyPy 的文档和社区积累不如 R 版完整版本迭代也不快。如果你是纯新手、没有 Python 基础R 版可能是更省心的选择如果你组里管线已经是 Python 为主mistyPy 能把流程统一起来省去跨语言传数据的麻烦。我自己当时是因为项目里所有预处理都在 scanpy 里做所以选了 mistyPy跑下来发现核心结果和论文里的逻辑完全对得上。3. 动手前的数据准备与环境配置3.1 Python 环境与依赖安装先把这个搞定能省一半的坑mistyPy 是个比较轻量的包依赖主要是 pandas、numpy、scikit-learn 这些常见库。我一般习惯用 conda 单独建一个环境避免和 base 环境打架。conda create -n misty python3.9 -y conda activate misty pip install mistypy pandas numpy scikit-learn如果你已经装了 anndata 和 scanpy那就一起装进去后面整理数据会用到pip install anndata scanpy matplotlib seaborn这里有个小经验不要图省事直接在 base 环境里装。mistyPy 在某些 Python 版本下会有依赖冲突单独环境可以随时删掉重来成本低很多。我最初在 Python 3.11 下 pip install mistypy 时遇到过依赖兼容问题换成 3.9 之后就顺畅了。如果你用 conda建议环境创建就用 3.9 或 3.10这是实测下来比较稳的组合。3.2 表达矩阵和坐标数据的格式整理规则其实很严格mistyPy 需要两个核心输入表达矩阵和坐标信息。我当时在这里卡了很久反复报错才发现是格式理解反了。表达矩阵要求是 DataFrame 格式一般情况下行是 spot空间点或细胞列是基因。这个设计和很多单细胞分析习惯基因行、细胞列不一样如果你是从 scanpy 的 AnnData 对象转出来要用adata.to_df()得到的就是 spot 为行、基因为列的结构。坐标数据是另一个 DataFrame行必须和表达矩阵的 index 完全一致列名必须包含x和y。很多同学从 scanpy 里取坐标的时候拿到的是一个 numpy 数组直接塞进去会报“找不到 x 列”之类的错误必须手动转成 DataFrame 并且把列名改好。下面是我常用的数据整理代码直接可以从 AnnData 对象生成 mistyPy 需要的数据import pandas as pd import anndata as ad # 读取空间转录组数据 adata ad.read_h5ad(visium_sample.h5ad) # 表达矩阵spot 为行基因为列 expr adata.to_df() # 基因名里可能有特殊字符最好统一替换一下 expr.columns [c.replace(-, _).replace(/, _) for c in expr.columns] # 坐标从 obsm 里取 spatial 坐标 coord_array adata.obsm[spatial] pos pd.DataFrame(coord_array, indexadata.obs_names, columns[x, y]) # 因为坐标可能不是整数转成 float 更保险 pos[x] pos[x].astype(float) pos[y] pos[y].astype(float) # 确认 index 一致 print(expr.shape, pos.shape) print(expr.index[:5]) print(pos.head())这里还有一个数据清洗建议空间转录组数据里低表达基因太多直接全量跑 misty 会非常慢而且意义不大。我一般会按平均表达量做一个过滤比如只保留平均表达 0.1 的基因把基因数压到几千个以内。# 过滤低表达基因减少计算量 expr expr.loc[:, expr.mean(axis0) 0.1] print(expr.shape)如果数据本身已经做过标准化比如 CPM、log 化直接用就行。没有做过的话建议至少做个简单的 library size 归一化否则 spot 之间的测序深度差异会严重干扰结果。4. 核心实操跑通一个完整的 mistyPy 分析4.1 最小可运行的代码从构建模型到拿到结果数据准备好之后运行 mistyPy 的代码其实不长。我当时用的是 mistypy 0.0.x 版本接口大概是下面这样。注意不同版本可能略有差异以官方仓库的最新文档为准。import misty as misty # 核心运行函数 misty_result misty.construct_misty_predictSpatial( expr, pos, exvarall, views{intraview: 10, juxtaview: 1, paraview: 5}, n_bs100, seed42, )这里的参数我逐个解释一下。exvarall表示把所有基因都作为解释变量。如果只想分析特定基因作为解释变量可以传入一个基因列表。views是一个字典定义每个视图用什么方式构建intraview: 10的意思是用 spot 自身表达谱中方差最大的前 10 个基因作为该视图的特征juxtaview: 1表示把物理相邻的 spot 的汇总表达作为特征paraview: 5表示把一定范围内这里用了 5 个距离层级的 spot 汇总表达作为特征。n_bs100是自助抽样bootstrap的次数用于估计系数的置信区间。这里我实际跑下来100 次是性价比比较高的选择。太少了结果不稳定太多了计算时间成倍增加。跑完之后结果会封装在一个对象里。想快速看结果可以用results misty_result.misty.results具体字段在不同版本里命名会有出入但大体上会包含这几块字段含义使用场景coefficients各视图内每个基因的回归系数找具体驱动基因importances基因重要性得分排序筛选关键基因gain添加某视图后的预测增益判断哪个空间尺度更重要R2方差分解结果整体判断空间上下文对目标基因的解释程度4.2 结果对象里到底有什么拆开看一遍以 R2 为例你最后会得到一个类似这样的表格结构每个目标基因一行几列分别是总 R2、intraview 解释比例、juxtaview 解释比例、paraview 解释比例。这里的 R2 越高说明目标基因的表达变异越能被空间上下文解释。我第一次看到某个结果里 R2 只有 0.15第一反应是模型做错了。后来反复验证发现这其实是正常的。在真实数据里基因表达的大部分变异是“噪声”和内在随机性能被空间结构解释的比例本来就不会特别高。真正值得关注的是视图之间的相对差异而不是 R2 的绝对值。再看 coefficients 和 importances。这两个字段能帮你找到“具体是哪些基因在影响这个 target gene”。比如你关注的基因是受体基因importances 排行榜前面正好出现了几个已知的配体基因那就可以作为配体-受体空间共定位的候选证据。当然这只能作为假说生成不能当因果结论后面还需要实验验证。4.3 视图参数选多少才合理经验值与实践逻辑视图参数的设置应该是 misty 分析里最需要动脑子的部分比代码本身重要得多。邻居数juxtaview 相关的 K决定“相邻”的范围。K 太小每个 spot 的邻居太少估计不稳定K 太大会把远处的信号也混进来局部分辨率下降。我在 Visium 数据上试过不同的 K一个经验区间是 5 到 20。具体选多少可以结合 spot 的间距和组织结构如果 spot 间距小、组织密度高可以选大一点如果稀疏选小一点。paraview 的范围选择更有讲究。它代表“更远处的信号”这个距离如果和你的生物学问题不匹配结果会非常苍白。比如你要研究肿瘤边界效应远处信号的范围设定在几十到几百微米范围内可能是合理的如果你研究的是组织分区的宏观差异范围就要更大。我自己的建议是不要一上来就全参数跑。先用小规模子集比如 500 个 spot、500 个基因试几组参数看结果是否合理再放全量数据跑。全量跑 5000 个基因可能要几十分钟到几个小时时间成本不值得浪费在试参数上。5. 结果的可视化与解读从数据表到故事5.1 哪些结果最值得先看图跑完 mistyPy你不会自动得到一张漂亮的图可视化基本要自己来。我一般先画三张图。第一张是 R2 分解图。把每个目标基因的 R2 分解画成堆叠条形图每个条形代表一个基因不同颜色代表 intraview、juxtaview、paraview 的贡献比例。这张图让你一眼看出来“哪些基因的表达变异主要是由空间邻域解释的”。第二张是预测值与真实值的散点图。选几个 R2 比较大、有代表性的基因把模型预测的表达量和实际表达量画成散点看拟合效果。第三张是把高 importance 基因映射回空间坐标。这样能直接看出这些基因在哪些区域活跃是否和你关注的解剖结构或细胞类型分布吻合。这一步可以用 matplotlib 的 scatter 函数来做颜色映射到基因表达非常简单。5.2 解读案例一组基因的三种空间行为模式我在课程复习时拿一个公共 Visium 数据集试跑了一轮印象很深的是同一个数据集里基因的空间行为模式有明显差异。有一类基因总 R2 很高其中 intraview 占大头。这说明它的表达变异主要和自身调控程序有关空间位置对它影响有限。另一类基因 juxtaview 占比很高比如一些细胞通讯相关的配体和受体它们很依赖邻近 spot 的状态。还有一类 paraview 占比高多为区域 marker 基因比如脑组织里不同解剖区域的标志基因它们的表达模式更多反映的是整个空间区域的分化状态而不是局部微环境。这三类基因对应完全不同的生物学机制也决定了后续分析方向。如果只在全局相关性矩阵里看这三类基因可能混成一团根本分不出来。5.3 最小可视化代码三行画一张堆叠条形图import matplotlib.pyplot as plt # 假设 r2_breakdown 是 mistyPy 输出的 R2 分解表格 # 列包括 target, intra, juxta, para r2_breakdown.plot( xtarget, kindbar, stackedTrue, figsize(12, 5), color[#4C72B0, #DD8452, #55A868], ) plt.ylabel(R2) plt.xticks(rotation90) plt.tight_layout() plt.show()堆叠条形图看整体分布如果要细看某几个基因可以单独取出来画横向条形图。空间映射图的代码也不复杂核心是让点的颜色对应表达量plt.figure(figsize(6, 6)) sc plt.scatter(pos[x], pos[y], cexpr[target_gene], cmapviridis, s5) plt.colorbar(sc, labelexpression) plt.title(target_gene spatial expression) plt.axis(off) plt.show()6. 常见问题与排错我踩过的坑你直接避开6.1 数据格式与维度不对齐这是 mistyPy 最常见的报错来源没有之一。常见的症状是报错信息里提示 index 对不上或者某个矩阵的维度不对。检查顺序一般是这样先用expr.shape和pos.shape确认行数一致再用expr.index.equals(pos.index)确认 index 完全一致最后确认表达矩阵是 spot 为行、基因为列。如果从 seurat 转过来的数据经常会出现基因行为、spot 为列的情况用.T转置一下就好。我当时还栽过一次坐标列名写了大写 X 和 Y结果报错说找不到 x 列。mistyPy 要求列名是小写的x和y这个细节不仔细看文档很容易漏。6.2 安装依赖与版本兼容问题mistypy 包比较小有时候会遇到 pip 安装时依赖解析卡住的情况。我的建议是用虚拟环境安装装不上就用 GitHub 源码安装pip install githttps://github.com/saezlab/mistyPy.git如果你本地的 scikit-learn 版本太新或太旧可能也会报一些奇怪的错误。固定用我上面推荐的环境组合能避开大部分问题。6.3 结果异常R2 全为 0、所有基因重要性都差不多如果跑出来的结果里 R2 全是 0先检查表达矩阵是不是没有过滤低表达基因。大量全零或接近全零的行会让模型无法学到有效信号。还有一个可能是数据还没有归一化spot 间的测序深度差异直接变成主导因素把真正的空间信号压在底下。如果所有基因的 importance 看起来都差不多、没有区分度大概率是视图参数设置不当比如邻居数太大把所有局部差异都平滑掉了。这时可以减小 juxtaview 的 K 值重新跑一次对比。6.4 运行速度太慢先跑子集再上全量mistyPy 的计算量会随 spot 数和基因数线性增长全量跑 10000 个 spot、10000 个基因可能要跑非常久。我的做法是先用随机抽样的 500 个 spot 和过滤后的 1000 个基因跑通流程确认参数没问题、结果不报错再全量跑。这就像先做菜谱实验确认了流程再做大锅菜避免两个小时跑完发现参数设错了又得重来。7. 从课程复习到实际课题还能往哪个方向扩展7.1 结合细胞类型比例做“邻居组成”分析空间转录组的 spot 通常是多个细胞的混合如果你有细胞类型反卷积结果比如 cell2location 或 SPOTlight 的输出可以把 spot 的细胞类型比例作为额外的协变量加进模型或者在拿到 misty 结果后把高 importance 基因对应的 spot 聚类再和细胞类型组成做关联。这一步能把“基因空间关系”升级为“细胞类型空间关系”回答的问题更接近临床和组织病理学的关注点。7.2 与配体受体分析和空间可变基因分析衔接mistyPy 的 importance 列表天然适合用来筛选空间配体-受体对。我在实际项目中就先用 mistyPy 找到某个信号通路受体基因的关键驱动因子再用 CellChat 或 NicheNet 验证这些因子里是否包含已知配体。两步可以形成互补misty 给出空间证据配体受体工具给出分子机制证据。7.3 适合解决的问题和不太适合的场景适合用 misty 解决的问题有三个特征一是研究对象明确是某个基因或某条通路二是它有明显的空间分布梯度或区域特异性三是你关心“周围环境”是否影响它。反过来如果你只是想知道两类细胞是否靠近直接做空间共定位统计或者邻域富集分析可能更简单直接不需要上 misty。我自己跑完一圈的体会是mistyPy 不是那种“导入数据、点一下按钮就出图”的工具它的价值在于迫使你去思考空间数据的尺度问题——你到底想解释哪个尺度上的变异邻居是几层远处是多远。把这些想清楚比纠结任何参数都重要这也是这门课复习下来我最受益的地方。最后分享一个小技巧跑 mistyPy 之前一定先把模型保存下来把关键中间结果缓存住。import pickle with open(misty_result.pkl, wb) as f: pickle.dump(misty_result, f)不然调一次参数重跑一次全量时间和耐心都会很快耗尽。这个习惯我现在每次跑空间分析都会保留折腾数据的时候能救回大把时间。
空间转录组分析：用mistyPy拆解基因表达的空间上下文
空间转录组数据拿到手之后最常被问到的一个问题就是某个基因的表达为什么会在空间上呈现梯度某个细胞群落为什么只在特定区域出现这些问题本质上都指向同一个分析目标——细胞的“空间关系”。我以前在课程里第一次接触 mistyR 的时候以为它只是又一个相关性分析工具后来真正用 Python 版本mistyPy跑完一轮数据才发现它解决的是普通相关系数完全回答不了的“多尺度拆解”问题。这篇文章就是把课程复习时整理的思路结合我自己跑数据的经验重新梳理了一遍从原理到代码到坑点都过一遍给打算做空间上下文分析的同学一个能直接上手的参考。1. 为什么要分析细胞的“空间关系”先搞清楚问题的边界1.1 拿到空间坐标之后数据分析不应该还停留在“找差异基因”做过单细胞转录组的人应该都有体会单细胞数据最大的遗憾是细胞在组织里的位置信息丢了。空间转录组技术比如 Visium、Slide-seq、Xenium 这类把“表达量”和“空间位置”对应起来数据维度一下子丰富了很多。但很多拿到空间数据的同学第一反应仍然是老一套找 marker、做聚类、跑差异表达。这些当然要做但空间数据最独特的信息——细胞在组织里的邻居是谁、处于什么微环境、周围信号的强弱如何——往往被忽略了。我在实际项目里就遇到过这种情况某个基因在整体表达量上没有明显差异但把它按空间位置画出来你会发现它只在肿瘤边缘的一个窄带里高表达而周围一圈细胞的配体基因恰好也在那里活跃。这种“位置决定的表达状态”如果只看全局表达均值根本看不出来。要回答这类问题就得把每个 spot或每个细胞的“空间背景”作为解释变量纳入模型看看它对目标基因的表达有没有贡献、贡献有多大。1.2 普通相关性分析为什么在这里不够用很多人第一个想法是直接算两个基因在空间上的相关系数不就行了吗比如 A 基因高表达的地方 B 基因也高表达就说明它们有空间共定位。这个思路在简单场景下能用但问题也不小。第一相关系数只描述“两者同步变化”无法区分因果关系和混杂因素。B 基因的高表达可能不是因为 A 基因而是因为第三个基因 C 在这个区域整体活性高。第二相关性分析没有“空间尺度”的概念。两个基因可能在相邻 spot 之间相互作用也可能在跨几个 spot 的较大范围内同步变化这两种关系的生物学含义完全不同。第三一个 spot 的基因表达从来不是被单个因子决定的而是自身调控、邻近细胞通讯、远处组织信号共同作用的结果。普通相关性分析给不出这种“多因子贡献分解”的答案。mistyR以及它的 Python 版本 mistyPy正是针对这个问题设计的它把目标基因的表达变异拆成几个空间尺度的贡献然后告诉你“自身基线解释了多少、邻域解释了多少、更大范围的空间背景解释了多少”。1.3 空间上下文spatial context到底是什么misty 的核心概念是“空间上下文”。简单理解每个 spot 的基因表达不仅由它自己的分子状态决定还受它周围一定范围内的其他 spot 影响。这个“周围”不是一个距离而是分层级的。我当时学这个的时候觉得特别像理解“人际关系”一个人现在的状态有一部分是自己长期以来的习惯自身基线有一部分是身边亲密朋友的影响物理相邻的 spot还有一部分是更大的社区氛围、远程协作或整个环境趋势的影响更大范围内的信号。这样分层之后你才能说清楚“到底是身边人把你带偏了还是整个大环境的问题”。2. mistyR 与 mistyPy 的多视图建模思想2.1 三种视图的直觉理解与数学表达misty 把空间上下文划分为三种视图view。intraview 代表 spot 自身的表达谱也就是作为基线的那部分juxtaview 代表物理上相邻的 spot 对该 spot 的影响通常定义为直接邻居例如距离最近的 K 个 spotparaview 代表更大范围内 spot 的贡献通过计算每个 spot 与周围一定距离内所有其他 spot 的平均表达来定义。用数学一点的表述就是对于目标基因 g 在 spot i 的表达 x_g,i模型写作x_g,i beta_intra * x_intra,i beta_juxta * x_juxta,i beta_para * x_para,i epsilon_i这里的 x_intra 是同一 spot 内其他基因的表达向量x_juxta 是相邻 spot 的表达汇总x_para 是远处 spot 的表达汇总。模型通过正则化线性回归可理解为加了惩罚项的线性拟合估计这些系数然后计算每个视图对方差的解释比例。这里有一个容易误解的点misty 并不是预测“某个基因在所有 spot 中的表达值”而是尝试用一个 spot 的“空间上下文”特征去解释该 spot 中目标基因表达的变异。所以输出的不是一张简单的相关性热图而是一套“方差分解表”。2.2 R2 分解、gain 与基因重要性三个最核心的输出第一次跑完 mistyPy我看着输出结果里的 R2、gain、importances 这几个字段有点懵后来才算彻底理清楚。R2 分解是最直观的它告诉你目标基因的表达变异中每个视图分别贡献了多少比例。比如 A 基因的 intraview 解释率是 20%juxtaview 是 35%paraview 是 10%。这意味着 A 基因的表达变异主要受邻近 spot 影响而不是自身基线决定的。gain 描述的是“增加某个视图后模型预测能力的提升幅度”。它衡量某个视图的增量价值。如果一个视图的 gain 很小说明去掉这个视图对预测几乎没影响。importance 则对应具体基因的重要性。它综合了回归系数和表达变异反映了“哪些基因的空间分布模式对目标基因最有解释力”。这比单纯看相关系数要更合理因为它是在控制其他因素之后计算的。2.3 Python 版 mistyPy 和 R 版 mistyR 怎么选mistyR 最早是 R 包发表在 2022 年前后社区用过的人相对多教程和 issue 也多。mistyPy 是 Python 实现核心算法思想一致包名在 PyPI 上叫 mistypyGitHub 仓库是 saezlab/mistyPy。对于习惯 Python 分析流程的人来说mistyPy 可以直接嵌入 scanpy 生态处理 AnnData 对象非常方便这是它最大的吸引力。但也要说实话mistyPy 的文档和社区积累不如 R 版完整版本迭代也不快。如果你是纯新手、没有 Python 基础R 版可能是更省心的选择如果你组里管线已经是 Python 为主mistyPy 能把流程统一起来省去跨语言传数据的麻烦。我自己当时是因为项目里所有预处理都在 scanpy 里做所以选了 mistyPy跑下来发现核心结果和论文里的逻辑完全对得上。3. 动手前的数据准备与环境配置3.1 Python 环境与依赖安装先把这个搞定能省一半的坑mistyPy 是个比较轻量的包依赖主要是 pandas、numpy、scikit-learn 这些常见库。我一般习惯用 conda 单独建一个环境避免和 base 环境打架。conda create -n misty python3.9 -y conda activate misty pip install mistypy pandas numpy scikit-learn如果你已经装了 anndata 和 scanpy那就一起装进去后面整理数据会用到pip install anndata scanpy matplotlib seaborn这里有个小经验不要图省事直接在 base 环境里装。mistyPy 在某些 Python 版本下会有依赖冲突单独环境可以随时删掉重来成本低很多。我最初在 Python 3.11 下 pip install mistypy 时遇到过依赖兼容问题换成 3.9 之后就顺畅了。如果你用 conda建议环境创建就用 3.9 或 3.10这是实测下来比较稳的组合。3.2 表达矩阵和坐标数据的格式整理规则其实很严格mistyPy 需要两个核心输入表达矩阵和坐标信息。我当时在这里卡了很久反复报错才发现是格式理解反了。表达矩阵要求是 DataFrame 格式一般情况下行是 spot空间点或细胞列是基因。这个设计和很多单细胞分析习惯基因行、细胞列不一样如果你是从 scanpy 的 AnnData 对象转出来要用adata.to_df()得到的就是 spot 为行、基因为列的结构。坐标数据是另一个 DataFrame行必须和表达矩阵的 index 完全一致列名必须包含x和y。很多同学从 scanpy 里取坐标的时候拿到的是一个 numpy 数组直接塞进去会报“找不到 x 列”之类的错误必须手动转成 DataFrame 并且把列名改好。下面是我常用的数据整理代码直接可以从 AnnData 对象生成 mistyPy 需要的数据import pandas as pd import anndata as ad # 读取空间转录组数据 adata ad.read_h5ad(visium_sample.h5ad) # 表达矩阵spot 为行基因为列 expr adata.to_df() # 基因名里可能有特殊字符最好统一替换一下 expr.columns [c.replace(-, _).replace(/, _) for c in expr.columns] # 坐标从 obsm 里取 spatial 坐标 coord_array adata.obsm[spatial] pos pd.DataFrame(coord_array, indexadata.obs_names, columns[x, y]) # 因为坐标可能不是整数转成 float 更保险 pos[x] pos[x].astype(float) pos[y] pos[y].astype(float) # 确认 index 一致 print(expr.shape, pos.shape) print(expr.index[:5]) print(pos.head())这里还有一个数据清洗建议空间转录组数据里低表达基因太多直接全量跑 misty 会非常慢而且意义不大。我一般会按平均表达量做一个过滤比如只保留平均表达 0.1 的基因把基因数压到几千个以内。# 过滤低表达基因减少计算量 expr expr.loc[:, expr.mean(axis0) 0.1] print(expr.shape)如果数据本身已经做过标准化比如 CPM、log 化直接用就行。没有做过的话建议至少做个简单的 library size 归一化否则 spot 之间的测序深度差异会严重干扰结果。4. 核心实操跑通一个完整的 mistyPy 分析4.1 最小可运行的代码从构建模型到拿到结果数据准备好之后运行 mistyPy 的代码其实不长。我当时用的是 mistypy 0.0.x 版本接口大概是下面这样。注意不同版本可能略有差异以官方仓库的最新文档为准。import misty as misty # 核心运行函数 misty_result misty.construct_misty_predictSpatial( expr, pos, exvarall, views{intraview: 10, juxtaview: 1, paraview: 5}, n_bs100, seed42, )这里的参数我逐个解释一下。exvarall表示把所有基因都作为解释变量。如果只想分析特定基因作为解释变量可以传入一个基因列表。views是一个字典定义每个视图用什么方式构建intraview: 10的意思是用 spot 自身表达谱中方差最大的前 10 个基因作为该视图的特征juxtaview: 1表示把物理相邻的 spot 的汇总表达作为特征paraview: 5表示把一定范围内这里用了 5 个距离层级的 spot 汇总表达作为特征。n_bs100是自助抽样bootstrap的次数用于估计系数的置信区间。这里我实际跑下来100 次是性价比比较高的选择。太少了结果不稳定太多了计算时间成倍增加。跑完之后结果会封装在一个对象里。想快速看结果可以用results misty_result.misty.results具体字段在不同版本里命名会有出入但大体上会包含这几块字段含义使用场景coefficients各视图内每个基因的回归系数找具体驱动基因importances基因重要性得分排序筛选关键基因gain添加某视图后的预测增益判断哪个空间尺度更重要R2方差分解结果整体判断空间上下文对目标基因的解释程度4.2 结果对象里到底有什么拆开看一遍以 R2 为例你最后会得到一个类似这样的表格结构每个目标基因一行几列分别是总 R2、intraview 解释比例、juxtaview 解释比例、paraview 解释比例。这里的 R2 越高说明目标基因的表达变异越能被空间上下文解释。我第一次看到某个结果里 R2 只有 0.15第一反应是模型做错了。后来反复验证发现这其实是正常的。在真实数据里基因表达的大部分变异是“噪声”和内在随机性能被空间结构解释的比例本来就不会特别高。真正值得关注的是视图之间的相对差异而不是 R2 的绝对值。再看 coefficients 和 importances。这两个字段能帮你找到“具体是哪些基因在影响这个 target gene”。比如你关注的基因是受体基因importances 排行榜前面正好出现了几个已知的配体基因那就可以作为配体-受体空间共定位的候选证据。当然这只能作为假说生成不能当因果结论后面还需要实验验证。4.3 视图参数选多少才合理经验值与实践逻辑视图参数的设置应该是 misty 分析里最需要动脑子的部分比代码本身重要得多。邻居数juxtaview 相关的 K决定“相邻”的范围。K 太小每个 spot 的邻居太少估计不稳定K 太大会把远处的信号也混进来局部分辨率下降。我在 Visium 数据上试过不同的 K一个经验区间是 5 到 20。具体选多少可以结合 spot 的间距和组织结构如果 spot 间距小、组织密度高可以选大一点如果稀疏选小一点。paraview 的范围选择更有讲究。它代表“更远处的信号”这个距离如果和你的生物学问题不匹配结果会非常苍白。比如你要研究肿瘤边界效应远处信号的范围设定在几十到几百微米范围内可能是合理的如果你研究的是组织分区的宏观差异范围就要更大。我自己的建议是不要一上来就全参数跑。先用小规模子集比如 500 个 spot、500 个基因试几组参数看结果是否合理再放全量数据跑。全量跑 5000 个基因可能要几十分钟到几个小时时间成本不值得浪费在试参数上。5. 结果的可视化与解读从数据表到故事5.1 哪些结果最值得先看图跑完 mistyPy你不会自动得到一张漂亮的图可视化基本要自己来。我一般先画三张图。第一张是 R2 分解图。把每个目标基因的 R2 分解画成堆叠条形图每个条形代表一个基因不同颜色代表 intraview、juxtaview、paraview 的贡献比例。这张图让你一眼看出来“哪些基因的表达变异主要是由空间邻域解释的”。第二张是预测值与真实值的散点图。选几个 R2 比较大、有代表性的基因把模型预测的表达量和实际表达量画成散点看拟合效果。第三张是把高 importance 基因映射回空间坐标。这样能直接看出这些基因在哪些区域活跃是否和你关注的解剖结构或细胞类型分布吻合。这一步可以用 matplotlib 的 scatter 函数来做颜色映射到基因表达非常简单。5.2 解读案例一组基因的三种空间行为模式我在课程复习时拿一个公共 Visium 数据集试跑了一轮印象很深的是同一个数据集里基因的空间行为模式有明显差异。有一类基因总 R2 很高其中 intraview 占大头。这说明它的表达变异主要和自身调控程序有关空间位置对它影响有限。另一类基因 juxtaview 占比很高比如一些细胞通讯相关的配体和受体它们很依赖邻近 spot 的状态。还有一类 paraview 占比高多为区域 marker 基因比如脑组织里不同解剖区域的标志基因它们的表达模式更多反映的是整个空间区域的分化状态而不是局部微环境。这三类基因对应完全不同的生物学机制也决定了后续分析方向。如果只在全局相关性矩阵里看这三类基因可能混成一团根本分不出来。5.3 最小可视化代码三行画一张堆叠条形图import matplotlib.pyplot as plt # 假设 r2_breakdown 是 mistyPy 输出的 R2 分解表格 # 列包括 target, intra, juxta, para r2_breakdown.plot( xtarget, kindbar, stackedTrue, figsize(12, 5), color[#4C72B0, #DD8452, #55A868], ) plt.ylabel(R2) plt.xticks(rotation90) plt.tight_layout() plt.show()堆叠条形图看整体分布如果要细看某几个基因可以单独取出来画横向条形图。空间映射图的代码也不复杂核心是让点的颜色对应表达量plt.figure(figsize(6, 6)) sc plt.scatter(pos[x], pos[y], cexpr[target_gene], cmapviridis, s5) plt.colorbar(sc, labelexpression) plt.title(target_gene spatial expression) plt.axis(off) plt.show()6. 常见问题与排错我踩过的坑你直接避开6.1 数据格式与维度不对齐这是 mistyPy 最常见的报错来源没有之一。常见的症状是报错信息里提示 index 对不上或者某个矩阵的维度不对。检查顺序一般是这样先用expr.shape和pos.shape确认行数一致再用expr.index.equals(pos.index)确认 index 完全一致最后确认表达矩阵是 spot 为行、基因为列。如果从 seurat 转过来的数据经常会出现基因行为、spot 为列的情况用.T转置一下就好。我当时还栽过一次坐标列名写了大写 X 和 Y结果报错说找不到 x 列。mistyPy 要求列名是小写的x和y这个细节不仔细看文档很容易漏。6.2 安装依赖与版本兼容问题mistypy 包比较小有时候会遇到 pip 安装时依赖解析卡住的情况。我的建议是用虚拟环境安装装不上就用 GitHub 源码安装pip install githttps://github.com/saezlab/mistyPy.git如果你本地的 scikit-learn 版本太新或太旧可能也会报一些奇怪的错误。固定用我上面推荐的环境组合能避开大部分问题。6.3 结果异常R2 全为 0、所有基因重要性都差不多如果跑出来的结果里 R2 全是 0先检查表达矩阵是不是没有过滤低表达基因。大量全零或接近全零的行会让模型无法学到有效信号。还有一个可能是数据还没有归一化spot 间的测序深度差异直接变成主导因素把真正的空间信号压在底下。如果所有基因的 importance 看起来都差不多、没有区分度大概率是视图参数设置不当比如邻居数太大把所有局部差异都平滑掉了。这时可以减小 juxtaview 的 K 值重新跑一次对比。6.4 运行速度太慢先跑子集再上全量mistyPy 的计算量会随 spot 数和基因数线性增长全量跑 10000 个 spot、10000 个基因可能要跑非常久。我的做法是先用随机抽样的 500 个 spot 和过滤后的 1000 个基因跑通流程确认参数没问题、结果不报错再全量跑。这就像先做菜谱实验确认了流程再做大锅菜避免两个小时跑完发现参数设错了又得重来。7. 从课程复习到实际课题还能往哪个方向扩展7.1 结合细胞类型比例做“邻居组成”分析空间转录组的 spot 通常是多个细胞的混合如果你有细胞类型反卷积结果比如 cell2location 或 SPOTlight 的输出可以把 spot 的细胞类型比例作为额外的协变量加进模型或者在拿到 misty 结果后把高 importance 基因对应的 spot 聚类再和细胞类型组成做关联。这一步能把“基因空间关系”升级为“细胞类型空间关系”回答的问题更接近临床和组织病理学的关注点。7.2 与配体受体分析和空间可变基因分析衔接mistyPy 的 importance 列表天然适合用来筛选空间配体-受体对。我在实际项目中就先用 mistyPy 找到某个信号通路受体基因的关键驱动因子再用 CellChat 或 NicheNet 验证这些因子里是否包含已知配体。两步可以形成互补misty 给出空间证据配体受体工具给出分子机制证据。7.3 适合解决的问题和不太适合的场景适合用 misty 解决的问题有三个特征一是研究对象明确是某个基因或某条通路二是它有明显的空间分布梯度或区域特异性三是你关心“周围环境”是否影响它。反过来如果你只是想知道两类细胞是否靠近直接做空间共定位统计或者邻域富集分析可能更简单直接不需要上 misty。我自己跑完一圈的体会是mistyPy 不是那种“导入数据、点一下按钮就出图”的工具它的价值在于迫使你去思考空间数据的尺度问题——你到底想解释哪个尺度上的变异邻居是几层远处是多远。把这些想清楚比纠结任何参数都重要这也是这门课复习下来我最受益的地方。最后分享一个小技巧跑 mistyPy 之前一定先把模型保存下来把关键中间结果缓存住。import pickle with open(misty_result.pkl, wb) as f: pickle.dump(misty_result, f)不然调一次参数重跑一次全量时间和耐心都会很快耗尽。这个习惯我现在每次跑空间分析都会保留折腾数据的时候能救回大把时间。