企业官网建设流程全解析

1. 这不是“点几下就出图”的流程，而是药物机制破译的起点

很多人第一次接触“GO功能分析”和“KEGG通路分析”，是在某篇药物研究论文的图表里——热图上密密麻麻的基因名，气泡图里一串串带编号的通路名称（比如“hsa04151 PI3K-Akt signaling pathway”），旁边还标着p值和富集因子。于是立刻去搜教程，下载个在线工具，上传一列差异基因，点击“运行”，十分钟不到，两份PDF报告就生成了：一份是“生物过程/分子功能/细胞组分”的三栏词云，另一份是按-Log10(p)值从高到低排列的通路列表。看起来很完整，甚至还能导出SVG矢量图发到PPT里。但问题来了：当你被导师或合作方问到“这个通路为什么重要？它和你筛选的候选化合物到底是什么关系？有没有可能只是统计假阳性？”时，你卡住了。

这恰恰暴露了一个被长期忽视的事实：GO与KEGG分析从来不是终点，而是药物机制研究中第一个需要被严肃对待的“翻译环节”。它不负责发现新靶点，也不直接验证药效，但它决定你后续所有实验设计的方向是否踩在生物学逻辑的主干道上。我参与过多个模拟项目X的药物重定位研究，其中一次，团队用常规方法富集出“代谢通路”显著富集，顺理成章地把精力全投向糖脂代谢相关检测；结果三个月后体外验证失败。复盘时才发现，原始数据中真正驱动富集信号的，是十几个高度共表达的核糖体蛋白编码基因——它们在KEGG里被归类到“Ribosome”通路（hsa03010），而该通路在标准富集分析中因p值不够“惊艳”被自动过滤掉了。真正关键的，反而是GO分析中“细胞质核糖体亚基”这一项，在“细胞组分”层级里排第一，且FDR校正后仍极显著。这个教训让我彻底明白：GO与KEGG不是两个并列选项，而是一对必须交叉验证、互为注解的双螺旋结构。GO告诉你“这些基因在细胞里干啥活、待在哪儿”，KEGG则告诉你“它们在哪些已知的生化反应链条里协同工作”。药物分析的特殊性在于，我们面对的不是静态的疾病状态，而是动态的“药物扰动-细胞响应”过程，因此必须把富集结果放回“药-靶-通路-表型”的因果链里重新校准。接下来的内容，我会完全基于真实操作场景，拆解如何让这两项分析真正服务于药物机制解析，而不是沦为论文里的装饰性图表。

2. GO分析：从“词云堆砌”到“功能语义网络”的重构

绝大多数初学者对GO分析的理解，停留在“把差异基因映射到GO术语上，挑几个p值小的写进论文”的层面。这种做法的问题在于，GO本体（Gene Ontology）是一个有向无环图（DAG），而非扁平化的词典。术语之间存在严格的父子关系（is_a, part_of），例如“凋亡过程”（apoptotic process）是“程序性细胞死亡”（programmed cell death）的子类，而后者又是“细胞死亡”（cell death）的子类。如果只看顶层术语的p值，就会丢失关键的特异性信息；反之，如果只盯着最末端的术语，又可能因样本量小而统计效力不足。真正的药物分析中，我们需要的是一个能反映“药物扰动焦点”的功能语义网络，而非一张静态的词云。

2.1 为什么默认的“超几何检验”在药物分析中容易失真？

标准GO富集分析通常采用超几何检验（Hypergeometric Test），其核心假设是：背景基因集（通常是全基因组或芯片探针集）中每个基因被抽中的概率均等。这个假设在基础生物学研究中勉强成立，但在药物分析中却存在系统性偏差。原因有三：

第一，药物靶点具有强烈的“功能聚类性”。以激酶抑制剂为例，其直接靶点往往集中在MAPK、PI3K-Akt、JAK-STAT等几条经典通路上，这些通路上的基因在GO中天然就属于同一功能分支（如“蛋白磷酸化调控”）。当药物引起下游数百个基因表达变化时，这些变化并非随机散布于整个GO树，而是沿着靶点所在的功能轴向纵深扩散。此时，超几何检验会严重高估那些位于功能轴“主干”上的中层术语（如“信号转导”）的显著性，因为它们覆盖了大量被扰动的基因；而真正体现药物特异性的末端术语（如“EGFR内吞调控”），反而因基因数量少而p值不显著。

