上一篇写的是R入门阶段最常用到的数据整理和小技巧留言区里问得最多的反而不是语法而是“为什么我的包装不上”“为什么同样的代码在别人电脑上能跑到我这里就报错”。这篇接着写我把最近半年在转录组分析、时间序列建模和日常编程题里频繁遇到的真实问题重新捋了一遍每一个都配了可以直接跑的代码和排查思路。适合已经会基础操作、但一碰到真实数据就卡壳的人也适合想看看R在不同领域里怎么落地的朋友。1. 环境是头号敌人装包失败的三种典型姿势先说环境因为这一块占了我收到的问题里将近一半。很多人写代码没问题打开项目第一件事就是装包结果包装不上后面的活儿全卡住。1.1 分清渠道BiocManager、devtools和install.packages不是一回事很多新手最纠结的一件事是同一个包为什么有的用install.packages(xxx)装不上换个方式就装上了。原因很简单R包有三个主要来源对应三种完全不同的安装方式CRAN上的正式发布包用install.packages()这类包比较稳定大多和统计、数据处理有关。Bioconductor上的生信包用BiocManager::install()转录组、单细胞、富集分析这类包基本都在这里。GitHub上的开发版包用devtools::install_github(用户名/仓库名)通常是作者还没正式发布的最新代码。热词里出现了“r biocmanager 必备”这个说法很准确。BiocManager的价值不只是换个安装入口它会自动根据你当前的R版本去匹配对应的Bioconductor版本然后把依赖包一并装好。手动用install.packages()去装生物信息包经常会出现“依赖不存在”的连环报错。如果你在国内网络环境建议提前配置镜像能省下大量等待时间。CRAN镜像和Bioconductor镜像是分开配的options(repos c(CRAN https://mirrors.tuna.tsinghua.edu.cn/CRAN/)) options(BioC_mirror https://mirrors.tuna.tsinghua.edu.cn/bioconductor)把这两行写进~/.RprofileR每次启动都会自动加载。实际测试下来清华源和中科大的源速度都比较稳。装包的本质是下载预编译二进制文件或源码并编译镜像源决定了下载速度选一个离你近的镜像比反复重试要有效得多。1.2 库路径权限那个“lib C:/Program Files/R/R-4.0.2/library”警告Windows上非常经典的一个问题用install.packages()安装包弹出类似这样的警告Warning in install.packages : lib C:/Program Files/R/R-4.0.2/library is not writable原因是R默认安装到了C:\Program Files目录下这个目录对普通用户是只读的而安装包需要往library目录里写入文件。解决办法有两条路第一条用管理员身份打开RStudio再重装。右键RStudio图标选“以管理员身份运行”装完再恢复正常启动。这个方法简单但不推荐长期用因为后面每次装包都要管理员权限很别扭。第二条配置个人库路径这才是更规范的做法。R允许把包装到用户自己的目录里优先级甚至比系统库更高.libPaths() # [1] C:/Program Files/R/R-4.0.2/library # 创建个人目录 dir.create(D:/Rlibs, showWarnings FALSE) .libPaths(c(D:/Rlibs, .libPaths()))然后把.libPaths(c(D:/Rlibs, .libPaths()))写进.Rprofile之后所有包都会默认装到D:/Rlibs再也不会碰到系统目录权限问题。install.packages()后面也可以跟lib参数指定目录但治标不治本不如一劳永逸。Linux上类似的问题表现为“permission denied”常见于用系统R时直接install.packages()包试图写入/usr/local/lib/R普通用户没有写权限。和Windows一样Either加sudo要么建一个个人目录加.libPaths()。更推荐后者因为sudo install.packages()会把包目录的owner都变成root后续维护很麻烦。顺带提一句热词里那个“r语言rgdal”是这个问题的典型延伸。rgdal这个包因为依赖GDAL系统库装起来一直很痛苦非要用需要先装好libgdal-dev、libproj-dev这些系统依赖。而且它已经停止维护、被CRAN归档了现在做空间分析首选是sf包同样的功能占用更少、安装也简单得多。1.3 00LOCK锁文件、并发安装和Docker数据卷权限“R锁”这个问题非常隐蔽。安装包的时候如果中途报错或者你同时开了多个R进程安装同一个包R会创建一个名为00LOCK-包名的目录来锁定安装位置防止多个进程写入冲突。但有时候另一个R进程卡住了或者安装中断锁文件没删掉后续再装就会报ERROR: failed to lock directory D:/Rlibs for modifying解决办法是在确认没有其他R进程还在运行的情况下直接删掉对应的00LOCK-包名目录即可。Windows下也可以调出任务管理器看有没有R.exe或Rscript.exe进程残留有就先结束再删锁。这个坑不难但一旦遇到会让人误以为包本身坏了排查大半天。Docker环境里还有个高频问题和R本身无关但常常被当成R报错来问permission denied while trying to connect to the Docker daemon socket at unix:///var/run/docker.sock这通常不是R的问题而是当前用户不在docker组里。解决方案是把用户加进docker组sudo usermod -aG docker $USER然后重新登录一下。如果是把宿主机的数据目录挂载进容器还经常出现容器内写不了文件的权限问题尤其当容器的用户ID是1000时应执行sudo chown -R 1000:1000 ./data把宿主机目录的所有权改到1000容器内就能正常读写。这个做法在跑R容器、Shiny容器的时候非常实用命令本身和R无关但能省掉很多莫名其妙的“写入失败”。2. 转录组表达量为什么FPKM要转TPM代码怎么走热词里有一条很典型“转录组测序fpkm值换算成tpm r语言步骤”。转录组分析里这是几乎绕不开的一步但很多人只记住了“要转”不知道为什么要转。2.1 FPKM的两处硬伤FPKM全称是Fragments Per Kilobase Million它先按基因长度归一化再按测序总深度归一化。问题是它把两步标准化反过来做导致一个数学缺陷同一个样本里所有基因的FPKM加总起来并不等于一个常数而且不同样本加总值差异很大。这意味着用FPKM去比较不同样本间同一个基因的表达量会被样本本身的测序深度波动干扰。TPM的全称是Transcripts Per Million它先按每个转录本的长度归一化再把样本内所有转录本的值求和后统一缩放成“每一百万条转录本里有多少条”。这样做的结果是每个样本的TPM总和恒定为1e6跨样本比较时天然消除了测序总深度差异。所以现在做样本间比较、做热图、做相关性分析几乎都用TPM而不是FPKM。一句话总结FPKM适合看单个基因在单个样本里的绝对表达TPM适合做跨样本的比较。差异分析时前者容易引入假阳性。2.2 R代码从FPKM换算TPM如果你手里已经有一份FPKM矩阵换算TPM的公式极其简单TPM_i FPKM_i / sum(FPKM) * 1e6R语言里用向量化写法几行就能完成df - data.frame( gene c(gene1, gene2, gene3, gene4), length c(1200, 3500, 800, 2000), fpkm c(12.5, 8.2, 20.1, 5.4) ) df$tpm - df$fpkm / sum(df$fpkm) * 1e6 head(df)更规范的做法是从原始counts矩阵直接计算。比如featureCounts出来的结果有每个基因的counts和长度这时候的计算过程是# 先求每千碱基的counts数再归一化到百万 df$rpkm - df$counts / (df$length / 1000) df$tpm - df$rpkm / sum(df$rpkm) * 1e6实际转录组项目里基因长度可以从GTF注释文件导出或者用GenomicFeatures::exonsBy()计算出每个基因的外显子总长度。这一点需要注意如果用不同的注释版本基因长度会有差异最终TPM也会略有不同文章里一定要写清楚用的是哪个GTF。很多人还会问要不要做log2转换。一般做热图或PCA前会做log2(tpm 1)加1是为了处理0值避免负无穷。这个转换仅仅是为了可视化或降维不影响后续差异分析里的统计方法选择。2.3 从表达量到α多样性归一化思想相通热词里连着出现了“α多样性r语言”这其实是生态学/微生物组那边常用的一套分析但底层的归一化思想和表达量标准化完全一致。α多样性描述的是单个样本内部的物种丰富度和均匀度常用的有Shannon指数、Simpson指数。R里最常用的是vegan包library(vegan) # 假设行是样本列是OTU/物种值是相对丰度或绝对丰度 otu_table - matrix( c(10, 20, 30, 5, 15, 40), nrow 2, byrow TRUE ) colnames(otu_table) - c(OTU1, OTU2, OTU3) # 自动计算Shannon、Simpson diversity(otu_table, index shannon) diversity(otu_table, index simpson)用的时候有一个关键点diversity()默认要求每个样本的测序深度一致或者已经做了归一化。如果你输入的是绝对丰度而没有抽平算出来的指数不是标准化的样本间不能直接比。通常先做vegan::rrarefy()抽平或者用decostand()做某种标准化。这正好呼应FPKM转TPM的逻辑——任何时候做跨样本比较第一步一定是归一化不然后面的指数没有任何可比性。3. 组间GO富集与t-SNE降维两个高频报错的完整排查单细胞和转录组分析里除了表达量换算还有两个几乎一定会碰到的坑一个是GO富集分析时报缺包另一个是t-SNE降维时出现保护栈溢出。3.1 GO富集分析clusterProfiler的组间工作流与getoptlong缺失组间GO富集分析R里的标准工具是clusterProfiler。它可以把你的差异基因列表在基因本体Gene Ontology的各个分类下做超几何检验判断哪些生物学过程被显著富集了。基本流程library(clusterProfiler) library(org.Hs.eg.db) # 人类基因注释库 gene_ids - c(TP53, EGFR, MYC, BRCA1, CDK2) ego - enrichGO( gene gene_ids, OrgDb org.Hs.eg.db, keyType SYMBOL, ont BP, # 可选 BP/MF/CC pAdjustMethod BH, pvalueCutoff 0.05, qvalueCutoff 0.05 ) # 查看结果 head(ego) dotplot(ego)组间比较可以分两路走一路是对每个组的差异基因分别跑enrichGO然后肉眼对比通路另一路是用compareCluster一次性比较多组# geneClusters列存储每个基因属于哪个比较组 comp - compareCluster( geneClusters ~ group, data df, fun enrichGO, OrgDb org.Hs.eg.db ) dotplot(comp)很多人照着教程跑到这里第一次就会遇到这个报错ERROR: there is no package called getoptlong这不是你的代码写错了而是clusterProfiler在后端绘图和富集可视化时依赖了getoptlong这个Bioconductor包。装依赖的时候如果网络抖动或者用了CRAN源去碰Bioconductor的包就会出现这种缺失。解决办法非常直接BiocManager::install(getoptlong) library(getoptlong)装完之后重新加载clusterProfiler一般就正常了。这类“缺了xx包”的问题正确的排查顺序是先BiocManager::install(缺少的包名)再看这个包是不是已经停止维护、需要换替代方案而不是去改自己的分析代码。曾经遇到一位朋友把整段富集代码改了三遍也没解决最后只是缺一个依赖包而已。3.2 protection stack overflow别急着加内存单细胞t-SNE分析里有一个报错特别吓人Error in protect(): protection stack overflow我第一次遇到时还以为是内存不够加到一个64G内存的机器上重跑依然报错。后来才明白这个“保护栈溢出”不是内存不足而是R内部用来跟踪对象引用的保护栈空间满了。常见触发场景是数据矩阵特别大、单次表达式特别深、某个函数递归调用层级过多。t-SNE这类算法在计算高维矩阵的相似度时很容易一边递归一边保护大量中间对象直接把保护栈撑爆。三种实际可用的解法拆长表达式。把一长串管道或嵌套写成中间变量分步执行减少单层保护的临时对象数量。提高R的保护栈上限。Linux下启动R时加参数R --max-ppsize500000或者在代码里用options(expressions 5000)但后者只影响循环展开和控制流相关设置。如果栈溢出是C语言层面的递归导致的还要先在shell里执行ulimit -s 65536再启动R给线程栈腾出更多空间。降低输入规模。t-SNE本质是先用PCA降维再计算高维相似度没必要保留全部两万多个基因。标准做法是先做PCA取前几十个主成分再跑t-SNElibrary(Rtsne) # 大矩阵建议先用稀疏矩阵存储 # mat_sp - Matrix::Matrix(mat, sparse TRUE) # 先log归一化再做PCA expr_log - log1p(expr_matrix) pca_res - prcomp(t(expr_log), scale. TRUE) # 取前30个主成分跑t-SNE set.seed(42) tsne_out - Rtsne(pca_res$x[, 1:30], perplexity 30)这样一改t-SNE输入维度从可能的上万降到30保护栈的压力指数级下降。实际测试中一个一千个细胞、两万个基因的表达矩阵做完PCA再跑t-SNE内存占用只有直接跑的四分之一左右而且速度更快、聚类效果通常也更稳定。单细胞项目里基因矩阵本身就是高度稀疏的直接转成dgCMatrix这种稀疏格式再参与计算比硬扛稠密矩阵要合理得多。4. SARIMA建模从白噪声到完整预测的实战顺序时间序列部分热词里出现了“sarima模型r语言”。SARIMA本质是ARIMA加了季节性差分适合月度、季度这类有固定周期的数据。很多人学的时候看了一堆理论真到自己上手反而不知道第一步做什么。这里给一个可以直接照抄的流程。4.1 什么时候必须用SARIMA而不是普通ARIMAARIMA模型的核心是自回归项、移动平均项和差分阶数它假设数据是平稳的或者经过d次差分后平稳。但ARIMA不处理周期性。比如月度销售数据每年12月冲高、1月回落这种周期规律ARIMA根本抓不住因为它只考虑前后相邻时间的依赖关系。SARIMA在ARIMA的基础上多了一组季节分量记作SARIMA(p,d,q)(P,D,Q)[m]其中m是季节周期长度。月度数据m12季度数据m4。判断方式很直接画一下时间序列图如果每年同一个月份总有相似的波峰或波谷就说明有季节性需要用SARIMA。4.2 建模步骤与R代码完整流程建议按下面顺序执行library(forecast) library(tseries) # 1. 构造时间序列对象 # 假设月度数据从2018年1月开始 sales - ts(data_vector, frequency 12, start c(2018, 1)) # 2. 平稳性检验 adf.test(sales) # p值小于0.05说明序列平稳否则需要差分 # 3. 自动判断差分阶数 ndiffs(sales) # 非季节差分阶数 nsdiffs(sales) # 季节差分阶数 # 4. 模型拟合 fit - auto.arima( sales, seasonal TRUE, stepwise FALSE, # 关闭逐步搜索结果更优 approximation FALSE # 关闭近似拟合结果更准 ) # 5. 残差诊断 checkresiduals(fit) # 残差应该是白噪声不能有明显的自相关 # 6. 预测未来12期 fc - forecast(fit, h 12) autoplot(fc)这套流程里最容易出错的是第2步。adf.test()原假设是“序列不平稳”所以p值小于0.05才说明平稳了。如果p值大于0.05就需要对原始序列做一阶或季节差分比如diff_sales - diff(sales, lag 12) # 季节差分 adf.test(diff_sales)有人读完auto.arima的输出会困惑为什么它选了个特别高或特别低的阶数auto.arima靠AIC、AICc、BIC这类信息准则在候选模型里挑一个最优的加了stepwise FALSE会搜索更大范围耗时更久但结果通常更可信。样本量不大的时候计算量完全可以接受。4.3 差分过度和周期误判是两个最常见的坑第一个常见问题是过度差分。差分是为了让序列平稳但差分次数不是越多越好。差分一次去除趋势差分两次可能引入原本不存在的伪自相关。普通一阶差分加一次季节差分已经覆盖了绝大多数情况如果ndiffs()返回2或更大就要警惕序列是不是已经有问题看看是不是有异常值或突变点。第二个常见问题是季节周期设置错误。月度数据频率是12季度数据频率是4但遇到周数据和日数据时很多人会搞混。判断周期最实用的办法是看谱图或者画ACF图找周期性峰值。ACF图在每个整数周期位置出现明显峰值比如第12、24、36条对齐基本就是月度序列的季节周期证据。我自己建模的习惯是模型出来后先不急着看预测值而是盯着残差图看。如果残差里还有明显的周期性说明季节分量没提取干净如果残差有趋势说明差分阶数不够。这一步做扎实了预测结果才敢拿到业务讨论里用。5. 一道字符串处理题拼数number的R解法热词里藏着一道非常经典的编程题“小 r 正在学习字符串处理。小 x 给了小 r 一个字符串 s。”加上后面能搜到的高频词“c. 拼数(number)”应该是同一类问题给你一堆数字让你把它们拼接成一个最大或最小的数。这道题我在不少社群里看到大家用Python、Java写反而R版本很少这里补充一下。5.1 题目描述与字典序陷阱题目大概是这样的给定若干个非负整数重新排列它们的顺序让它们拼接起来组成的数字最大。例如[3, 30, 34, 5, 9]最大结果应该是9534330。最容易掉进去的坑是直接把数字转成字符串按字典序降序排序。如果这样做3和30会排成30在3前面因为字符串比较时3和30先比首位3相等然后30就比3“长”所以30 3。但实际拼接时330 330而303 303前者更大。所以字典序排序直接跪了。正确的比较方式是比较两个字符串a和b的两种拼接结果。若ab ba则a应排在b前面否则反过来。这个比较规则满足全序关系直接作为排序比较器即可。5.2 自定义比较器的R实现R语言里的sort()函数不支持传入自定义比较器这是这类题在R里稍显麻烦的地方。如果数字数量很少可以用gtools::permutations()列出所有排列取拼接结果最大的那个library(gtools) nums - c(3, 30, 34, 5, 9) perm - permutations(length(nums), length(nums), nums) candidates - apply(perm, 1, function(x) paste0(x, collapse )) max(candidates)但排列复杂度是阶乘级一超过10个数字就爆了。更通用的办法是自己写一个基于比较的排序比如冒泡排序把自定义比较器传进去nums - c(3, 30, 34, 5, 9) strs - as.character(nums) # 自定义比较若 a 拼 b 大于 b 拼 a则 a 应该排前面 greater_str - function(a, b) { paste0(a, b) paste0(b, a) } # 冒泡排序按 greater_str 规则降序 bubble_sort_str - function(x, cmp) { n - length(x) if (n 1) return(x) for (i in seq_len(n - 1)) { for (j in seq_len(n - i)) { if (!cmp(x[j], x[j 1])) { tmp - x[j] x[j] - x[j 1] x[j 1] - tmp } } } x } sorted_strs - bubble_sort_str(strs, greater_str) paste0(sorted_strs, collapse ) # [1] 9534330冒泡排序时间复杂度是O(n²)数据量大时不划算但作为R的练习很直观。数据量更大时可以把自定义比较器包装成一个辅助映射转成R能原生排序的形式或者在Rcpp里写C比较器一个晚上搞定。R的向量化处理强在矩阵计算对这类两两比较的排序场景并不擅长认清这个边界反而能少踩坑。5.3 顺手补充一个开发效率快捷键写完这个自定义函数在RStudio里想给每个函数生成标准的文档说明不要手动敲。光标放在函数名的上一行或开头按CtrlShiftRRStudio会自动插入roxygen2的注释骨架# Title # # param x 参数说明 # return 返回值说明 # export # # examples这套注释配合devtools::document()可以直接生成R包文档。热词里有“devtools的ctrl加r”说的应该就是这件事。用惯了之后写函数顺手就把文档补了后面再用?函数名查看项目可维护性上了不止一个台阶。6. 两则向量化小题大象喝水与数值比较的R写法日常咨询里也会收到一些偏“算法入门”的小题。热词里有“一只大象口渴了要喝20升水才能解渴但现在只有一个深h厘米、底面半径r厘米的小桶”还有“r(ab)小于等于”这类判断式。这两个题都很简单但正好能看出R新手和老手写法的差距。6.1 圆柱体积大象喝20升水一桶够不够基础版本水桶是圆柱体底面半径r厘米高h厘米。体积公式是pi * r^2 * h得到的是立方厘米。1升等于1000立方厘米所以20升等于20000立方厘米。新手写法通常是一个循环r - 10 h - 30 vol_cm3 - pi * r^2 * h vol_liters - vol_cm3 / 1000 vol_liters 20但其实一旦要多试几组桶的参数就该立刻切换成向量化写法df - data.frame(r c(10, 15, 20, 25), h c(30, 40, 50, 60)) df$vol_liters - pi * df$r^2 * df$h / 1000 df$enough - df$vol_liters 20 dfR的向量化不是“推荐的优化技巧”而是这个语言的根本设计。数据框的一整列就是一个向量pi * df$r^2 * df$h / 1000会同时处理所有行。不要把Python里那种逐行循环的习惯带进来那会让代码又慢又啰嗦。6.2 从“判断(ab)是否小于等于C”看浮点陷阱热词里有个“r(ab)小于等于”大概是判断(a b) c这类条件。看起来简单实际有浮点数精度问题。比如a - 0.1 b - 0.2 (a b) 0.3 # [1] FALSE原因在于0.1、0.2和0.3在二进制浮点里都不能被精确表示0.1 0.2算出来是0.30000000000000004大于0.3所以返回FALSE。这在统计计算里非常要命因为很多过滤逻辑都依赖这种比较。稳妥的做法是加一个可容忍的误差范围is_less_equal - function(a, b, c, tol 1e-9) { (a b) (c tol) } is_less_equal(0.1, 0.2, 0.3) # [1] TRUER语言里的all.equal()提供了类似功能near()函数出自dplyr也是为浮点比较设计的。写数字判断时保留一个误差参数是好习惯尤其是做p值阈值过滤、表达量阈值筛选时很多看似“莫名其妙多一个”或“少一个”的样本都是浮点比较边界问题造成的。最后说点个人体会。R的报错信息经常不是病根“permission denied”“library路径”“不存在某个程辑包”这类问题八成出在环境而不是代码。我现在每开一个新项目固定会先执行一遍.libPaths()、sessionInfo()和镜像源检查把环境底子打好再跑分析。这个习惯帮我省掉了大量反复试错的时间。这个“R编程示例”系列后面打算写可视化和R包封装如果你在跑上面任何一段代码时遇到新的报错欢迎把报错原样留下来我下篇一起拆。