做生物信息的人大概都经历过这个时刻手里拿着一批差异基因或者课题组筛出来的候选靶点每个单独看都有文献支持但放在一起到底谁在上下游、谁是核心节点、哪些蛋白天然抱团单靠Excel拉表格根本讲不清楚。我第一次正经跑蛋白互作网络就是在这种场景下逼出来的一批来自转录组和蛋白质组交叉验证的基因数量说多不多说少不少十几个二十个发文章或者给导师汇报的时候总不能只扔一张火山图吧。当时最先想到的就是STRING数据库配合R语言做清洗、构图和拓扑分析。这套组合的好处很直接STRING把全世界已发表和被数据库收录的蛋白互作证据集中在一起还给出一个综合打分R语言则把导出的数据变成可以统计、筛选、重复操作的分析流程。对于需要处理多组学结果、做药物靶点筛选、或者纯想给某个通路画一张“关系地图”的人来说这个是性价比最高的入门方案。这篇文章我会把我常用的整个流程完整走一遍包括网页端参数怎么设置、TSV文件怎么导入R、igraph怎么画网络、以及我在这个过程里踩过的坑。代码我尽量贴全你可以直接复制改成自己的基因列表。1. 为什么蛋白互作网络分析都在用STRING1.1 STRING到底干了什么事从“谁和谁认识”到“证据加权”很多人都听说过“蛋白互作网络”这个说法但真正上手时会发现难点不在于画图而在于怎么判断两个人之间的事靠不靠谱。生物体内大多数蛋白不是单打独斗的它们组成复合体、形成信号通路、响应外界刺激但文献里这一条互作可能来自酵母双杂交那一条来自免疫共沉淀另外一些可能只是高通量筛选里的一个候选可靠性天差地别。STRINGSearch Tool for the Retrieval of Interacting Genes/Proteins做的事情就是把不同来源的互作证据收集起来覆盖上万种物种、数千万个蛋白然后针对每一对蛋白给出一个综合分数。这个分数不是简单的“有”或者“没有”而是把多个维度的证据加权融合在一起。你拿到手的不只是一堆蛋白的名字而是一套带可信度标签的关系数据。所以做蛋白互作网络分析第一步不是急着画图而是先想清楚我手中这批基因在STRING里能够匹配到的有效互作到底有多少这些互作的证据来源是否足够支持我的结论。另外STRING检索时支持多种标识符基因Symbol、Ensembl ID、UniProt ID都能识别。我自己常用的基因列表通常来自转录组差异表达分析输出的是基因Symbol直接粘贴到STRING搜索框就能用。不过这里有一个小坑如果你用的是旧式别名比如有些人习惯把TP53BP1写成53BP1把H2AFX写成H2AXSTRING虽然内置同义词映射但未必每个别名都能命中。所以最好提前用HGNC的标准Symbol去跑匹配率会明显提高这一点后面我还会单独说。1.2 从数据源看懂打分机制实验、数据库、文本挖掘与共表达STRING里每条互作的综合分数由几个证据通道组合而成。网页导出的full TSV文件里会明确列出neighborhood_on_chromosome、experiments、databases、textmining等列名。很多人看到这些列名就晕其实每个名词对应一类证据来源experiments来自实验验证的物理互作或生化互作比如co-IP、pull-down、酵母双杂交、X射线晶体结构等。这类证据通常被认为可信度较高。databases来自第三方通路数据库或蛋白复合物数据库收录的互作比如Reactome、KEGG等。textmining通过文本挖掘从文献摘要和全文中提取的共现关系。这类证据覆盖面广但也最容易出现“看起来相关、实际并间接”的情况。coexpression基于基因在不同条件下的表达量相关性推断出的潜在关联比如大量RNA-seq样本中两个基因共表达这类证据偏向功能关联而非物理结合。neighborhood_on_chromosome、fusion、cooccurrence这些在细菌、古菌基因组里比较常见对人类基因的贡献相对较小。我在实际操作中有一个习惯如果最终目的是找核心Hub基因、做药物靶点筛选我会把textmining通道的贡献单独看一下。默认网络里如果textmining占比过高得到的边可能会比较“软”因为两篇文献里提到两个基因并不代表两个蛋白一定在同一个复合物里工作。此时可以在STRING网页端调整证据通道的显示或者在R里读入full TSV后只保留那些experiments或databases分数较高的边构建子网络。后面代码部分我会给出具体做法。2. STRING网页端实操从基因列表到TSV导出2.1 关键参数怎么选物种、网络类型和置信度阈值进入STRING主页后最常用的入口是“Search”下的“Multiple proteins”。贴入基因列表前第一步不是急着点Search而是先确认物种。这个整不好后面全是白干。我见过不少初学者把人类基因贴进去默认物种却是某种模式生物最后匹配率惨不忍睹。STRING支持按物种名搜索输入Homo sapiens即可。如果你分析的是小鼠、斑马鱼、水稻或者某种细菌一定记得先改成对应物种。接下来是网络类型选择。STRING里有“full STRING network”和“physical subnetwork”两个大的方向。前者包含所有类型的互作证据包括功能关联和共表达后者限定为物理互作比如直接结合、形成复合体。如果研究目标是信号通路、调控机制full STRING network更常用如果你想找的是一条明确的蛋白复合物核心链physical subnetwork更干净。我的建议是第一轮跑full network看全局情况然后再用physical subnetwork做敏感性验证。置信度阈值是另一个高频设置项。STRING默认的minimum required interaction score是0.400也就是medium confidence。这个值到底选多少取决于你的网络规模和目的0.150low基本什么都留着网络会非常密集适合做探索性观察但不太适合直接下结论。0.400medium默认值兼顾覆盖率和可信度用于常规差异基因网络分析比较合适。0.700high边少而精留下来的大多是高置信互作适合寻找核心复合物或验证特定通路。0.900highest非常严格通常只剩实验证据很强的少数互作。我平时跑一批差异基因会先用0.400出一张全貌图再看0.700的收敛情况。如果两个阈值下核心Hub基因基本一致那结论就相对稳如果差异巨大就得回头检查基因列表质量比如是否混入了大量表达量低、注释不清晰的基因。2.2 一次真实导出过程从输入基因到下载TSV我用一个和DNA损伤修复相关的基因列表做例子实际跑一遍ATM、ATR、BRCA1、BRCA2、TP53、CHEK1、CHEK2、RAD51、TP53BP1、H2AFX、PALB2、BARD1、MRE11、EXO1。这14个基因在DNA损伤应答里都有明确功能用来演示网络分析很合适。在STRING搜索页粘贴这串基因物种选Homo sapiens点击Search网站会跳转到结果页面。这里你会看到一张默认的网络图图上每个节点是一个蛋白点与点之间的连线就是预测或验证过的互作关系。颜色和线型代表证据类型把鼠标悬停在边上可以查看具体来源。这个页面还可以调整“meaning of network edges”默认是confidence线粗代表综合分数高也可以切换成experiments、databases、textmining看不同证据通道的分布。确认网络形态合理后点击页面顶部的“Exports”按钮在下载区域找到“as simple tabular text output (TSV)”下载得到的就是标准互作表。文件名一般是string_interactions.tsv。这个文件就是R语言分析的原材料。如果你需要更丰富的信息比如每个证据通道的独立分数可以选择下载“full TSV”。两种我都用过日常分析其实simple TSV就够但如果你想在R里做“只保留实验证据”这种精细过滤full TSV更有用。下载后别急着关网页。我强烈建议顺手截几张图尤其是“Summary”页面的统计信息和网络图原图写报告或者组会汇报时可以直接用。另外STRING的“Analysis”页面会给出GO富集、KEGG通路、UniProt关键词等一系列富集结果这些后面作为网络分析的补充讨论非常有价值导出成表格收好。3. R语言读数据与网络构建代码逐段拆解3.1 数据读入正确处理注释行和列名收到string_interactions.tsv之后第一步是在R里读进来。这个文件本身是以Tab分隔的文本但有个小陷阱文件开头有一行以#开头的注释记录版本信息和网络链接。如果直接用read.delim读注释行可能会被当成表头导致后面的列名全是错的。我最早吃这个亏时names()输出X.STRING.network...整个人愣了一下才反应过来。解决办法是给read.delim加上comment.char#参数让它忽略#开头的内容。读取代码不算复杂library(igraph) # 读取STRING导出的TSV文件 ppi - read.delim(string_interactions.tsv, sep \t, header TRUE, comment.char #, stringsAsFactors FALSE) # 查看数据结构 names(ppi) head(ppi)读进来以后建议先看一眼列名。如果下载的是simple TSV正常会包含node1、node2、combined_score等列如果是full TSV列名会更多包含neighborhood_on_chromosome、experiments、databases、textmining等。不同版本的STRING导出格式偶尔会有细微差异所以先用names()确认是最稳妥的。另外如果文件里某列末尾有奇怪的字符尤其是Windows系统下载的文件偶尔会在行尾多出\r我会加一步gsub清洗把\r去掉再读。这一步虽然不属于标准流程但确实能避免不少莫名其妙的对不上。3.2 构建igraph网络对象节点、边、权重数据读进来以后下一步是把表格转成网络图对象。这里最常用的是igraph包它对图论分析的封装非常成熟性能也够用。核心代码就一行# 用前两列node1和node2构建无向网络 g - graph_from_data_frame(ppi[, c(node1, node2)], directed FALSE) # 把combined_score作为边的权重 E(g)$weight - ppi$combined_score这里有一个细节STRING导出的表格里可能包含重复边比如同一对蛋白在多条证据中都出现或者存在node1和node1自己连着一条自环的情况。这些会对后续的度、介数、聚类系数等拓扑指标造成干扰所以构图前我会先做一次清理g - simplify(g, remove.multiple TRUE, remove.loops TRUE)经过这步处理后网络里每个节点和每条边都是唯一。之后看一眼网络规模对整体情况有个数vcount(g) ecount(g)14个输入基因的列表实际网络可能不止14个节点因为STRING默认允许在结果中加入少量与输入蛋白直接互作的额外节点。这是它的设计逻辑你的基因列表之间如果互作关系很少网络会显得很空加入第一层邻居有助于补充上下文。但如果你做定量分析比如统计网络中的关键节点就需要明确要不要让这些额外节点留在里面。想只保留自己输入的基因时可以用induced_subgraph做子网筛选代码我放在3.3节。3.3 子网筛选与证据通道过滤让网络更贴近你的生物学问题实际科研场景里网络不是越全越好。假设你的基因列表本身只有14个基因但STRING返回了一张40个节点的网络你直接拿去统计度排名靠前的可能全是那些“被添加进来”的非目标蛋白这会干扰你对核心功能的判断。解决办法是两类第一类筛选回输入基因列表形成的子网第二类按证据通道和分数过滤边。先看子网筛选。假设我的输入基因列表是my_genes - c(ATM,ATR,BRCA1,BRCA2,TP53,CHEK1,CHEK2, RAD51,TP53BP1,H2AFX,PALB2,BARD1,MRE11,EXO1) g_sub - induced_subgraph(g, which(V(g)$name %in% my_genes)) vcount(g_sub) ecount(g_sub)这段代码的关键是which(V(g)$name %in% my_genes)它会找到原始网络里那些名字在输入列表中的节点位置然后诱导出这些节点以及它们之间的边。注意这里只会保留输入基因彼此之间的互作如果某些输入基因在STRING里和任何其他输入基因都没有边它们会成为孤立点。孤立点本身也是信息说明这批基因之间的直接互作证据确实薄弱。再看证据通道过滤。如果你下载的是full TSV可以用类似下面这种方式构建一个只包含实验证据分数达到0.700的边网络# 过滤实验证据分数 0.7 ppi_exp - ppi[ppi$experiments 0.7, c(node1, node2)] g_exp - graph_from_data_frame(ppi_exp, directed FALSE) g_exp - simplify(g_exp, remove.multiple TRUE, remove.loops TRUE)这种过滤方式对寻找物理互作核心非常有帮助。文本挖掘带来的边往往有很多假阳性的嫌疑而实验证据的边更接近真实的复合物和直接结合关系。我自己在汇报时会准备两套图一套全证据版本用来展示整体格局一套experiments过滤版本用来锁定高置信核心两张图配合着讲逻辑会清楚很多。4. 网络可视化与核心节点挖掘4.1 让图能看出信息节点大小、颜色、布局和保存可视化不是单纯的画图而是要把网络里的结构信息体现出来。igraph默认画法很朴素直接plot会导致节点大小一样、颜色一样、标签挤成一团根本看不出谁是核心。我会做三件事节点大小映射到度、边的粗细映射到combined_score、颜色映射到社区模块。# 设置随机种子保证每次布局一致 set.seed(42) # 节点度 deg - degree(g_sub) # 用Fruchterman-Reingold布局适合中小规模网络 plot(g_sub, vertex.size deg * 3 5, vertex.color #7FB3D5, vertex.frame.color #2E86C1, vertex.label.cex 0.75, vertex.label.color gray20, edge.width E(g_sub)$weight * 3, layout layout_with_fr(g_sub), main DNA损伤修复相关蛋白互作网络)如果你觉得基础plot排版不够美观可以用ggraph包配合ggplot2的风格来画但igraph的灵活性和上手难度对于第一轮探索其实更友好。我的建议是先用igraph快速出图看结构和Hub节点等确定好展示版本后再考虑用ggraph或Cytoscape精修。布局算法也值得说两句。igraph里常见的layout_with_fr是力导向布局节点之间像有弹簧一样互相排斥边多的节点被拉向中心比较适合展示尺度在100个节点以内的网络。节点数量一旦上到几百FR布局会非常慢建议改用layout_with_kk或者先用社区发现算法分组后再布局。遇到几百上千个节点的大图说实话折线图或者Cytoscape的yFiles布局更好使igraph并不是万能的。导出图片时我一般用pdf设备矢量图放进论文里不会糊pdf(ppi_network.pdf, width 10, height 8) plot(g_sub, vertex.size deg * 3 5, vertex.color #7FB3D5, vertex.frame.color #2E86C1, vertex.label.cex 0.75, vertex.label.color gray20, edge.width E(g_sub)$weight * 3, layout layout_with_fr(g_sub)) dev.off()以我的经验png这类位图格式在屏幕上看着不错但投稿时编辑经常要求矢量图。PDF是最通用的选择。如果你在RStudio里画还可以直接点击Export菜单选择PDF或PNG但那样不好自动化推荐还是用代码控制输出。4.2 Hub基因和关键模块怎么找度、介数、社区发现网络图看个大概之后量化分析必须跟上。最基础的指标是度也就是每个节点连接了多少个邻居。度排名靠前的基因通常被称为Hub基因。在疾病机制或药物靶点研究里Hub基因往往提示该蛋白在网络中处于中心位置调控的潜在影响面更大。# 按度排序取前10 top10_deg - sort(degree(g_sub), decreasing TRUE)[1:10] print(top10_deg)光看度还不够介数中心性betweenness是另一个我常看的指标。它衡量的是在所有节点对的最短路径中经过某个节点的比例。直观理解就是一个人虽然不是朋友圈里加好友最多的人但所有消息传递都要经过他中转那他在信息流通中的角色也非常关键。放在蛋白互作网络里高介数的蛋白往往是不同模块之间的连接者可能是信号通路的转换枢纽。btw - betweenness(g_sub) top10_btw - sort(btw, decreasing TRUE)[1:10] print(top10_btw)把度和介数的Top基因列表拿过来对比你会发现有些基因两个榜单都在这些就是核心中的核心有些只在度榜上说明它连接很多但比较局部有些只在介数榜上说明它是一个重要的跨界节点。模块检测同样很有价值。社区发现算法可以把网络划分成几个内部联系紧密、彼此联系较少的子模块。我常用的是walktrap和fast_greedy# 社区发现 cl - cluster_walktrap(g_sub, weights E(g_sub)$weight) membership(cl) modularity(cl)把社区分配结果映射到节点颜色能在图上直观地看到哪些蛋白抱团往往是同一个生物学通路或复合物的成员。比如在DNA损伤修复的例子中你可能看到一条由ATM、BRCA1、BRCA2、PALB2、BARD1组成的支路和另一条由CHEK1、CHEK2、TP53组成的支路。这种模块划分结果和STRING自带的富集分析结合起来就能对“这批基因在哪些通路协同工作”给出比较完整的回答。4.3 导出给Cytoscape继续加工一种常规交接方式igraph能做的统计和基础可视化已经很强但你真要把网络图打磨成发文章级别的效果Cytoscape往往是绕不开的工具。igraph和Cytoscape之间的交接我推荐用GraphML格式它能保留节点属性、边权重、社区分组等元数据。# 把社区分组写入节点属性 V(g_sub)$module - membership(cl) # 导出GraphML write_graph(g_sub, ppi_network.graphml, format graphml)在Cytoscape里导入GraphML后可以用“Style”面板按module给节点着色按degree设置节点大小用NetworkAnalyzer插件重新计算一次拓扑参数还能用yFiles Layout里的“Organic”布局让整张图更好看。很多人觉得Cytoscape上手门槛高其实流程并不复杂关键是先把igraph这边的数据分析做完Cytoscape只负责呈现和微调这样两边的工作量都不大。我个人习惯是R里面出结果表Cytoscape出最终展示图。分析过程中的反复试错在R里完成因为脚本化操作可以随时重跑最终封面图用Cytoscape慢慢调因为交互式调整到底还是方便。这个配合方式用顺了之后处理几十批基因列表都不慌。5. 高频报错与实操心得5.1 常见问题和排查速查表字符串这个流程跑的次数多了我自己积累了一份问题速查表列举几种最常踩的坑和对应的解决办法现象可能原因解决方法读入TSV后列名全是X.node1之类注释行被当成表头重新读取加上comment.char#参数基因匹配率只有30%-50%基因名不规范、别名未被识别先用HGNC统一为标准Symbol再查一下物种选择是否正确STRING识别出的节点和输入基因数差别很大默认加入了额外的互作伙伴在结果页设置里降低“additional nodes”或者在R里用induced_subgraph筛选网络特别稀疏几乎没有边阈值过高或输入基因之间本身缺乏直接互作证据降到0.400或改用full STRING network增加证据通道网络特别密所有节点连成一团阈值过低或textmining证据占比过高提高阈值到0.700在设置里关闭textmining通道或按experiments列过滤R提示找不到对象node1列名里带了其他前缀先用names()检查列名再用colnames()重命名write.graph导出的GraphML在Cytoscape里打不开节点属性类型不兼容比如NA值转换前用V(g)$attr[is.na(V(g)$attr)] - NA 清理缺失值第2行“基因名未识别”真的是高频问题。STRING虽然做了同义词映射但你在文献里看到的很多名称其实是非官方写法。比如“53BP1”对应的HGNC标准Symbol是“TP53BP1”“H2AX”对应“H2AFX”“p21”对应“CDKN1A”。如果一开始没有统一成标准Symbol匹配率会很难看而且别人复现你数据时也容易对不上。所以我现在拿到任何基因列表的第一件事就是去HGNC网站上做一次批量标准化或者在Ensembl BioMart里转换这个习惯帮我省掉了至少一半的检查时间。第5行“网络过密”也值得特意说一下。很多做药物靶点筛选的人喜欢追求“多条边”的富集感觉得网络越密说明这批基因越相关。但网络过密往往意味着无法区分关键节点所有基因互相连着Hub分析等于失效。所以密度过大时第一反应不应该是接受而是检查是不是爬进了一堆低质量textmining边。5.2 几个我的习惯性操作除了那张速查表还有几个习惯我每次都会执行。第一个是跑完网络后一定会回到原始基因列表重新确认网络里出现的节点名和输入基因的对应关系。STRING会自动帮你做ID映射但偶尔会遇到一个输入Symbol跳转成了另一个种属的同名同源基因尤其是小鼠研究里一个不小心的物种混淆会导致整张网络错得离谱。第二个习惯是保存每一步的关键对象。R的RDS格式保存rds对象非常方便saveRDS(g_sub, g_sub.rds) saveRDS(cl, clusters.rds)这样即使脚本中途出错或者过了两周回来想继续分析也不用从头再跑一遍。第三个习惯是在汇报结果时永远同时报告网络在0.400和0.700两个阈值下的节点数、边数、Hub基因列表。这不只是应对审稿人的稳健性检验也是帮自己确认结论是否对参数敏感。如果两个阈值下关键Hub基因高度一致那这个结果就很有说服力如果完全变样那就得思考这批基因之间的互作证据本身是否足够强。6. 进阶扩展当课题不止于一张网络图6.1 把STRING结果接进富集分析网络分析只是开始真正发文章时审稿人往往会问这些关键Hub基因富集在哪些通路STRING网页端自带的Analysis页面可以直接给出一批GO和KEGG富集结果但如果你想在R里做更系统、可重复的富集分析可以用clusterProfiler。这种方式的操作逻辑是先把网络中的核心基因列表提取出来再用Ensembl或org.Hs.eg.db做ID转换最后跑go富集和KEGG通路分析。R代码如下# 提取Hub基因按度数取前30 hub_genes - names(sort(degree(g_sub), decreasing TRUE)[1:30]) writeLines(hub_genes, hub_genes.txt)拿到Hub基因列表后你可以把它粘贴到clusterProfiler的enrichGO或enrichKEGG接口里结合STRING的富集结果做交叉验证。如果网络分析提示某个模块以BRCA1/PALB2/BARD1为核心而KEGG富集又显示同源重组通路显著富集那这两个结果就互相印证了结论会扎实很多。6.2 用STRINGdb包在R里直接拉数据除了网页端导出TSV还有另一个路线值得知道Bioconductor的STRINGdb包。它可以让你在R环境里直接调用STRING的API查询互作、获取富集结果、甚至可以给基因打分。# BiocManager::install(STRINGdb) library(STRINGdb) string_db - STRINGdb$new(version 12.0, species 9606, score_threshold 400)这里要说明的是STRINGdb的学习成本和维护状态都不如TSVigraph这套组合来得稳定而且网络请求受数据库服务器状态影响跑大批量数据时速度不一定理想。如果你只是处理十几个到两三百个基因我仍然建议网页导出TSV简单可靠运行过程中也能随时去网页上检查数据。STRINGdb更适合那些需要把STRING整合进自动化pipeline、每天跑不同基因列表的场景。6.3 把方法论沉淀成自己的模板流程最后分享一点方法论层面的心得。蛋白互作网络分析本质上是一个从基因列表到关系洞察的加工过程。它的输入是一串基因名输出却不是一张图而是一套结构化的判断哪些蛋白是核心枢纽哪些模块对应哪些生物学功能这些模块之间如何衔接。为了让结果可复现我强烈建议把整个分析沉淀成一套固定的脚本模板。我自己电脑上有一个叫做ppi_pipeline.R的脚本开头是读TSV中间是构建网络和输出统计指标末尾是出图。每次拿到新的基因列表只需要改文件路径和my_genes向量运行一遍就能得到基础结果。这种模板化的习惯比每次临时敲代码高效得多也能避免重复踩坑。顺便说一句经常有人问我为什么不用Python做这个。Python的networkx确实也能实现类似功能但R环境里和DESeq2、clusterProfiler这些组学分析工具的衔接更顺畅我个人觉得R在生物信息这个场景下更顺手。如果你刚接触蛋白互作网络或者正卡在某一步跑不通建议你按这篇文章的流程走一遍文件的读取、参数设置、代码运行整个周期大概半小时。真正难的其实不是跑代码而是跑完之后能不能把网络结果和你的生物学问题对应上。第一步先把自己手头那批基因列表扔进去跑通再慢慢地去调整阈值、挖掘模块你会发现这一套工具组合能做的事情远比想象中多。