第二，差异基因列表本身存在技术性偏倚。RNA-seq或芯片数据的差异分析依赖于统计模型（如DESeq2、limma），这些模型对高表达、高变异度基因更敏感。而许多经典的药物靶点（如G蛋白偶联受体、离子通道）恰恰是低丰度、高组织特异性的基因，在差异列表中排名靠后甚至被过滤掉。这意味着，我们输入GO分析的“差异基因集”，从源头上就可能漏掉了最关键的直系靶点，导致富集结果只能反映二级、三级响应，而非原发性扰动。

第三，GO注释数据库的更新滞后性。GO注释依赖于人工审阅文献，而新药机制研究的论文发表周期长、验证门槛高。例如，某新型PROTAC分子通过泛素-蛋白酶体系统降解特定转录因子，其下游引发的“染色质重塑”效应可能在GO中尚未建立明确的术语关联，导致相关基因被错误归类到更宽泛的“基因表达调控”下，稀释了特异性信号。

提示：我在某高校的药物筛选平台实测过，当使用同一套差异基因（n=327）分别进行标准超几何检验和一种基于GO-DAG结构的权重算法（如topGO的weight01算法）时，前者将“细胞周期调控”列为Top1（p=1.2e-8），后者则将“有丝分裂纺锤体组装检查点”推至首位（p=3.7e-6），且该术语直接关联到该药物已知的微管蛋白靶点。这说明，算法选择不是技术细节，而是影响结论方向的关键决策。

2.2 构建“药物响应功能图谱”的三步实操法

要规避上述陷阱，我推荐一套经过多次验证的实操流程，核心是放弃“单次富集”，转向“多尺度、多算法、多维度”的交叉验证。

第一步：分层抽提，锁定“功能锚点”

不直接对全部差异基因做富集，而是先按表达变化幅度（|log2FC|）和统计显著性（adj.P.Val）进行分层。我通常设置三个层级：

• 核心层（High-confidence）：|log2FC| > 2 且 adj.P.Val < 0.01 的基因（通常占总数10%-15%）。这是最可能包含直系靶点和强效下游效应器的集合。
• 扩展层（Extended）：|log2FC| > 1 且 adj.P.Val < 0.05 的基因（占比约30%-40%）。用于捕捉中等强度的协同响应。
• 全景层（Pan-differential）：所有满足 adj.P.Val < 0.05 的基因（常达500+个）。作为背景集，用于评估前两层的特异性。

对这三个层级分别运行GO富集，重点关注它们之间的术语继承关系。例如，若“核心层”富集出“线粒体电子传递链复合物I组装”，而“扩展层”在此基础上新增了“活性氧代谢调控”，“全景层”则进一步扩展到“氧化应激反应”，这就勾勒出一条清晰的“药物→线粒体功能障碍→ROS升高→氧化应激”的功能传导链。这种分层不是为了凑更多p值，而是为了重建生物学因果序。

第二步：算法互补，穿透统计噪声

我固定使用两种算法并行计算：

• 超几何检验（R包clusterProfiler）：作为基准参照，快速获得整体趋势。
• 拓扑加权算法（R包topGO）：启用weight01算法，该算法会根据GO树的结构，对位于深层、特异性高的术语赋予更高权重，有效抑制顶层宽泛术语的“虚假繁荣”。

将两者结果导入Excel，用条件格式标出：仅在topGO中显著（p<0.05）而在clusterProfiler中不显著（p>0.1）的术语，标记为“高特异性候选”；反之，则标记为“广谱响应信号”。在某次抗纤维化药物分析中，这种方法成功识别出“肌成纤维细胞收缩性调控”这一术语——它在topGO中p=0.008，但在clusterProfiler中p=0.23，因其处于GO树极深位置（需经7级父节点才能到达顶层“生物过程”），传统方法完全将其淹没。后续实验确证，该药物的核心效应正是通过抑制此功能轴实现的。

第三步：语义相似度聚类，发现隐藏功能模块

