这篇是ArchR学习记录系列的第十二篇。前面几篇已经把数据从Fastq一路走到Peak Matrix、聚类、marker peak到了这一步十有八九会遇到同一个问题眼前的细胞群已经有了身份可这群细胞的开放状态到底是谁在背后调控ArchR里解决这个问题最常用的模块就是chromVAR——它能把每个细胞、每个转录因子基序对应的染色质可及性偏差变成一个连续的活性分数而不是简单统计一堆峰里有没有某个motif。这篇文章就把chromVAR的原理、ArchR里的标准跑法、可视化套路和报错坑一次讲透适合刚跑完单细胞ATAC-seq、想继续往下游做转录因子推断的同学。1. chromVAR到底在算什么先搞懂偏差分数1.1 从“峰里有没有motif”到“TF活不活跃”先说动机。单细胞ATAC-seq给你的本质是一个“峰×细胞”的稀疏计数矩阵Peak Matrix。每一格代表某个细胞里某个染色质峰上有多少可及性信号。这个矩阵能告诉你哪些区域开着但没法直接告诉你哪个转录因子在调控。转录因子发挥作用的位置通常是对应基序motif所在的那一小段序列而基序可以用PWM来描述。ArchR的addMotifAnnotations做的工作就是把每个峰的基因组序列拿去做PWM匹配生成“这个峰是否含有某个TF的基序”的0/1标记看起来只是给Peak Matrix加了一层注释。如果按传统思路接下来很多人会干一件事统计每个簇里有多少峰含有某个motif然后算个比例或做个超几何检验。这个做法不是不行但它在单细胞层面上很吃亏。因为不同峰的GC含量、平均可及性、以及细胞本身的开放程度差异很大简单统计会把技术偏差当成生物学信号。chromVAR的核心思路正好用在刀刃上它不回答“这个motif在这个簇里富不富集”而是回答“在扣除峰的可及性背景后这个细胞里带这个motif的峰是否比预期更开放”。这个更开放的信号才是转录因子活性更合理的代理指标。1.2 偏差分数和z-score先认识这两个数字chromVAR输出两个你一定会碰到的量deviations和z。deviations是原始偏差分数直观理解是“观察到的信号”和“预期信号”的比值取对数。每个TF motif会对应一组“匹配峰”每个细胞在这组峰里的真实信号是observed然后chromVAR会用一组背景峰做预测得到expected。如果某细胞里GATA1 motif相关峰的实际可及性明显高于这类峰应有的平均水平deviations就是正的反过来就是负的。z是在deviations基础上做的标准化相当于把偏差分数除以标准差便于跨motif比较。ArchR里跑完chromVAR后MotifDeviations矩阵里同时存了这两个assay我实际使用中几乎都拿z做后续的可视化和比较。用个生活例子班里批作文不能直接拿字数比谁写得好因为有人字体大、有人字间距大。chromVAR做的就是先按字体大小和排版分行也就是GC含量和平均可及性把同水平的作文分到一组再比谁实际写的内容更超出预期。这个“超出预期”的幅度就是偏差分数。1.3 与简单motif富集分析的区别为了帮大家避开分析思路上的混淆我把普通motif富集和chromVAR偏差分析放在一张表里对比维度普通motif富集chromVAR偏差分析输入差异峰/某簇的峰集合峰×细胞计数矩阵 motif注释核心问题某组峰里是否富集某motif某细胞/某群中某TF活性是否高于预期输出p值、富集倍数每个细胞每个motif的连续偏差/标准化z分数分析单位通常是亚群单个细胞可聚合到群受到GC/开放度影响需要额外校正内置背景峰校正适合场景组间差异峰的motif构成单细胞水平的TF活性推断这两者不是替代关系而是视角不同。做HOMER式motif富集时你已经在问“这群特异性开放的峰里面有没有某个TF的偏好位置”而chromVAR是在问“这群细胞里某个TF的调控信号是否真的比背景活跃”。很多项目需要两种一起上先用chromVAR挑候选TF再用传统富集佐证生物学含义。2. ArchR里跑chromVAR的完整流程2.1 前置检查PeakMatrix不齐后面全是空谈在跑addMotifAnnotations和computeDeviations之前先确认项目里已经有PeakMatrix。ArchR的chromVAR流程是基于峰集和细胞计数矩阵的如果你只做了addReproduciblePeakSet但忘了addPeakMatrix后面一定会报“PeakMatrix not found”之类的错。我习惯先跑一句getAvailableMatrices(proj)正常应该能看到PeakMatrix。如果不在补齐proj - addPeakMatrix(proj, force TRUE)这步会按细胞重新量化峰上的片段数是后续一切的基础。另外提醒一句之前做过的过滤很关键doublet、低质量细胞必须提前清理干净。chromVAR的背景峰模型是全局统计混入大量低质量细胞会把背景拉偏最后所有TF的z分数都失真。2.2 addMotifAnnotations把PWM打到全基因组峰上前置条件就绪后加motif注释proj - addMotifAnnotations( ArchRProj proj, motifSet JASPAR2020, name Motif, force TRUE )motifSet参数默认常用JASPAR2020也有JASPAR2018、cisBP等选择。ArchR内部会用TFBSTools做PWM扫描匹配到每个峰后生成一张二进制MotifMatrix记录“该峰是否包含该motif”。这一步产生的行名一般是Motif_1这样的编号TF名字放在元数据里需要时可以通过rowData查回来。首次运行JASPAR相关注释时需要联网下载对应的数据库。如果网络不稳或处在离线环境建议提前手动安装对应包或者准备好自己的PFM/PWM列表。自定义motif集合并不少见特别是做非模式物种或特定TF家族时可以用library(TFBSTools) pfmL - readRDS(my_pfms.rds) # PFMatrixList 或 PWMatrixList proj - addMotifAnnotations(proj, motifSet pfmL, name MyMotif)自定义列表时一个容易踩的坑是ID重复。如果两个PFM的ID相同ArchR构建MotifMatrix时会产生冲突后续rowData(devMat)$name会乱掉。我建议在导入前先用uniquify或name()自查一遍保证每条记录的唯一性。跑完注释后可以用getFeatures(proj, useMatrix MotifMatrix)快速看一眼是否生成成功。提示force TRUE是为了强制重跑注释。如果你换过PeakSet或者之前注释时用的数据库版本不对一定要加这个参数否则ArchR可能直接返回旧结果。2.3 背景峰的选择为什么每个motif都需要一组“对照组”chromVAR最核心的设计是背景峰校正。光有MotifMatrix还不够你还需要算出每个motif匹配峰的“对照组”也就是背景峰。bgdPeaks - getBackgroundPeaks(ArchRProj proj, n 200)getBackgroundPeaks会为每个motif的峰集合挑选一批在GC含量、平均可及性上相似的峰作为背景。n是重采样次数ArchR默认通常取200表示每个背景集合被重复构建200次。n越大背景估计越稳但耗时和内存也会同步上升。数据量特别大的时候有人会把n降到100我一般是在正式分析时保留200做快速探索时临时降到100先跑通。新版ArchR里也有addBgdPeaks这种写法本质是一样的但API在不同版本间有差异。如果你的代码报“Background peaks not found”或“bgdPeaksmust be non-empty”先查一下当前ArchR版本的函数签名。别在一棵树上吊死换个正确函数通常马上好。一个容易被忽略的细节背景峰是基于当前PeakSet和细胞集生成的。如果你后来改了PeakSet、过滤了更多细胞必须重新生成背景峰不能复用旧对象。否则对照背景与真实峰集不匹配z分数会被系统性污染。2.4 computeDeviations计算细节与资源规划motif注释和背景峰都有了正式计算偏差proj - computeDeviations( ArchRProj proj, backgroundPeaks bgdPeaks, force TRUE ) devMat - getMatrixFromProject(proj, useMatrix MotifDeviations)computeDeviations封装了chromVAR的核心逻辑ArchR会先校验输入矩阵的类型确保是整数型计数然后对每个细胞每个motif计算观察值、预期值和偏差分数。最终devMat是一个SummarizedExperiment对象含deviations和z两个assay行对应motif条目列对应细胞。这里必须提醒computeDeviations是整套流程里最吃CPU的一步。我跑过一个四万多细胞的PBMC项目motif数量在800左右单机八线程大概跑了四十分钟到一小时。如果细胞数超过十万建议把它放到后台Rscript里跑前端不要挂着IDE防止内存和会话状态互相干扰。nohup Rscript run_chromVAR.R run_log.txt 21 跑之前还可以用addArchRThreads(threads 8)把线程数拉起来但线程不是越高越好。每开一个线程就会多一份内存占用如果机器只有32GB内存建议线程数控制在4到6个否则可能跑到一半因为内存不足被系统杀掉。3. 拿到偏差分数后可视化与群体差异分析3.1 plotDeviations先看全貌再谈细节computeDeviations跑完后第一批图通常用plotDeviations出。它可以按细胞或按簇展示一个或多个motif的偏差变化p - plotDeviations( proj, name c(GATA1, CEBPA, RUNX1), pal paletteContinuous(solarExtra) )不过有一个小坑ArchR的motif名在MotifMatrix里通常带后缀编号比如GATA1_1、GATA1_2。同一TF可能对应多个PWM条目所以直接传GATA1有时候会匹配不上。稳妥做法是先从rowData(devMat)$name里查出对应的完整条目名再传给plotDeviations。plotDeviations返回的是一张热图每行是一个motif条目每列是一个细胞或一组细胞。如果按簇分组颜色越偏红说明该簇里这个TF的偏差越高。我一般不会只画一两个TF而是先画一批候选TF看整体图谱确认哪些在簇间有明显分隔再去细看信号。3.2 plotVarDev高变异TF才是真正值得追的候选全局热图适合看已知候选筛新TF还是得靠变异性排序p - plotVarDev(proj, plotName VarDev)plotVarDev输出一个点图横轴是每个motif纵轴是变异程度。变异性越高说明这个TF在细胞间活性差异越大越可能是驱动细胞状态分化的因素。实际操作中从最右侧那一撮点里挑Top 20或Top 50基本上能把项目里值得关注的TF圈出来。这个方法很适合用来讲故事你不需要对800个motif做一遍差异检验只用plotVarDev把注意力锁在少部分高变异条目上再进入下游验证。我个人觉得这比一上来就做全矩阵聚类更高效因为全矩阵聚类会被大量表达模式相似的“管家motif”干扰。3.3 把TF活性画到UMAP上偏差分数是连续值最好的展示方式之一就是迭到UMAP上。前提是项目里已经做了UMAP并且需要先计算插补权重proj - addImputeWeights(proj) p - plotEmbedding( proj, colorBy MotifMatrix, name GATA1_1:score, imputeWeights getImputeWeights(proj) )addImputeWeights会基于KNN对全矩阵做插补跑起来也挺耗时但只需跑一次之后画任意TF都可以复用。colorBy MotifMatrix表示从MotifMatrix里取特征后面的GATA1_1:score表示取该motif的插补后偏差分数。注意这里不是取MotifDeviations而是MotifMatrix中扩展出来的score特征ArchR内部会把偏差分数包装成可用特征。如果你发现在plotEmbedding里传name始终无效可以先运行getFeatures(proj, useMatrix MotifMatrix)把可用的特征名打印出来再复制完整的条目名进去省得再猜。UMAP迭加图的好处是信息量很大一张图上能看到某个TF的活性集中在哪个细胞亚群、是否有连续梯度、是否存在少见的高活性细胞团。这些信息在热图里容易丢失。3.4 组间比较一行代码背后的统计推断可视化只能直观判断真正要说“群体A里TF X活性显著高于群体B”还是得做统计比较。ArchR本身的差异分析功能更多面向peak和基因表达对偏差分数我一般直接把devMat取出来自己算devMat - getMatrixFromProject(proj, useMatrix MotifDeviations) z - assays(devMat)$z groupA - which(proj$Clusters %in% c(C1, C2)) groupB - which(proj$Clusters %in% c(C5)) # 以某一行motif为例做wilcoxon检验 motifIdx - which(rowData(devMat)$name GATA1) p - wilcox.test(z[motifIdx, groupA], z[motifIdx, groupB])$p.value这样单个TF遍历跑一遍不难但涉及几十个TF时注意多重检验校正。我通常会用p.adjust(ps, method BH)做一个简单的FDR控制。另一个更稳妥的做法是先用伪bulk的方式把同簇的细胞平均或求和得到伪bulk偏差分数再做组间比较避免把单细胞的不独立性问题带进检验。这个思路未必写进ArchR教程但实际在做课题汇报时会经得起审稿人追问。4. 高频报错、解读误区和调优经验4.1 高频报错速查表很多报错不是基因问题而是流程顺序和参数版本问题。我把自己反复遭遇过的几类整理成表便于快速定位报错/现象常见原因解决办法PeakMatrix not found没做addPeakMatrix或矩阵被覆盖补齐addPeakMatrix(proj, forceTRUE)Error in ... JASPAR2020数据库包未安装或网络失败手动安装对应数据库包或改为本地PFM文件seqlevels不匹配PeakSet的染色体命名和BSgenome不一致用seqlevelsStyle统一命名或指定正确的BSgenomeBackground peaks not found版本API差异或背景峰变空检查getBackgroundPeaks返回的对象重新生成computeDeviations长时间无输出线程/内存设置不合理降低线程数转后台Rscript运行z全为NA或大量NA某些簇细胞太少、峰覆盖太稀疏过滤过小簇或尝试提高背景峰稳健性最让我头疼的是染色体命名不一致如果PeakSet是chr1而BSgenome是1addMotifAnnotations里做序列提取时会失败或静默丢行。一定要在所有步骤前统一风格尤其做非人类样本时别忘了对应物种的BSgenome。4.2 解读chromVAR结果时必须记住的三条红线第一z分数高不代表TF一定结合了DNA。它只说明“带这个motif的开放峰活性高于预期”中间还有可及性到结合的距离。想坐实结合需要ATAC footprints、ChIP-seq或突变验证。写文章或作报告时建议使用“推断活性”而不是“结合活性”这类表述。第二同一TF的多个PWM条目不要随便合并。GATA1_1和GATA1_2可能是不同长度或不同来源的PWM匹配到的峰集合不一样直接合并会把噪声带进来。如果你确实想合并得按“共享峰的信号汇总”来做而不是简单把两个z相加。第三单细胞的z值非常吵。不要根据单个细胞的偏差分数下结论也不要拿一个细胞的z值把细胞硬分成“有活性和无活性”两组。正确姿势是先看群水平或UMAP插补结果用群体统计量给结论背书。4.3 实测调优从半天到两小时如果你跑大样本总觉得chromVAR慢我的建议按优先级排列。第一先过滤再计算。把所有可疑细胞、低质量和极小簇先清掉背景峰会干净很多计算量也会明显下降。第二getBackgroundPeaks中的n可以探索式降到100但正式结果最好保留一份200的版本做稳健性验证。第三把线程数和内存预算先算好我倾向于一个线程至少给3GB内存开8个线程就需要24GB以上可用内存不然只是看起来在用多线程。第四中间结果一定要保存saveRDS(devMat, devMat.rds) saveArchRProject(proj)devMat本身是一个SummarizedExperiment后续所有可视化、统计比较都能基于它重来不需要每次都重新跑一遍computeDeviations。保存好这个对象你后面就不必反复加载整个ArchRProject能省不少内存。5. 一些个人体会把这条系列写到第十二篇chromVAR是我觉得单细胞ATAC-seq分析里“投入产出比”很高的一个模块。它不需要额外造数据也不需要复杂的模型调试只要流程顺序正确出来的结果在生物学上往往能直接指向一批转录因子候选给下游实验设计提供线索。我现在的标准流程是先plotVarDev扫Top 20高变异TF再用plotDeviations看这些TF在簇间的整体分布挑出感兴趣的几个上UMAP验证空间模式最后用z矩阵做组间统计比较。这个顺序能尽量避免“看到什么说什么”的偶然性也能让p值检验建立在可视化预筛之上而不是反过来硬凑结果。还有一个很实用的习惯如果项目里有多个样本或多个条件的ATAC数据不要一上来就在全局所有细胞上跑一套chromVAR先按大群或者条件分层跑比较同一细胞类型在不同条件下的TF活性差异往往比全局偏差分数更容易讲清楚生物学问题。毕竟chromVAR算的是相对活性背景变了z值的含义也会跟着变。 
长河君 · 2026/10/8 17:55:00 · 阅读 1.2k 单细胞ATAC-seq转录因子活性推断：chromVAR原理与ArchR实操指南
这篇是ArchR学习记录系列的第十二篇。前面几篇已经把数据从Fastq一路走到Peak Matrix、聚类、marker peak到了这一步十有八九会遇到同一个问题眼前的细胞群已经有了身份可这群细胞的开放状态到底是谁在背后调控ArchR里解决这个问题最常用的模块就是chromVAR——它能把每个细胞、每个转录因子基序对应的染色质可及性偏差变成一个连续的活性分数而不是简单统计一堆峰里有没有某个motif。这篇文章就把chromVAR的原理、ArchR里的标准跑法、可视化套路和报错坑一次讲透适合刚跑完单细胞ATAC-seq、想继续往下游做转录因子推断的同学。1. chromVAR到底在算什么先搞懂偏差分数1.1 从“峰里有没有motif”到“TF活不活跃”先说动机。单细胞ATAC-seq给你的本质是一个“峰×细胞”的稀疏计数矩阵Peak Matrix。每一格代表某个细胞里某个染色质峰上有多少可及性信号。这个矩阵能告诉你哪些区域开着但没法直接告诉你哪个转录因子在调控。转录因子发挥作用的位置通常是对应基序motif所在的那一小段序列而基序可以用PWM来描述。ArchR的addMotifAnnotations做的工作就是把每个峰的基因组序列拿去做PWM匹配生成“这个峰是否含有某个TF的基序”的0/1标记看起来只是给Peak Matrix加了一层注释。如果按传统思路接下来很多人会干一件事统计每个簇里有多少峰含有某个motif然后算个比例或做个超几何检验。这个做法不是不行但它在单细胞层面上很吃亏。因为不同峰的GC含量、平均可及性、以及细胞本身的开放程度差异很大简单统计会把技术偏差当成生物学信号。chromVAR的核心思路正好用在刀刃上它不回答“这个motif在这个簇里富不富集”而是回答“在扣除峰的可及性背景后这个细胞里带这个motif的峰是否比预期更开放”。这个更开放的信号才是转录因子活性更合理的代理指标。1.2 偏差分数和z-score先认识这两个数字chromVAR输出两个你一定会碰到的量deviations和z。deviations是原始偏差分数直观理解是“观察到的信号”和“预期信号”的比值取对数。每个TF motif会对应一组“匹配峰”每个细胞在这组峰里的真实信号是observed然后chromVAR会用一组背景峰做预测得到expected。如果某细胞里GATA1 motif相关峰的实际可及性明显高于这类峰应有的平均水平deviations就是正的反过来就是负的。z是在deviations基础上做的标准化相当于把偏差分数除以标准差便于跨motif比较。ArchR里跑完chromVAR后MotifDeviations矩阵里同时存了这两个assay我实际使用中几乎都拿z做后续的可视化和比较。用个生活例子班里批作文不能直接拿字数比谁写得好因为有人字体大、有人字间距大。chromVAR做的就是先按字体大小和排版分行也就是GC含量和平均可及性把同水平的作文分到一组再比谁实际写的内容更超出预期。这个“超出预期”的幅度就是偏差分数。1.3 与简单motif富集分析的区别为了帮大家避开分析思路上的混淆我把普通motif富集和chromVAR偏差分析放在一张表里对比维度普通motif富集chromVAR偏差分析输入差异峰/某簇的峰集合峰×细胞计数矩阵 motif注释核心问题某组峰里是否富集某motif某细胞/某群中某TF活性是否高于预期输出p值、富集倍数每个细胞每个motif的连续偏差/标准化z分数分析单位通常是亚群单个细胞可聚合到群受到GC/开放度影响需要额外校正内置背景峰校正适合场景组间差异峰的motif构成单细胞水平的TF活性推断这两者不是替代关系而是视角不同。做HOMER式motif富集时你已经在问“这群特异性开放的峰里面有没有某个TF的偏好位置”而chromVAR是在问“这群细胞里某个TF的调控信号是否真的比背景活跃”。很多项目需要两种一起上先用chromVAR挑候选TF再用传统富集佐证生物学含义。2. ArchR里跑chromVAR的完整流程2.1 前置检查PeakMatrix不齐后面全是空谈在跑addMotifAnnotations和computeDeviations之前先确认项目里已经有PeakMatrix。ArchR的chromVAR流程是基于峰集和细胞计数矩阵的如果你只做了addReproduciblePeakSet但忘了addPeakMatrix后面一定会报“PeakMatrix not found”之类的错。我习惯先跑一句getAvailableMatrices(proj)正常应该能看到PeakMatrix。如果不在补齐proj - addPeakMatrix(proj, force TRUE)这步会按细胞重新量化峰上的片段数是后续一切的基础。另外提醒一句之前做过的过滤很关键doublet、低质量细胞必须提前清理干净。chromVAR的背景峰模型是全局统计混入大量低质量细胞会把背景拉偏最后所有TF的z分数都失真。2.2 addMotifAnnotations把PWM打到全基因组峰上前置条件就绪后加motif注释proj - addMotifAnnotations( ArchRProj proj, motifSet JASPAR2020, name Motif, force TRUE )motifSet参数默认常用JASPAR2020也有JASPAR2018、cisBP等选择。ArchR内部会用TFBSTools做PWM扫描匹配到每个峰后生成一张二进制MotifMatrix记录“该峰是否包含该motif”。这一步产生的行名一般是Motif_1这样的编号TF名字放在元数据里需要时可以通过rowData查回来。首次运行JASPAR相关注释时需要联网下载对应的数据库。如果网络不稳或处在离线环境建议提前手动安装对应包或者准备好自己的PFM/PWM列表。自定义motif集合并不少见特别是做非模式物种或特定TF家族时可以用library(TFBSTools) pfmL - readRDS(my_pfms.rds) # PFMatrixList 或 PWMatrixList proj - addMotifAnnotations(proj, motifSet pfmL, name MyMotif)自定义列表时一个容易踩的坑是ID重复。如果两个PFM的ID相同ArchR构建MotifMatrix时会产生冲突后续rowData(devMat)$name会乱掉。我建议在导入前先用uniquify或name()自查一遍保证每条记录的唯一性。跑完注释后可以用getFeatures(proj, useMatrix MotifMatrix)快速看一眼是否生成成功。提示force TRUE是为了强制重跑注释。如果你换过PeakSet或者之前注释时用的数据库版本不对一定要加这个参数否则ArchR可能直接返回旧结果。2.3 背景峰的选择为什么每个motif都需要一组“对照组”chromVAR最核心的设计是背景峰校正。光有MotifMatrix还不够你还需要算出每个motif匹配峰的“对照组”也就是背景峰。bgdPeaks - getBackgroundPeaks(ArchRProj proj, n 200)getBackgroundPeaks会为每个motif的峰集合挑选一批在GC含量、平均可及性上相似的峰作为背景。n是重采样次数ArchR默认通常取200表示每个背景集合被重复构建200次。n越大背景估计越稳但耗时和内存也会同步上升。数据量特别大的时候有人会把n降到100我一般是在正式分析时保留200做快速探索时临时降到100先跑通。新版ArchR里也有addBgdPeaks这种写法本质是一样的但API在不同版本间有差异。如果你的代码报“Background peaks not found”或“bgdPeaksmust be non-empty”先查一下当前ArchR版本的函数签名。别在一棵树上吊死换个正确函数通常马上好。一个容易被忽略的细节背景峰是基于当前PeakSet和细胞集生成的。如果你后来改了PeakSet、过滤了更多细胞必须重新生成背景峰不能复用旧对象。否则对照背景与真实峰集不匹配z分数会被系统性污染。2.4 computeDeviations计算细节与资源规划motif注释和背景峰都有了正式计算偏差proj - computeDeviations( ArchRProj proj, backgroundPeaks bgdPeaks, force TRUE ) devMat - getMatrixFromProject(proj, useMatrix MotifDeviations)computeDeviations封装了chromVAR的核心逻辑ArchR会先校验输入矩阵的类型确保是整数型计数然后对每个细胞每个motif计算观察值、预期值和偏差分数。最终devMat是一个SummarizedExperiment对象含deviations和z两个assay行对应motif条目列对应细胞。这里必须提醒computeDeviations是整套流程里最吃CPU的一步。我跑过一个四万多细胞的PBMC项目motif数量在800左右单机八线程大概跑了四十分钟到一小时。如果细胞数超过十万建议把它放到后台Rscript里跑前端不要挂着IDE防止内存和会话状态互相干扰。nohup Rscript run_chromVAR.R run_log.txt 21 跑之前还可以用addArchRThreads(threads 8)把线程数拉起来但线程不是越高越好。每开一个线程就会多一份内存占用如果机器只有32GB内存建议线程数控制在4到6个否则可能跑到一半因为内存不足被系统杀掉。3. 拿到偏差分数后可视化与群体差异分析3.1 plotDeviations先看全貌再谈细节computeDeviations跑完后第一批图通常用plotDeviations出。它可以按细胞或按簇展示一个或多个motif的偏差变化p - plotDeviations( proj, name c(GATA1, CEBPA, RUNX1), pal paletteContinuous(solarExtra) )不过有一个小坑ArchR的motif名在MotifMatrix里通常带后缀编号比如GATA1_1、GATA1_2。同一TF可能对应多个PWM条目所以直接传GATA1有时候会匹配不上。稳妥做法是先从rowData(devMat)$name里查出对应的完整条目名再传给plotDeviations。plotDeviations返回的是一张热图每行是一个motif条目每列是一个细胞或一组细胞。如果按簇分组颜色越偏红说明该簇里这个TF的偏差越高。我一般不会只画一两个TF而是先画一批候选TF看整体图谱确认哪些在簇间有明显分隔再去细看信号。3.2 plotVarDev高变异TF才是真正值得追的候选全局热图适合看已知候选筛新TF还是得靠变异性排序p - plotVarDev(proj, plotName VarDev)plotVarDev输出一个点图横轴是每个motif纵轴是变异程度。变异性越高说明这个TF在细胞间活性差异越大越可能是驱动细胞状态分化的因素。实际操作中从最右侧那一撮点里挑Top 20或Top 50基本上能把项目里值得关注的TF圈出来。这个方法很适合用来讲故事你不需要对800个motif做一遍差异检验只用plotVarDev把注意力锁在少部分高变异条目上再进入下游验证。我个人觉得这比一上来就做全矩阵聚类更高效因为全矩阵聚类会被大量表达模式相似的“管家motif”干扰。3.3 把TF活性画到UMAP上偏差分数是连续值最好的展示方式之一就是迭到UMAP上。前提是项目里已经做了UMAP并且需要先计算插补权重proj - addImputeWeights(proj) p - plotEmbedding( proj, colorBy MotifMatrix, name GATA1_1:score, imputeWeights getImputeWeights(proj) )addImputeWeights会基于KNN对全矩阵做插补跑起来也挺耗时但只需跑一次之后画任意TF都可以复用。colorBy MotifMatrix表示从MotifMatrix里取特征后面的GATA1_1:score表示取该motif的插补后偏差分数。注意这里不是取MotifDeviations而是MotifMatrix中扩展出来的score特征ArchR内部会把偏差分数包装成可用特征。如果你发现在plotEmbedding里传name始终无效可以先运行getFeatures(proj, useMatrix MotifMatrix)把可用的特征名打印出来再复制完整的条目名进去省得再猜。UMAP迭加图的好处是信息量很大一张图上能看到某个TF的活性集中在哪个细胞亚群、是否有连续梯度、是否存在少见的高活性细胞团。这些信息在热图里容易丢失。3.4 组间比较一行代码背后的统计推断可视化只能直观判断真正要说“群体A里TF X活性显著高于群体B”还是得做统计比较。ArchR本身的差异分析功能更多面向peak和基因表达对偏差分数我一般直接把devMat取出来自己算devMat - getMatrixFromProject(proj, useMatrix MotifDeviations) z - assays(devMat)$z groupA - which(proj$Clusters %in% c(C1, C2)) groupB - which(proj$Clusters %in% c(C5)) # 以某一行motif为例做wilcoxon检验 motifIdx - which(rowData(devMat)$name GATA1) p - wilcox.test(z[motifIdx, groupA], z[motifIdx, groupB])$p.value这样单个TF遍历跑一遍不难但涉及几十个TF时注意多重检验校正。我通常会用p.adjust(ps, method BH)做一个简单的FDR控制。另一个更稳妥的做法是先用伪bulk的方式把同簇的细胞平均或求和得到伪bulk偏差分数再做组间比较避免把单细胞的不独立性问题带进检验。这个思路未必写进ArchR教程但实际在做课题汇报时会经得起审稿人追问。4. 高频报错、解读误区和调优经验4.1 高频报错速查表很多报错不是基因问题而是流程顺序和参数版本问题。我把自己反复遭遇过的几类整理成表便于快速定位报错/现象常见原因解决办法PeakMatrix not found没做addPeakMatrix或矩阵被覆盖补齐addPeakMatrix(proj, forceTRUE)Error in ... JASPAR2020数据库包未安装或网络失败手动安装对应数据库包或改为本地PFM文件seqlevels不匹配PeakSet的染色体命名和BSgenome不一致用seqlevelsStyle统一命名或指定正确的BSgenomeBackground peaks not found版本API差异或背景峰变空检查getBackgroundPeaks返回的对象重新生成computeDeviations长时间无输出线程/内存设置不合理降低线程数转后台Rscript运行z全为NA或大量NA某些簇细胞太少、峰覆盖太稀疏过滤过小簇或尝试提高背景峰稳健性最让我头疼的是染色体命名不一致如果PeakSet是chr1而BSgenome是1addMotifAnnotations里做序列提取时会失败或静默丢行。一定要在所有步骤前统一风格尤其做非人类样本时别忘了对应物种的BSgenome。4.2 解读chromVAR结果时必须记住的三条红线第一z分数高不代表TF一定结合了DNA。它只说明“带这个motif的开放峰活性高于预期”中间还有可及性到结合的距离。想坐实结合需要ATAC footprints、ChIP-seq或突变验证。写文章或作报告时建议使用“推断活性”而不是“结合活性”这类表述。第二同一TF的多个PWM条目不要随便合并。GATA1_1和GATA1_2可能是不同长度或不同来源的PWM匹配到的峰集合不一样直接合并会把噪声带进来。如果你确实想合并得按“共享峰的信号汇总”来做而不是简单把两个z相加。第三单细胞的z值非常吵。不要根据单个细胞的偏差分数下结论也不要拿一个细胞的z值把细胞硬分成“有活性和无活性”两组。正确姿势是先看群水平或UMAP插补结果用群体统计量给结论背书。4.3 实测调优从半天到两小时如果你跑大样本总觉得chromVAR慢我的建议按优先级排列。第一先过滤再计算。把所有可疑细胞、低质量和极小簇先清掉背景峰会干净很多计算量也会明显下降。第二getBackgroundPeaks中的n可以探索式降到100但正式结果最好保留一份200的版本做稳健性验证。第三把线程数和内存预算先算好我倾向于一个线程至少给3GB内存开8个线程就需要24GB以上可用内存不然只是看起来在用多线程。第四中间结果一定要保存saveRDS(devMat, devMat.rds) saveArchRProject(proj)devMat本身是一个SummarizedExperiment后续所有可视化、统计比较都能基于它重来不需要每次都重新跑一遍computeDeviations。保存好这个对象你后面就不必反复加载整个ArchRProject能省不少内存。5. 一些个人体会把这条系列写到第十二篇chromVAR是我觉得单细胞ATAC-seq分析里“投入产出比”很高的一个模块。它不需要额外造数据也不需要复杂的模型调试只要流程顺序正确出来的结果在生物学上往往能直接指向一批转录因子候选给下游实验设计提供线索。我现在的标准流程是先plotVarDev扫Top 20高变异TF再用plotDeviations看这些TF在簇间的整体分布挑出感兴趣的几个上UMAP验证空间模式最后用z矩阵做组间统计比较。这个顺序能尽量避免“看到什么说什么”的偶然性也能让p值检验建立在可视化预筛之上而不是反过来硬凑结果。还有一个很实用的习惯如果项目里有多个样本或多个条件的ATAC数据不要一上来就在全局所有细胞上跑一套chromVAR先按大群或者条件分层跑比较同一细胞类型在不同条件下的TF活性差异往往比全局偏差分数更容易讲清楚生物学问题。毕竟chromVAR算的是相对活性背景变了z值的含义也会跟着变。