先交代一下背景。我做基因组非编码元件相关分析也有一段时间了日常打交道最多的就是ATAC-seq peak、增强子、CTCF结合位点这一类不落在基因外显子上的区域。最近在做一个增强子集合的功能注释时又被问到“非编码元件的GO富集到底用rGREAT还是chipseeker”这个问题其实挺有代表性的几乎每个做非编码调控元件注释的人都会卡在这一步。简单说我的结论是如果你是拿着一堆peak、增强子、绝缘子这类非编码元件去做功能富集优先用rGREAT如果你手里是落在启动子、外显子、内含子这类基因结构上的位点chipseeker会更顺手。但要说清楚为什么就得把两个工具对“元件-基因关联”这件事的处理逻辑讲明白。1. 为什么非编码元件的GO富集这么容易翻车1.1 问题本质元件离基因太远不知道算谁的做GO富集的第一步是把基因组上的区域映射到基因上。对编码区的变异或差异表达基因来说这一步毫无争议基因就在那里转录本结构清清楚楚。但非编码元件不一样增强子可能离最近的基因启动子有几十万甚至上百万碱基CTCF结合位点经常落在基因间区的“荒郊野岭”这时候就会面临一个尴尬问题这个元件到底算哪个基因的把元件注释给距离最近的基因还是注释给被它调控的基因如果是后者怎么知道它调控谁我用一个生活化的类比解释一下chipseeker的逻辑像“查户口”它看你这个位置是不是在某个人家里如果在院子里、客厅里、卧室里就记在这个人名下如果站在两栋楼之间的空地上它就标记成“无主地”。rGREAT的逻辑像“划学区”它不管你是不是住在房子里只要你在某所学校的招生范围内就把你算作这个学校的学生。非编码元件恰恰就是那批“站在空地上的人”所以rGREAT这种学区制天然就更适合。1.2 富集背景的选择不要忽略这一步除了元件到基因的映射富集分析的另一个关键点是背景基因集怎么选。GO富集本质上是在做超几何检验或费希尔精确检验看在你的目标元件关联的基因里某个功能的基因是否显著多于随机预期。这个“随机预期”就取决于背景基因集。如果背景选的是全部编码基因结果往往会被那些基因数目特别多的功能类别主导比如各种代谢过程、基础细胞过程这些结果虽然p值好看但跟你的非编码元件可能根本没有生物学关系。更合理的做法是在rGREAT中把背景限制在与你的元件有类似可及性特征或类似基因组分布的基因集合上或者在chipseeker中自定义背景基因集比如用表达量匹配的基因或相同染色体分布的基因。这一步容易被忽略但对结果的影响非常大我后面会详细讲参数设置。2. chipseeker和rGREAT的注释逻辑与适用边界2.1 chipseeker的核心思路区域落在哪里就注释到哪个基因chipseeker是Bioconductor生态里处理peak注释的老牌工具它的核心函数annotatePeak会先把你的peak比对到基因模型的各个功能部件上包括promoter、5‘ UTR、3’ UTR、exon、intron、下游区域和基因间区。它会给出一个非常细的注释结果比如“这个peak落在某某基因的第二个内含子里”或者“落在某某基因启动子下游2kb处”。这个逻辑对落在基因结构内部或其附近的peak非常友好注释结果可以直接用于后续的GO分析因为每个peak都明确关联到了一个具体的基因。但问题是对非编码元件来说大量peak都落在Intergenic区域chipseeker里虽然也有genicAnnotation参数可以调整但默认情况下intergenic的peak会直接被归入一个单独的类别不会关联到任何基因然后这些peak就被排除在后续富集之外。有人可能会说那我调一下参数把intergenic的peak距离最近的基因关联上不就行了chipseeker确实提供了一些选项但它的底层逻辑还是基于基因结构注释的线性距离不会考虑染色质构象、增强子-启动子互作这类三维基因组信息本质上还是在“查户口”的框架里打转。这意味着即使你强行把intergenic peak注释给最近基因结果也可能是有偏的因为最近的基因不一定是它真正调控的基因。2.2 rGREAT的核心思路用一个“调控域”把元件圈给基因rGREAT是Bioconductor上对原版GREAT算法的重新实现它处理非编码元件的方式和chipseeker完全不同。rGREAT会为每个基因定义一个调控域这个调控域分两部分一部分是basal regulatory domain从基因的TSS出发上游和下游各延伸一段距离默认是上下游各5kb这个区域是基因的基础调控区另一部分是extension domain在basal domain的基础上继续向上下游延伸默认向上游延伸1000kb、向下游延伸1000kb但会止步于最近的其他基因的basal domain边界。这个设计的妙处在于rGREAT把一个非编码元件关联到基因的方式从“点是否落在基因结构内”变成了“点是否落在基因的调控领地内”。一个增强子哪怕离基因50kb远只要它落在该基因的extension domain范围内它就会被算作这个基因的关联元件。而且一个元件可以被关联到多个基因因为它的调控域可能与多个基因的调控域重叠。这种多重关联关系更符合增强子可以同时调控多个基因这一生物学事实。我在实际使用中的体会是rGREAT注释非编码元件时确实很少出现“无家可归”的情况每个元件几乎都能关联到至少一个基因而且通常能关联到多个基因。但这里也藏着一个隐患如果调控域范围设得太大关联到的基因会非常多做GO富集时背景噪声会变大显著富集的结果可能会被稀释。2.3 怎么选看你的数据落在哪里现在可以把选型逻辑说清楚了。如果你的peak集合主要是启动子区域、外显子、内含子这类明确落在基因结构内的位点比如做H3K27ac的启动子区域分析或者做ChIP-seq且peak大多靠近TSS那么chipseeker够用而且注释分辨率更高可以直接看到peak落在哪个具体的基因结构部件上。但如果你的peak集合是增强子、开放染色质区域、CTCF位点、或者各种远端调控元件这些元件大概率有大半落在基因间区那么chipseeker会把大量peak归为intergenic然后束手无策rGREAT则可以通过调控域把这些远端元件跟基因关联起来得到更完整的富集结果。另外还有个实际体验上的差别。rGREAT的富集分析是直接在注释阶段完成的调用great()函数后就能直接得到GO和通路富集结果chipseeker则分两步走先annotatePeak再用注释后的基因去做clusterProfiler的富集。对只想快速看结果的人来说rGREAT明显更省事。3. rGREAT实操从peak集合到GO富集结果的全流程3.1 准备输入文件GRanges对象是标准格式rGREAT的使用流程非常简洁最关键的是把你的元件集合转成GRanges对象。我一般这样读入数据。假设有一个BED格式的增强子peak文件前几列是chr、start、endlibrary(GenomicRanges) library(rGREAT) peaks - read.table(enhancer_peaks.bed, sep \t, header FALSE, stringsAsFactors FALSE) colnames(peaks)[1:3] - c(chr, start, end) gr_peaks - GRanges(seqnames peaks$chr, ranges IRanges(start peaks$start, end peaks$end))这里要注意一个细节BED格式的坐标是0-based的而GRanges对象的坐标是1-based的。rGREAT官方推荐的做法是BED文件读取后start列要加1即IRanges(start peaks$start 1, end peaks$end)。如果坐标不偏移虽然大部分情况下结果偏差不大但遇到边界情况时peak可能会错位一个碱基导致注释结果轻微改变。我习惯在数据读入时就修正避免后面排查问题浪费时间。3.2 运行great()函数并设置合适的参数GRanges对象准备好后直接调用great()函数即可。这里有几个参数非常关键直接决定注释质量和富集结果。job - great(gr_peaks, gene_sets GO Molecular Function, background genome, mode basalPlusExtension, min_bases 1000, max_bases 1000000)第一个参数gene_sets对应的不是rGREAT自己的注释而是从MSigDB或GO数据库里取的功能基因集。可以选“GO Molecular Function”、“GO Biological Process”、“GO Cellular Component”也可以选“MSigDB H”、“MSigDB Oncogenic”等。我一般建议多用几个gene_sets做交叉验证比如同时跑GO Biological Process和GO Cellular Component因为非编码元件调控的功能往往不只体现在生物学过程上还可能影响细胞组分中蛋白复合物的形成。第二个重要参数是background选“genome”时背景基因是所有编码基因的所有基因选“data”时背景基因是所有出现在你的输入元件关联结果中的基因。我强烈建议除非你明确知道自己在做什么否则默认选“genome”就好。选“data”会把背景基因集限制得特别小富集检验的统计功效会大幅下降本来显著的GO term可能全部变成不显著。然后是mode参数默认是basalPlusExtension这个策略前面已经介绍过。如果你只想用basal domain注释可以设成basal但非编码元件尤其是增强子大部分不在TSS上下游5kb范围内用basal模式会让很多元件关联不到基因。反过来如果你把max_bases设得特别大比如500kb甚至1000kb以上关联的基因数量会激增但背景噪声也会变大。我自己的经验是对增强子类元件max_bases设置在500kb到1000kb之间比较合理对CTCF位点这种可能参与染色质结构组织的元件这个范围也可以接受。3.3 查看和导出富集结果运行完great()后可以用getEnrichmentTables()把富集表拿出来也可以直接用plotRegionGeneAssociationGraphs()做可视化。我一般会先把富集表完整导出再做后续筛选。tbls - getEnrichmentTables(job, download FALSE, min_region_gene_ratio 0.05, min_gene_set_size 10, max_gene_set_size 2000)这里有几个筛选参数值得说明。min_region_gene_ratio表示一个GO term关联到的输入元件基因数占输入元件总关联基因数的比例低于这个比例会被过滤。这个filter很重要它确保富集结果不是由少数几个元件异常关联到的基因撑起来的。min_gene_set_size和max_gene_set_size过滤掉那些基因数目过少或过多的GO term基因数目过少的term没有统计效力基因数目过多则可能太宽泛结果没有生物学聚焦度。取回富集表后通常会按p-value排序再看fold enrichment和region_gene_ratio。我习惯额外看一个指标expected和observed的比值。有时候某个GO term的p值很显著但observed的关联基因数其实很少这种结果往往是被个别极端元件影响的生物学意义有限。反过来那些observed关联基因数多、富集倍数适中、p值稳健的结果才是值得深入挖掘的。3.4 富集结果可视化条形图与基因组关联图rGREAT自带的一些图比如根据富集表画的条形图或者点图一般够用。但我更常用的是plotRegionGeneAssociationGraphs()它能展示你的元件跟关联基因之间的对应关系在汇报数据时需要用来支撑“这些增强子可能调控哪些基因”的结论。plotRegionGeneAssociationGraphs(job)这张图会把你输入的元件按基因组位置展开画出基因、TSS、元件之间的关联线条。如果元件集中在某些基因附近线条会非常密集说明这些区域可能是一个调控热点。我在做增强子功能注释时经常先跑一遍这个图看看整体注释质量如何。如果关联线条乱成一团且元件关联到上千个基因那大概率是max_bases设得太大需要调小后重新跑。3.5 实操中的候选基因导出除了富集表rGREAT还能拿到每个元件关联了哪些基因这一步在后续做motif分析或构建调控网络时非常有用。assoc_genes - getRegionGeneAssociations(job)这个函数返回一个GRangesList每个元件对应一组关联基因。我通常会把这个结果转成data.frame保存下来然后根据元件ID去查找它关联的基因名。在做增强子-基因关联时我一般还会把rGREAT的关联结果和Hi-C数据或eQTL数据做交叉验证。毕竟rGREAT的关联是基于线性距离的虽然比chipseeker合理很多但仍在二维基因组范畴。如果同一个增强子在Hi-C数据和rGREAT结果中都与同一个基因关联那我对这个调控关系的信心会大很多。4. chipseeker实操什么时候用、怎么配合clusterProfiler4.1 annotatePeak的用法和输出解读顺手把chipseeker的流程也交代清楚。使用场景很明确你的peak落在基因结构附近或内部时chipseeker能给出更精细的结构注释。代码很简洁library(ChIPseeker) library(TxDb.Hsapiens.UCSC.hg38.knownGene) library(clusterProfiler) library(org.Hs.eg.db) txdb - TxDb.Hsapiens.UCSC.hg38.knownGene peakAnno - annotatePeak(gr_peaks, TxDb txdb, annoDb org.Hs.eg.db, level gene, overlap original)这里有一个容易忽略的点level参数可以设成gene或transcript。默认按gene注释时多个转录本共享同一个基因一个peak如果落在某个基因的任意转录本上就算注释到这个基因。设成transcript时注释粒度更细但同一个peak可能会因为落到多个转录本上而产生重复注释后续统计时要先去重。我一般用gene级别分析非编码元件时gene级别的稳定性和可解释性都更好。4.2 用注释后的基因集做GO富集拿到注释结果后提取出被注释到的基因ID然后交给clusterProfiler跑富集就行了。gene_ids - unique(peakAnnoanno$geneId) ego - enrichGO(gene gene_ids, OrgDb org.Hs.eg.db, keyType ENTREZID, ont BP, pAdjustMethod BH, qvalueCutoff 0.05)这里必须强调直接用unique(peakAnnoanno$geneId)取基因ID是个偷懒但常见的方式问题在于如果某个基因同时被100个peak注释到它在列表里只出现一次而另一个基因只被1个peak注释到也只出现一次。在富集检验时两者的权重完全一样。这其实隐含了一个假设即每个被注释基因都有同等重要性但实际情况往往是被更多peak注释到的基因更值得关注。如果想要更精细的加权富集可以按peak计数对基因加权或者用ngs数据做rank-based的富集分析。不过大多数场景下这种简单的unique处理是可以接受的。chipseeker的另一个常见坑是如果gr_peaks里夹杂了大量intergenic区域的peak这些peak在annotatePeak结果里会落在annotation列为Intergenic的条目中对应的geneId列会是NA。如果不加过滤就提取geneId做富集NA会被直接忽略相当于这些intergenic peak被静默丢弃了。这会导致一个错觉你跑完富集发现结果还挺好但实际上你只分析了不到一半的peak。做这类分析时我强烈建议先统计各注释类别的分布看看Intergenic占比是多少。如果超过一半就该认真考虑rGREAT了。4.3 可视化peak注释分布图chipseeker自带的可视化函数相当好用plotAnnoPie()和plotAnnoBar()能非常直观地展示peak在基因组功能区域的分布比例。我通常在分析开始时先跑这两个图快速判断数据特性。plotAnnoPie(peakAnno)如果饼图里“Intergenic”占了大半而你的研究对象确实是非编码元件这就是一个明确信号说明chipseeker不适合这次分析。反过来如果你的peak主要落在Promoter、Exon、Intron里那chipseeker的注释结果就很有价值。5. 两种工具对比选型决策表与组合使用思路5.1 核心差异对照表我用一张表把两个工具的核心差异列清楚方便在团队协作或写方法部分时直接引用。对比维度rGREATchipseeker核心注释逻辑调控域basal extension基因结构部件promoter、exon、intron等非编码远端元件处理通过调控域关联到基因适合远端增强子大部分非编码远端元件被归为Intergenic无基因关联一个元件关联基因数可以关联多个基因通常关联一个基因基因背景设置genome或data需手动提取注释基因配合clusterProfiler富集功能直接输出GO/通路富集表需要额外走clusterProfiler流程适合的场景增强子、ATAC-seq、CTCF等非编码元件启动子、外显子、内含子等基因结构相关peak对结构注释的依赖低不需要精确的转录本结构高依赖TxDb的转录本结构注释这张表基本能回答“什么时候用哪个”的问题。需要强调的是这两个工具不是互斥关系在同一个项目里完全可以先用chipseeker做结构注释看看分布情况再用rGREAT做功能富集。5.2 组合使用先chipseeker看结构再rGREAT做富集我在真实的增强子注释项目里经常采用组合方案。第一步先用chipseeker的annotatePeak对全部peak做结构注释看看落在启动子、外显子、内含子和基因间区的比例。这一步的价值在于了解数据的基本构成比如如果只有20%的peak落在基因结构内那剩下80%的intergenic peak必须用rGREAT来覆盖。第二步对全部peak直接用rGREAT做GO富集获得完整的调控功能的图景。第三步如果某个GO term特别有兴趣再回到chipseeker的注释结果里看与该GO term相关的基因具体落在哪些结构区域进一步细化候选增强子-基因调控对。这个流程的好处是每个工具都用在了它擅长的位置而且最终结果可以互相验证。比如rGREAT富集到某个神经发育相关的GO term然后chipseeker的注释结果显示这些基因主要被intergenic区域的peak关联这个信息本身就很有生物学意义说明这些调控可能来自远端元件而非启动子附近的变化。6. 常见报错与高发问题排查实录6.1 rGREAT查询失败或返回空结果rGREAT的great()函数在运行时背后会请求GREAT的API服务网络不稳定或GREAT服务更新时经常会出现连不上或返回空结果的情况。这个问题特别常见我自己就遇到过好几次尤其是跑大批量任务时。排查思路分几步。第一步先确认网络连通性GREAT API是外部服务如果你在防火墙受限的环境里跑大概率连不上。第二步检查是不是参数设置导致的结果全被过滤了比如min_region_gene_ratio设得偏高把一些弱关联的GO term全过滤掉了。第三步尝试把输入元件集合拆分成小批量比如每次跑5000个peak分批次运行后再合并结果。GREAT原始网站对大输入有长度限制rGREAT在这方面虽然有所优化但太大的集合仍然会遇到超时或查询失败。如果GREAT服务长时间不可用还有一个备选方案用great()函数时设置use_gi TRUE它会使用Bioconductor自带的GO注释来做区域关联不再依赖外部API。我实测下来use_gi TRUE的结果与标准GREAT查询在主要富集通路上基本一致但在一些具体的GO term细节上会有差异可以作为紧急情况下的替代方案。6.2 chipseeker的TxDb物种版本不匹配chipseeker依赖TxDb对象提供基因结构注释而TxDb的版本必须跟你的参考基因组版本匹配。hg19的peak用了hg38的TxDb.hsapiens.UCSC.hg38.knownGenechromosome命名可能对不上比如chr1与1混用导致所有peak都注释不上。我遇到过一种更隐晦的情况参考基因组版本装对了但TxDb对象是从AnnotationHub下载的没有仔细核对assembly日期导致TxDb里的基因模型与用到的参考基因组之间有细微差异。解决办法很简单一方面是统一从Ensembl或UCSC下载参考序列和TxDb另一方面是在分析开始前检查seqlevels是否一致比如seqlevels(gr_peaks)和seqlevels(txdb)做一次比较。这类小检查能省掉后面大量排查时间。6.3 rGREAT关联基因过多导致富集结果失真这是一个方法论层面的坑。rGREAT把max_bases设成1Mb后一个增强子可能关联到几十个甚至上百个基因这时GO富集的统计检验会在一个很大的关联网络里进行理论上能检测到的显著term会减少因为背景关联基因多了同样的observed值显著性会下降。如果你发现rGREAT跑出来的富集结果几乎没有显著项但直觉上你的元件集合应该有明确功能倾向不妨回头看看getRegionGeneAssociations()返回的关联基因数分布。如果中位关联基因数超过50建议逐步缩小max_bases到500kb、200kb再重新跑。对绝大多数增强子研究来说200kb到500kb是一个比较平衡的范围。不过要注意这个调整也要结合你的生物学问题来定如果研究的是超级增强子这类跨越大范围调控的元件太小的范围反而会漏掉真实调控关系。6.4 GO富集结果里全是“看门基因”功能类别做非编码元件GO富集时最让人沮丧的结果就是一堆“细胞粘附”、“代谢过程”等大而化之的GO term排在最前面。这通常说明背景基因集有问题。在使用chipseekerclusterProfiler时如果背景基因是全部注释基因而你的元件关联到的基因偏向某些特定染色体区域则富集结果会被染色体上的基因密度特征主导。一个常用的修正策略是计算你的元件关联基因在每个染色体上的分布再模拟相同染色体分布、相同数量规模的随机基因集合反复多次后用模拟集合作为背景。目前Bioconductor上已有类似的包可以做这个事但rGREAT的“genome”背景其实已经考虑到了基因组结构的影响所以在这个问题上rGREAT会比chipseeker更让人省心。7. 进阶经验非编码元件功能注释的几条心得如果想在“rGREAT还是chipseeker”之外更进一步我还有几条实操心得供各位做非编码元件功能注释时参考。第一GO富集的输入元件应该是经过严格质量控制的非冗余peak集合而不是原始的全部peak。我见过不少分析直接用MACS2输出的原始peak做富集里面包含大量低置信度的宽峰或背景噪声peak这些噪声peak在rGREAT的调控域关联机制下会把很多无关基因拉入关联集稀释真正的信号。建议先做IDR分析保留高置信度peak或者至少用peak的q-value做个阈值筛选再拿去富集。第二如果有组蛋白修饰或染色质可及性数据富集前先按元件信号强度分层。比如把H3K27ac信号最强的增强子单独拿出来做一组GO富集把信号较弱的增强子放在另一组两组结果对比往往能看出功能差异。强弱增强子在功能富集上的差别可能比增强子与随机区域的差别更具生物学意义。第三rGREAT的关联结果只能说明“统计上的潜在调控关系”不代表真实的物理调控。增强子与启动子的互作还受三维基因组结构影响所以我在项目里一般会把rGREAT结果当作候选列表然后再用Hi-C、eQTL或CRISPR验证来缩小范围。这一点在写文章时尤其重要避免把预测结果当作实验结论来描述。我个人的习惯是启动子相关peak用chipseeker远端非编码元件用rGREAT遇到模块化的调控项目就两种都用、互相补充。分析工具的选择不是站队问题而是看数据形态和要回答的科学问题。如果你手头正好有一批增强子peak不知道怎么做功能注释把坐标整理成GRanges对象跑一遍rGREAT再对照本文提到的参数设置和坑点调整几轮基本就能拿到可用的GO富集结果了。后续想深入的话还可以把rGREAT关联到的基因和差异表达基因取交集看这些非编码元件是不是真的和特定通路或表型相关这一套组合拳做下来结论的可靠性会高很多。