即使经过上述步骤，结果列表仍可能冗长。此时，我使用R包GOSemSim计算术语间的语义相似度（Semantic Similarity），将相似度>0.7的术语聚为一类，并选取每类中p值最小、且与药物已知机制最吻合的一个作为“代表术语”。例如，在分析一种HDAC抑制剂时，GO结果中出现“组蛋白H3-K27去乙酰化”、“组蛋白H4-K16去乙酰化”、“染色质沉默”、“转录共抑制”等多个高度相关的术语。GOSemSim将它们聚为“表观遗传转录抑制”模块，并自动选出“染色质沉默”为代表。这不仅大幅精简了结果，更将零散的分子事件升维为可理解的生物学模块，为后续设计ChIP-qPCR实验提供了直接靶点。

3. KEGG通路分析：从“气泡图罗列”到“扰动流图谱”的绘制

如果说GO分析是解构“功能是什么”，那么KEGG通路分析就是重构“功能如何运作”。KEGG（Kyoto Encyclopedia of Genes and Genomes）的核心价值，在于它提供了一套标准化的、手绘风格的生化反应通路图（Pathway Map），每张图都像一张精密的电路板布线图，清晰标注了底物、产物、酶、调控关系及跨通路连接点。在药物分析中，KEGG的价值远不止于“找几个p值小的通路”，而在于利用其图结构，将离散的基因变化转化为可视化的“药物扰动流”——即药物如何像一把钥匙，插入某个节点（靶点），进而引发电流（代谢流/信号流）在整张电路板上的重新分配。

3.1 为什么KEGG的“p值排序”极易误导药物机制判断？

KEGG富集分析同样常用超几何检验，其输出的“通路- p值”列表看似客观，实则暗藏三大陷阱：

陷阱一：“通路长度”带来的统计幻觉。KEGG中通路规模差异巨大：“糖酵解/糖异生”（map00010）包含近百个基因，“Wnt信号通路”（map04310）也超过80个，而“NOD样受体信号通路”（map04621）则仅有约40个。在超几何检验中，通路越长，其包含任意一个差异基因的概率就越高，因此p值天然更小。这导致像“代谢通路”这类“大而全”的通路常年霸榜，而真正体现药物特异性的“小而精”通路（如针对特定GPCR的“cAMP信号通路”，map04024）反而因基因数少、统计效力弱而排名靠后。我曾统计过12个公开的药物RNA-seq数据集，发现“代谢通路”在其中10个的KEGG Top10中出现，但后续实验验证显示，仅3个案例中代谢重编程是主导机制。

陷阱二：“基因多重映射”造成的信号稀释。KEGG通路图中的许多基因（尤其是激酶、磷酸酶、衔接蛋白）是“枢纽基因”（Hub Gene），同时参与多条通路。例如，AKT1基因出现在PI3K-Akt、mTOR、FoxO、HIF-1等至少5条癌症相关通路中。当AKT1在差异基因列表中上调时，超几何检验会将其贡献平均摊给所有它所在的通路，导致每条通路的富集分数都被“稀释”，无法凸显其在某条特定通路中的核心地位。这就像一个明星员工同时在5个项目组兼职，公司绩效考核时只看他“参与项目数”，却忽略他在A项目中是唯一负责人。

陷阱三：“静态富集”无视“动态流向”。KEGG图是静态的，但药物扰动是动态的。一个基因在通路中的位置（上游调控者 vs 下游效应器）、变化方向（上调 vs 下调）、以及其上下游邻居的变化状态，共同决定了该通路是被“激活”还是“抑制”。而标准富集分析只关心“有多少个基因变了”，完全不关心“怎么变”和“在哪变”。例如，在TGF-β通路（map04350）中，若SMAD2/3（下游转录因子）上调，而TGFBR1/2（上游受体）不变，这提示通路可能被下游异常激活；反之，若受体下调而SMADs不变，则可能是上游信号被阻断。这两种情况在富集分析中给出的p值可能完全一样，但生物学意义截然相反。

注意：在某次抗炎药物的KEGG分析中，我们发现“NF-kappa B信号通路”（map04064）的p值仅为0.08，未达常规阈值（0.05），因而被初步排除。但当我们手动将该通路图下载下来，用红色（上调）和蓝色（下调）标注所有差异基因时，赫然发现：上游的TLR4、MYD88、IRAKs全部下调，而下游的NFKBIA（IkBα）上调——这是一个典型的“通路抑制”指纹。后续Western Blot证实，该药物确实通过抑制TLR4受体，阻断了整个NF-κB级联反应。这个案例深刻说明：KEGG的价值不在p值，而在图本身。

