R语言做科研绘图做到后面你会发现二维图形能解决90%的问题但就是有那么10%的场景非3D不可。这个系列更新到第98期今天专门聊聊R语言里的3D形状绘制——从3D散点图、3D PCA投影到响应曲面再到怎么把三维图从屏幕上搬进论文。不只是贴代码我会把选型思路、参数踩坑、导出经验一并讲清楚。适合被“3D效果”吸引过来但还不确定自己需不需要3D图的科研党尤其是做转录组、微生物多样性、富集分析这类组学分析的朋友。先说明白一件事3D图不是万能药用对了是画龙点睛用错了就是给读者添堵。1. 先泼一盆冷水这些科研场景真的需要3D图吗1.1 3D图的真正价值结构展示与交互探索先说个反直觉的结论3D图在科研里的作用从来不是“精确展示数值”而是“快速传达分布结构和整体趋势”。举个典型例子。你做PCA降维如果前两个主成分只能解释60%的方差样本在二维平面上挤成一团组间差异看不清楚。这时候把第三主成分拿出来三个轴一起转往往就能发现原来被压扁在二维平面里的分组规律。这种场景二维图替代不了。再比如设计实验时考察两个连续因子的交互效应。二维等高线图虽然能表达但响应曲面的“山峰”和“低谷”在三维视角下给人的冲击力完全不同。你给导师或者合作者展示的时候一个可以旋转的曲面比一张静态等高线图直观得多。还有一种场景是几何结构展示比如分子构象、晶体结构、空间坐标数据。这种本身就是三维数据不用3D反而说不清楚。1.2 那些“伪3D需求”转录组、α多样性、富集分析有意思的是我最近收到的咨询里大量做转录组和微生物多样性分析的朋友会在FPKM换算成TPM之后或者算完α多样性指数之后来问“能不能用3D图画一下”。我的回答通常是大概率不需要。转录组样本间的差异关系优先用二维PCA、PCoA、聚类热图就能讲清楚。α多样性指标本质上是组间比较的分布图箱线图、小提琴图比任何3D散点都更精确。GO富集分析的结果气泡图、条形图是标准配置强行三维化反而丢掉了P值和富集因子的层级信息。这些场景有个共同点你关心的本质是分组差异是否显著而不是样本在三维空间里的几何结构。3D散点图在这种需求下只会让读者花更多时间去“找点”而不是看结论。所以当有人问我要不要学3D绘图我第一句永远是先拿二维图试试如果二维图已经把结论讲清楚了3D就是个伪需求。1.3 先判断场景再选技术路线如果确认要上3D优先判断属于下面哪一类动态探索类数据量大、分组多、需要旋转观察适合rgl或者plotly这类交互式工具。出版展示类论文插图、组会汇报静态图适合plot3D或者scatterplot3d这类偏静态渲染的工具。几何结构类分子、晶体、空间点云rgl是最顺手的因为它自带OpenGL渲染和材质控制。我的建议是探索阶段用交互沉淀结论时再画静态图。千万别屏幕上一转就截图转到最后连自己都不知道视角是从哪来的。2. 选型篇rgl、plot3D、scatterplot3d、plotly哪个才是你的菜2.1 四个包各自的门派和脾气R语言里能做3D图形的包不少但真正常用、文档全、社区活跃的主要是这四个。rgl走的是OpenGL交互路线。它画出来的图可以在窗口里用鼠标拖动旋转滚轮缩放命令风格跟基础绘图比较接近但很多函数的名字不一样。它的强项在于动态探索和几何对象构建——你可以画球体、椭球、凸包、光照表面还能组合出复杂的场景。缺点是出版级矢量图导出能力有限复杂场景经常导出失败。plot3D是静态出版向的代表。它从传统的persp和scatterplot3d思路发展而来提供了一整套统一命名的函数比如scatter3D、surf3D、ribbon3D、image3D。它是通过计算好的投影矩阵把三维坐标画在二维平面上不依赖OpenGL所以在任何环境里都能跑非常稳。输出走的是标准R绘图设备PDF、TIFF、PNG都跟普通二维图一样没有任何额外依赖。缺点是不能交互。scatterplot3d是我看过的代码量最小但效果最稳定的包。它专门做一件事画3D散点图。语法极其简单几分钟就能出图。缺点也明显——只适合散点曲面、椭球、光照这些功能基本没有。plotly虽然本身不是R原生的但R语言的plotly包可以直接把R数据转成交互的Web图形。它输出的是HTML/JavaScript可以嵌入Shiny应用或者R Markdown报告。如果你要给合作者发一个能自己转动的图plotly是最方便的选择。代价是语法跟R基础绘图差异较大有时数据需要转换成plotly的格式静态导出还要单独配置orca工具。特性rglplot3Dscatterplot3dplotly交互能力强无无强出版级输出一般高中高一般环境依赖需要OpenGL无无无学习成本中高中极低中典型场景探索、演示、几何结构论文插图快速散点图Web、Shiny2.2 我的组合策略我自己的习惯是看情况混用。拿到一批新数据先用rgl快速撒点转几圈找视角确认哪些维度的组合有看头。找到合适的投影角度之后再用plot3D重新画一张出版向的静态图。这样既享受了交互探索的效率又保证最后放进论文的图清晰稳定。如果只是内部快速看一眼分布不打算放进正式报告scatterplot3d就够了一行代码完事。如果是要做Web端的数据汇报plotly直接上配合Shiny还能让合作者自己选变量。这里顺便说一句很多人一上来就装rgl结果在虚拟机或者远程桌面上各种报错然后转头说R不能画3D。其实不是R的问题是环境问题。遇到这种情况要么回物理机跑要么换plot3D这种不依赖OpenGL的方案。记住一句话——工具服务于场景不是场景服务于工具。3. 实测一用rgl画一幅能转动的3D散点图和3D PCA3.1 先搞定rgl的OpenGL运行环境装rgl只需要一条命令install.packages(rgl)但装完不等于能跑。rgl依赖系统的OpenGL支持Windows物理机一般没问题macOS通常也正常最容易出问题的是虚拟机、远程桌面、云服务器这些没有完整图形栈的环境。安装完成后我建议先跑一个最小验证library(rgl) open3d() points3d(1, 1, 1, col red, size 8)如果弹出一个窗口并且能看到一个红点说明环境没问题。如果直接报错或者黑屏大概率是OpenGL不支持去检查显卡驱动或者干脆换plot3D方案。3.2 从iris开始最简3D散点图rgl画3D散点图的核心函数是plot3d参数跟基础绘图的plot很相似library(rgl) plot3d(iris$Sepal.Length, iris$Sepal.Width, iris$Petal.Length, col as.numeric(iris$Species), size 8, type s)这里的type s表示用球体表示点比默认的平面圆点更有体积感。col用as.numeric把因子转成数字rgl会自动映射到调色板。跑完这条命令你应该能看到一个可以拖动的三维散点图。鼠标操作要记牢左键拖拽是旋转滚轮是缩放右键拖拽是平移。这三个操作决定你能不能找到一个好视角。3.3 加分组椭球让类群差异一眼可见散点图只是第一步很多场景下需要在图上叠加统计结构。最常见的就是为每个组画95%置信椭球。rgl自带的ellipse3d函数就是干这个的scores3d - prcomp(iris[, 1:4], scale. TRUE)$x[, 1:3] plot3d(scores3d, col as.numeric(iris$Species), size 8, type s) for (sp in levels(iris$Species)) { idx - which(iris$Species sp) clr - c(red, green3, blue)[which(levels(iris$Species) sp)] ell - ellipse3d(cov(scores3d[idx, ]), centre colMeans(scores3d[idx, ]), level 0.95) shade3d(ell, col clr, alpha 0.2) }ellipse3d的核心逻辑是根据协方差矩阵和中心点生成一个椭球网格level 0.95控制椭球大小。shade3d负责渲染alpha 0.2让椭球半透明这样能看到内部的散点。加完椭球之后组间分离情况一目了然。3.4 3D PCA实战从表达矩阵到三主成分投影真正面对组学数据时我们通常是对表达矩阵或者丰度矩阵做PCA。这里模拟一个类似转录组FPKM/TPM数据结构的矩阵set.seed(123) n_genes - 2000 n_samples - 30 group - factor(rep(c(Control, Treat_A, Treat_B), each 10)) expr - matrix(rnbinom(n_genes * n_samples, mu 100, size 5), nrow n_genes, ncol n_samples) colnames(expr) - paste0(S, 1:n_samples) # 模拟分组差异 expr[1:200, group Treat_A] - expr[1:200, group Treat_A] * 3 expr[201:400, group Treat_B] - expr[201:400, group Treat_B] * 4 # 常规分析流程过滤、标准化、log转换、PCA keep - rowSums(expr 1) 5 expr_filt - expr[keep, ] expr_log - log2(expr_filt 1) pca - prcomp(t(expr_log), scale. TRUE) scores - pca$x[, 1:3] pca_var - round(100 * summary(pca)$importance[2, 1:3], 1) library(rgl) plot3d(scores, col c(darkorange, steelblue, darkgreen)[as.numeric(group)], size 10, type s, xlab paste0(PC1 (, pca_var[1], %)), ylab paste0(PC2 (, pca_var[2], %)), zlab paste0(PC3 (, pca_var[3], %))) # 给每个样本添加标签 text3d(scores, texts rownames(scores), adj c(0.5, 0.5), cex 0.7, family sans)注意几个细节。PCA之前一定要对表达矩阵做过滤和标准化里面的log转换参数要根据数据分布调整plot3d的xlab、ylab、zlab直接写主成分方差占比比只写“PC1”有信息量得多text3d用来标注样本名adj参数控制文字对齐避免文字和点重叠。加上分组椭球之后三个处理组的关系在三维空间里会非常清楚。3.5 组会演示利器录一段旋转动画rgl有一个很方便的功能可以把视角旋转过程录制成视频或者GIF。比如想做一个持续6秒、每秒10转的旋转动画movie3d(spin3d(axis c(0, 0, 1), rpm 10), duration 6, dir ., convert TRUE)这段代码会让图绕Z轴匀速旋转并把过程逐帧保存然后拼接成动画文件。注意movie3d依赖外部工具Windows下通常需要安装ImageMagick或者FFmpeg并配置好系统路径。如果不想折腾可以直接按住鼠标手动旋转一遍用录屏软件记录效果一样。动画非常适合组会汇报动态旋转的过程会让导师快速理解样本在三维空间的分布比放三张不同角度的静态图省事得多。4. 实测二用plot3D画响应曲面从网格数据到出版级插图4.1 曲面图的前提网格数据是怎么来的第二类高频需求是响应曲面图常用于研究两个连续因子对响应变量的交互影响。很多新手直接用原始观测点去画曲面结果画出来是一堆乱线。原因在于曲面图本质上是z f(x, y)的函数图像x和y必须是规则的网格坐标z必须是对应网格上的函数值矩阵。实际操作中网格数据有三种来源实验设计模拟生成、回归模型预测、插值计算。最常见也最规范的做法是用原始数据拟合一个模型然后在理论网格上计算预测值。4.2 模拟实验数据并拟合响应面这里用一个简单的二次响应面例子演示。假设两个因子x和y响应变量z满足z 0.5x² 0.3y² 0.4xy 随机噪声首先构造模拟数据library(plot3D) set.seed(42) x_obs - runif(30, -2, 2) y_obs - runif(30, -2, 2) z_obs - 0.5 * x_obs^2 0.3 * y_obs^2 0.4 * x_obs * y_obs rnorm(30, 0, 0.2) fit - lm(z_obs ~ x_obs y_obs I(x_obs^2) I(y_obs^2) x_obs:y_obs)这里的lm拟合实际上就是在做响应面回归。真实科研里x和y可以是温度、pH、浓度这类实际因子z是某个检测指标。接下来构造网格并计算模型预测值grid_x - seq(-2, 2, length.out 50) grid_y - seq(-2, 2, length.out 50) grid_dat - expand.grid(x grid_x, y grid_y) grid_dat$z - predict(fit, newdata data.frame( x_obs grid_dat$x, y_obs grid_dat$y )) z_mat - matrix(grid_dat$z, nrow length(grid_x), ncol length(grid_y))注意expand.grid把所有x坐标和y坐标两两组合生成完整的网格。predict在网格上给出模型预测值。z_mat必须被转成矩阵列数对应y的个数行数对应x的个数否则surf3D会报错。4.3 surf3D出图曲面、原始点、投影的一次成型plot3D的surf3D函数可以直接画网格曲面。先把曲面画出来再把实测点叠加到同一坐标系上surf3D(x grid_x, y grid_y, z z_mat, col lightblue, theta -30, phi 20, xlab Factor X, ylab Factor Y, zlab Response, border grey40, lwd 0.3)theta和phi控制观察视角theta是水平旋转角phi是仰角。负值表示观察方向在另一个象限。再叠加原始观测点scatter3D(x_obs, y_obs, z_obs, add TRUE, type p, pch 19, col black, cex 1.2)add TRUE是plot3D体系里直接把图形叠到当前画布上的关键参数。原始点画成黑色实心圆点贴在曲面上读者可以直观看到观测值落在曲面附近的程度。4.4 底层等高线投影哪个区域响应最高曲面图的一个常见增强操作是在底部平面上叠加等高线投影相当于给曲面加了一个“海拔地图”。plot3D里可以用contour3D或者直接在底面上画线实现scatter3D(x_obs, y_obs, z_obs, colkey FALSE, col black, pch 19, theta -30, phi 20, xlab Factor X, ylab Factor Y, zlab Response) surf3D(x grid_x, y grid_y, z z_mat, col heat.colors(100, alpha 0.8), add TRUE, border NA) contour3D(x grid_x, y grid_y, z z_mat, colvar z_mat, col black, add TRUE, zlim c(min(z_mat), min(z_mat) 0.01), levels seq(0.4, 1.6, by 0.2))注意contour3D里的zlim被我设成接近曲面的最低点这样等高线就会被压到底面附近。levels参数控制画哪些等高线这里每隔0.2画一条。配合colkey参数控制颜色图例的显示这张图放到论文里既有曲面又有等高线很多期刊会喜欢这种组合尤其是做工艺优化、环境因子交互分析的研究。5. 决定3D图“能不能看”的细节颜色、光照与视角5.1 颜色不是装饰是数据编码3D图的颜色比二维图更容易让人产生误读。原因很简单第三维本身已经叠加了深度信息如果颜色再乱来读者的视觉系统就过载了。连续的第三维变量比如曲面高度或者PC3得分建议用连续的色标。viridis包是我这几年最推荐的它对灰度打印友好而且对色盲读者也相对友好library(viridis) plot3d(scores, col viridis(100)[cut(scores[, 3], 100)], ...)离散分组变量要选区分度高的颜色。注意别直接用红配绿很多读者和审稿人是红绿色盲。用色盲安全的调色板更稳妥。如果曲面图用默认的rainbow色虽然颜色丰富但色相的跳跃会让人误以为存在边界而实际上数据是平滑变化的。这是很多新手图“看起来花哨但经不起推敲”的原因。5.2 光照和材质曲面有没有立体感就看这里同样的曲面光照参数不同立体感天差地别。rgl里对材质和光源的控制非常精细library(rgl) open3d() surface3d(grid_x, grid_y, z_mat, col steelblue, lit TRUE, specular white, shininess 80, alpha 0.9) light3d(theta -30, phi 30, viewpoint TRUE) bg3d(white)lit TRUE开启光照计算specular white给材质加高光shininess越大高光越集中。light3d用来设置光源位置viewpoint TRUE表示光源跟随视角变化这样旋转图形时光照方向始终对着观察者不会出现转到背面就全黑的尴尬。曲面图如果不开光照默认是均匀着色看起来就像一块彩色塑料板。开了光照之后山峰和谷底会因为受光角度不同产生明暗变化立体感立刻增强。这个细节对响应曲面类图形的展示效果提升非常明显。5.3 视角找到一个信息量最大的角度3D图形的最大坑是“看不见”某些角度下一个组完全被另一个组挡在后面你以为样本缺失了其实只是视角问题。rgl可以用rgl.viewpoint精确控制视角rgl.viewpoint(theta -30, phi 20, fov 60, zoom 0.8)plot3D里则是通过theta和phi参数实现同样效果。我自己的经验是旋转图形时先找一个能看到所有组散布范围的角度确保三个坐标轴的刻度都可见然后小幅旋转微调直到没有明显的遮挡。判断标准很简单——在最终视角下每个组的点都至少能看到一半以上。有个小技巧如果遇到严重遮挡可以适当降低点的尺寸或者提高透明度的区分度。3D图被遮挡是常态制图者的职责是选一个遮挡最小的角度而不是把所有点都清清楚楚地画出来。5.4 标签与图例别让读者猜坐标3D图的坐标轴标签比二维图更容易被忽略。很多人画完rgl图轴标签没加、单位没写、图例没有读者根本看不懂坐标代表什么。rgl里加标签用text3dtext3d(max(scores[, 1]), 0, 0, PC1 (41.2%), family sans, cex 1.2)注意rgl对中文字体的支持比较弱直接用中文经常显示成方框。如果必须用中文可以尝试family sans但更稳妥的做法是英文标签或者先在R里设置好全局字体。图例建议用legend3d放在三维场景里或者直接导出后在图片编辑软件里叠加——后者更灵活不容易遮挡数据点。plot3D这边就省心很多scatter3D自带colkey参数控制图例xlab/ylab/zlab直接标注坐标轴出版向的需求基本都考虑了。6. 论文插图导出与典型报错修复实录6.1 先想清楚最终输出形式画3D图之前先想清楚这张图最终出现在哪里是投期刊的静态插图还是组会上的动态演示还是给合作者的网页链接。这个决定直接影响到导出方式。期刊插图通常需要高分辨率或者矢量图。矢量图好处是缩放不糊但3D场景有大量渐变、透明度、光照计算导出矢量格式很容易出现兼容性问题。我个人的折中方案是用高分辨率PNG/TIFF导出配合后期软件补充标注绝大多数期刊都能接受。动态演示则简单得多录屏成视频或者直接用rgl的movie3d生成动画文件又或者导出成WebGL的HTML页面让合作者自己转动。6.2 rgl的导出快照、矢量、网页三选一rgl导出静态图最直接的方式是快照rgl.viewpoint(theta -30, phi 20, zoom 0.8) rgl.snapshot(3d_pca.png, fmt png, width 2000, height 1500)先固定视角再快照这个顺序千万别反。很多人旋转完图直接快照导出来角度不对只能重画。rgl.snapshot支持高分辨率2000像素级别的文件放进论文完全够用。如果确实需要矢量图rgl有一个rgl.postscript函数但限制很大只支持点线类的简单场景。复杂一点的光照曲面或者椭球导出会丢失元素甚至报错。所以我不太推荐把rgl的复杂场景导出成PDF/EPS。想要交互式结果可以用rgl的WebGL导出功能writeWebGL(dir webgl_figure, width 800, height 600)生成的HTML文件可以用浏览器打开读者可以自由旋转缩放特别适合补充材料或者个人主页展示。6.3 plot3D的导出最稳的传统R绘图路线plot3D没有交互需求导出方式跟普通二维图完全一致这是它作为出版方案的最大优势pdf(response_surface.pdf, width 6, height 6) surf3D(grid_x, grid_y, z_mat, col lightblue, theta -30, phi 20) dev.off()如果要投期刊很多编辑部要求600dpi以上的TIFFtiff(response_surface.tiff, width 6, height 6, units in, res 600, compression lzw) surf3D(grid_x, grid_y, z_mat, col lightblue, theta -30, phi 20) dev.off()lzw压缩可以在不损失质量的前提下减小文件体积这个细节对投稿系统很友好。6.4 踩坑一OpenGL不可用导致rgl黑屏这个坑我在前面提过但值得单独展开。rgl在虚拟机、远程桌面、云服务器上最常见的报错是“OpenGL support not detected”或者打开窗口后一片黑。原因虚拟化环境通常没有完整的GPU图形栈OpenGL上下文创建失败。排查顺序先跑open3d()和points3d(1,1,1)这种最简命令确认是不是环境问题检查显卡驱动Windows物理机更新驱动后基本能解决虚拟机用户建议开启3D加速功能如果实在解决不了换plot3D方案。plot3D是纯软件计算投影不依赖OpenGL在任何环境都能稳定出图6.5 踩坑二不存在叫“getoptlong”这个名字的程辑包这个报错在最近的热搜里出现得非常频繁很多人加载绘图扩展包时被卡住。先说明getoptlong本身不是一个3D绘图包它是一个Bioconductor生态里的命令行参数解析包但很多组学和绘图相关的Bioconductor包在依赖链上需要它。报错原因你安装的某个包声明依赖getoptlong但你只装了主包没装依赖运行library或者调用功能时就会报“不存在叫‘getoptlong’这个名字的程辑包”。正确解法if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(getoptlong)这里有个常见的歧路用install.packages(getoptlong)去装会报“package ‘getoptlong’ is not available”因为CRAN里根本没有这个包必须走BiocManager。装完这个依赖之后再重新加载当初报错的主包问题基本就解决了。这类“依赖地狱”在R生态里很常见。遇到缺失包的报错我的经验是先别急着到处搜答案用library命令确认整个依赖链上缺了哪些然后根据包来源CRAN还是Bioconductor选择正确的安装通道。6.6 踩坑三与四中文标签、颜色发黑中文标签变方框rgl对中文字体支持弱text3d默认字体不含中文字形。稳妥方案是英文标签或者先设置好全局字体。如果是plot3D在PDF设备里输出中文建议用支持嵌入字体的Cairo设备或者干脆保持英文。曲面颜色发黑这个现象常常出现在rgl绘制曲面时。原因一般是光照参数没有正确配置曲面材质在计算光照时没有得到足够的环境光。解决方法是设置环境光系数或者增加辅助光源material3d(ambient gray30, specular white, shininess 60) light3d(theta 0, phi 0, viewpoint TRUE)颜色发黑还有一个容易被忽略的原因曲面数据矩阵里包含大量NA值。NA区域在渲染时会显示为黑色空洞乍一看像光照问题。先检查z_mat是否完整再检查图形设置。最后想再多说一句3D绘图本质上是一个“用复杂度换信息量”的博弈。三维图形天然比二维难读它适合的是那些二维确实装不下的结构信息。这几年我在帮不同课题组审图、改图的过程里最大的体会就是一张好的3D图不是炫技不是为了酷而是让审稿人或者读者在5秒内get到你想让他们看到的结构特征。所以不管用什么包、什么参数先问自己一句这张图需要第三维吗需要的话就把第三维做到极致不需要就安心用二维。这个思路想清楚了今天聊的这些技术细节才有意义。