计算病理学图像与基因表达预测：数据预处理与线性注意力实践

做计算病理学相关项目的人看到这个标题应该会会心一笑——“digital profiling of gene expression from histology images”说白了就是给定一张HE染色的病理切片让模型把这张图“读”成一份基因表达谱。这里面的核心手段是linearized attention线性注意力而“数据来源和预处理”恰恰是整个项目里最不性感但最能决定成败的一段。这篇文章我就专门拆这一部分数据从哪来、怎么筛、怎么清洗、怎么把组织学图像和基因表达矩阵对齐以及我在实际操作中踩过的那些坑。适合看这篇文章的人我猜大致有三类一是准备在这个方向做研究但还没摸清公开数据门路的学生二是想把这个流程接到自己课题里的临床或生物背景研究者三是已经在跑模型但总觉得验证集指标虚高、训练不稳定、想回头检查数据处理环节的工程向同学。无论哪一种我都尽量把话说得具体、说得直接。1. 项目整体设计与思路拆解1.1 这任务到底在预测什么先说清楚问题的边界。histology images和gene expression之间的关系本质上是一个跨模态映射问题一张20倍镜下的HE切片里面包含细胞形态、组织排列、核浆比、间质比例等等信息这些形态学特征背后是基因表达程序在驱动。所以模型的输入是一张或一组组织图像块tile输出的是一组连续的基因表达值比如数百到数千个基因的TPM或log-transformed表达量。这个思路有价值是因为它把传统的bulk RNA-seq和空间信息联系起来了。拿一张存档切片不需要再跑测序就能推断出大致的转录组状态在临床样本尤其有吸引力——很多历史样本只有FFPE切片不可能回头再做RNA-seq但形态学信息一直都在。再加上空间转录组数据比如10x Visium出来之后这类图像到表达谱的建模又有了新的验证和训练途径。从任务类型上看这是一个典型的回归任务图像编码器把tile序列映射成patch embedding经过注意力层做全局聚合最后接一个回归头去预测基因表达。也有一部分工作把它做成分类/排序任务比如预测某个通路活性高低但本质上还是依赖同一套数据管线。1.2 为什么线性注意力是关键技术为什么这个方向会引入linearized attention而不是直接用ViT、ResNet或者常规Transformer原因特别实际一张WSI全切片图像切成256×256的tile一张切片往往能产生3000到10000个tile。如果把这些tile当成token序列输入标准Transformer自注意力的复杂度是O(n²)。n5000意味着每次前向要做2500万次两两交互这对显存和算力都是灾难更别提训练时还要叠加反向传播。线性注意力解决的就是这个问题。它的基本思路是把注意力矩阵的计算线性化将复杂度降到O(n)级别或者至少接近线性。常见的实现方式包括基于核函数近似的随机特征映射比如Performer的FAVOR、线性点积变换、或者一些稀疏/分块注意力的变体。你可以这样理解普通注意力是让每个人和会场里每个人都分别握手才能算出谁跟谁该多交流线性注意力则先让所有信息汇到一个“黑板”上每个人只需看黑板就能获得全局信息。这个类比不完全精确但方向是对的。不过要提醒一句线性注意力并不是白送的收益。很多线性化方案在长序列上确实省显存但在表达能力和收敛稳定性上会和完整softmax attention有差距需要配合位置编码、窗口注意力或更好的归一化策略才能拉回来。所以在我们项目里线性注意力是骨架但真正决定模型学没学到东西的还是数据侧的质量——也就是这篇文章要讲的内容。数据不对再强的注意力也白搭。2. 数据来源深度解析2.1 TCGA为什么它是这个项目的第一选择做histology images gene expression的配对数据绕不过TCGAThe Cancer Genome Atlas。TCGA的最大优势是同一个病人同时具备三类数据HE染色的WSI、同一组织块的RNA-seq表达谱、详细的临床元数据。这三样东西放在一起直接构成了监督学习所需的配对样本。从GDC Data Portal拿数据的路径大概是这样的进入portal.gdc.cancer.gov在Repository页面按case筛选左边栏的Data Category里选“Tissue Slide Image”拿到WSI文件选“Transcriptome Profiling”里的“Gene Expression Quantification”拿到RNA-seq数据。筛选完成后把文件加入cart下载manifest文件然后通过GDC的客户端工具批量拉取gdc-client download -m manifest.txt -d raw_files/这个manifest文件是核心索引。里面每一行包含文件ID、文件名、MD5哈希等字段。下载完成后一定校验MD5TCGA的文件偶尔会因为网络原因下载不全解压时才发现损坏就浪费时间了。项目里实际用到的TCGA癌种常见的有TCGA-BRCA乳腺癌、TCGA-LUAD和TCGA-LUSC肺癌、TCGA-HNSC头颈癌、TCGA-COAD结直肠癌。选择哪个癌种取决于你对疾病本身和HE形态特征的熟悉程度以及样本量是否足够。我的建议是第一步先选一个样本量大、切片质量相对稳定的癌种跑通全流程比如HNSC或BRCA不要一开始就搞多癌种泛化。值得一提的还有TCGA的biospecimen文件。这个文件记录了每个样本的组织来源比如“01”表示原发肿瘤“11”表示实体正常组织。如果你要做肿瘤-正常对比或者做正常组织基线这些字段是必须的。后面我会详细讲barcode怎么解析。2.2 GTEx与空间转录组补充数据怎么用TCGA虽然覆盖面广但它几乎全是肿瘤样本正常组织样本量极少且分布不均。如果你的模型希望学习正常组织形态到表达谱的映射或者希望模型的表征不要让肿瘤和正常组织差异淹没其他生物学信号就需要考虑补充GTExGenotype-Tissue Expression的数据。GTEx提供了多组织正常样本的bulk RNA-seq和对应的病理数据。不过要注意GTEx的样本形态学图像和TCGA的WSI格式不完全一致处理路径要单独适配。另一条补充路线是使用空间转录组公开数据最典型的就是10x Visium。Visium数据给了你一张HE染色图像和一张spot级的基因表达矩阵每个spot大约覆盖5到10个细胞比tile粒度细得多。虽然Visium样本量少但它提供的空间配对关系是天然的训练/验证资源。很多团队的做法是先在TCGA大样本上做预训练再用Visium做细粒度验证或迁移检查模型预测的空间表达模式是否和组织分布一致。2.3 下载后的数据组织与版本管理数据下载下来之后千万别急着开搞。先做一个清晰、可复现的目录规划我自己的标准结构是这样project_root/ ├── raw_files/ # TCGA下载的原始SVS文件和表达矩阵 ├── processing/ # 预处理中间产物 │ ├── tiles/ # 切好的图像块 │ ├── tissue_mask/ # 组织掩膜 │ └── expression/ # 清洗后的表达矩阵 ├── splits/ # train/val/test划分表 ├── logs/ # 预处理日志 └── scripts/ # 可复现的处理脚本这里有一个容易被低估的点表达矩阵的ID版本问题。TCGA的RNA-seq数据有多种量化版本比如STAR - Counts、STAR - TPM等而且不同时间下载的矩阵基因注释版本可能不同Gencode v22、v36等。分析之前一定要统一基因ID版本否则下游合并特征和模型输出时会出现大量ID对不上的情况。唯一可靠的路径是下载原始count数据后用同版本的Gencode注释重新做一次TPM标准化或者至少确认所有样本用的是同一个注释版本。另外强烈建议把元数据单独存一份csv包含sample_barcode、patient_id、slide_id、file_path、expression_file_path、组织类型编码。这份文件是后续所有对齐和划分的地基。地基没打好后面每个步骤都会出幺蛾子。3. 预处理全流程3.1 HE切片清洗与tile切割HE切片到达手里时通常是SVS、NDPI、TIFF这类金字塔格式OpenSlide库是处理它们的事实标准。直接用OpenSlide读取全图Level 0不仅内存吃不消也没有必要。标准做法是读取一个合适的缩放级别比如目标放大倍率对应的Level然后在这个级别上切tile。切tile之前必须做组织区域分割。不要试图在全白背景上切图那样只会把大量无信息tile送进模型拖慢训练速度还可能导致模型学到“背景偏白就预测成某种组织”的虚假模式。我自己常用的方法很简单在低倍率下读一张下采样图转成灰度用Otsu大津法做一个全局阈值分割再做几次形态学开闭运算去掉零散噪点生成tissue mask然后把这个mask映射到目标Level只保留mask覆盖率达到阈值的tile。import openslide import numpy as np import cv2 slide openslide.OpenSlide(TCGA-XX-XXXX-01A-01-TS1.svs) # 读低倍率大图用于组织分割 level_thumbnail slide.get_thumbnail((512, 512)) thumb_rgb np.array(level_thumbnail.convert(RGB)) thumb_gray cv2.cvtColor(thumb_rgb, cv2.COLOR_RGB2GRAY) _, tissue_mask_low cv2.threshold(thumb_gray, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU) tissue_mask_low cv2.morphologyEx(tissue_mask_low, cv2.MORPH_CLOSE, np.ones((11, 11), np.uint8))常见参数是tile256×256或512×512level对应20倍或10倍放大。什么时候用20倍当肿瘤异质性高、想保留细胞核形态细节时。什么时候用10倍当你想扩大感受野、让每个tile包含更多组织结构时。老实说这个选择跟任务和模型结构都有关我通常先切20倍如果模型训练指标偏低再考虑混合倍率。从切割逻辑上还有一个细节tile之间需不需要overlap如果后续要做强分类或回归overlap能够缓解tile边界断裂问题但它会显著增加数据量而且如果按tile划分训练集和验证集overlap区域会造成严重的数据泄漏。我建议训练初期不设overlap先用干净的低冗余数据跑通基线。切完的tile还需要过滤。常见的过滤规则有三条tissue覆盖率低于5%的丢弃Laplacian方差过低图像太模糊的丢弃饱和度过高且大面积无结构纹理的丢弃。前两条最常用。每丢掉一批tile记录丢弃原因和数量这些统计信息是你排查数据问题的第一手材料。3.2 染色归一化Reinhard与Macenko怎么选HE切片最大的“特色”之一就是颜色不稳定。不同医院、不同染色批次、不同扫描仪产出的切片色相分布差异非常大。如果不做染色归一化模型很容易学到染色批次特征而不是组织形态特征。这一点在跨数据中心测试时会直接表现为性能暴跌。处理染色归一化主流方法有两个Reinhard方法和Macenko方法。Reinhard方法的思路非常简单——把图像从RGB转到LAB颜色空间然后把每个通道的均值和标准差匹配到一个目标图像上。它实现容易速度快但它假设色彩差异仅由全局亮度和饱和度差异引起对HE这种双通道染色的情况不够物理。Macenko方法则是基于染色分离的思路HE染色本质上由苏木精和伊红两种染料的叠加构成如果能在光学密度空间里估计出染色浓度矩阵就能把每张切片的染色浓度重新调整到标准模板的水平。这个更接近HE染色的机理在大多数情况下效果更好。我自己在项目里优先使用Macenko。给你一个简单的方法示例这里用staintools库import staintools # 选择一个“标准”的参考切片 target staintools.read_image(reference_tile.png) target staintools.LuminosityStandardizer.standardize(target) normalizer staintools.StainNormalizer(methodmacenko) normalizer.fit(target) # 对切好的tile逐张归一化 tile staintools.read_image(tile.png) tile_std staintools.LuminosityStandardizer.standardize(tile) tile_norm normalizer.transform(tile_std)但这里有几个必须注意的坑。第一个坑是参考图像的选择直接影响全部数据的颜色分布一定要选一张染色典型、组织区域丰富、没有大面积空白或折叠的切片。第二个坑是对于染色极浅或者裁剪到角落里几乎没有组织的tileMacenko的染色矩阵估计可能失败输出出现异常伪影。第三个坑是一些FFPE老切片的染色已经完全退化任何归一化方法都救不回来这类tile直接过滤掉不要勉强。3.3 基因表达矩阵的清洗与标准化基因表达侧的处理没有图像那么花哨但有它自己的讲究。TCGA下载的RNA-seq数据有count、TPM、FPKM等多种格式。对于回归任务我推荐使用TPM然后做log1p变换log2(TPM1)更常见因为基因表达的原始值跨度极大几个高表达基因比如一些核糖体蛋白基因会把Loss主导不转换基本没法稳定训练。标准化之前先做基因过滤。一般规则是在所有样本中表达量过低比如TPM 1且表达样本比例过低的基因直接去掉。这里没有绝对阈值建议结合任务需求来裁——如果项目不关心低表达基因你可以保留平均TPM前5000的基因这能有效降维并减少噪声。基因矩阵还需要处理批次效应。TCGA样本跨越多个测序中心和三/四种测序平台不同平台之间的表达定量存在系统性偏移。处理batch效应有两种策略一种是在数据层面用ComBat、limma的removeBatchEffect这类方法做校正另一种是把批次信息作为模型输入的一部分。前者简单直接但可能抹掉部分生物学差异后者保留完整信号但引入额外变量。我目前倾向在回归任务里先用limma做轻量校正同时把测序平台作为配对的metadata记录方便后续做分层评估。最后是标签预处理。如果我们直接预测数千个基因的TPM每个标签的分布差异很大有的基因方差小、表达稳定有的基因只有少数样本高表达不利于稳定训练。一个实用做法是对每个基因的表达值做z-score标准化或者做rank-gauss变换把标签映射到近正态分布。训练结束后再把标准化参数存下来方便把预测值还原回原始表达单位。3.4 空间转录组的预处理可选但推荐如果你手头有Visium数据预处理流程要比TCGA复杂一些。Visium原始输出是SpaceRanger生成的一个HDF5格式的feature-barcode矩阵filtered_feature_bc_matrix.h5和一个包含每个spot像素坐标的tissue_positions_list.csv。图像侧SpaceRanger还会输出一张对齐过的HE图像spot坐标可以和图像坐标匹配。做迁移或验证时一个标准的做法是以spot为中心切一个适当大小的图像块比如224×224或299×299对应的标签就是该spot的基因表达值。注意Visium的spot看起来就像二维格子spot间距大约100μm但组织边缘的spot常常落到空腔内需要根据SpaceRanger输出的in_tissue字段过滤。Visium数据的另一个便利是它天然包含空间结构可以直接用预测得到的空间表达模式来评估模型——例如某个marker基因在形态学上对应的区域是否被模型正确预测出来。这一层验证是TCGA bulk数据提供不了的。4. 常见问题与排查技巧实录4.1 数据泄漏与划分分割的坑这是我见过发生频率最高、后果最隐蔽的问题。很多人先在tile层面切训练集和验证集一张切片前半部分进训练集后半部分进验证集。表面上看数据没有重叠但是同一张切片生成的tile之间有极强的相关性——它们来自同一个组织、同一染色批次、同一扫描仪。模型很容易记住这种“主场特征”导致验证集指标虚高。等模型换到外部数据集时性能立刻跳水。正确做法是按病人或至少按slide划分数据保证同一个病人的所有tile全部落在同一个集合里。在TCGA这种每个患者通常只有一张WSI的数据集里按slide划分和按病人划分基本等价但在有多张切片的少数样本上要注意患者维度的去重。另外TCGA有部分患者的正常组织样本与肿瘤组织样本来自同一个人如果同时拿来做训练也必须按patient_id做split防止信息跨样本泄漏。4.2 显存爆炸与dataloader效率即便用了linearized attentiontile数量也还是会很快逼近显存极限。很多实验不得已要用batch size1或2这是常态。这时候要解决的问题不再是模型太大而是GPU利用率太低。建议把预处理后的tile缓存成本地文件而不是在训练时实时读SVS切片。你可以将所有归一化后的tile存成HDF5或LPZ格式并提前附加所属slide_id和坐标信息。训练时用多进程Dataloader做数据读取避免CPU读盘成为瓶颈。# splits示例按slide_id划分 train_slides [TCGA-XX-0001, TCGA-XX-0002] val_slides [TCGA-XX-0003] # 对tile索引表做划分 tile_df.set_index(slide_id, inplaceTrue) train_df tile_df.loc[train_slides] val_df tile_df.loc[val_slides]如果显存仍然吃紧还有一个很实用的技巧梯度累积。把batch size设为1累积16或32步后再做一次参数更新等效于batch size 16或32。虽然BatchNorm之类的层在小batch下表现会波动但配合LayerNorm为主的Transformer架构问题不大。4.3 染色差异导致的批次特征污染如何判断染色归一化到底有没有起作用一个快速诊断方式切tile时顺便记录来源切片ID然后对tile做embedding再用聚类或PCA可视化看tile是否按切片来源聚集。如果聚类结果和癌种、组织类型完全无关而是按切片ID形成明显分群就说明颜色特征已经泄漏进表征了。这种问题在融合多个数据集、比如TCGA不同癌种一起训练时尤其严重。我建议在归一化后做一次statistical check统计每个切片的平均RGB值或HSV色相分布看看各批次的中心是否已经对齐。如果还有明显偏移就要考虑更激进的归一化参数或直方图匹配。4.4 元数据与barcode对齐错误TCGA的样本barcode是一长串字段比如TCGA-A1-A0SB-01A-01R-A099-07。其中第4-15个字符即TCGA-A1-A0SB-01A代表sampleA1-A0SB代表患者01代表组织来源01原发肿瘤11正常A代表分析版本。在合并WSI和表达谱时要用sample级别的barcode做join而不是用patient级别。如果只对到patient同一患者的多块切片和多份RNA-seq会被错误交叉配对。还有一个常见的坑同一患者可能存在多个aliquot的RNA-seq数据例如同一患者做了二次测序。遇到重复样本我的建议是保留质量更高的aliquot记录去重规则并在日志里写明。这些细节看似琐碎但会在模型复现和结果审计时救你命。5. 实操经验小结最后说几句真心话。这个项目做到后面你会发现模型结构只是冰山一角真正影响最终表现的是数据管线里一个个“小决定”。同样一套linearized attention架构数据清洗严格与否结果可能差出一大截。我个人花在预处理上的时间和模型调参时间差不多是一半对一半。每次实验前先把数据处理脚本固化下来记录每一个过滤参数和样本数量变化这样后面无论调整模型还是写论文方法部分都会很从容。再分享一个小技巧在做tile切割时把每个tile对应的原切片ID、坐标、组织区域覆盖率、模糊得分全部落盘成一个parquet或者csv。这份清单看似不起眼但它是你做样本筛选、可视化归因、模型错误分析的基础。没有它你只能靠猜去定位问题出在哪一步。这个方向后续其实还能扩展很多从单癌种扩展到多癌种泛化、引入更多形态学特征做多模态融合、或者把线性注意力模块换成更高效的长序列变体都有得做。但无论怎么扩展数据底座始终是根。希望这篇文章能帮你把那层最容易出错的地基铺稳。