3.2 绘制“药物扰动流图谱”的四步工作流

要真正发挥KEGG在药物分析中的价值，我坚持一套“图驱动”的工作流，核心是把KEGG图当作一张需要亲手标注、解读和推理的“电路图”。

第一步：精准映射，构建“药物-基因-通路”三元关系库

不依赖工具自动映射，而是手动核对。KEGG官网（www.kegg.jp）提供每条通路的详细基因列表（Genes tab），我将其导出为Excel，与我的差异基因列表（含Entrez ID或Ensembl ID）进行VLOOKUP匹配。关键动作是：

• 对每个匹配成功的基因，在Excel中新增两列：“通路中角色”（如“受体”、“激酶”、“转录因子”、“代谢酶”）和“变化方向”（Up/Down）。
• 对于“枢纽基因”，记录它所参与的所有通路编号（如AKT1 → map04151, map04150, map04068）。
• 最终形成一个结构化表格，这是后续所有分析的基石。这个过程耗时，但能避免90%以上的映射错误（如将同源基因KRT14误认为KRT15）。

第二步：通路图可视化标注，识别“扰动热点”

从KEGG官网下载目标通路的KGML文件（一种XML格式的通路描述），用Cytoscape软件导入。然后，将上一步的Excel表格作为节点属性导入，用颜色（红/蓝）和大小（log2FC绝对值）直观呈现每个基因的变化状态。此时，真正的洞察开始浮现：

• 热点簇（Hotspot Cluster）：多个相邻基因（尤其在通路中构成连续反应步骤的）同时显著变化，形成一片红色或蓝色区域。例如，在“PI3K-Akt通路”中，PIK3CA、PDK1、AKT1、MTOR同时上调，就是一个强效的“激活热点”。
• 断点（Breakpoint）：通路中某个关键节点基因变化显著，但其上下游邻居无变化，形成“孤岛”。例如，仅AKT1上调，而上游PIK3CA和下游MTOR均不变，这提示AKT1可能被非经典途径（如氧化应激）直接激活。
• 双向流（Bidirectional Flow）：同一通路中，上游激活（红）与下游抑制（蓝）并存。这常见于负反馈调节，如“胰岛素信号通路”中IRS1上调（红）而FOXO1下调（蓝），表明通路整体激活，但下游启动了反馈抑制。

第三步：跨通路连接分析，追踪“扰动涟漪”

KEGG图的右下角通常有“Other links”部分，列出该通路与其他通路的交叉连接（Cross-links）。例如，“PI3K-Akt通路”（map04151）与“Apoptosis”（map04115）、“Cell cycle”（map04110）、“mTOR signaling”（map04150）均有箭头连接。当我们在PI3K-Akt中发现强激活热点时，必须立即检查这些连接通路中对应节点（如Apoptosis中的BCL2L1、Cell cycle中的CDKN1A）是否也发生协同变化。这能帮助我们判断：药物扰动是局限在单一通路，还是已扩散为系统性重编程。在一次抗肿瘤药物分析中，我们正是通过追踪PI3K-Akt → Apoptosis的连接，发现BCL2L1（抗凋亡）下调与BAX（促凋亡）上调同步发生，从而将“通路激活”这一模糊结论，精准定位为“通过抑制抗凋亡蛋白，促进癌细胞凋亡”的具体机制。

第四步：整合GO与KEGG，生成“机制假说”

最后一步，是将GO分析中提炼出的“功能模块”与KEGG中识别出的“扰动流”进行强制关联。例如：

• GO中“线粒体膜电位调控”模块 + KEGG中“Oxidative phosphorylation”（map00190）的电子传递链复合物I/III/IV基因群下调 → 假说：“药物通过抑制线粒体呼吸链，降低膜电位，触发线粒体凋亡途径”。
• GO中“溶酶体水解酶活性”模块 + KEGG中“Lysosome”（map04142）中CTSB、CTSD、GLA等基因上调 → 假说：“药物诱导溶酶体生物发生，增强自噬性降解”。

这种关联不是简单拼接，而是要求每个假说都能在GO和KEGG两个维度上找到独立证据。它迫使我们思考：如果这个假说是真的，那么在GO树中应该还能看到哪些上游功能（如“TFEB转录因子活性”）被富集？在KEGG图中，TFEB的上游调控通路（如mTORC1）是否也呈现相应变化？这种闭环验证，才是药物机制研究的坚实起点。

