前一阵子处理一批百创智造平台的联合测序数据空间转录组和空间ATAC各有一张组织切片图像看起来都是同一块组织的相邻切片可两张图放进同一个坐标系就是错位、旋转、缩放还不一样。我盯着它们看了好一阵子甚至怀疑制片时切片方向就反了。把整个过程整理成两篇笔记之后这篇专门说第二种切片对齐思路组织轮廓形状配准法。这个方法解决的问题和上一篇一样把空间ATAC和空间转录组的spot坐标统一到同一个组织图像坐标系上让后续的peak-gene联合分析能以空间位置为基础展开。它更适合样本量大、手动选点容易疲劳的项目阶段也是我目前批量处理数据时的默认对齐方案。如果你刚接触空间多组学对齐环节常常是第一个劝退点希望这篇能帮你少走一些弯路。1. 为什么切片对齐要准备第二套方案1.1 手动选点法在真实数据里的明显短板上一篇说的方法一核心是人工在两张图上找公共解剖结构点比如腺体边界、血管分叉、组织边缘拐角然后用这些点求解一个变换矩阵。这个方法听上去简单但真正拿几张百创智造平台导出的原图试一次你就会发现几个麻烦。第一两张图来自相邻切片染色方式不同分辨率也不同空间转录组多用HE或荧光染色空间ATAC的组织成像往往更粗糙肉眼能稳定辨认的公共特征其实非常少。组织如果恰好是比较均质的肝、肾或脑皮层没有明显边界结构选点本身就有很大歧义。第二旋转和缩放耦合在一起普通人坐在屏幕前很难判断到底转了多少度、放大了多少倍选出来的点经常互相矛盾解出来的矩阵往往只对选中的点有效对其他区域完全不对。第三误差无法量化换一个人来选结果差别很大。样本一旦上了规模比如一个项目十几张切片手动选点基本是劝退流程。所以我在第二批数据开始处理之前就决定必须引入可复现的自动对齐方案。不是手动选点不好用而是它更适合单个样本的精细验证不适合作为批量流程的第一道对齐基线。1.2 轮廓配准把“点对点”换成“形状对形状”方法二的思路实际上非常朴素既然两张图的组织来自相邻切片那它们的整体轮廓形状应当高度相似只是摆放位置、方向和缩放比例不一样。我只需要从两张图上分别分割出组织区域提取轮廓然后找一个相似变换让两个轮廓尽量重合。你可以这样理解两张图像相当于同一页纸的两份复印件只不过复印件之间发生了旋转、平移和缩放。我们不需要逐字检查复印件上的文字只需要把纸张边缘轮廓对齐内容自然就大致叠上了。组织切片轮廓承担的就是“纸张边缘”的角色。这个思路最大的好处是不依赖手动选点的主观判断。组织掩膜一旦分割出来后续的质心、主轴方向、尺度信息全部由算法自动计算整个流程可以被量化、可复现、可批量跑。它是“形状对形状”的整体配准而不是“点对点”的局部拟合所以对个别局部形变也有更好的容忍度。1.3 为什么用相似变换而不是整个仿射变换做配准时有一个很容易踩的设计分叉到底用相似变换旋转、缩放、平移共4个自由度还是用完整仿射变换6个自由度我的建议是第一层对齐只用相似变换。相邻切片在理想情况下只存在摆放位置差异组织本身没有经历剪切变形所以一个相似变换已经足够描述两个坐标系之间的主要差异。仿射变换虽然能拟合得更细腻但它多出来的自由度很容易吃掉图像噪声导致变换矩阵在未选中的区域严重外推。说白了自由度越多就越容易过拟合而这种过拟合在坐标变换里是非常隐蔽的。如果某个样本确实存在明显的局部拉伸正确做法是先跑完相似变换把整体对齐误差压下来再评估是否需要局部弹性修正。我在第6节会专门讲这个因为你一定会遇到类似的样本。2. 动手之前先把两个平台的数据坐标体系摸清楚2.1 输出文件里面的坐标到底是什么百创智造平台的空间转录组产物和空间ATAC产物通常各自导出一套坐标文件两套文件里的横纵坐标字段很容易被当成同一套体系直接用实际上它们的原点和单位很可能完全不同。转录组这边一般会有一个类似barcode位置表的文件每个捕获spot对应一组图像像素坐标常见字段是barcode、x、y原点通常在组织图像的左上角x轴向右y轴向下。ATAC这边输出的坐标则可能是spot在低倍捕获图上的像素坐标也有可能以bin坐标或者格网索引的形式给出甚至文件名里还带一个分辨率字段不仔细看metadata根本发现不了。这里有一件事必须一开始就确认坐标字段的顺序。有的文件里第一列是col第二列是row有的反过来读进DataFrame之后不做可视化直接算变换矩阵你的结果可能整个镜像。我的习惯是拿到任何坐标文件先打印前五行再画一个散点图看分布形状是否与组织图像轮廓一致确认之后再进入后续步骤。2.2 单位统一像素、微米、bin 的换算对齐前最容易被忽略却影响最大的一步是做好单位归一。两张图分别来自转录组和ATAC两个产品线图像的物理像素尺寸经常不一样ATAC坐标系里一个单位甚至可能对应的是bin编号而不是像素。我一般这样处理先从平台导出的参数里找到每张图像的像素尺寸单位通常是微米/像素。如果图像像素尺寸字段缺失就用已知参考换算比如文件中标注了组织实际宽度或者图上带有标尺条。随后把所有坐标系统一到同一个参考分辨率下通常是高清转录组图像使用的微米坐标。换算公式本身很简单# 统一到转录组图像的像素坐标系 # 假设参考图像每像素对应 pixel_size_ref 微米 x_px x_um / pixel_size_ref y_px y_um / pixel_size_ref但要注意如果ATAC坐标已经是bin编号那得先知道每个bin对应的物理长度是500bp还是1万个碱基的单位以及这个格网与图像像素的换算关系。忽视这一步的典型症状是轮廓形状看起来完全一致但缩放倍率差了几十倍对齐后ATAC点全部堆在图像的一个角落里。2.3 图像文件本身也要先处理大组织切片的HE扫描图动辄上万像素乘上万像素直接用原始分辨率读进内存做配准非常吃力。我的做法是先看平台有没有提供降采样图如果没有就用tifffile分块读取或者把整图缩放到一个方便操作的尺度比如把最长边压到5000像素以内。这里有一个关键细节在降采样图像上做的所有对齐计算最终都要通过缩放因子还原到原始像素坐标系。如果你在配准时用的是缩略图而ATAC原始坐标是基于真实物理微米的那就先转成微米再基于微米算变换随后把变换结果落到降采样图上做可视化两个层面互不干扰。我踩过一次直接在缩略图上算完变换就保存坐标的坑导致下游全部错位最后花了一整天才排查出来。方向问题也要在预处理阶段处理掉。tif文件里可能包含方向元数据但很多软件读图时并不会自动应用。最简单可靠的检查方法就是画一个初始叠加图把ATAC组织轮廓按当前坐标直接叠到转录组组织轮廓上如果两者形状一致但镜像存在就需要在预处理阶段翻转图像或者变换坐标中的对应轴向。3. 方法二实操轮廓配准加相似变换3.1 从两张图上分别提取组织掩膜整个方法二的第一个实操步骤是把组织区域从背景里分割出来。百创智造平台给的图大部分是白背景加染色组织用OTSU阈值分割通常就能拿到不错的初始掩膜。我写了一个通用的掩膜提取函数核心流程是读图转灰度高斯模糊OTSU二值化形态学闭运算取最大连通域。闭运算这一步非常重要因为染色组织里经常有细小的空隙和杂点不闭合会让轮廓变成一堆碎片。import cv2 import numpy as np def extract_tissue_mask(image_path, close_kernel_size50): img cv2.imread(image_path, cv2.IMREAD_GRAYSCALE) img_blur cv2.GaussianBlur(img, (5, 5), 0) _, bin_img cv2.threshold(img_blur, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU) kernel np.ones((close_kernel_size, close_kernel_size), np.uint8) closed cv2.morphologyEx(bin_img, cv2.MORPH_CLOSE, kernel) # 取最大连通域排除墨迹、标注文字等杂散区域 num_labels, labels, stats, _ cv2.connectedComponentsWithStats(closed) if num_labels 1: largest 1 np.argmax(stats[1:, cv2.CC_STAT_AREA]) mask (labels largest).astype(np.uint8) * 255 else: mask closed return mask如果组织切片本身不是一个连续整体比如肠道截面或者坏死脱落的组织最大连通域策略会丢掉重要信息。这时候可以把面积大于某个阈值例如最大面积20%的多个连通域都保留下来再合并成一张掩膜。图像上有黑色墨迹时墨迹也会被OTSU划成前景但只要它不是组织区域的一部分用连通域面积过滤就能基本排除。3.2 用形状矩算质心、主方向和缩放得到掩膜之后下一步就是计算初始变换参数。我用的是图像矩也就是形状矩它由掩膜自动计算出来不需要任何人工标注。一阶矩给质心二阶中心矩给主轴方向面积比给缩放比例。质心的计算在OpenCV里很直接。主轴方向来自二阶中心矩构成的协方差矩阵其特征向量指向形状延伸的主要方向。说得直白一点把组织形状想象成一根不太规则的胡萝卜主轴就是这根胡萝卜横躺的方向。两张图上胡萝卜摆放角度不同这个角度差就是旋转参数。缩放比直接用两个掩膜面积的平方根比值。面积开方之后得到的是特征长度用它作为缩放比例比直接用面积比值更稳定因为面积是平方量级受小毛刺影响大。下面这段代码展示了完整流程def initial_pose_from_masks(mask_ref, mask_mov): # 参考掩膜转录组 M_ref cv2.moments(mask_ref) cx_ref M_ref[m10] / M_ref[m00] cy_ref M_ref[m01] / M_ref[m00] mu20_ref M_ref[mu20] / M_ref[m00] mu02_ref M_ref[mu02] / M_ref[m00] angle_ref 0.5 * np.arctan2(2 * M_ref[mu11], mu20_ref - mu02_ref) # 待变换掩膜ATAC M_mov cv2.moments(mask_mov) cx_mov M_mov[m10] / M_mov[m00] cy_mov M_mov[m01] / M_mov[m00] mu20_mov M_mov[mu20] / M_mov[m00] mu02_mov M_mov[mu02] / M_mov[m00] angle_mov 0.5 * np.arctan2(2 * M_mov[mu11], mu20_mov - mu02_mov) scale np.sqrt(M_ref[m00] / M_mov[m00]) angle angle_ref - angle_mov tx cx_ref - scale * (cx_mov * np.cos(angle) - cy_mov * np.sin(angle)) ty cy_ref - scale * (cx_mov * np.sin(angle) cy_mov * np.cos(angle)) return scale, angle, tx, ty这里有一个非常隐蔽的坑主轴方向存在180度歧义。胡萝卜横躺时两端谁是头谁是尾二阶矩无法区分所以算出来的旋转角可能正好差了180度。我在代码里会自动跑两个候选角度一个是算出来的angle另一个是angle加π分别计算轮廓重叠率保留重叠率更高的那个。这一步看起来愚蠢却帮我省下了很多人工核对的时间。3.3 精配阶段让轮廓严丝合缝初始位姿通常已经能把轮廓叠个七八成但还达不到可以直接用于下游的精度。原因在于质心和主轴对形状变化的响应不够灵敏组织边缘的小突起就会轻微改变主轴方向。精配我一般分两步走。第一步是网格搜索微调在初始旋转角附近加减5度按0.2度步长遍历同时小幅搜索x、y平移量每一步计算变换后掩膜与参考掩膜的重叠面积记录最大重叠率对应的参数。这个做法简单可靠跑起来也不慢。第二步是用OpenCV的ECC图像配准做细化它对输入图片的噪声比较鲁棒不需要标注关键点直接优化掩膜图像之间的相似度。需要注意ECC的输入要求是单通道浮点图掩膜需要先转成float32并归一化到0到1。def refine_with_ecc(mask_ref, mask_mov, warp_matrix, max_iters50): ref_f cv2.normalize(mask_ref.astype(np.float32), None, 0, 1, cv2.NORM_MINMAX) mov_f cv2.normalize(mask_mov.astype(np.float32), None, 0, 1, cv2.NORM_MINMAX) # 用相似变换的2x3矩阵作为初始值 warp_matrix warp_matrix.astype(np.float32) criteria (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, max_iters, 1e-5) cc, warp_matrix cv2.findTransformECC( ref_f, mov_f, warp_matrix, cv2.MOTION_AFFINE, criteria, None, 5 ) return warp_matrix这个阶段别追求过度拟合。ECC迭代次数不要开太大50次以内足够。我之前一个样本把迭代次数调到500结果把噪声也一起配进去了局部区域看着比50次更整齐但实际变换矩阵严重畸变换到其他区域全错了。3.4 把ATAC spot坐标变换到目标坐标系配准完成后所有参数都收敛到一个2×3的变换矩阵上接下来就是纯粹的坐标变换。这一步要把每个ATAC spot的原始坐标通过矩阵映射到转录组图像坐标系。实际写代码时要注意ATAC spot坐标文件里可能还有质量分、峰名等字段变换过程中必须保留标签信息后续要做join和关联分析全指望这些标签。建议输出成固定格式的CSV同时把变换矩阵存成JSON方便以后复核或者复现同一个参数。def transform_spot_coordinates(df, M, src_cols(x, y)): xy df[list(src_cols)].to_numpy(dtypenp.float64) # 齐次坐标转换 ones np.ones((xy.shape[0], 1)) xy_homo np.hstack([xy, ones]) xy_new xy_homo M.T # 2x3矩阵作用于齐次坐标 df df.copy() df[x_aligned] xy_new[:, 0] df[y_aligned] xy_new[:, 1] return df # M 为2x3相似变换矩阵 aligned_atac transform_spot_coordinates(atac_df, M) aligned_atac.to_csv(atac_aligned_coords.csv, indexFalse)变换完第一件事不是做下游分析而是画overlay检查图。把转录组组织轮廓作为底图把变换后的ATAC spot按坐标散点画上去肉眼确认形状重合后再进入验收流程。4. 别急着往下游先做三关验收4.1 第一关轮廓叠加的可视化检查自动对齐跑完我从来不会直接信任输出坐标哪怕是算法自己算出来的。第一关永远是画一张图一边是配准前的ATAC spot叠加在转录组图像上一边是配准后的ATAC spot叠加图人工肉眼确认。这一关看起来原始但筛掉了相当大比例的方向错误。我就遇到过主轴180度歧义没被程序识别、重叠率计算又因为掩膜噪声产生误判的情况最后靠人眼一眼看出组织轮廓的尖端朝向反了。程序可以辅助判断但最终决定必须由人来做。如果对比图显示组织轮廓整体重叠得很好只有个别spot零星飘在外面通常是正常现象因为捕获spot本来就会有一些落在组织边缘外。如果出现系统性偏移比如上方多了一块、下方缺了一块那就说明变换矩阵可能带了剪切分量或者参考图像和ATAC图像的裁剪范围不一致。4.2 第二关边缘距离的定量指标肉眼检查通过之后我会用定量指标给配准质量打分。方法很简单分别在参考掩膜和变换后掩膜上提取边缘轮廓点然后对每个参考边缘点用KDTree找最近邻的移动边缘点统计距离分布。from scipy.spatial import KDTree def edge_distance_metric(mask_ref, mask_mov): ref_pts np.argwhere(mask_ref 0) mov_pts np.argwhere(mask_mov 0) tree KDTree(mov_pts) dist, _ tree.query(ref_pts) return np.median(dist), np.percentile(dist, 95)我的验收标准很简单表里列一下指标建议阈值说明边缘点对中位距离小于50微米大约覆盖2到3个spot的直径边缘点对95%分位距离小于150微米防止局部大残差被中位数掩盖ATAC spot落入组织掩膜比例大于90%衡量坐标变换后数据有效利用率如果指标不过不要急着调参数先画出距离热图看看误差集中在哪些区域。热图能快速区分是整体旋转问题、缩放问题还是局部形变问题这决定了后续调整方向。4.3 第三关生物学特征交叉验证前面两关都是纯几何层面的检验第三关属于生物学层面的交叉验证也是我最喜欢的一关。它的逻辑很朴素如果对齐正确那么已知在空间转录组上表达高的区域对应的ATAC可及性也应该在空间上高度相关。我会挑一个组织里空间分布比较有特色的标志基因比如在某个解剖区域特异表达的基因同时看对应的ATAC峰信号。把两者分别映射到对齐后的坐标系上叠在一起观察。如果两个信号热点重合说明对齐方向正确没有隐藏的镜像或旋转残留。这一步虽然不能给出严格的误差数字但它在真实生物学场景下的说服力非常强。我之前调试一个脑组织样本几何指标全部通过结果在用标志基因验证时发现两套信号一个偏前一个偏后最后查明是参考图像降采样时插值方式导致了几十微米的系统性偏移重新恢复原始分辨率后再对齐才彻底解决。5. 从对齐到联合一个可落地的多组学数据组合流程5.1 网格bin划分多组学矩阵切片坐标对齐只是手段最终目的还是做空间ATAC和空间转录组的联合分析。最直接的做法是把对齐后的坐标空间划分成规则格网以格子为单位进行多组学数据汇总。bin的边长取值要综合考虑两类数据的分辨率。空间转录组spot之间通常相隔几十到一两百微米空间ATAC的捕获分辨率也可能类似所以我一般先用500微米左右的bin试跑再根据spot密度下调到200到250微米。bin太大会丢掉空间细节bin太小则大量格子只覆盖到单组学数据后续相关性分析不稳定。import pandas as pd import numpy as np def assign_bin(df, coord_col, bin_size, start0): edges np.arange(start, coord_col.max() bin_size, bin_size) return pd.cut(coord_col, edges, labelsFalse) atac_aligned[bin_x] assign_bin(atac_aligned[x_aligned], 500) atac_aligned[bin_y] assign_bin(atac_aligned[y_aligned], 500) rna[bin_x] assign_bin(rna[x_px], 500) rna[bin_y] assign_bin(rna[y_px], 500)随后按bin分组汇总。RNA这边按spot的UMI计数分bin加总ATAC这边按spot的峰可及性或peak count分bin加总。最终得到两份矩阵RNA表达矩阵和ATAC峰可及性矩阵它们共享同一套行索引也就是空间bin。这份bin矩阵多了一个之前没有的维度空间位置。后面无论做聚类还是做相关性分析都可以结合坐标进行可视化和空间注释比单纯把两类细胞谱系数据拼接在一起更符合组织真实状态。5.2 peak-gene空间共定位与联合聚类有了bin层面的多组学矩阵我优先建议做两件事。第一件事是peak-gene空间共定位分析对每个bin内的基因表达量和对应区域的ATAC可及性计算相关性筛选出在空间上共变的peak与基因组合。这种空间共定位关系可以作为增强子、启动子相互作用推断的辅助证据。但是这里要克制一些空间共定位只能说明两件事发生在相近的位置不代表调控因果关系。我在报告里都会加一句这类关联必须结合序列信息、转录因子基序和实验验证来解释。第二件事是联合聚类。最简单的路线是把RNA矩阵和ATAC矩阵分别做标准化然后拼接成一个联合矩阵再跑常规的聚类算法。更讲究一点可以用多组学因子分析类的方法把两种数据各自的共享因子和特有因子分离出来。我不建议一上来就追求复杂模型先把bin矩阵和相关矩阵建好再根据数据量决定是否升级。6. 实操中踩过的坑与排错记录6.1 组织掩膜分割失败全是杂点、空洞或墨迹我遇到最多的问题就出在掩膜分割这一步。某些样本组织坏死区域多或者染色偏浅OTSU阈值会把大片组织划成背景掩膜出现严重的内部空洞。这时候算主轴和质心都会受到干扰因为我们其实是在拿一块千疮百孔的橡皮泥做形状分析。解决办法是强化形态学闭运算把核设大一点。另外不要只依赖图像分割还可以用平台导出的spot坐标做凸包把组织区域整体框出来再用凸包掩膜修正分割结果。毕竟spot本身已经标记了组织覆盖范围用它们来兜底掩膜比自己纯靠图像处理稳定得多。6.2 切片方向翻转上下翻、左右翻的自动判断相邻切片制样时两张切片可能直接在贴片方向上一个正一个反导致形状轮廓完全一致但左右镜像。形状矩和相似变换的旋转只能处理平面内的转动处理不了镜像翻转如果不先检测后面整个流程全废。我现在的做法是在配准前增加一个枚举步骤把ATAC掩膜和坐标分别做四种变换不翻、左右翻、上下翻、上下左右都翻然后逐一计算与参考掩膜的重叠率把重叠率最高的那个作为预处理结果。四种变换跑完不过几秒钟却能把方向错误从根源上排除。这个逻辑同样适用于后期排错当你发现轮廓怎么看都不对劲时先把这四种情况跑一遍。6.3 坐标映射时的Scale错误metadata字段别想当然坐标系换算里最隐蔽的坑是误以为两个平台输出文件的像素尺寸一致。我遇到过一份数据转录组图像像素尺寸是0.5微米/像素ATAC图谱像素尺寸是2微米/像素相差了4倍。最开始的代码里我用了一个默认值跑出来就是典型的“形状对但大小差4倍”症状一眼就能看出来但如果你不做可视化这个错误会一直潜伏到下游分析。更麻烦的是有时候平台会在不同文件中使用不同单位同一个字段名一个文件存像素另一个文件存微米。我的经验是每次处理新样本前都先用文件里的已知标尺或者组织实际尺寸做一次单点验证确认坐标读出来之后数量级合理再跑全套流程。6.4 相邻切片的局部形变相似变换的精度边界相似变换假设组织整体是刚性的但相邻切片在展片、贴片和染色过程中会产生局部形变。有的样本边缘基本能对齐中间区域却整体偏移了50到100微米。这已经超过了相似变换能表达的范围强行用全局参数去拟合结果就是中间对齐了边缘又歪掉。遇到这种情况我不会急着上非线性配准。第一步先确认不是翻转和缩放问题第二步检查数据是否分成多个不连续的组织块第三步才考虑局部修正。局部修正可以用基于径向基的薄板样条变换但要控制正则化强度避免组织内部出现扭曲。一个可信的准则是非线性修正要保留组织的整体结构只处理那些全局变换无法解释的局部残差。做了这么多轮对齐之后我的习惯是给每张切片都留一份transform_matrix.json里面记录变换矩阵、图像缩放因子、像素尺寸和掩膜参数。后续重分析或者换下游工具时直接复用这个矩阵不用重新配一次。如果你想批量处理样本建议先自动跑完全流程再抽样20%的样本做人工复核兼顾效率和可靠性。方法二和方法一并不是互斥的我现在的标准流程是先用轮廓配准自动初对齐再针对置信度低的样本手动微调两套方案配合起来基本能覆盖空间ATAC与空间转录组联合分析里绝大多数对齐场景。 
长河君 · 2026/10/22 9:09:54 · 阅读 1.2k 空间多组学切片对齐：组织轮廓配准法详解
前一阵子处理一批百创智造平台的联合测序数据空间转录组和空间ATAC各有一张组织切片图像看起来都是同一块组织的相邻切片可两张图放进同一个坐标系就是错位、旋转、缩放还不一样。我盯着它们看了好一阵子甚至怀疑制片时切片方向就反了。把整个过程整理成两篇笔记之后这篇专门说第二种切片对齐思路组织轮廓形状配准法。这个方法解决的问题和上一篇一样把空间ATAC和空间转录组的spot坐标统一到同一个组织图像坐标系上让后续的peak-gene联合分析能以空间位置为基础展开。它更适合样本量大、手动选点容易疲劳的项目阶段也是我目前批量处理数据时的默认对齐方案。如果你刚接触空间多组学对齐环节常常是第一个劝退点希望这篇能帮你少走一些弯路。1. 为什么切片对齐要准备第二套方案1.1 手动选点法在真实数据里的明显短板上一篇说的方法一核心是人工在两张图上找公共解剖结构点比如腺体边界、血管分叉、组织边缘拐角然后用这些点求解一个变换矩阵。这个方法听上去简单但真正拿几张百创智造平台导出的原图试一次你就会发现几个麻烦。第一两张图来自相邻切片染色方式不同分辨率也不同空间转录组多用HE或荧光染色空间ATAC的组织成像往往更粗糙肉眼能稳定辨认的公共特征其实非常少。组织如果恰好是比较均质的肝、肾或脑皮层没有明显边界结构选点本身就有很大歧义。第二旋转和缩放耦合在一起普通人坐在屏幕前很难判断到底转了多少度、放大了多少倍选出来的点经常互相矛盾解出来的矩阵往往只对选中的点有效对其他区域完全不对。第三误差无法量化换一个人来选结果差别很大。样本一旦上了规模比如一个项目十几张切片手动选点基本是劝退流程。所以我在第二批数据开始处理之前就决定必须引入可复现的自动对齐方案。不是手动选点不好用而是它更适合单个样本的精细验证不适合作为批量流程的第一道对齐基线。1.2 轮廓配准把“点对点”换成“形状对形状”方法二的思路实际上非常朴素既然两张图的组织来自相邻切片那它们的整体轮廓形状应当高度相似只是摆放位置、方向和缩放比例不一样。我只需要从两张图上分别分割出组织区域提取轮廓然后找一个相似变换让两个轮廓尽量重合。你可以这样理解两张图像相当于同一页纸的两份复印件只不过复印件之间发生了旋转、平移和缩放。我们不需要逐字检查复印件上的文字只需要把纸张边缘轮廓对齐内容自然就大致叠上了。组织切片轮廓承担的就是“纸张边缘”的角色。这个思路最大的好处是不依赖手动选点的主观判断。组织掩膜一旦分割出来后续的质心、主轴方向、尺度信息全部由算法自动计算整个流程可以被量化、可复现、可批量跑。它是“形状对形状”的整体配准而不是“点对点”的局部拟合所以对个别局部形变也有更好的容忍度。1.3 为什么用相似变换而不是整个仿射变换做配准时有一个很容易踩的设计分叉到底用相似变换旋转、缩放、平移共4个自由度还是用完整仿射变换6个自由度我的建议是第一层对齐只用相似变换。相邻切片在理想情况下只存在摆放位置差异组织本身没有经历剪切变形所以一个相似变换已经足够描述两个坐标系之间的主要差异。仿射变换虽然能拟合得更细腻但它多出来的自由度很容易吃掉图像噪声导致变换矩阵在未选中的区域严重外推。说白了自由度越多就越容易过拟合而这种过拟合在坐标变换里是非常隐蔽的。如果某个样本确实存在明显的局部拉伸正确做法是先跑完相似变换把整体对齐误差压下来再评估是否需要局部弹性修正。我在第6节会专门讲这个因为你一定会遇到类似的样本。2. 动手之前先把两个平台的数据坐标体系摸清楚2.1 输出文件里面的坐标到底是什么百创智造平台的空间转录组产物和空间ATAC产物通常各自导出一套坐标文件两套文件里的横纵坐标字段很容易被当成同一套体系直接用实际上它们的原点和单位很可能完全不同。转录组这边一般会有一个类似barcode位置表的文件每个捕获spot对应一组图像像素坐标常见字段是barcode、x、y原点通常在组织图像的左上角x轴向右y轴向下。ATAC这边输出的坐标则可能是spot在低倍捕获图上的像素坐标也有可能以bin坐标或者格网索引的形式给出甚至文件名里还带一个分辨率字段不仔细看metadata根本发现不了。这里有一件事必须一开始就确认坐标字段的顺序。有的文件里第一列是col第二列是row有的反过来读进DataFrame之后不做可视化直接算变换矩阵你的结果可能整个镜像。我的习惯是拿到任何坐标文件先打印前五行再画一个散点图看分布形状是否与组织图像轮廓一致确认之后再进入后续步骤。2.2 单位统一像素、微米、bin 的换算对齐前最容易被忽略却影响最大的一步是做好单位归一。两张图分别来自转录组和ATAC两个产品线图像的物理像素尺寸经常不一样ATAC坐标系里一个单位甚至可能对应的是bin编号而不是像素。我一般这样处理先从平台导出的参数里找到每张图像的像素尺寸单位通常是微米/像素。如果图像像素尺寸字段缺失就用已知参考换算比如文件中标注了组织实际宽度或者图上带有标尺条。随后把所有坐标系统一到同一个参考分辨率下通常是高清转录组图像使用的微米坐标。换算公式本身很简单# 统一到转录组图像的像素坐标系 # 假设参考图像每像素对应 pixel_size_ref 微米 x_px x_um / pixel_size_ref y_px y_um / pixel_size_ref但要注意如果ATAC坐标已经是bin编号那得先知道每个bin对应的物理长度是500bp还是1万个碱基的单位以及这个格网与图像像素的换算关系。忽视这一步的典型症状是轮廓形状看起来完全一致但缩放倍率差了几十倍对齐后ATAC点全部堆在图像的一个角落里。2.3 图像文件本身也要先处理大组织切片的HE扫描图动辄上万像素乘上万像素直接用原始分辨率读进内存做配准非常吃力。我的做法是先看平台有没有提供降采样图如果没有就用tifffile分块读取或者把整图缩放到一个方便操作的尺度比如把最长边压到5000像素以内。这里有一个关键细节在降采样图像上做的所有对齐计算最终都要通过缩放因子还原到原始像素坐标系。如果你在配准时用的是缩略图而ATAC原始坐标是基于真实物理微米的那就先转成微米再基于微米算变换随后把变换结果落到降采样图上做可视化两个层面互不干扰。我踩过一次直接在缩略图上算完变换就保存坐标的坑导致下游全部错位最后花了一整天才排查出来。方向问题也要在预处理阶段处理掉。tif文件里可能包含方向元数据但很多软件读图时并不会自动应用。最简单可靠的检查方法就是画一个初始叠加图把ATAC组织轮廓按当前坐标直接叠到转录组组织轮廓上如果两者形状一致但镜像存在就需要在预处理阶段翻转图像或者变换坐标中的对应轴向。3. 方法二实操轮廓配准加相似变换3.1 从两张图上分别提取组织掩膜整个方法二的第一个实操步骤是把组织区域从背景里分割出来。百创智造平台给的图大部分是白背景加染色组织用OTSU阈值分割通常就能拿到不错的初始掩膜。我写了一个通用的掩膜提取函数核心流程是读图转灰度高斯模糊OTSU二值化形态学闭运算取最大连通域。闭运算这一步非常重要因为染色组织里经常有细小的空隙和杂点不闭合会让轮廓变成一堆碎片。import cv2 import numpy as np def extract_tissue_mask(image_path, close_kernel_size50): img cv2.imread(image_path, cv2.IMREAD_GRAYSCALE) img_blur cv2.GaussianBlur(img, (5, 5), 0) _, bin_img cv2.threshold(img_blur, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU) kernel np.ones((close_kernel_size, close_kernel_size), np.uint8) closed cv2.morphologyEx(bin_img, cv2.MORPH_CLOSE, kernel) # 取最大连通域排除墨迹、标注文字等杂散区域 num_labels, labels, stats, _ cv2.connectedComponentsWithStats(closed) if num_labels 1: largest 1 np.argmax(stats[1:, cv2.CC_STAT_AREA]) mask (labels largest).astype(np.uint8) * 255 else: mask closed return mask如果组织切片本身不是一个连续整体比如肠道截面或者坏死脱落的组织最大连通域策略会丢掉重要信息。这时候可以把面积大于某个阈值例如最大面积20%的多个连通域都保留下来再合并成一张掩膜。图像上有黑色墨迹时墨迹也会被OTSU划成前景但只要它不是组织区域的一部分用连通域面积过滤就能基本排除。3.2 用形状矩算质心、主方向和缩放得到掩膜之后下一步就是计算初始变换参数。我用的是图像矩也就是形状矩它由掩膜自动计算出来不需要任何人工标注。一阶矩给质心二阶中心矩给主轴方向面积比给缩放比例。质心的计算在OpenCV里很直接。主轴方向来自二阶中心矩构成的协方差矩阵其特征向量指向形状延伸的主要方向。说得直白一点把组织形状想象成一根不太规则的胡萝卜主轴就是这根胡萝卜横躺的方向。两张图上胡萝卜摆放角度不同这个角度差就是旋转参数。缩放比直接用两个掩膜面积的平方根比值。面积开方之后得到的是特征长度用它作为缩放比例比直接用面积比值更稳定因为面积是平方量级受小毛刺影响大。下面这段代码展示了完整流程def initial_pose_from_masks(mask_ref, mask_mov): # 参考掩膜转录组 M_ref cv2.moments(mask_ref) cx_ref M_ref[m10] / M_ref[m00] cy_ref M_ref[m01] / M_ref[m00] mu20_ref M_ref[mu20] / M_ref[m00] mu02_ref M_ref[mu02] / M_ref[m00] angle_ref 0.5 * np.arctan2(2 * M_ref[mu11], mu20_ref - mu02_ref) # 待变换掩膜ATAC M_mov cv2.moments(mask_mov) cx_mov M_mov[m10] / M_mov[m00] cy_mov M_mov[m01] / M_mov[m00] mu20_mov M_mov[mu20] / M_mov[m00] mu02_mov M_mov[mu02] / M_mov[m00] angle_mov 0.5 * np.arctan2(2 * M_mov[mu11], mu20_mov - mu02_mov) scale np.sqrt(M_ref[m00] / M_mov[m00]) angle angle_ref - angle_mov tx cx_ref - scale * (cx_mov * np.cos(angle) - cy_mov * np.sin(angle)) ty cy_ref - scale * (cx_mov * np.sin(angle) cy_mov * np.cos(angle)) return scale, angle, tx, ty这里有一个非常隐蔽的坑主轴方向存在180度歧义。胡萝卜横躺时两端谁是头谁是尾二阶矩无法区分所以算出来的旋转角可能正好差了180度。我在代码里会自动跑两个候选角度一个是算出来的angle另一个是angle加π分别计算轮廓重叠率保留重叠率更高的那个。这一步看起来愚蠢却帮我省下了很多人工核对的时间。3.3 精配阶段让轮廓严丝合缝初始位姿通常已经能把轮廓叠个七八成但还达不到可以直接用于下游的精度。原因在于质心和主轴对形状变化的响应不够灵敏组织边缘的小突起就会轻微改变主轴方向。精配我一般分两步走。第一步是网格搜索微调在初始旋转角附近加减5度按0.2度步长遍历同时小幅搜索x、y平移量每一步计算变换后掩膜与参考掩膜的重叠面积记录最大重叠率对应的参数。这个做法简单可靠跑起来也不慢。第二步是用OpenCV的ECC图像配准做细化它对输入图片的噪声比较鲁棒不需要标注关键点直接优化掩膜图像之间的相似度。需要注意ECC的输入要求是单通道浮点图掩膜需要先转成float32并归一化到0到1。def refine_with_ecc(mask_ref, mask_mov, warp_matrix, max_iters50): ref_f cv2.normalize(mask_ref.astype(np.float32), None, 0, 1, cv2.NORM_MINMAX) mov_f cv2.normalize(mask_mov.astype(np.float32), None, 0, 1, cv2.NORM_MINMAX) # 用相似变换的2x3矩阵作为初始值 warp_matrix warp_matrix.astype(np.float32) criteria (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, max_iters, 1e-5) cc, warp_matrix cv2.findTransformECC( ref_f, mov_f, warp_matrix, cv2.MOTION_AFFINE, criteria, None, 5 ) return warp_matrix这个阶段别追求过度拟合。ECC迭代次数不要开太大50次以内足够。我之前一个样本把迭代次数调到500结果把噪声也一起配进去了局部区域看着比50次更整齐但实际变换矩阵严重畸变换到其他区域全错了。3.4 把ATAC spot坐标变换到目标坐标系配准完成后所有参数都收敛到一个2×3的变换矩阵上接下来就是纯粹的坐标变换。这一步要把每个ATAC spot的原始坐标通过矩阵映射到转录组图像坐标系。实际写代码时要注意ATAC spot坐标文件里可能还有质量分、峰名等字段变换过程中必须保留标签信息后续要做join和关联分析全指望这些标签。建议输出成固定格式的CSV同时把变换矩阵存成JSON方便以后复核或者复现同一个参数。def transform_spot_coordinates(df, M, src_cols(x, y)): xy df[list(src_cols)].to_numpy(dtypenp.float64) # 齐次坐标转换 ones np.ones((xy.shape[0], 1)) xy_homo np.hstack([xy, ones]) xy_new xy_homo M.T # 2x3矩阵作用于齐次坐标 df df.copy() df[x_aligned] xy_new[:, 0] df[y_aligned] xy_new[:, 1] return df # M 为2x3相似变换矩阵 aligned_atac transform_spot_coordinates(atac_df, M) aligned_atac.to_csv(atac_aligned_coords.csv, indexFalse)变换完第一件事不是做下游分析而是画overlay检查图。把转录组组织轮廓作为底图把变换后的ATAC spot按坐标散点画上去肉眼确认形状重合后再进入验收流程。4. 别急着往下游先做三关验收4.1 第一关轮廓叠加的可视化检查自动对齐跑完我从来不会直接信任输出坐标哪怕是算法自己算出来的。第一关永远是画一张图一边是配准前的ATAC spot叠加在转录组图像上一边是配准后的ATAC spot叠加图人工肉眼确认。这一关看起来原始但筛掉了相当大比例的方向错误。我就遇到过主轴180度歧义没被程序识别、重叠率计算又因为掩膜噪声产生误判的情况最后靠人眼一眼看出组织轮廓的尖端朝向反了。程序可以辅助判断但最终决定必须由人来做。如果对比图显示组织轮廓整体重叠得很好只有个别spot零星飘在外面通常是正常现象因为捕获spot本来就会有一些落在组织边缘外。如果出现系统性偏移比如上方多了一块、下方缺了一块那就说明变换矩阵可能带了剪切分量或者参考图像和ATAC图像的裁剪范围不一致。4.2 第二关边缘距离的定量指标肉眼检查通过之后我会用定量指标给配准质量打分。方法很简单分别在参考掩膜和变换后掩膜上提取边缘轮廓点然后对每个参考边缘点用KDTree找最近邻的移动边缘点统计距离分布。from scipy.spatial import KDTree def edge_distance_metric(mask_ref, mask_mov): ref_pts np.argwhere(mask_ref 0) mov_pts np.argwhere(mask_mov 0) tree KDTree(mov_pts) dist, _ tree.query(ref_pts) return np.median(dist), np.percentile(dist, 95)我的验收标准很简单表里列一下指标建议阈值说明边缘点对中位距离小于50微米大约覆盖2到3个spot的直径边缘点对95%分位距离小于150微米防止局部大残差被中位数掩盖ATAC spot落入组织掩膜比例大于90%衡量坐标变换后数据有效利用率如果指标不过不要急着调参数先画出距离热图看看误差集中在哪些区域。热图能快速区分是整体旋转问题、缩放问题还是局部形变问题这决定了后续调整方向。4.3 第三关生物学特征交叉验证前面两关都是纯几何层面的检验第三关属于生物学层面的交叉验证也是我最喜欢的一关。它的逻辑很朴素如果对齐正确那么已知在空间转录组上表达高的区域对应的ATAC可及性也应该在空间上高度相关。我会挑一个组织里空间分布比较有特色的标志基因比如在某个解剖区域特异表达的基因同时看对应的ATAC峰信号。把两者分别映射到对齐后的坐标系上叠在一起观察。如果两个信号热点重合说明对齐方向正确没有隐藏的镜像或旋转残留。这一步虽然不能给出严格的误差数字但它在真实生物学场景下的说服力非常强。我之前调试一个脑组织样本几何指标全部通过结果在用标志基因验证时发现两套信号一个偏前一个偏后最后查明是参考图像降采样时插值方式导致了几十微米的系统性偏移重新恢复原始分辨率后再对齐才彻底解决。5. 从对齐到联合一个可落地的多组学数据组合流程5.1 网格bin划分多组学矩阵切片坐标对齐只是手段最终目的还是做空间ATAC和空间转录组的联合分析。最直接的做法是把对齐后的坐标空间划分成规则格网以格子为单位进行多组学数据汇总。bin的边长取值要综合考虑两类数据的分辨率。空间转录组spot之间通常相隔几十到一两百微米空间ATAC的捕获分辨率也可能类似所以我一般先用500微米左右的bin试跑再根据spot密度下调到200到250微米。bin太大会丢掉空间细节bin太小则大量格子只覆盖到单组学数据后续相关性分析不稳定。import pandas as pd import numpy as np def assign_bin(df, coord_col, bin_size, start0): edges np.arange(start, coord_col.max() bin_size, bin_size) return pd.cut(coord_col, edges, labelsFalse) atac_aligned[bin_x] assign_bin(atac_aligned[x_aligned], 500) atac_aligned[bin_y] assign_bin(atac_aligned[y_aligned], 500) rna[bin_x] assign_bin(rna[x_px], 500) rna[bin_y] assign_bin(rna[y_px], 500)随后按bin分组汇总。RNA这边按spot的UMI计数分bin加总ATAC这边按spot的峰可及性或peak count分bin加总。最终得到两份矩阵RNA表达矩阵和ATAC峰可及性矩阵它们共享同一套行索引也就是空间bin。这份bin矩阵多了一个之前没有的维度空间位置。后面无论做聚类还是做相关性分析都可以结合坐标进行可视化和空间注释比单纯把两类细胞谱系数据拼接在一起更符合组织真实状态。5.2 peak-gene空间共定位与联合聚类有了bin层面的多组学矩阵我优先建议做两件事。第一件事是peak-gene空间共定位分析对每个bin内的基因表达量和对应区域的ATAC可及性计算相关性筛选出在空间上共变的peak与基因组合。这种空间共定位关系可以作为增强子、启动子相互作用推断的辅助证据。但是这里要克制一些空间共定位只能说明两件事发生在相近的位置不代表调控因果关系。我在报告里都会加一句这类关联必须结合序列信息、转录因子基序和实验验证来解释。第二件事是联合聚类。最简单的路线是把RNA矩阵和ATAC矩阵分别做标准化然后拼接成一个联合矩阵再跑常规的聚类算法。更讲究一点可以用多组学因子分析类的方法把两种数据各自的共享因子和特有因子分离出来。我不建议一上来就追求复杂模型先把bin矩阵和相关矩阵建好再根据数据量决定是否升级。6. 实操中踩过的坑与排错记录6.1 组织掩膜分割失败全是杂点、空洞或墨迹我遇到最多的问题就出在掩膜分割这一步。某些样本组织坏死区域多或者染色偏浅OTSU阈值会把大片组织划成背景掩膜出现严重的内部空洞。这时候算主轴和质心都会受到干扰因为我们其实是在拿一块千疮百孔的橡皮泥做形状分析。解决办法是强化形态学闭运算把核设大一点。另外不要只依赖图像分割还可以用平台导出的spot坐标做凸包把组织区域整体框出来再用凸包掩膜修正分割结果。毕竟spot本身已经标记了组织覆盖范围用它们来兜底掩膜比自己纯靠图像处理稳定得多。6.2 切片方向翻转上下翻、左右翻的自动判断相邻切片制样时两张切片可能直接在贴片方向上一个正一个反导致形状轮廓完全一致但左右镜像。形状矩和相似变换的旋转只能处理平面内的转动处理不了镜像翻转如果不先检测后面整个流程全废。我现在的做法是在配准前增加一个枚举步骤把ATAC掩膜和坐标分别做四种变换不翻、左右翻、上下翻、上下左右都翻然后逐一计算与参考掩膜的重叠率把重叠率最高的那个作为预处理结果。四种变换跑完不过几秒钟却能把方向错误从根源上排除。这个逻辑同样适用于后期排错当你发现轮廓怎么看都不对劲时先把这四种情况跑一遍。6.3 坐标映射时的Scale错误metadata字段别想当然坐标系换算里最隐蔽的坑是误以为两个平台输出文件的像素尺寸一致。我遇到过一份数据转录组图像像素尺寸是0.5微米/像素ATAC图谱像素尺寸是2微米/像素相差了4倍。最开始的代码里我用了一个默认值跑出来就是典型的“形状对但大小差4倍”症状一眼就能看出来但如果你不做可视化这个错误会一直潜伏到下游分析。更麻烦的是有时候平台会在不同文件中使用不同单位同一个字段名一个文件存像素另一个文件存微米。我的经验是每次处理新样本前都先用文件里的已知标尺或者组织实际尺寸做一次单点验证确认坐标读出来之后数量级合理再跑全套流程。6.4 相邻切片的局部形变相似变换的精度边界相似变换假设组织整体是刚性的但相邻切片在展片、贴片和染色过程中会产生局部形变。有的样本边缘基本能对齐中间区域却整体偏移了50到100微米。这已经超过了相似变换能表达的范围强行用全局参数去拟合结果就是中间对齐了边缘又歪掉。遇到这种情况我不会急着上非线性配准。第一步先确认不是翻转和缩放问题第二步检查数据是否分成多个不连续的组织块第三步才考虑局部修正。局部修正可以用基于径向基的薄板样条变换但要控制正则化强度避免组织内部出现扭曲。一个可信的准则是非线性修正要保留组织的整体结构只处理那些全局变换无法解释的局部残差。做了这么多轮对齐之后我的习惯是给每张切片都留一份transform_matrix.json里面记录变换矩阵、图像缩放因子、像素尺寸和掩膜参数。后续重分析或者换下游工具时直接复用这个矩阵不用重新配一次。如果你想批量处理样本建议先自动跑完全流程再抽样20%的样本做人工复核兼顾效率和可靠性。方法二和方法一并不是互斥的我现在的标准流程是先用轮廓配准自动初对齐再针对置信度低的样本手动微调两套方案配合起来基本能覆盖空间ATAC与空间转录组联合分析里绝大多数对齐场景。