简介多重免疫荧光技术正在成为组织微环境研究的关键手段但单细胞级的空间关系分析长期缺乏高效工具。细胞邻域统计、空间遍历顺序优化等问题本质上可抽象为图论中的路径规划。哈密尔顿路径作为图遍历的基础理论提供了一种不重复访问所有节点的最优顺序方案。结合QuPath的细胞检测与通道荧光定量能力将分割结果导出为结构化数据后利用MATLAB构建邻接图并求解近似哈密尔顿路径可以显著提升空间指标计算的稳定性和可复现性。这种跨平台配合适用于荧光切片中的细胞邻域分析、组织区域划分及多重标记的通道优化尤其适合处理万级细胞规模的空间组学数据。本文从实际工程角度拆解该开源工作流的部署逻辑、关键参数与常见坑点为病理空间分析人员提供一套可直接落地的参考方案。1. 多重荧光分析还要碰哈密尔顿这套系统开源工具包把两件事连起来了做病理切片空间分析的人迟早会卡在同一个问题上QuPath 导出几百列单细胞数据之后下一步怎么算细胞邻域、怎么规划荧光通道的顺序却找不到现成脚本。我最早看到这套资源里同时出现 matlab 哈密尔顿代码和 QuPath_multiplex 时第一反应是是不是打包错了。跑通之后才确认它就是要把 QuPath 的细胞识别结果交给 MATLAB用哈密尔顿路径算法去优化遍历顺序和空间指标计算。对做多重免疫荧光、空间组学、细胞邻域统计的从业者来说这套工具能省掉从零写图论代码的时间而且完全开源可以按项目需求改。下面拆一下它的功能、使用流程和容易翻车的地方。2. 从单标记到多重荧光哈密尔顿路径在组织分析里到底解决什么问题2.1 多重荧光数据不是多张图而是一张图的多个通道叠加多重荧光染色得到的切片本质上是在同一组织位置上采集了多个荧光通道的信号。QuPath_multiplex 这类工作流的核心是把荧光信号拆成单细胞级别的量化数据每个细胞核的位置坐标、面积、核周信号强度、以及每个通道对应的表达量。常见做法是先在 QuPath 里做细胞分割比如使用 watershed 算法识别细胞核再把不同通道的像素强度取平均输出成 CSV 或 Parquet 表格。这张表格里每一行是一个细胞每一列是一个通道强度。很多人直接拿去做聚类但忽略了空间关系。细胞之间的物理距离、哪些细胞紧挨着、哪些细胞在某个标记阳性细胞的邻域范围内这些都需要你把坐标数据构造成图结构而哈密尔顿路径正是处理这类图遍历问题的一个基础工具。2.2 哈密尔顿路径的核心遍历每个节点且只经过一次哈密尔顿路径Hamiltonian path是指一条不重复地经过图中所有节点的路径。如果终点又回到起点就叫哈密尔顿回路。在实际的多重荧光分析里这个路径不一定指物理上的移动路线更多时候是一种空间遍历策略。比如你想沿着组织切片扫描一遍所有细胞或者给每个细胞按空间顺序编号让相邻细胞的编号也尽可能相邻这时候哈密尔顿路径就能提供一个理论上的最优顺序。不过需要说明一般的组织切片细胞图都不是规整网格完全求解哈密尔顿路径是 NP 完全问题细胞数量一多回溯法直接不可解。所以这套资源里的 MATLAB 代码实际用的是贪婪算法加 2-opt 局部优化只求近似解。它的价值不在于数学上的最优而在于给出一个稳定的、可复现的空间遍历顺序让你后续的计算不依赖随机扫描顺序。方法复杂度适用场景暴力回溯O(n!)节点数 15DP 状态压缩O(2^n * n^2)节点数 20最近邻贪心 2-optO(n^2)~O(n^3)节点数 1000~10000模拟退火参数敏感节点数 5000这套工具里默认走的是最近邻贪心加 2-opt我一般会先跑一遍默认参数再根据坐标分布调整输入点顺序。2.3 为什么选 MATLAB 而不是 Python 或 JavaQuPath 本身是基于 Java 的生态里也有 Groovy 脚本。但这份资源把求解部分放在 MATLAB核心原因是矩阵运算和批量坐标处理在 MATLAB 里写起来更直接。再有就是很多病理分析团队的老代码都是 MATLAB能直接接进现有 pipeline。如果你不喜欢 MATLAB也可以用 Python 照着思路改写但那就等于丢掉了一份能开箱即用的后端。3. 部署和执行QuPath 脚本与 MATLAB 代码的配合方式3.1 运行环境与文件依赖我建议直接在 Windows 或 macOS 上装 QuPath 0.4.x 和 MATLAB R2020b 及以上版本。项目里主要有两个目录QuPath_scripts 放 Groovy 脚本Matlab_codes 放哈密尔顿求解和指标计算函数。目录结构大致这样QuPath_multiplex/ ├─ QuPath_scripts/ │ ├─ script1_export_all_tables.groovy │ └─ script2_export_detection_with_channels.groovy ├─ Matlab_codes/ │ ├─ main_hamiltonian_pipeline.m │ ├─ build_cell_graph.m │ ├─ hamiltonian_path_approx.m │ └─ options.m └─ sample_data/ └─ example_cells.csvQuPath 侧需要先安装好扩展或者直接把脚本拖到 QuPath 的 Script Editor 里运行。MATLAB 侧不需要额外的工具箱只用到了统计和优化里的基础函数这大大降低了部署门槛。3.2 QuPath 侧细胞检测与通道强度提取第一步不是跑 MATLAB而是在 QuPath 里保证导出的字段是完整的。这里我一般会先惯例跑一遍下面的 Groovy 脚本它会遍历当前项目里的所有图像检测 DAPI 通道的细胞核然后输出包含每个通道平均强度的 CSV。import qupath.lib.objects.PathDetectionObjects import qupath.lib.analysis.stats.HaralickFeatures def imageData getCurrentImageData() def classifier new qupath.opencv.tools.OpenCVDensityMap() def detections getDetectionObjects() // 检查是否有检测对象没有则提示先做细胞分割 if (detections.isEmpty()) { print No detections found, run cell detection first return } // 收集每个检测对象的特征坐标、面积、通道均值 def rows [] detections.each { detection - def roi detection.getROI() def centroid roi.getCentroidX() , roi.getCentroidY() def area roi.getArea() def intensities [] def channels imageData.getServer().getMetadata().getChannels() for (int i 0; i channels.size(); i) { def stats qupath.lib.analysis.stats.StatisticsHelper.getChannelStatistics(detection, i) intensities stats.getMean() } rows [centroid, area] intensities } // 写入 CSV 到项目目录 def path buildFilePath(PROJECT_BASE_DIR, cell_data.csv) def csv new File(path) csv.withPrintWriter { out - out.println(centroid_x,centroid_y,area, channels.collect { it.getName() }.join(,)) rows.each { r - out.println(r.join(,)) } } print Exported: path这段脚本的核心逻辑是遍历当前图像的检测对象取出每个细胞核的重心坐标和面积再循环每个通道用StatisticsHelper.getChannelStatistics获取该检测对象在对应通道的平均强度。注意它假设你已经跑过细胞分割如果没跑检测对象列表是空的脚本会提前退出。参数说明detection.getROI()返回的是 ImageJ ROIType 的区域对象getCentroidX()和getCentroidY()给出的是当前图像像素坐标系下的坐标。channels.size()对应图像元数据里的通道数量通常是荧光标记数加上 DAPI。导出的 CSV 第一行是表头后面每行是一个细胞。3.3 MATLAB 侧读 CSV 和构建邻接矩阵拿到 CSV 后MATLAB 代码的主入口通常是main_hamiltonian_pipeline.m。它会先读表再根据细胞的坐标计算两两距离。这里有个关键参数邻域半径。不是所有细胞都需要连边通常设置一个半径阈值比如 30 像素只有在这个半径内的细胞才构成图中的边。这样能大大减少哈密尔顿路径搜索的规模。% main_hamiltonian_pipeline.m opt load_options(options.m); % 读取配置 % 读取 QuPath 导出的 CSV T readtable(opt.csv_path); x T.centroid_x; y T.centroid_y; % 构建邻接图返回稀疏距离矩阵 [cell_graph, dist] build_cell_graph(x, y, opt.edge_radius); % 求解哈密尔顿路径近似解 [path_order, path_dist] hamiltonian_path_approx(cell_graph, x, y, opt.method); % 输出路径顺序表保存为 CSV output_order table(path_order, x(path_order), y(path_order), ... VariableNames, {order, cell_x, cell_y}); writetable(output_order, opt.output_path);build_cell_graph内部用的是rangesearch函数找出每个细胞邻域内的邻居索引再计算对应的距离。hamiltonian_path_approx默认采用最近邻贪婪法从距离质心最近的细胞开始每次找当前路径端点最近的未访问节点。如果这个节点与端点之间的距离大于edge_radius它仍然会强制连接因为哈密尔顿路径要求必须是连通图否则无法遍历所有细胞。参数说明opt.edge_radius是邻域半径单位与坐标一致。如果设得太小图可能不连通设得太大计算量会上升。我建议先让 MATLAB 打印出连通分量数确认图是连通的再跑路径求解。opt.method可选greedy或greedy2opt后者会额外做几轮 2-opt 翻转优化把路径总长度缩短约 5%~15%。4. 动手跑通一个示例从切片到哈密尔顿路径输出4.1 准备示例数据和参数预设这套资源里带一个sample_data/example_cells.csv内容大概是 200~500 个模拟细胞坐标和 5 个通道强度。如果你没有完整的 QuPath 项目也可以直接用这个示例数据验证 MATLAB 流程。我的习惯是先把这个 CSV 复制到一个独立目录避免污染原项目。打开options.m需要改这几个地方opt.csv_path sample_data/example_cells.csv; opt.edge_radius 50; % 正方形像素邻域半径 opt.method greedy2opt; % 近似解方法 opt.output_path output/path_order.csv; opt.visualize true; % 是否画路径图其中visualize设为 true 时MATLAB 会画一条路径连线图方便你直观看到遍历顺序。第一次跑建议打开因为能立刻发现问题比如图不连通或坐标轴比例不对。4.2 运行 QuPath 导出脚本如果你有自己的切片如果有自己的荧光切片步骤是打开 QuPath创建项目并添加图像先在“细胞检测”面板里选 DAPI 通道设好阈值和核分割参数运行检测。检测完成后把前面 3.2 节的 Groovy 脚本贴到 Script Editor改成自己的导出路径运行。它会把所有检测对象输出成一个 CSV。注意QuPath 的坐标是像素坐标不是物理微米坐标。如果你在 MATLAB 里计算距离需要知道像素尺寸。通常用默认的 0.5 µm/pixel但这必须与图像元数据一致否则后续的距离指标全部偏差。4.3 在 MATLAB 里跑通完整流程在 MATLAB 当前目录切换到Matlab_codes运行main_hamiltonian_pipeline。屏幕会打印每个阶段的进度最后生成output/path_order.csv。这个 CSV 只有三列order路径序号、cell_x、cell_y表示按哈密尔顿路径遍历时每个细胞的访问顺序和坐标。之后你可以加载这个 CSV做一个最简单的验证绘制路径线观察是否有长距离跳跃。正常情况相邻节点距离应该集中在edge_radius附近如果出现跨越整个视野的长线说明图不连通或者贪婪算法在某个节点被迫强制跳跃了。4.4 把路径结果转成空间指标哈密尔顿路径本身不是终点它还能用来计算一组空间统计量。比如沿路径计算相邻细胞之间的表达量差值可以反映组织区域中某个标记信号的空间连续性或者把路径按固定步长切段计算每段区域内细胞类型的多样性。% compute_spatial_metrics.m T readtable(output/path_order.csv); expr readtable(sample_data/example_cells.csv); % 按路径顺序排列表达量 sorted_intensity expr.channel_3(T.order); % 计算相邻细胞表达量的绝对差分 diff_series abs(diff(sorted_intensity)); % 输出均值和标准差 mean_diff mean(diff_series); std_diff std(diff_series);这段代码的作用是把表达式强度按哈密尔顿路径顺序重排然后计算相邻细胞之间的表达变化幅度。如果你把同一批细胞按随机顺序计算得到的mean_diff会显著不同。用它来判断某个标记是否在组织里呈斑块状分布比纯粹的聚类结果更直观。5. 常见问题与避坑记录从翻车到跑通的五个血泪点5.1 现象QuPath 导出的 CSV 通道顺序对不上第一次我拿到同事导出的表发现第二列写的是 DAPI但实际读回来强度最高的是别的通道。后来核对发现QuPath 的通道列表不总是按照荧光滤光片的顺序排列有些时候是自动识别的顺序跟图像元数据里的通道名不一致。原因QuPath 图像服务器读取通道信息时可能按波长排序或者按文件内嵌顺序不是你在拼图时定义的顺序。解决在 Groovy 脚本里显式获取channels.size()后不要直接按索引写死列名而是把channels.collect { it.getName() }打印出来确认一遍。我现在的做法是先把通道名输出到一个单独 log再和 CSV 表头比对确保每个通道列的语义正确。这一步没做后面的所有通道分析都是错的。5.2 现象MATLAB 内存直接爆掉读 50 万细胞卡死用整套流程跑真实的整张切片时细胞数有几十万rangesearch如果传的是全量坐标矩阵瞬间要构建 50 万乘 50 万的稠密距离矩阵内存直接溢出。原因build_cell_graph里如果先用pdist2算所有点对距离那 O(n^2) 矩阵就把内存吃光了。解决先用 k-d tree 范围搜索即 MATLAB 的rangesearch或knnsearch返回邻域索引后再构建稀疏矩阵。我现在会在build_cell_graph.m里检查如果numel(x) 20000强制使用块式范围搜索每次只处理一个子块避免一次性生成全距离矩阵。5.3 现象哈密尔顿路径求解慢到怀疑人生跑 2 万细胞贪心加 2-opt 花了半小时还没结束一直在 2-opt 循环里出不来。原因默认的 2-opt 实现会反复扫描所有边对当路径节点数达到数万时每轮改进都是 O(n^2)。解决把 2-opt 的最大迭代轮数限制在 5~10 轮并且当连续三次未改进就提前终止。另外可以把edge_radius调小让图更稀疏贪心寻找最近邻时搜索范围更小。如果细胞数超过 5 万我一般直接只用贪婪法不要 2-opt因为最终空间指标差异小于 2%没必要耗时间。5.4 现象路径图在 MATLAB 里拉伸变形坐标不是正方形画出来的路径连线图横轴和纵轴单位不一致看起来像被压扁了。原因MATLAB 的 figure 默认坐标比例不是等比的而 QuPath 的像素不一定是正方形常见于各向异性的扫描仪。解决在画图前检查metadata里的像素尺寸并在 MATLAB 里用axis equal强制坐标比例一致。如果图像的像素宽高比不是 1:1要在读坐标时乘以对应的校正系数。否则计算出的路径长度会多出虚假分量。5.5 现象mein 脚本在别人机器上报错找不到函数换了一台电脑跑MATLAB 提示找不到rangesearch但代码是对的。原因rangesearch需要统计和机器学习工具箱基础 MATLAB 安装没有。另外某个版本里rangesearch的返回值格式变了。解决在代码开头用exists(rangesearch,file)检测如果没有就回退到分块距离计算。我后来干脆在build_cell_graph里自己实现了基于 k-d tree 的邻域搜索彻底不依赖工具箱这样移植到没有额外工具箱的机器上也能跑。6. 进阶用法用已验证的哈密尔顿路径结果反向校正 QuPath 检测质量跑通路径只是第一步更实用的是把 MATLAB 算出来的路径顺序作为空间依据反过来检查 QuPath 检测结果里有没有漏检或重复检出的细胞。比如你可以按路径顺序计算每一步的步长正常情况下相邻步长应该落在现细胞核半径附近。如果某一步突然出现超长跳跃那个位置大概率有个细胞没被检测到。具体做法是在 MATLAB 里找出所有大于阈值比如 5 倍中位步长的路径段把两端细胞的坐标输出成suspicious_pairs.csv然后回到 QuPath用这个表中的坐标定位到原图上肉眼检查是不是漏了一个核或两个核贴太近被合并了。我经常用这个办法快速筛查批量导出的数据比盯着切片来回扫效率高得多。验证路径正确性也有一个习惯把路径长度和随机顺序的总长度做对比。在同一组细胞上用randperm生成 20 次随机排列计算每个排列相邻点距离之和再和哈密尔顿路径的距离之和比。如果哈密尔顿结果比中位数还长说明求解算法可能没收敛或图短路太重要检查参数。从那以后我每次跑新批次的切片都强制走一遍检查 CSV 表头、检测连通性、对比随机长度这三步再进入正式分析。这套流程忙下来虽然多花十几分钟但能省下后续分析翻车的几天时间。希望帮到你。本文还有配套的精品资源点击获取