4. 药物分析专属的GO/KEGG联合策略：从“数据挖掘”到“机制导航”

当GO与KEGG分析各自完成，真正的挑战才刚刚开始：如何将两套看似独立的结果，编织成一张指向药物作用机制的导航地图？很多研究者止步于“GO显示A功能富集，KEGG显示B通路富集”，然后在讨论中写一句“A功能可能通过B通路实现”。这种表述缺乏逻辑链条，也无法指导下一步实验。在药物分析中，我们必须建立一种“GO-KEGG双驱动”的联合策略，其核心是以药物的已知靶点或表型为圆心，让GO和KEGG结果围绕这个圆心旋转、聚焦、相互印证。

4.1 “靶点中心化”分析：让所有数据回归药物本源

药物分析的最大优势，是我们通常知道它的“起点”——即已知的靶点（Target）或主要表型（Phenotype）。这个已知信息，是破解GO/KEGG数据迷雾的唯一罗盘。我称之为“靶点中心化”策略，操作分为三步：

第一步：锚定“靶点功能坐标”

无论靶点是激酶、GPCR、离子通道还是核受体，都必须在GO和KEGG中为其精确定位。例如：

• 靶点为EGFR（表皮生长因子受体）：在GO中，其核心功能是“受体酪氨酸激酶活性”（GO:0004713）、“EGFR信号通路”（GO:0007173）；在KEGG中，它位于“EGFR tyrosine kinase inhibitor resistance”（map05225）和“PI3K-Akt signaling pathway”（map04151）的最上游。
• 靶点为HDAC（组蛋白去乙酰化酶）：GO中为“组蛋白去乙酰化酶活性”（GO:0004407）、“染色质修饰”（GO:0016569）；KEGG中则在“Histone modification”（map04650）和“Transcriptional misregulation in cancer”（map05202）中。

这个“坐标”不是静态的，而是动态的参考系。所有后续的GO/KEGG富集结果，都要问一句：“这个术语/通路，距离我的靶点坐标有多远？是直接下游（1步），还是间接下游（3步以上）？”

第二步：构建“距离-显著性”二维矩阵

将GO和KEGG结果整理成一个二维表格，X轴为“与靶点的生物学距离”（Distance），Y轴为“统计显著性”（-Log10(p)）。距离的定义如下：

• 0步（Direct）：靶点基因本身，或其直接互作蛋白（如EGFR的配体EGF、适配蛋白GRB2）。
• 1步（Immediate）：靶点的直接底物或下游转录因子（如EGFR磷酸化的STAT3、AKT；HDAC去乙酰化的H3K27）。
• 2步（Proximal）：1步分子的直接靶基因或效应器（如STAT3调控的BCL2L1；H3K27ac修饰的MYC启动子）。
• 3+步（Distal）：更远端的、常表现为细胞表型的基因（如凋亡执行者CASP3、细胞周期蛋白CCND1）。

然后，将每个富集到的GO术语或KEGG通路，根据其涉及的基因，估算其在该距离层级上的“重心”。例如，一个通路若包含3个0步基因、5个1步基因、2个2步基因，则其重心距离为1.2。将所有结果填入矩阵，我们会得到一张类似雷达图的分布。真正的药物机制信号，往往密集出现在“距离1-2步、显著性高”的象限。而那些“距离0步但p值不显著”的结果，恰恰提示：药物可能并未直接改变靶点表达，而是通过翻译后修饰（如磷酸化）或亚细胞定位来起效——这正是需要后续WB或IF实验验证的关键线索。

第三步：表型反向验证，锁定“功能-表型”因果链

