单细胞转录组的数据分析做到一定阶段,大家一定会撞上同一个瓶颈:矩阵拿在手里,几万个细胞、两万多个基因,PCA跑完还是分不开群。我印象很深的一次,是某份肝病样本的建模任务里,两群已知功能不同的细胞在PCA图上黏成一团,换了特征选择也拆不开。后来把深度学习引入降维步骤,用自编码器压缩后的隐向量重新聚类,边界一下就出来了。这篇文章是“生物深度学习”系列的第二篇,我们把焦点收窄到单细胞RNA测序(scRNA-seq)数据上,从预处理、网络结构、损失函数到下游聚类验证,完整走一遍实操链路。适合已经认得Python、跑过简单神经网络,但还没真正把深度学习用在测序数据上的朋友。内容会比较多,建议边看边动手。
1. 为什么第二篇首选单细胞转录组:高维稀疏数据给的三个难题
1.1 单细胞表达矩阵长什么样
单细胞转录组测序的核心产出,是一张细胞乘以基因的表达矩阵。细胞数从几千到几百万不等,基因数通常在两万以上。矩阵里的数值代表每个细胞中某个基因检测到的转录本数量,也就是表达量。听起来和普通表格没区别,但这份数据有三个让传统机器学习头疼的特点。
第一是高维。两万个基因意味着每个样本(细胞)落在两万个特征空间里,远高于图像等常见数据。第二是稀疏。表达矩阵里零值的占比经常在90%以上,大量基因只在一小部分细胞中有检测信号。第三是技术噪声。不同样本的测序深度、建库日期、操作批次都会带来系统性偏差,生物学差异和技术差异混在一起,很难直接分辨。
1.2 为什么PCA、t-SNE这些传统手段会失效
PCA的数学本质是在原始特征空间中寻找方差最大的线性方向。问题是细胞状态之间的转变,比如发育过程里一个谱系分化为另一谱系,往往呈现非线性关系,多基因之间是复杂协作而非简单线性叠加。PCA只能捕捉线性结构,遇到“开关式表达”“梯度渐变”这类模式时,压缩后的主成分会丢失大量关键组织关系。
t-SNE和UMAP虽然能做非线性映射,但它们通常作用在原始表达矩阵或线性降维结果上。原始维度太高的时候,距离计算会被成千上万个近乎不表达的基因稀释掉,几个真正的关键基因表达差异被“淹没”在零值里。这时候不管怎么挑困惑度参数,可视化的团块结构都很难反映真实的细胞身份。
1.3 自编码器给了另一种“压缩-重建”范式
自编码器的思路很直接:用编码器把高维表达谱压成低维隐向量,再用解码器从隐向量重建原来的表达谱。训练目标是让重建结果尽可能接近原始输入。如果一台网络能够用几十个数值重建几万个基因的表达模式,那说明这几十个数值确实抓住了细胞的核心特征。
这个思路和PCA的最大区别在于非线性。编码器和解码器之间用带非线性激活函数的隐藏层连接,所以可以拟合基因之间复杂的协作关系。更棒的是,扩张结构很灵活,可以把样本来源、测序深度这类技术变量作为额外输入,让模型在压缩时主动排除批次信息。这也是单细胞领域越来越多研究把自编码器当作降维主力工具的原因。
2. 数据进网络前,预处理这关躲不掉
2.1 细胞质量过滤:垃圾进则垃圾出
很多人觉得深度模型能自动处理噪声,其实在测序数据上这是个危险误区。一个破裂的细胞,线粒体基因占比会异常高,核基因信号严重丢失;一个空液滴,可能只捕获到环境RNA。这些不合格的细胞如果混进训练集,网络会把破损状态当作有意义的生物学状态学进去,污染整个隐空间。
常规做法是设置三道门槛:每个细胞至少检测到一定数量的基因(例如500个),总转录本数达到某个下限(例如2000),线粒体基因占比低于20%。三条规则可以过滤掉大部分离群细胞。过滤阈值需要根据数据质量调整,不要死套。判断标准很简单:过滤后的矩阵中,绝大多数细胞的检测基因数落在合理区间,且没有明显的尾部离群。
2.2 标准化和取对数,到底在修正什么
测序深度的意思是,一个细胞里所有转录本的总数。两个细胞可能生物学上完全一样,但一个被测出5万条转录本,另一个只测出2千条,直接比较某个基因的计数根本没有意义。因此第一步要按深度归一化,通常把每个基因的值除以该细胞总转录本数,再乘以固定系数(比如10000),得到每万个转录本中的比例。
但这还不够。单个基因的表达量波动范围很大,从几个count到几千个count,直接作为网络输入会让数值大的基因主导损失函数。取对数(log1p)可以把动态范围压到适合神经网络处理的区间。处理结束后基因表达量大致呈连续分布,MSE这类重建损失才比较合理。注意,是否需要进一步做z-score标准化,视模型情况而定。自编码器的输入层可以直接吃log后的值,标准化反而可能放大零值之间的差异。
2.3 挑高变基因:不是数据越多越好
两万个基因全部喂给网络不是不行,但效率很低且噪声很大。大量基因在细胞间几乎是常量,它们不仅帮不上区分细胞类型的忙,还会稀释信号。标准做法是先筛选高变基因(HVGs),通常取2000到5000个。具体选多少取决于数据复杂度——已知细胞类型多,可以留5000;只看两三种细胞状态,2000足够。
筛选原理是找到表达量在细胞间波动最大的一批基因,常用方差稳定转换后的离散度排序。这一步本质是多做一次特征工程,即使深度学习声称能自动学特征,把维度和噪声降下来,训练稳定性和下游聚类效果都会明显提升。
2.4 训练集和验证集的划分很容易被忽略
无监督自编码器训练时,很多人把所有细胞一股脑全丢进去训练。这没有致命问题,但你会失去判断过拟合的依据。我的习惯是随机留出10%的细胞作为验证集,训练时只更新训练集的梯度,每隔若干轮在验证集上算一次重建损失。如果验证集损失开始反弹,说明网络在记忆训练集噪声,该早停或者调正则化了。
另一个容易踩的坑是:不要在降维完成之后,为了展示效果又对整个矩阵重新做一次标准化。前后预处理不一致会导致隐空间向量坐标体系错位,下游聚类的分辨率、样本间可比性都会受影响。数据预处理完之后,后续所有环节都应该锁定在同一份处理好的矩阵上。
3. 网络结构怎么选:欠完备自编码器和变分自编码器的实际表现
3.1 先把自编码器拆开看
一个最基本的自编码器包含两部分。编码器输入维度等于基因数(比如3000),经过一两个隐藏层压缩,到达瓶颈层,这层维度就是隐向量维度,常用32或64。解码器从瓶颈层出发,经过对称的隐藏层,输出维度重新变回3000。损失函数衡量输出和输入的差距,通过反向传播同时更新两块模块的权重。
瓶颈维度是一个核心超参数。太小,信息塞不下,细胞细节丢失,聚类结果粗糙;太大,几乎不压缩,网络退化成记忆输入,隐空间没有泛化能力。我在不同数据集上试下来,32到128是比较稳妥的搜索区间。可以先固定64跑一版,再根据结果上下调整。
3.2 欠完备自编码器的典型翻车点
欠完备自编码器指隐维度严格小于输入维度,这是最基础的配置。它最大的问题不是结构,而是缺少隐空间约束。标准自编码器训练时只惩罚重建误差,对隐向量的分布没有要求。结果就是训练到后期,隐向量可能分布在多个互不连接的区域,细胞类型之间出现奇怪的断层。
另一个翻车点出现在网络层级太深、参数太多时。充足拟合能力的网络完全可以记住训练集中每个细胞的原始表达模式,验证损失却很高。更隐蔽的问题是,某些稀有细胞类型因为数量少,重建损失占比小,网络干脆放弃学习它们,分解键在损失函数设计上才能改善,这在下一节会专门谈。
3.3 变分自编码器为什么更贴合生物数据
变分自编码器(VAE)在隐空间上多了一个约束:假定隐向量来自一个连续的高斯分布。编码器不只输出一个向量,而是输出均值和方差,隐向量从该分布中采样得到。损失函数中增加KL散度项,逼迫隐空间分布接近标准正态。这个约束听上去复杂,效果却非常实用。
隐空间连续意味着,两个相似细胞在隐空间中的距离是平滑变化的。细胞分化轨迹、状态梯度这类生物学连续过程,在连续隐空间里能够得到更真实的刻画。而且采样过程带了一点随机性,让网络不能死记每个细胞的表达谱,等于自带正则化,对dropout噪声更鲁棒。目前很多单细胞深度降维框架都基于VAE思想,加入批次变量输入后效果很稳。
3.4 几个超参数去哪抄
网络层数不用贪多,2到3层隐藏层足够。单细胞表达模式复杂度有限,过深网络带来的提升很小,反而让训练不稳定。激活函数隐藏层用ReLU或LeakyReLU都行,输出层千万不要接ReLU,否则会把重建结果的下限锁在零,导致所有零值预测都是零,损失函数失去梯度意义。
Dropout建议设置在0.2到0.5之间,一开始可以取0.3。它能有效防止网络记住单细胞样本的稀疏噪声。批归一化或层归一化可以加在隐藏层之间,但要注意,推理阶段和训练阶段的归一化参数使用差异,不要在保存模型后混淆。我的建议是新手先把VAE跑通,再考虑加这些细节。
4. 损失函数与训练策略:重建不是越小越好
4.1 MSE对count数据的问题出在哪
MSE损失假设误差服从高斯分布,这在取对数后的表达谱上还算凑合,但面对原始count数据就很不合适。count数据是离散的,且方差常常大于均值(过离散),直接套MSE会把大量损失集中在高表达基因上,低表达但可能关键的基因完全得不到惩罚。
更符合数据特性的做法是使用负二项分布损失。让网络的输出层预测两个参数:均值和分散度,用负二项分布的对数似然衡量重建概率。这个损失函数能天然处理过离散和零膨胀问题,训练出的隐空间对稀有兴趣的基因保持敏感。缺点是实现稍复杂,需要自己推导稳定版本的损失函数表达。如果是初学者,我建议先保留log1p之后用MSE,至少能跑通链路;吃到苦头后,再决定要不要换成负二项分布。
4.2 KL散度权重不是固定不变的
VAE的总损失是重建损失加上一项KL散度,通常写作重建损失加β乘以KL项。β等于1是标准VAE,β大于1是β-VAE,即更强调隐空间的正则化。这个β直接影响最终隐空间长什么样。
β偏大,隐空间分布非常平滑,聚类之间界限模糊;β偏小,隐空间几乎没有约束,容易过拟合出多个孤岛。从我的经验看,先用β=1跑一版,观察聚类结果。如果细胞群挤在一起分不开,把β降到0.1到0.5试试;如果出现大量离群点或孤立团,把β往上调到2左右。这属于低成本高收益的调参策略。
4.3 训练过程的具体参数建议
优化器直接选Adam,学习率从1e-3开始,训练中如果损失下降放缓就衰减到1e-4。batch size建议128或256,设得太大会让GPU显存爆炸,设得太小则损失振荡严重。训练轮数不要拍脑袋固定,用early stopping来判断,监控验证集上的重建损失,连续10轮不下降就停止。
早停之后,再回头看看训练轮数和损失曲线。正常情况训练损失和验证损失应该同步下降,间距不大。如果训练损失已经降到很低但验证损失停滞,大概率模型容量过大,这时候优先减小隐藏层宽度,而不是加Dropout。因为Dropout能一定程度上缓解过拟合,但修改数据表达范围才是治本。
4.4 无监督学习里的“过拟合”怎么判断
无监督模型没有标签准确率可用,容易被忽视。我在项目里最常用的判断指标有两个。第一个是验证重建损失是否回升。第二个是隐空间形态:把隐向量可视化之后,正常训练会出现连续或明暗渐变的流形结构,过拟合时会出现大量离散的孤立点,每个点仿佛在“背诵”某一个训练样本。
还有一个很有效的检查做法:用同一份数据训练三次,每次改变随机种子,然后对比同一个细胞在三次实验中的隐向量坐标。如果三个坐标相距很远,说明模型不稳定,问题多半出在过拟合或学习率过高。反复做这个稳定性检查,能帮你提前发现很多隐蔽故障。
5. 从隐空间到细胞类型:聚类、可视化和marker验证
5.1 聚类:把隐向量交给社区发现算法
拿到隐向量后,别忘了原始的自编码器降维任务其实到这里还没结束。我们会把隐向量看作每个细胞的新特征表达,构建细胞之间的近邻图。方法是对每个细胞找它的k个最近邻居(比如k取15到30),用这些连接关系建成一张图,再跑社区发现算法如Leiden来划分模块。
这里经常有人犯一个错误:重建损失已经很小,隐空间似乎很好,但直接用原始表达矩阵聚类反而效果更好这类情况也发生过。原因是隐向量弱化了细胞间细微差距,某些罕见亚群可能被压缩到同一区域。所以在聚类之前,先看隐空间可视化,确认预期的分化关系大致存在,再调聚类算法的分辨率。
5.2 UMAP可视化:展示工具的边界要清楚
UMAP可以基于隐向量继续降维到二维平面,用于观察细胞亚群的分布。这里的n_neighbors和min_dist两个参数影响巨大:n_neighbors越大,展示的全局结构越清晰;min_dist越小,团块越紧凑。调整它们会让同一份隐向量呈现出完全不同的“视觉结论”。
必须强调,UMAP的二维坐标不能当作聚类依据。它的目的是展示而不是发现。我见过有人直接对UMAP坐标再做一次聚类,把坐标上的聚集点当成细胞类型,这完全是用展示工具的投影结果代替真实数据结构,很容易得到不可靠结论。正确顺序永远是:隐向量先聚类,再基于聚类结果指定颜色映射到UMAP图上。
5.3 marker基因验证:确认分群不是算法自嗨
聚类得到的数字编号没有生物学含义,需要去回找每个细胞群的marker基因。常规流程是对每个群与其他所有细胞做差异表达分析,列出表达量变化倍数高且显著的基因;再拿这些候选基因去对照已知的细胞类型marker表。比如免疫细胞里T细胞通常高表达CD3家族,B细胞高表达MS4A1和CD79A,髓系细胞高表达LYZ和CD14。
看到marker匹配良好,才能确认这个群对应了什么细胞类型。如果某个群的候选marker全是网上搜不到、注释不到功能的基因,多半是聚类过度细分造成的假群。这时应当调低聚类分辨率,让小的破碎群合并到邻近大群中,而不是强行保留并解释它。
6. 踩坑实录:批次效应、稀疏性处理与可复现性
6.1 批次效应没控制住,训练白费
批次效应指的是不同样本、不同实验日期、不同测序芯片带来的技术性表达差异。如果训练数据里包含多个批次,而模型没有显式接收批次信息,隐空间里很容易出现按批次来源分团的结构。最典型的临床样本场景是:本来要区分的肿瘤细胞和正常细胞没分开,反而是不同病人的细胞被完美分开——这种情况一看就知道模型学到的是技术差异。
对策不是事后指责数据质量,而是把批次信息作为额外输入加入网络。在VAE中,可以让encoder接收一个批次标签的嵌入向量,训练时告诉模型哪些细胞来自同一批次;或者用对抗训练的方式,让判别器无法从隐向量中猜出批次来源。需要权衡的是,批次校正过强会把真实的生物学差异也抹掉。建议处理前后做对比:校正后如果已知的marker基因差异仍然存在,就说明保留住了生物信号。
6.2 稀疏矩阵的“零”比表面更复杂
表达矩阵里的零值有两种来源:基因确实验不到,属于生物学上的不表达;或者基因实际有表达但测序没捕获到,这叫dropout事件。这两种情况混在一起,不能用简单的方式处理。直接把稀疏矩阵转成稠密数组去喂网络,会瞬间撑爆内存。两万个基因乘以十万个细胞,稠密矩阵是十几GB规模,普通机器根本扛不住。
实操中应当保留稀疏矩阵格式,按批量读取数据,训练循环里每次取一个批次并转为稠密后进入网络。输出层的设计也要配合稀疏性。如果使用MSE,隐层的输出直接预测log后的表达值即可;如果使用负二项分布损失,输出的是均值参数,而不是硬预测某个基因“有无表达”。训练过程中也可以对每批次内的零值加权降低,让网络把更多注意力放在非零信号上。
6.3 可复现性:同一份代码为什么不同机器结果不一样
深度学习的随机性在单细胞数据分析里是个大问题。不同机器、不同显卡、甚至同一台机器的不同线程调度,都会让隐向量产生细微差异。对于生物分析,这种差异如果传导到聚类,可能改变细胞类型的判定结果,完全不可接受。
解决办法是在代码开头统一固定随机种子,并关闭框架中的自动优化搜索。还要注意多线程数据加载:固定num_workers并给每个worker设定独立的随机种子。模型训练过程中定期保存隐向量和超参数记录。数据版本也要锁定,预处理后的矩阵计算一个哈希值,避免某天别人用了一版更新过的数据复现不出你的结论。
6.4 一次真实的挫败经历,也是最有价值的一课
某次做某份正常胰腺组织的单细胞数据,我信心满满地搭了一个VAE,训练损失曲线很漂亮,重建损失降到了足够低。结果把隐向量拿去做聚类,发现最大的两个群分别是“第一次测序的细胞”和“第二次补测的细胞”。再细查,数据合并时基因符号大小写没有统一,部分基因被当成两个不同的特征写进了矩阵。
那次排错花了整整一个下午。最终的修复不是调任何神经网络参数,而是回到预处理环节,统一基因命名格式并重新跑完整套流程。修复之后,VAE给出的细胞分群与已知marker高度吻合。这件事给我一个很深的体会:深度学习模型在测序数据上跑得再顺,也只是整条链路的最后一环;数据前期的标准、合并、命名、质量控制,哪一个都比网络结构更可能决定成败。
最后再分享一个我长期养成的习惯:正式训练全量数据之前,先随机抽取几千个细胞,快速跑一个简化版VAE,把隐空间可视化出来。肉眼确认大体的分群拓扑没问题,再上全量数据和完整训练流程。这个习惯帮我省掉了大量在深夜重新训练的时间。单细胞数据加深度学习的组合现在仍然有很多细节要靠人肉试错,但把每一步的原理和常见坑都弄清楚,至少能把试错周期从“几天”压缩到“几小时”。这些经验,希望对正在做类似工作的你同样有用。