简介QuPath_Multiplex是QuPath项目的多路复用扩展代码包QuPath本身是一款面向数字病理学的开源整张载玻片图像分析软件由贝尔法斯特女王大学团队维护支持IHC与HE分析工作流、细胞分割、组织微阵列解阵列、交互式机器学习及分层对象数据模型并可通过脚本或插件与MATLAB、ImageJ等外部工具集成。代码包共700个文件主体为466个Java源码文件另含Gradle构建脚本、CSS/HTML界面资源、jar依赖库、txt说明文档及多类许可文件压缩包大小约69.43MB目录结构完整适合生物医学图像分析研究者、Java开发者和拟扩展QuPath的算法工程师学习参考。该资源已有1852人浏览学习。借助这套源码可快速理解QuPath的分层数据模型、算法插件机制与工程构建流程为深度定制病理图像分析工具提供直接参考。1. QuPath_Multiplex 到底是什么后面的 matlab 哈密尔顿代码解决什么问题做多重免疫荧光mIF图像分析时QuPath_Multiplex 这类工作流负责把一张多通道全片扫描图像变成一张单细胞测量表它输出的不是好看的染色图而是每个细胞对象的坐标、面积、通道平均强度和表型。很多团队在这里停住只会统计阳性率忽略细胞之间怎么排布。把测量表导入 MATLAB 做图论分析是补齐空间维度的一条低成本路径。如何把散点状细胞按空间访问顺序串起来对应图论里的哈密尔顿路径问题经过每个细胞恰好一次、不重复不遗漏。matlab 哈密尔顿代码就是解决这类问题的从 QuPath 的 multiplex 单细胞表建图求一条空间路径再沿路径提取连续表型段等特征。这个组合适合已经跑通 QuPath 细胞检测、想用空间排列特征做组间对比或临床相关性研究的开发者也适合对 Python/R 生态不熟、但仍依赖 MATLAB 做后续统计的团队。下文从 QuPath 怎么导出最省心的测量表开始一直讲到路径求解的代码、参数边界和常见翻车场景。2. 从 QuPath 拿 multiplex 单细胞表导出脚本与字段取舍QuPath 的 multiplex 分析流程可以简化为导入多光谱图像、拆通道、细胞检测、表型分类、测量导出。后面 MATLAB 端能不能跑通取决于这一步导出表中的坐标与强度测量是否干净。很多问题不是 MATLAB 算错而是源头表里混进了像素单位、通道串扰和边界伪影。2.1 导入和通道拆分决定强度测量命名的第一步常见做法是把多光谱平台输出的 OME-TIFF 直接拖进 QuPath。导入后先确认图像类型为 Fluorescence并检查 Channels 面板里的通道名。不要一上来就点“Estimate stain vectors”对 mIF 数据它估计的向量可能和你实际使用的荧光素不一致导致拆出来的通道强度有串扰。如果你有单染对照切片就把光谱向量手动填到染色向量里。通道名尽量写成 CD8-Cy3、CD4-TRITC 这种格式而不是 Channel 1 / Channel 2后续脚本解析强度测量项时会直观很多避坑也能少走一步。操作推荐设置原因图像导入OME-TIFFFluorescence 类型保留多通道元数据通道命名标记物-荧光素强度测量项名可读染色向量单染对照优先避免通道串扰进入强度表像素校准按扫描仪 metadata 设置坐标转换为微米的前提2.2 细胞检测与分类参数直接影响后续连通度细胞检测使用 QuPath 预置的 Cell detection关键参数有 Requested pixel size、Background radius、Nucleus diameter。我一般先按 0.5 µm/px 的图像重采样来跑背景半径取 8 µm核直径取 12 µm然后再根据当前切片微调。Cell detection的阈值会让导出表里的 Cell: Area 和强度分布发生很大变化而后面 MATLAB 建图又依赖坐标和面积过滤所以这一步值得来回试几个参数。如果你不打算在 QuPath 里做机器学习分类可以在 Positive cell detection 里为每个标记物设阈值。注意一旦设了阈值导出的 classification 字段会变成你写的阳性/阴性标签如果后续想在 MATLAB 里重新划阈值就不要在 QuPath 里做二次分类否则导出的分类字段会掩盖真实强度分布。导出表必须包含通道平均强度且保留未分类细胞便于 MATLAB 建图前的独立筛选。2.3 用 Groovy 脚本一次导出测量表下面这段脚本会遍历当前图像里的所有细胞对象输出像素坐标、微米坐标、细胞面积、所有以 mean intensity 结尾的通道强度以及表型字符串。import qupath.lib.scripting.QP def imageData QP.getCurrentImageData() def cal imageData.getServer().getPixelCalibration() def cells QP.getCellObjects() if (cells.isEmpty()) { print(No cell objects found) return } // 取所有以 mean intensity 结尾的测量项 def intensityCols cells.get(0).getMeasurementList().getMeasurementNames() .findAll { it.endsWith(mean intensity) } def sb new StringBuilder() sb.append(cell_id,centroid_x_px,centroid_y_px,centroid_x_um,centroid_y_um,area_um2,) sb.append(intensityCols.join(,)) sb.append(,classification\n) cells.eachWithIndex { c, idx - def roi c.getROI() double xUm roi.getCentroidX() * cal.getPixelWidthMicrons() double yUm roi.getCentroidY() * cal.getPixelHeightMicrons() double area c.getMeasurementList().getMeasurementValue(Cell: Area) def row [idx, roi.getCentroidX(), roi.getCentroidY(), xUm, yUm, area] intensityCols.each { name - row c.getMeasurementList().getMeasurementValue(name) } row (c.getPathClass() ? c.getPathClass().toString() : unclassified) sb.append(row.join(,) \n) } new File(E:/analysis/phenotype_table.csv).text sb.toString()逻辑说明脚本先把所有关于mean intensity的测量列名抓出来这样后面如果改了通道名脚本不用大改然后逐细胞读取 ROI 质心和校准系数分别保存像素坐标和微米坐标。classification取自 PathClass如果细胞没有分类就写成 unclassified避免空值把 MATLAB 的 readtable 搞乱。参数说明缺口路径要改成你自己的实际保存路径。cal.getPixelWidthMicrons()返回的是单个像素的物理宽度如果 QuPath 没校准它可能是 1这时微米坐标等于像素坐标后文会专门讲怎么发现这个问题。导出的列含义单位注意cell_id循环序号无只做行标识centroid_x_px / centroid_y_px质心坐标像素方便以后回看原图centroid_x_um / centroid_y_um质心坐标微米MATLAB 建图主键area_um2细胞面积平方微米过滤碎片各通道 mean intensity通道平均荧光强度强度单位建图前可复核classification表型分类字符串无可能含 unclassified3. MATLAB 建邻接图从散点到可解哈密尔顿路径的空间结构拿到单细胞表以后第一件事不是跑哈密尔顿路径而是先把细胞点变成一张无向带权图。只有带上空间邻接关系路径问题才有意义。如果直接拿所有细胞做完全图哈密尔顿路径会退化成按编号排序的折线基本没有生物学解释价值。3.1 读取数据并做坐标/强度检查用readtable读入 CSV 后先删除坐标缺失行再根据面积过滤碎片细胞。这一步能避免把染色折叠处的小碎核当成真实细胞。T readtable(phenotype_table.csv, TextType, string); % 删除坐标缺失行避免后面建图出现 NaN T rmmissing(T, DataVariables, {centroid_x_um, centroid_y_um}); % 如果 QuPath 未正确设置像素校准坐标范围会偏小 xrange range(T.centroid_x_um); yrange range(T.centroid_y_um); if xrange 200 yrange 200 warning(坐标范围感觉像像素而不是微米请回到 QuPath 重新校准); end % 过滤碎片细胞面积太小的大多是边缘伪影或细胞碎片 validArea T.area_um2 5; T T(validArea, :);逻辑说明rmissing只针对坐标列表型列里出现 unclassified 字符串不会被当成缺失值删掉。range检查是一种粗筛一张 1 mm × 1 mm 的组织切片如果以微米为单位坐标范围通常是几百到几千如果范围只有几十多半是像素单位。参数说明面积阈值 5 µm² 不是固定值要根据你的扫描分辨率和细胞大小调整。如果用的核分割是 10 µm 直径最小也能到 70 µm² 左右所以 5 µm² 只负责过滤碎片。3.2 半径图构建与主连通分量保留邻接图最常见的构建方法是半径图距离小于阈值的两个细胞连一条边边权就是欧氏距离。我一般先用 30 µm 作为初始半径再根据连通率调整。R 30; % 微米先按组织内细胞平均间距估计 pts [T.centroid_x_um, T.centroid_y_um]; [idx, dcell] rangesearch(pts, pts, R); s []; t []; w []; for i 1:height(T) for k 1:length(idx{i}) j idx{i}(k); if j i % 只保留上三角避免重复边 s(end1) i; t(end1) j; w(end1) dcell{i}(k); end end end A sparse(s, t, w, height(T), height(T)); A max(A, A); % 对称得到无向带权邻接矩阵 G graph(A); G.Nodes.OriginalIdx (1:height(T)); % 只保留最大连通分量保证后续路径存在 [bin, binsize] conncomp(G); [~, mainId] max(binsize); mainIdx find(bin mainId); G subgraph(G, mainIdx); % 提取主连通分量的稠密距离矩阵 A_main full(adjacency(G, weighted));逻辑说明rangesearch返回每个半径范围内邻居的序号和距离循环里用j i只保留上三角再对称化防止后面 graph 统计边数时重复。conncomp会给出所有连通分量的编号最大分量才是我们要分析的组织主体。参数说明R 从 30 µm 起步不代表一定要用 30。判断标准是主连通分量占比如果max(binsize) / sum(binsize)低于 70%通常说明半径太小或者切片里有大量空白坏死区。另一种情况是平均度数过高到几十说明半径偏大路径会失去局部约束。把半径看作灵敏度参数而不是固定常数。3.3 节点属性挂载与建图策略取舍构建图之后把表型字符串和数值编码挂到节点上这样后面沿路径算连续同类段时不需要再查表。origIdx G.Nodes.OriginalIdx; classVals T.classification(origIdx); phenoNames unique(classVals); G.Nodes.Phenotype classVals; G.Nodes.PhenoNum zeros(height(G), 1); for p 1:length(phenoNames) G.Nodes.PhenoNum(classVals phenoNames(p)) p; end逻辑说明subgraph会把节点重编号所以必须提前保存OriginalIdx。表型字符串转成数值编码是为了在哈密尔顿路径输出后用向量计算连续段。建图策略不止半径图一种这里列一个对比建图方式边规则适用场景主要问题半径图距离小于 R 的细胞连边与邻域作用尺度一致R 需要调k 近邻图每个细胞连最近 k 个细胞密度不均的组织空间尺度容易被密度污染Delaunay/RNG基于三角剖分细胞间物理边界会生成穿过空白区的长边完全图任意两点都连边无空间约束哈密尔顿路径失去意义本文采用半径图是因为 30 µm 邻域在 mIF 里通常能覆盖细胞间直接接触和旁分泌尺度比 kNN 更容易向审稿人解释。4. 哈密尔顿路径近似求解贪心最近邻与 2-opt 的落地代码图建好后求一条经过所有细胞节点一次的路径就是哈密尔顿路径问题。精确求解在大规模图上不可行但细胞图像有空间嵌入属性用贪心初始化加 2-opt 局部优化就能拿到可解释的近似路径。4.1 为什么用哈密尔顿路径而不是聚类边界聚类只能给出“哪些细胞靠得近”哈密尔顿路径给出的是全局顺序沿着这条顺序走完所有细胞总移动距离尽量小。如果某个表型的细胞在组织里形成条带或簇它们会在这条路径上连成一段较长的连续序列如果表型分布随机路径上同类别的连续段就会很短。所以这个方法的输出不是聚类标签而是路径本身。后续可以统计连续段长度做置换检验判断空间聚集是否显著。4.2 近似哈密尔顿路径函数下面这个函数接收稠密距离矩阵 A先做贪心最近邻再做 2-opt 反转优化function path approx_hamiltonian_path_2opt(A, startIdx) % A: 稠密距离矩阵A(i,j)0 表示 i 和 j 之间无直接边 % startIdx: 起始细胞在当前矩阵中的节点编号 N size(A, 1); D A; D(1:N1:end) 0; % 对角线归零 D(D 0) inf; % 无直接边设为无穷大 % 第一步贪心最近邻 visited false(N, 1); path zeros(N, 1); cur startIdx; visited(cur) true; path(1) cur; validLen N; for step 2:N rowD D(cur, :); rowD(visited) inf; [minD, nextIdx] min(rowD); if isinf(minD) validLen step - 1; % 图不连通提前结束 break; end path(step) nextIdx; visited(nextIdx) true; cur nextIdx; end path path(1:validLen); len validLen; % 第二步2-opt 局部优化 if len 3 return; end totalDist sum(D(sub2ind([N, N], path(1:end-1), path(2:end)))); maxIter 300; improved true; iter 0; while improved iter maxIter improved false; for i 1:len-2 for j i2:len oldLen D(path(i), path(i1)) D(path(j-1), path(j)); newLen D(path(i), path(j)) D(path(j-1), path(i1)); if oldLen newLen 1e-6 path(i1:j-1) flip(path(i1:j-1)); totalDist totalDist - oldLen newLen; improved true; end end end iter iter 1; end end逻辑说明贪心步骤从startIdx出发每次从当前节点的邻居里选一个最近且未访问的节点直到走完所有节点如果遇到不连通图会中途退出。2-opt 步骤尝试把路径中间的一段整体反转如果反转后路径总长更短就保留反转结果。由于新连接边必须已经存在于图中D(path(i), path(j))为无穷大的情况不会触发交换所以不会引入非法边。参数说明startIdx可以任意选我一般选最大连通分量的第一个节点maxIter 300是迭代上限防止死循环。1e-6是浮点比较的容差用来忽略毫米级数值波动。对 1000 个节点的图这组参数在普通台式机上几秒到几十秒能跑完如果节点数上万建议先做分块或采样后面避坑章节会说。4.3 路径上连续表型段的显著性检验拿到路径序列后最直接的空间特征是最长连续同类表型段的长度。用随机置换的零模型做检验能回答“这个连续段是不是碰巧排出来的”。phenoSeq G.Nodes.PhenoNum(path); maxReal maxRunLength(phenoSeq); nPerm 1000; nullMax zeros(nPerm, 1); for k 1:nPerm shuffled phenoSeq(randperm(length(phenoSeq))); nullMax(k) maxRunLength(shuffled); end pval mean(nullMax maxReal); fprintf(最长同类表型连续段 %d, 置换检验 p %.3f\n, maxReal, pval); function r maxRunLength(v) r 1; cur 1; for i 2:length(v) if v(i) v(i-1) cur cur 1; r max(r, cur); else cur 1; end end end逻辑说明randperm打乱的是表型序列而不是节点坐标所以零模型保留路径本身的空间约束只破坏表型与空间位置的对应关系。如果真实maxRunLength在 1000 次随机排列里排到了前 5%说明这类细胞确实更倾向于在路径上连续出现。参数说明nPerm 1000是最低可接受的置换次数如果最终 p 值在 0.05 附近建议加到 5000 或 10000 再下结论。多标记物同时检验时需要做 FDR 校正否则运气成分会被系统性放大。5. 避坑QuPath MATLAB 邻域分析常见的 4 个翻车点这一章专门讲我在这条链路里踩过的四个典型问题。每个都按现象、原因、解决来写方便你在自己的数据上直接定位。5.1 坐标单位没统一距离阈值全盘失效现象R 设 30 µm结果图里边的数量要么少到几乎不连通要么平均度数高到像完全图路径总长度大得离谱。原因Groovy 脚本里乘了getPixelWidthMicrons()但 QuPath 的像素校准在导入某些格式时并没有生效返回的是 1所以微米坐标其实等于像素坐标。这时 30 µm 的半径实际是 30 像素尺度完全错误。解决在 QuPath 里先打开图像的像素校准面板确认每个方向的物理尺寸导出后在 MATLAB 里打印max(T.centroid_x_um)如果数值和切片标尺对不上就在导出脚本里手动传入已知校准系数不要在 MATLAB 里再乘一次。坐标系统一这步宁可慢不能省。5.2 强度批次效应让表型阈值失效路径上的类别分布不可比现象同一套代码在第二张切片上跑出的邻域图看起来合理但连续段长度从 12 掉到 3反复检查代码也找不到原因。原因荧光淬灭、曝光时间、组织厚度差异会让相同阳性信号的通道强度整体偏移。如果你在 QuPath 里用固定阈值做 Positive detection第二张切片里大量真实阳性细胞被分成了阴性表型字符串本身就是错的。解决不要在 QuPath 里把阈值写死。导出原始通道平均强度在 MATLAB 里按样本分位数或混合高斯模型把每个通道重新分成低/中/高三档。哈密尔顿路径的建图只用坐标对强度不敏感但表型标签一错后面的连续段检验就是白算。5.3 哈密尔顿回路和哈密尔顿路径混用导致路径一直找不到现象期望得到一条首尾相接的闭环调用回路函数或自己写闭合检验时经常报错但改成开放路径后很快就能跑通。原因回路要求最后一格回到起点条件比路径强得多。组织边缘本身就是断开的半径图很难天然存在完整回路尤其是切片边界不规则的时候路径能通过所有细胞但首尾之间却没有一条满足规则的边。解决默认按开放路径处理除非是类器官、球状体这类边界紧凑且在视野内完整的样本才考虑闭环。即使要用闭环也应先算开放路径再额外检查首尾边是否存在于邻接矩阵中而不是一上来就强行求回路。5.4 边缘伪影和碎片细胞把路径拉成一条假直线现象路径上出现一长串面积很小的点连续段检验结果显示 p 值极显著但看着原图那些点都在组织折叠或切痕上。原因QuPath 的细胞检测把荧光边缘、折叠区域和细胞碎片都识别成了对象。这些假细胞在空间上常沿折叠线排列天然会在路径上形成长的连续序列但它们不表达真实的邻域关系。解决在 MATLAB 建图前用area_um2下限过滤同时增加一个 DAPI 或核通道强度下限保留下来的细胞必须有真实核信号。如果切片里有明确坏死区建议在 QuPath 里用Tissue Detection注释做掩膜导出前就把坏死区排除。这一步能避免大量“路径显著但生物学无关”的伪发现。6. 把路径结果可视化回 QuPath用 overlay 检查哈密尔顿路径是否跑偏我的习惯是在批量分析之前先导出一条 ROI 的路径把它叠在原图上肉眼扫一遍。所谓“跑偏”是指路径上相邻两点跨过明显大块空白区而图构建时这些空白区本来不应该产生边。如果出现了说明最大连通分量里混进了较远的碎块或者半径阈值太大。如果你在导出表里保留了像素坐标列从 MATLAB 里把路径节点转成像素坐标写到 CSV% mainIdx 是 G 中各节点在 T 里的原始行号 px T.centroid_x_px(mainIdx(path)); py T.centroid_y_px(mainIdx(path)); writematrix([px, py], path_overlay.csv);然后在 QuPath 的脚本编辑器里用类似下面的方式创建一条折线对象。不同版本的 QuPath 对 ROI 构造函数的名称略有差别运行前先查一下当前版本的 API 说明。// 读取 MATLAB 写出的像素坐标 def rows new File(/path/to/path_overlay.csv).readLines().drop(1) def xPts rows.collect { Double.parseDouble(it.split(,)[0]) } def yPts rows.collect { Double.parseDouble(it.split(,)[1]) } def roi new qupath.lib.roi.PolylineROI(xPts as double[], yPts as double[], 1.0) def pathObj new qupath.lib.objects.PathObject(roi) QP.getCurrentImageData().getHierarchy().addObject(pathObj)把折线对象加进 Hierarchy 后原图上会直接显示这条路径。我一般还会让 MATLAB 同时输出每个路径节点序号这样能原地点开某个细胞看它在路径上的第几格确认不是一条乱跳的折线。这种肉眼验证比任何数值指标都更能暴露建图参数的问题。我自己在处理这类图的时候还有一个更保守的习惯先选一块 100×100 像素的小 ROI把整条链路的参数都调通再整张切片批量跑。Path 分析一旦进入批量阶段参数错误会被复制到所有样本里后面的统计补救基本是后悔药。宁可慢一点也别跳过 overlay 检查。希望帮到你。本文还有配套的精品资源点击获取