最终，所有分析必须落脚到可测量的表型。例如，如果我们的药物预期效果是“抑制肿瘤细胞增殖”，那么GO/KEGG分析中，所有指向“细胞周期进程”（GO:0022402）、“DNA复制”（GO:0006260）、“E2F靶点”（KEGG map04110）的结果，都必须与实际的CCK-8、EdU或流式细胞周期检测结果进行比对。我坚持一个铁律：任何GO/KEGG富集结果，若不能在至少一项功能表型实验中得到定量支持，就必须视为待验证假说，不得写入论文结论。在某次抗纤维化药物项目中，GO分析强烈富集“细胞外基质组织”（GO:0030198），KEGG也显示“ECM-receptor interaction”（map04512）显著，但体外胶原分泌检测却无明显变化。深入排查发现，差异基因中富集的是“基质金属蛋白酶”（MMPs）家族，而非胶原合成基因（COL1A1/COL3A1），这提示药物的作用可能是“降解已有ECM”，而非“抑制新ECM合成”。这一修正，直接改变了后续动物模型的评价指标设计。

4.2 实战避坑：药物分析中GO/KEGG的五大高频雷区

基于多年处理各类药物数据的经验，我总结出五个极易被忽视、却足以颠覆结论的“雷区”，每一个都附有真实案例和解决方案：

雷区一：“背景基因集”选错，全盘皆输

最常见的错误，是直接使用KEGG或GO官方提供的“人类全基因组”作为背景集。但药物分析中，我们的样本是特定细胞系或组织（如肝癌HepG2细胞、人原代肺成纤维细胞）。这些细胞类型中，大量基因根本不表达（RNA-seq TPM=0），却仍被计入背景集，严重扭曲了超几何检验的概率模型。正确做法是：用该样本自身的“所有检测到的基因”（Detected Genes）作为背景集。在DESeq2中，这对应rowSums(counts(dds) > 0) > 1的基因；在limma中，则是avereps后!is.na(exprs)的基因。我曾帮某实验室复现一个失败的分析，仅将背景集从“全基因组”（~20,000基因）切换为“HepG2检测基因”（~12,000基因），就使“胆固醇代谢通路”的p值从0.15骤降至0.003，后续实验证实该药物确有降脂效应。

雷区二：“同源基因”混淆，张冠李戴

KEGG和GO数据库中，许多基因条目是“同源基因组”（Ortholog Group），如“KRT14/KRT15/KRT16/KRT17”在KEGG中可能共用一个KO号（K07222）。当你的差异分析鉴定出KRT14上调，而KEGG映射工具自动将其归入该KO组时，整个KO组都会被计为“变化”。但KRT15在该样本中可能根本未表达。解决方案：永远使用基因的精确ID（Entrez ID或Ensembl ID）进行映射，拒绝任何“组映射”。在R中，用bitr()函数（clusterProfiler）时，务必指定fromType = "ENSEMBL"或"ENTREZID"，而非"SYMBOL"，因为Symbol在不同物种或数据库版本中易歧义。

雷区三：“时间点”缺失，动态变静态

药物处理是一个时间过程。2小时、24小时、72小时的RNA-seq数据，反映的是截然不同的响应阶段：早期是信号转导，中期是转录重编程，晚期是表型执行。若将多个时间点的数据混合分析，GO/KEGG结果必然是“大杂烩”，失去时序逻辑。正确策略：对每个时间点单独进行GO/KEGG分析，然后用“术语一致性”（Term Consistency）来评估稳健性。即，一个术语若在2h、24h、72h三个时间点均显著富集（p<0.05），则其可靠性远高于仅在单一时间点显著的术语。在分析一种新型抗病毒药物时，我们发现“干扰素信号通路”仅在6h显著，而“细胞周期阻滞”在24h和72h均显著，这完美对应了其“早期激活先天免疫、晚期抑制病毒复制所需宿主细胞分裂”的双重机制。

雷区四：“剂量梯度”忽略，线性变非线性

药物效应常呈U型或倒U型曲线。低剂量可能无效，中剂量最佳，高剂量反而因毒性导致非特异性基因变化。若只分析一个剂量点，很可能捕获的是毒性反应而非治疗效应。解决方案：至少设置3个剂量点（如IC10, IC50, IC90），对每个剂量点单独分析，然后寻找“剂量依赖性富集”（Dose-dependent Enrichment）。即，某个GO术语的富集分数（如ES）随剂量增加而单调上升或下降。这比单点p值更能反映真实的药理学关系。在某次G蛋白偶联受体激动剂分析中，“cAMP代谢”通路在IC10时无变化，IC50时轻度上调，IC90时剧烈上调，且其上调幅度与细胞内cAMP浓度测定值高度线性相关（R²=0.98），这为机制提供了无可辩驳的生化证据。

