做单细胞转录组分析的人每天打交道最多的图大概就是UMAP和热图。UMAP负责给你“讲故事”的轮廓而热图则负责把“基因表达差异”这件事摊开了给人看。但这个热图在单细胞数据面前往往非常“水土不服”——动不动就是上万个基因、几万个细胞真按传统的画法糊在一张图上基因名看不见分组规律也出不来最后只能截取一小块区域局部看整体信息几乎全部丢失。这篇是《单细胞基因可视化之热图的根本改造》系列的第二篇。上一篇聊了基础思路和简单调优这篇我直接分享一套我目前一直在用的完整改造方案重点解决三个问题基因怎么筛、颜色怎么映射、细胞/聚类怎么排序最终输出的是一张既能在组会讲清楚结论、又能直接放进论文里的“三维”热图。换句话说这次不是换个配色、调个字体那种小打小闹而是把热图的整个生成逻辑重新捋了一遍。1. 改造思路先想清楚一张热图到底要承担什么任务1.1 传统热图在单细胞数据中为何常常失效很多人画单细胞热图习惯直接把表达矩阵丢给pheatmap或者ComplexHeatmap跑完才发现几个很头疼的问题。第一矩阵太大了。单细胞数据动辄几千上万个基因、上万个细胞画出来的PDF文件可能有几百兆打开要卡半天。即便勉强渲染出来了基因名密密麻麻叠在一起完全无法阅读更别说从中发现什么规律。第二单细胞表达矩阵极度稀疏。绝大多数基因在绝大多数细胞里压根没有表达表达值为0的比例通常能到90%以上。如果保留0值去画热图整个图大面积都是一个颜色的底色真实的信息反而被淹没在“无表达”的海洋里。如果直接过滤掉0值又等于丢掉了一个重要信息——某个基因在哪些细胞里沉默、在哪些细胞里激活这本身就是一种规律。第三颜色的映射逻辑不对。传统热图对连续变量做得比较稳妥但单细胞数据的表达值分布极度偏态少数高表达基因的数字能到几千甚至几万而大多数表达值集中在个位数。如果不做截断或缩放颜色的渐变几乎会被极端值霸占低表达的差异完全拉不开。第四基因没有分层。大多数热图是“一个基因一行”地平铺但单细胞分析里基因之间是有结构的——某些基因簇构成一个模块共同在某个细胞亚群中激活。如果每一行都是平权摆放这种模块结构就体现不出来热图就退化成了一张散点大杂烩。1.2 改造要落地的四个设计原则我这次改造遵循了四条原则后面所有步骤都是围绕它们展开的。第一空间换信息。既然所有细胞画上去看不到规律那就换一种思路列方向不再是“单个细胞”而是“细胞聚类”。每个聚类先做内部基因表达均值再画热图整张图的横向规模一下子从几万降到十几个信息不丢但噪声少了很多。单个细胞的信息可以另外画一张小倍图来补充而不是强行塞进同一张图。第二基因必须筛选。画热图不是把所有基因都搬上去。筛基因的核心逻辑不是“表达量高”而是“能区分细胞/聚类状态的基因”。我习惯于用“高变基因 各聚类marker基因”的并集前者保证覆盖全转录组的变异性后者保证关键的细胞身份基因一定在场。这个组合能兼顾探索性和解释性。第三颜色必须受控。对表达值做行方向基因方向的z-score标准化然后强制截断在 ±2.5 这个范围再映射到颜色。这样处理之后颜色的变化完全由“聚类间差异”驱动而不是被某个基因在某个细胞里的绝对高表达绑架。零值和惰性值单独处理不参与渐变映射用灰色或特定记号标识。第四条也是最重要的一条聚类排序不能随便来。细胞列的顺序、基因行的顺序这两者的排列逻辑决定了整张图能不能“读出规律”。我自己习惯的做法是先对矩阵做层次聚类聚类的距离度量用correlation聚类算法用ward.D2然后再结合分群结果做手动微调。这样出来的顺序既反映了数据本身的结构又照顾到了解释的方便性。理清楚这些原则接下来的所有参数调整、代码细节其实都是在为这四条原则服务。2. 核心细节解析基因筛选、颜色映射与零值处理2.1 基因筛选的第一原则高变不等于有用很多教程教你先用Seurat::FindVariableFeatures()挑2000个高变基因然后直接画热图。这是一个可行的起点但单靠高变基因当热图的行有一个非常明显的问题。高变基因反映的是“在细胞间表达变异大”但它既包括真实的细胞类型差异也包括那些随机噪声大、或者只在少数细胞内非特异波动的基因。换句话说高变基因集合里常常混杂着不少“和分群结果无关”的基因。画到热图里这些基因的条纹一团乱会严重干扰看图视线把真正关键的分类基因的条纹湮没掉。我的做法是高变基因和marker基因联合分两步。第一步用FindVariableFeatures取特征基因通常取3000个。第二步用各个聚类群的marker基因来“补充”。具体说跑完FindAllMarkers之后对每个cluster取avg_log2FC最高的前20~30个基因所有cluster的marker取并集。然后和3000个高变基因再取并集。这里有一个细节marker基因中如果有表达量极低、只在零星几个细胞里有信号的基因建议过滤掉。判断标准可以看pct.1如果一个marker在目标cluster里也只有不到10%的细胞表达那它作为marker的可靠性就要打问号了画出来也多半是一坨淡色块看着热闹实则没有信息量。如果样本分组复杂、一次要画几组对照也可以进一步对基因数量做约束。我自己一般控制在120~200个基因之间。太少了缺少全局感太多了图又回到拥挤状态。值得一提的是这个数量跟聚类数有关系10个聚类和5个聚类marker并集的大小差异是很大的所以这个数量不是死数是一个经验区间。2.2 颜色映射的数学细节为什么要做z-score截断热图的颜色本质上是一个从数值到色标的映射函数。如果不做标准化表达量数值范围是0到几千甚至上万映射函数只能把大部分区间分配给高表达的少数基因低表达区间的颜色过渡极其粗糙。标准做法是行方向z-score。对每个基因每行计算所有细胞或者所有聚类表达值的均值和标准差然后(x - mean) / sd。这样每行基因的均值是0标准差是1颜色渐变就可以以0为中间点向两边展开。但z-score之后还有一个坑个别基因的z-score可能达到5甚至10以上。如果不做截断这些极端值仍会霸占色标两端低到中等差异的基因还是看不出颜色差异。所以必须做截断我一般用 ±2.5极端的情况用 ±3。截断后再映射颜色大部分基因的色差才能拉开。颜色渐变的方向上我习惯用“蓝-白-红”不连续映射函数相当于RColorBrewer的RdBu反向。下调用蓝色上调用红色零差异是白色。这样一张图扫过去上调基因和下调基因一眼就能分清。如果你喜欢viridis那种连续色阶也可以但从生物学阅读习惯来说“红上调、蓝下调”是最自然的认知习惯。需要注意颜色映射这件事必须分清楚“展示型”和“分析型”。分析阶段可以用viridis等连续色阶看趋势最终汇报或文章里原则上选择分档明确、色觉友好、打印出来也清晰的方案不要选那些需要在屏幕上反光才能看清的颜色。2.3 零表达值的处理既是噪声也是信息单细胞数据零表达比例高这是一个绕不开的问题。传统的连续色阶会把0值映射为蓝色或某种冷色调然后真正低表达、中等表达的细胞也都挤在类似冷色调区域导致整张图的下半部分全是颜色污染。我的处理策略是将基因表达层和“是否检测到”信息拆成两层来看。在做聚类水平热图的时候这个问题的严重性会小很多因为求均值之后一个聚类只要有部分细胞表达了某基因平均表达值就会大于0。但如果画单细胞水平的热图零值占比大就必须考虑掩码方案。我的做法是额外构建一个布尔矩阵——表达量等于0的位置标记为NA非0的位置保留原值。然后在Heatmap函数里设置na_col参数为浅灰色或白色这样零表达区域就不会参与颜色映射而是以独立的背景色呈现。这种处理的直接效果是基因如果有表达标的颜色是有意义的没表达的区域则清清楚楚地告诉你“这个基因在这里没启动”而不是假装是一个很低的连续值。这层处理还有一个进阶用法。在最终展示的图中我会在热图主体之外额外交代一个“百分表达率”图按聚类展示每个基因在有表达细胞中的阳性比例。把这个信息放到主热图的右侧或者下方读者能同时看到“表达有多强”和“有多少细胞在表达”两个维度也就是我标题里说的“三个视角的热图”——强度、广度、差异方向。3. 实操过程与核心环节实现从Seurat对象到一张可发表的图3.1 数据准备矩阵、基因筛选、聚类因子所有的可视化改造都离不开一个结构清晰、经过质控的Seurat对象。假设你已经有了一个seu对象已经跑过NormalizeData、FindVariableFeatures、ScaleData、RunPCA、RunUMAP、FindClusters。第一步是提取表达矩阵。这里有一个容易被忽略的细节如果你已经执行过ScaleData那么GetAssayData提取到的数据已经是scale后的数据。画热图时我通常重新提取原始归一化数据自己做后续处理因为ScaleData对全部基因强制定中心化而且数值会被scale.max截断再用来做热图的z-score会失真。# 原始归一化矩阵行是基因列是细胞 expr_mat - GetAssayData(seu, assay RNA, layer data) expr_mat - as.matrix(expr_mat)第二步筛选基因set.seed(42) hvgs - VariableFeatures(seu, nfeatures 3000) markers - FindAllMarkers(seu, only.pos TRUE, min.pct 0.25, logfc.threshold 0.5) top_markers - markers %% group_by(cluster) %% top_n(25, wt avg_log2FC) %% pull(gene) %% unique() genes_use - unique(c(hvgs, top_markers)) genes_use - intersect(genes_use, rownames(expr_mat))这里set.seed很重要因为有些步骤会涉及随机抽样设置种子可以保证结果可复现。第三步把细胞按照聚类信息准备好cell_info - data.frame( cell colnames(seu), cluster as.character(Idents(seu)), stringsAsFactors FALSE )3.2 构建聚类均值矩阵大规模单细胞热图的关键压缩这是整个改造最核心的一步。全细胞热图如果你坚持要画也不是不行但那是一张“细节图”需要在几十万像素的大画布上慢慢看。大多数场景下我们要的是“规律图”——聚类与基因之间的关联。先做聚类均值矩阵cluster_level - lapply(unique(cell_info$cluster), function(cl) { cells_in - cell_info$cell[cell_info$cluster cl] Matrix::rowMeans(expr_mat[genes_use, cells_in, drop FALSE]) }) cluster_mat - do.call(cbind, cluster_level) colnames(cluster_mat) - unique(cell_info$cluster)这个矩阵的行是筛选出来的基因列是聚类。每个数值是某个聚类中某个基因的平均表达量。然后做行方向的z-scorescale_rows - function(x) { rm - rowMeans(x) rs - apply(x, 1, sd) rs[rs 0] - 1 (x - rm) / rs } mat_scaled - scale_rows(cluster_mat) mat_scaled[mat_scaled 2.5] - 2.5 mat_scaled[mat_scaled -2.5] - -2.5这个步骤里面有一个极其重要的细节sd为0的基因必须把分母换成1否则会出现NaN。这些基因通常是那些在样本中表达恒定、或完全不表达的基因如果不处理后面聚类会直接报错或者断层。3.3 三层面的排序与注释让热图可解释的关键一步热图的列聚类顺序我不建议直接用自动聚类的输出。而是先跑一个层次聚类再根据聚类数和分群做调整。hc_col - hclust(dist(t(mat_scaled)), method ward.D2) col_order - hc_col$order这里用t(mat_scaled)做列聚类是因为我们要对聚类列之间的距离做衡量。距离度量默认是欧氏距离实际效果还可以但如果你希望更稳健一点可以换成as.dist(1 - cor(mat_scaled))来做相关性距离。两种我都试过correlation距离的聚类结果更稳定不容易被单个基因的噪声影响欧氏距离对表达尺度更敏感有时候会把表达丰度相近而模式不同的聚类强硬归到一起。行基因的排序同样重要。按层次聚类结果排序后可以手动再把几个关键的“明星基因”调整到显眼位置必要时用row_split把基因按模块分隔开。这个分隔动作在ComplexHeatmap里很简单在语义上却很重要——它把“基因模块”的概念实体化了读者能直接看到一条一条的功能模块竖切信号。注释信息我一般加三层顶部放聚类标签左侧放基因类别注释marker基因、高变基因等右侧放每个基因在每个聚类中的“表达率”。top_anno - HeatmapAnnotation( cluster colnames(mat_scaled), col list(cluster cluster_colors), annotation_height unit(6, mm) )注意聚类颜色必须和UMAP图中用的颜色保持一致。这是组会和三方沟通时最容易被质问的问题为什么热图的聚类颜色和UMAP对应不上这个很讨厌务必一开始就统一。3.4 主热图、表达率副图、富集条目的三联拼装这是我做这个热图改造之后形成的一个固定套路。我把它叫做“三联拼装”实际使用的ComplexHeatmap用户可以理解成多个Heatmap对象拼在同一个画布上。第一联主热图就是上面生成的聚类均值z-score热图。第二联表达率热图。同样的基因和聚类组合数值换成该聚类中表达该基因的细胞百分比。注意这里的数值范围是0到1不需要再走z-score直接用连续色阶比如白到深紫或白到深绿就行。pct_mat - lapply(unique(cell_info$cluster), function(cl) { cells_in - cell_info$cell[cell_info$cluster cl] expr_sub - expr_mat[genes_use, cells_in, drop FALSE] rowMeans(expr_sub 0) }) pct_mat - do.call(cbind, pct_mat) colnames(pct_mat) - colnames(mat_scaled)第三联富集条目面板。对每一行的基因在其所属的模块里标注对应的富集条目比如GO term或KEGG通路或者至少把marker基因的来源标注出来。这一联可以用文本注释或板块色块实现。实际做法是根据基因模块的划分结果把模块对应的富集条目写在右侧如果一个模块有多条注释就选最显著的2~3条。拼装代码大致长这样ht_list - Heatmap( mat_scaled, name z-score, col col_fun, cluster_rows hc_row, cluster_columns hc_col, show_row_names TRUE, row_names_gp gpar(fontsize 6), top_annotation top_anno, column_split col_split_factor, border TRUE ) Heatmap( as.matrix(pct_mat), name pct, col c(white, darkgreen), cluster_rows FALSE, cluster_columns FALSE, show_row_names FALSE, width unit(2, cm), border TRUE ) rowAnnotation( module module_anno, col list(module module_colors), width unit(2, cm) )画完以后别忘了出矢量图pdf(heatmap_three_dimension.pdf, width 12, height 10) draw(ht_list) dev.off()这套流程跑通后从数据处理到出图基本稳定在15分钟以内而且图的信息量非常大组会上放出来别人通常第一反应是“这个数据结构很清楚”。4. 常见问题与排查技巧实录4.1 z-score之后出现NaN或者整行全空这是最常踩的坑。原因基本都出在sd 0的基因上。这类基因通常是“在所有聚类中表达量完全一样”的基因它们对可视化没有任何贡献直接删掉就好。我的经验是在取交集之后顺便加一步filter把标准差为0的基因剔除掉避免后面所有步骤连环报错。4.2 聚类顺序每次跑都不一样层次聚类的顺序有时候会因为平局处理、数据微小变化导致分支顺序变化。如果希望结果稳定一个是在聚类前固定随机种子另一个是保存一次聚类结果后续直接复用hc_col$order和hc_row$order不要每次重新跑。我在项目里一般都会把hc_col、hc_row这两个对象保存成RDS方便调试阶段反复出图而不改变版面结构。4.3 颜色映射被几个高表达基因绑架如果某几个基因的表达量特别高其他基因全部挤在色标的一端基本可以判断是没有做z-score截断。把mat_scaled[mat_scaled 2.5] - 2.5这步加上去就好了。如果做完截断后还是偏看看是不是行列方向搞反了——z-score必须是按行基因做不能按列做。4.4 想做单细胞水平大图但矩阵太稀疏、太卡如果真的需要看单个细胞级别的热图我的建议是换一个作图思路不要用完整的热图库去渲染几万乘几万的矩阵而是先做降采样从一个聚类中随机抽50~100个细胞控制总量在一万以内再画。否则就算出图了文件巨大AI或者PDF软件打开都费劲细节也完全看不清。展示单细胞水平的 “原始样貌” 和展示聚类水平的 “规律模式”这两件事分开来处理是最省心的分工。4.5 富集条目不匹配或者注释位置跑偏富集条目的位置偏移大多数是因为行数不匹配。因为做了基因筛选和去重基因数量和富集分析输入的背景基因数量经常对不上。我的解决方法是所有注释信息在筛选之后统一重新排列用match()对齐基因名而不是直接相信cbind之后的天然顺序。这一点在做rowAnnotation时特别容易出问题建议养成用match(genes_use, rownames(expr_mat))保底的习惯。4.6 聚类数量和颜色映射不匹配导致美术灾难有些聚类数一多默认配色就会非常相近打印出来几乎无法区分。建议提前备好一套超过20种颜色的自定义色板按顺序取。同时把列注释中的聚类标签和UMAP中的标签颜色严格绑定不同图之间保持一致这一点在投稿时会被审稿人看得非常仔细。5. 从“能看”到“能讲”目前我在用的完整流程清单最后把整套流程再串一遍方便你直接对着做。这套流程我目前已经用在了两个单细胞项目里稳定性和解释力都经住了组会和外部合作的检验。第一步确定分析目标。是展示细胞类型之间的差异还是展示某个处理条件引起的变化。目标不同基因筛选的侧重会不一样前者偏重marker基因后者可能更适合差异表达基因。第二步整理数据对象。确认Seurat对象的聚类、注释信息是否准确这一步需要结合UMAP、marker列表、先验知识综合判断。聚类注释错误后面画得再漂亮也是错的。第三步筛选基因。FindVariableFeatures取3000个叠加各个聚类的top marker基因对标准差为0的基因做剔除。第四步构建聚类均值矩阵并z-score处理。记住截断 ±2.5这是核心。第五步聚类和排序。列用相关性距离加ward.D2行用同样的方式聚类再根据基因模块做切分和手动调整。第六步组装注释和副图。列注释放聚类标签和样本来源行注释放基因模块副图放表达率热图最后用富集条目收尾。第七步导图输出。先出PDF或SVG矢量图再根据目标期刊或汇报屏幕的需求转成300dpi以上的PNG。第八步反馈迭代。把图发给合作者看收集“这块看不懂”“这个颜色什么意思”之类的反馈回到第三步调整基因列表或注释方式。坦白说热图可视化这个事真正难的不是运行代码而是清楚自己想展示什么信息。单细胞数据的信息密度太高一张图不可能塞下所有内容。改造的本质就是有选择地舍弃一部分信息突出另一部分信息。我踩了很多坑之后最大的体会是一张好的单细胞热图并不是把所有数据都画出来而是让看的人只需要三秒钟就能读出你要讲的结论——哪个基因在哪个细胞类群中特异表达、趋势是什么、可信度如何。能做到这一步这张图就算是改造成功了。