雷区五：“对照选择”不当，噪音变信号

最致命的错误，是使用“未处理细胞”作为唯一对照。但药物溶剂（如DMSO）、处理操作（如换液、离心）本身就会引起基因表达变化。正确的对照必须是“溶剂对照”（Vehicle Control），即除不含药物外，其他所有条件（DMSO浓度、处理时间、操作步骤）与药物组完全一致。我见过太多案例，其中“氧化应激反应”通路的富集，最终被证实完全源于DMSO溶剂本身，而非药物。因此，在差异分析前，务必先对“药物组 vs 溶剂对照组”和“溶剂对照组 vs 未处理组”分别进行分析，只有前者显著而后者不显著的基因，才真正属于药物特异性响应。

5. 从分析到落地：一份可直接复用的药物GO/KEGG分析工作清单

纸上得来终觉浅，绝知此事要躬行。再精妙的理论，若不能转化为可执行、可复现的操作步骤，就只是空中楼阁。基于前述所有原则和经验，我为你梳理了一份详尽的、开箱即用的药物GO/KEGG分析工作清单。它不是抽象的方法论，而是我每天在电脑前实际敲下的命令、打开的网页、填写的参数。你可以把它打印出来，贴在显示器边框上，逐项打钩完成。

5.1 数据准备与质控：确保输入“干净无污染”

这是整个分析的地基，90%的后续问题都源于此步的疏忽。

• [ ] 原始数据来源确认：明确你的数据是来自RNA-seq（FASTQ）、微阵列（CEL）还是qPCR（CT值）。不同来源，预处理流程天壤之别。本文默认以RNA-seq的count矩阵为起点。
• [ ] 差异分析复核：不直接信任他人提供的差异基因列表。自己用DESeq2（推荐）或limma-voom重跑一遍。关键参数：
  • DESeq2：DESeqDataSetFromMatrix(countData, colData, design = ~ condition)；dds <- DESeq(dds)；res <- results(dds, contrast=c("condition","Drug","Control"))；res <- lfcShrink(dds, coef="condition_Drug_vs_Control", type="apeglm")。必须使用apeglm进行LFC收缩，它能极大减少低表达基因的假阳性。
  • 过滤标准：res <- res[!is.na(res$padj) & res$padj < 0.05 & abs(res$log2FoldChange) > 1, ]。注意，这里用padj（BH校正）而非pvalue，且|log2FC| > 1是药物分析的硬性门槛，低于此值的微小变化，生物学意义存疑。
• [ ] 背景基因集构建：如前所述，绝不使用“全基因组”。在DESeq2对象dds中执行：background_genes <- rownames(dds)[rowSums(counts(dds) > 0) > 1]。将background_genes保存为background_gene_list.txt，后续所有富集分析均以此为准。
• [ ] ID统一转换：将差异基因的ID（无论是Ensembl、Symbol还是RefSeq）统一转换为Entrez ID，因为GO/KEGG数据库最稳定。使用R包org.Hs.eg.db：

  library(org.Hs.eg.db) eg_ids <- mapIds(org.Hs.eg.db, keys=row.names(res), column="ENTREZID", keytype="ENSEMBL") # 过滤掉转换失败的NA eg_ids <- eg_ids[!is.na(eg_ids)] write.table(eg_ids, "diff_entrez_ids.txt", quote=F, row.names=F, col.names=F)

5.2 GO分析实操：分层、多算法、语义聚类

• [ ] 分层基因集生成：基于res对象，用以下代码生成三个txt文件：

  # 核心层 (|log2FC|>2 & padj<0.01) core_genes <- rownames(res)[abs(res$log2FoldChange) > 2 & res$padj < 0.01] # 扩展层 (|log2FC|>1 & padj<0.05) extended_genes <- rownames(res)[abs(res$log2FoldChange) > 1 & res$padj < 0.05] # 全景层 (padj<0.05) pan_genes <- rownames(res)[res$padj < 0.05] # 转换为Entrez ID并保存 write.table(mapIds(org.Hs.eg.db, keys=core_genes, column="ENTREZID", keytype="ENSEMBL"), "core_entrez.txt", ...) # 同理生成 extended_entrez.txt 和 pan_entrez.txt

• [ ] 双算法富集运行：
  • clusterProfiler（超几何）：

    library(clusterProfiler) ego