Constructingmodule maps for integrated analysis of heterogeneous biological networks
DavidAmar and Ron Shamir*
BlavatnikSchool of Computer Science, Tel Aviv University, Tel Aviv 69978, Israel
Received September 24,2013; Revised December 17, 2013; Accepted January 10, 2014
(译者注:在本文中,译者尽量使用恰当的术语进行翻译,但是苦于专业知识短缺,大部分都保留了原文术语。译文中标记为红色的部分是需要重点学习的地方,原因有三:1,未能找到合适的术语进行翻译,术语理解存在分歧。2,原文关键部分,写出了关键的技术等。3,原文过于复杂,在翻译时可能存在歧义。其中大部分标记的部分,译者都已经附上原文作为参考。最后,译者仍建议对照原文进行阅读。 最后,本文未上传原文,需要对照观看的可以自行搜索下载,原文图片也未上传,请在原文中观看)
摘要:
针对大规模不均匀的omic data的综合分析的方法亟需提高。本文提出了一个基于网络的方法来解决这个问题。对于给定的两个网络,分别代表不同类型的基因相互作用,本文建立了一个模块(基因)链接网络,在第一个网络中强连通的,链接表示第二个网络中的强连通的模块链接。本文提出的新型算法在三个独特的领域,经过模拟和真实数据的验证,都比过去的方法更加有效。第一个,通过分析在酵母中蛋白—蛋白相互作用和负遗传作用,我们发现在蛋白复合物之间存在上位关系。第二个,我们分析在想酵母中蛋白—蛋白相互作用、DNA损伤具体的积极基因的相互作用和在蛋白复合物之间的显示功能的重新连接,显示出了DNA损害反应的新型机制。最后,使用肺癌患者的非小细胞中的吸血虫,我们分析全局差分工表达网络和疾病相关的差分共表达以及识别免疫激活过程的两个模块之间的相关性急剧下降网络,在可能的微型RNA控制下。我们的研究论证modulemaps 在分析不均匀的高通量(大数据)omicdata是一个非常有用强大的工。
Introduction:
生物网络可以对生物系统提供一个综述,帮助更好滴理解生物系统、基因功能和其他分子化合物。在一些系统中,它已经被用来进行基因相互作用、基因功能和疾病基因的关联的预测。
在这些网络中,节点代表分子化合物,边代表相关性。例如,在PPI(蛋白相互作用)网络中,节点表示蛋白质,边表示物理相互作用。在GI(基因相互作用)网络中,节点表示基因,边表示双敲除扰动下生物体的健康,主要包含两种边:减轻(阳性)GIs和加重GIs。在减轻GIs中,也别叫做阳性GIs,生物体的健康在双敲除扰动之后比基于单敲除扰动所做的期望要更好。在加重GIs中,也叫做阴性GIs,正好相反。在基因共表达网络中,节点表示基因,边对两个基因之间的表现相关性进行数值评分(即为有权图)。在基因差分相关(DC)网络中,边对两个样本集合之间基因对关联的差异进行评分。随着生物网络使用和普及(网络类型),对大数据进行挖掘的计算方法变得尤为重要。
计算方法利用多个网络所得到的结果往往比只分析单个网络所得结果要更好。比如,PPI网络和基因共表达网络的综合分析用来探测共表达而且在PPI网络中链接的基因集合。这样的分析方法比标准的聚类算法更加优秀,并被成功用于基因功能预测。阳性和阴性GIs数据被用于探测基因组内和基因组之间的异位显性。在阳性GIs往往发生在补偿途径之间,阴性GIs往往发生在途径(pathway or complexes)之内的前提下,GIs数据的分析用于将功能性基因模块之间的上位关系(epistaticrelations)显示为网络图形。在PPI模块网络中增加连通性约束之后,效果取得显著提升已被报道。针对多个网络建立一个summary map的能力可以进行识别被发现模块之间的关联分析,因此,相比于针对单个网络的标准聚类方法,提高了结果的解释能力(interpretability)。
建立在针对特定一些网络组(对,pair)的优先研究的基础上,我们介绍并研究建立一个针对两个生物网络H和G的summary map的基本问题,在两个网络H和G中所有的节点都相应代表相同的基因或蛋白质,而且两个网路中的边各自表示一种独特的关系。Summary map中的节点是网络H中强连通的基因集合,并且集合成对地被边链接起来。其中的边表示了在网络G中强连通的两个基因集合。目的是通过优化一个特定的目标函数,找到网络H中的基因模块,同时找出通过网络G找出模块之间的相互关系,我们称这个计算问题为“module map problem”。
针对这个问题迄今为止大多数算法用于找到among pathways的上位相互作用的summarymap。Kelley and Ideker 提出了一个方法,基于图中局部搜索来发现关联模块对。Ulitsky et al.使用一组网络H作为初始点,然后通过合并模块来提升效果。一个近似的方法最近被出用于分析基因共表达和基因差分相关(DC)网络。用于发现基因组的对这些网络中的节点分析,或多或少地与一类生物体是相关的。虽然之前用于解决“module map problem”的算法被证实是有价值的,但是针对这个问题的完全分析、这些算法在不同情况下的有点和缺陷的分析也是必要的。
在大多数情况下,找到最优的modulemap的问题的计算复杂度是NP难度的,因为它包含了H的聚类作为一个子问题。因此,启发式方法就被使用。这类算法通常包含两个阶段。我们称第一个阶段为启动阶段(“initiators”):算法找到一个初始解,也许其中包含了许多小模块。第二个阶段称为改进阶段(”improvers”):算法通过预先定义的目标函数来改进初始解。使用不同的两个阶段的组合,可以形成各种各样的算法。
在此,我们研究新奇的和现存的initiators和improvers两个阶段。我们发现一个基于带有统计学上全局提升改进措施的maximal bicliques in G 的新启动阶段(“initiators”),效果总是不差于(consistentlybetter or equal to)基于综合的和真实的几种类型。我们称这个resulting algorithm 为ModMap。我们将这个方法应用到三种生物情景下的实验数据中:1,使用酵母PPIs和阴性GIs,我们发现蛋白复合物之间的上位关系(epistatic relations)。2,使用酵母PPIs和特定DNA损伤阳性GIs,我们发现和DNA损伤应答相关的蛋白复合物中的新链接的形成。3,使用针对肺癌组织非小细胞(NSCLC)的基因共表达谱(geneexpression profiles)的DC分析,我们确认在免疫活化过程和发现特定疾病(disease-specific)的微型RNA(microRNAs)之间特定疾病相关损失。
材料和方法 Materials and methods
Module map problem 的定义
该问题的输入是两个网络H=(V,EH,WH)和G=(V,EG,WG),定义在相同的顶点集合上。这些网络可以是有权重的,也可以没有。目标是找到一个可以概述(summary)两个网络的module map。Module map是一个图F=(M,L),其中M是不相交的节点集合的集合,不相交的点集合叫做module,M=(M1,…..,Mp),Mi⊆V, Mi∩Mj = ∅,并且L是module对的一个集合{(U1,V1),….,(Up,Vp)},其中每一个Ui和Vi都在M中。这些对被称为map links(图中的边)。另外,每一个module必须至少连接到不少于一个其他module上面。简单来讲,我们的目标就是找到一个modulemap,其中每个module相当于H的一个重子图(heavysubgraph),并且每条边表示两个module在G中的重型偶子图(heavybipartite subgraph)。关于重子图(heavysubgraph)的正式的概念介绍将会在后面引入。图片1D展示了一个玩具例子,包含两个无权图和他们的modulemap。
之前有关建立modulemap的算法在定义目标函数和链接时有所差异。DICER算法每次找到一堆链接的modules。如果他们之间的权重和WG足够高,则认为这一对modules是有边链接的。我们称这个方法为“local”,因为它每次只能找到一个module对。Ulitsky算法旨在最大化“global score”,也就是H中所有module内的score加上G中所有边(links)的总和。除了累加global score之外,modules之间的边(links)只有在通过统计数据测试后,才能建立。我们称第二个方法为”global”。两个方法同事用来确定modules和links。
图片1D展示了局部法和全局法的不同。假设在两个图中边的权重(edgeweight)是1,不存在的边(non-edge)权重为-1,局部法使用的阈值为0,两个modules之间WG的总和高于0时则意味着产生一个link。在两个方法中,modules都是H中高密度的节点群集(clusters of nodes)。按照这两个方法,module 1 将链接到module2上,它的局部分是4(8条edge,4条non-edge),全局分析中链接的P值(P-value)<0.05,并且modulepair的总分值为13(modulescore 6+3 + link score 4)。Module 2和module3之间的权重WG的总和是-4(10edges and 14 non-edges),所以局部法拒绝了链接。但是,全局法依然要链接module2 和3:链接值P-value很大(P=0.039),并且增加这个链接后,将把全局map的分值提升至24,[13for the (1,2) pair +15 for module 3–4 for the (2,3) link]。这个列子阐述了全局法在稀疏图上的优势,在这样的图中,大的modules并不像期望的那样紧密地分布在一起。
算法 Algorithm
我们进行了系统研究并改进了一些的双相算法,用来做modulemap detection,用于发现一个初始解(可能包含很多小的modules),然后进行改进。我们第一阶段为启动阶段,第二个阶段为改进阶段。为简单起见,我们描述算法,假定所有正向权重的边都被认为是很重要的。对于无权图,我们假设边的权重为1,non-edge的权重为-1,对于有权图,所有的边(点对)都有权重,所以不存在non-edges。
启动阶段 Initiators
我们测验五种不同的启动阶段的算法:1),DICER(24),每次发现一对链接的modules。2),图H的分层次聚类(25),它发现一组modules。3),节点增加贪心算法来找到H中的modules。4),DICERk,DICER算法的一个变种,将module的最小size限制到k.5),使用穷举法来列举出图G中最全面的bicliques,然后使用DICER方法来找到剩余的modules。我们称第五种算法为MBC-DICER,对于所有的启动详情请看追加文本和追加的图片S1。每一个启动阶段(initiator)都会创立一个module集合,但是在被使用聚类算法建立起来的map中,modules不一定linked。
改进阶段improver
局部改进法通过增加单个节点到module中或合并两个modulemap links 来扩展module。这种方法的一个缺陷就是它不能初始解中不存在的新module。另一个缺陷就是他不能合并一个module,如果他的两个部分链接到两个不同而且彼此没有链接到一起的modules.请看补充的图片S2。接下来,会介绍全局改进法,可以经常用来克服以上两个问题。
我们的全局改进法基于(17)中的过程。让M={M1,…,Mn}代表一组彼此不相交的点集合(集合可以是单个点,而且和任何其他集合不相链接)。给定集合(U,V),U,VϵM,并且xϵU, x和V的链接的重要程度(thesignificance of the linkage of x with V)通常使用Wilcoxon rank-sum test来比较x与V中的边权重WG 和 x与所有不在V的边权重来计算的。这样所有U和V中节点的P-value值都被计算得到,并且混合使用了Stoufer’smethod。如果最终的P-value值p(U,V)至多是α,然后U和V则被map中的一个link链接起来。让L={(U1,V1),…,(Up,Vp)}代表links的结果集合。
全局score的解决方案是Mi内所有边的权重WH的和 加上 linked node sets 之间的边权重WG的和:
如果全局score增加并且新链接可以通过significancetest,则这个改进阶段合并一对nodesets(two modules or a module and a single gene)。考虑到一次合并将会创建一个新的moduleY,为了计算新的全局score,对于在M中的其他所有的moduleZ需要重新计算p(Y,Z)。这个过程使用贪心、迭代(greedily,iteratively)的方式进行,产生最好提升结果的合并方案将会被进行,直到没有可以提升全局score的合并方案为止。
我们改良之前提到的算法是的如下所提到的大规模图能够进行快速分析。首先,当计算p(U,V),我们考虑G’(G的无权重版)中的links。我们使用一个超几何的测试方法来鉴定一个节点是否有足够数量的边在G’的oppositeset中(例如,从一个节点vϵV到集合U),合并前所有节点的P-values都使用Fisher’smethod。如果resulting value ≤α,集合U和V将会链接在一起。这样的测试方法更快,而且对于G,它将提供和Wilcoxon test 一样质量(equalquality)的maps。带权重的测试,比如Wilcoxontest,并非总是适用于基因模块(genemodule)之间的linkage检测。例如,在DC图中,一个强link必须包含很多阳性边,但是Wilcoxontest却只关心edgescore 的排名。
第二,我们使用另一个参数β>>α,如果在一些情况下(at some point)两个集合之间可能存在的link的P-value至少是β,我们则称这两个集合是“anti-linked”。在一些原始的算法中(original algorithm),当考虑合并两个集合U和V到W中时,W和其他每一个集合Y可能存在的links都必须计算。但是,如果U和Y是“anti-linked”,或者V和Y是“anti-linked”,则我们将W和Y也标记为“anti-linked”,避免需要进一步考虑可能存在的链接(W,Y)。在实践中,我们使用α=0.005来处理类似在(17)中的酵母数据(yeast data),并且测试基因表达数据(gene expression data)的一些情况(补充文本中有说明)。这就使得在实践中在保证了结果质量无损失(without loss of solutionquality)的情况下,实现了两倍甚至更多的加速处理。
模拟(simulations)
我们最开始建立两个空的500个node的图H和G,然后加入边来创建一个完美的modulemap,其中modules都是H中的一个小集合(cliques),links都是G中的bicliques。这个modulemap的拓扑结构(M,L)是一个|M|=6的随机树(randomtree)。然后我们加入twoH-cliques and two G-bicliques 到图中来代表并非该modulemap中一部分的附加的‘decoy’structures。Clique,biclique和modules的size都是在10-20的范围内随机选择的,服从均匀分布和点集合不相交(uniformdistribution and disjoint node sets)。称生成的边集合分别为EH*和EG*。最后,我们通过引入随机噪声来改良这些图:G和H中的每一条边都以概率P被删除,每一个non-edge都以概率P被edge替换掉。所有这些步骤都是独立完成的。对于创建有权图,相同的步骤被使用,但是所有可能的边都被在最终的H和G中被代替:如果(u,v)在EH*或EG*,则W(u,v) is sampled from N(1,σ),否则是N(-1,σ)。我们也测试了这个方式,用1000个点的图,10到20个modules,5 decoys(cliques he bicliques)。
酵母中蛋白-蛋白相互作用和阴性基因相互作用的分析nalysisof negative genetic interactions and protein-protein interactions in yeast
PPIs和阴性GIs都是从BIOGRID中下载的。这些网络被用于检测蛋白复合物之间的上位关系(epistaticrelations among protein complexes)。PPI网络被作为H使用,GI网络被当做G使用(详情看附表S1)。
Analysis of DNA damage-specificgenetic interactions data
我们使用(21)中的数据,在418基因中所有成对的GIs都被检测,也使用(31)中的数据,其中也检测了2022个基因(gene)和55个querygene之间的GIs。‘DNAdamage-specific positive GI’ 的定义中,在untreatedcells中S<0,在treatedcells中S>0.5而且对于不同的GI的P-value<0.01。这分析使用(21)中数据产生了840个interactions,使用(31)中数据产生了1677interactions。此外,我们定义阳性GI为稳定(‘stable’),如果它满足:在untreatedcells 和DNA damage cells中都有S>1.5。这分析使用(21)中数据产生了491个interactions,使用(31)中数据产生了3139interactions。这是因为大多数GIs的实验装置都是不能直接相互比较的。
计算微分相关评分Calculatingdifferential correlation scores
给定一个包含一些主题的基因表达谱的训练数据,我们使用(24)中的统计方法(statisticalmethod)来计算每一对基因的CC(consistentcorrelation)和DC score。首先,DCscore 是使用样品的真实labels 来计算的。然后,score将会被转化成为对数似然比(LLR,log-likelihoodratio) score,通过对比最初的DCscore和用相同数据但是带有随机打乱的标签计算得到的score来得到。因此,阳性LLRscore使用DC中的显著变化来标记基因对。在真实DC变化中优先概率被设置到只有相关变化不少于0.4的时候,相关变化(correlationchange)才会有一个阳性的LLR score。相比于(24),这个方法保证了一种相似的但是有点更加严格的acceptance threshold。See Supplementary Text foradditional information。
GO and microRNAenrichment analysis
我们使用TANGO(32)用于基因本体论分子功能(geneontology molecular function )和modules的生化过程改进分析(biologicalprocess enrichment analysis),FAME(33)用于microRNAenrichment analysis。两个工具(TANGO和FAME)都可以作为EXPANDER软件的部分可以使用。当一组modules被分析时,我们使用错误发现率FDR,以参数q=0.05进行多路测试(multipletesting)的矫正。改进分析的背景设置是只有一部分基因在网络network中,并不是生物体所有基因都在。以防网络中的GOterms 的比例过高,这一步过滤减少了一些偏见bias。
Networkvisualization
使用Cytoscape来完成。
实用性Availability
在网站http://acgt.cs.tau.ac.il/modmap/上可以免费得到工具ModMap,只需要一行命令就可以运行。
结果 results
模拟结果Simulations
我们首先在综合图(syntheticgraph)H和G上测试了几种不同的算法。从一个完美的perfectmodule map 开始,我们首先加入H中的cliques和G中的bicliques来表示附加的结构,附加结构不是map的一部分,并且引入随机噪声到edges。为了生成稀疏图和密集图,我们测试了一个大范围的噪声参数σ和p分别在有权和无权的情况下。呈现出的结果是每一个map都拥有500nodes和6 modules的图(for graphs with500 nodes and six modules per map)。我们同样检测了更大的带有相似结果的图(详情见S3和S4)。
对于每一个p和σ,我们测试了10组initiator and improver,使用10组随机数据集。我们使用Jaccard coefficient 来评测了产生的解决方案的质量,通过比较reportedmodules 和 the known modules。无权和有权modules的结果分别显示在图片1A和1B中。只有四种在每一次模拟中平均效果最好的算法才被展示出来。附表S2 包含了所有组合的结果。Local improvement 算法甚至在噪声更少的数据上都没能达到perfect scores。相反,跟在global improver 后的MBC-DICER 和 DICER5 算法在没有统计噪声干扰时,可以达到perfect Jaccard scores。MBC-DICER算法具有很好的鲁棒性(robust),甚至当噪声水平高达无权module中p=0.15和有权module中σ=1.2。一个在噪声水平达到p=0.15时,在带有1000 nodes和10 modules的无权图上,各种算法结果的比较在图1C中有所显示。虽然图更大,但是效果依然很好。使用改进阶段(improvers)比只使用initiator solution 更加有效,尤其是DICER算法的变体。带有全局improver的MBC-DICER算法达到了最好的效果(0.87)。有趣的是,local improver对于其他所有算法都比global improver更好(比如,0.71 versus 0.59 for DICER5 )。这可能是因为MBC-DICER的initiator阶段发现了robust fully connected modules,这在高噪声水平为global improver提供了更好的起始点(starting point)。经过不同的k value测试,我们最终选定了5,确定了DICER5的算法(更多看S4)。另外,我们比较了使用global improver 的超几何测试(hypergeometric Test)和 Wilcoxon rank-sum test,这都是先前研究中使用的算法。我们的结果显示两个方法效果几乎一致,但是hypergeometric Test有更快的速度(更多看附加文本)。总的来说,结果显示带有global improver 的MBC-DICER算法在有权和无权图上都可以产生最好结果。我们称这个算法为ModMap,并且从现在开始使用它作为默认的算法(the algorithm ofchoice)。
酵母但白蛋白相互作用和阴性基因相互作用数据 Yeastprotein-protein interaction and negative genetic interaction data
我们是使用来自BIOGRID (30) 的PPIs和阴性GIs来检测蛋白复合物之间的上位关系(epistatic relations)。只有那些具有两种类型相互作用的基因被使用。总的而言,网络(network)包含3979基因,45456PPIs,和76237 阴性GIs(相互作用列在附表S1)。这个数量的基因和edges比先前研究的要更大。比如,(22)包含1460 基因,(17)包含743 基因。因此,我们的网络具有提供给酵母相互作用的强大潜能,并且考虑到了不同算法的综合性能测试(our networks have the potential to provide a broader overview ofthe yeast interactome and allow for a comprehensive performance testing of thedifferent algorithms)。
就像之间研究中完成的那些,我们使用它们的统计数据和模块的功能特征化(17,22)来评测解决方案(we evaluated solutions by theirstatistics and the functional characterization of the modules)。统计数据主要包括modules的数量、covered genes的数量和module的最大size。我们使用TANGO(34)来评测模块功能(module functional)的改进,并且记录发现的GO terms的数量,提升的module的比例和module map中links的两个module都被提升的比例(带有相同或不同的function),称之为提升链接(“enrichedlinks”)。Enriched links代表了已知的生物terms之间的紧密GIs(Enriched links represent dense GIs among known biological terms.)。
所有的解决方案的统计数据在附表S3中显示。Global 比 local improver 具有明显的优越性。与global improver 形成对照,我们记录了至少100 modules,涵盖了800-1000 genes,local improver 只找到了 2-28 modules,涵盖了15-192 genes。除了DICER算法,所有solutions的结果都是相似的而且高质量的。ModMap是最好的,按照的是enriched modules 是 87%, enriched links是80%。总的而言,在综合功能详尽和质量上(combining functionalcomprehensiveness and quality),ModMap的map是最好。我们也比较了基于Collins(37)数据的其他有权算法来进行GIs data analysis(22,36)。详情看附加TXT。我们的实验结果显示ModMap算法可以产生高质量的maps,相比现在的有权方法有很多的提升。
图2显示了一份使用ModMap建立的map,其中links局限于P < 10E-50(详情看附表S4-S6)。其中每一个节点代表了一个module,边表示map links。显示的map中所有的module都被提升了0.05FDR,每个module至少有一个GO term。节点的标签显示了veil最大程度提升的term。三个重要的额hub被标记成绿色: Rpd3L complex (14 genes, P=4.35E-38),Swr1 complex (13 genes, P=1.08E-35) 和the mediator complex(17 genes, P=4.89E-43)。Rpd3L 和 Swr1复合物是chromatin (核染色质)相关的,并且之间都被在一个基于基因的研究中被注释为(annotated)GIs的hubs。Bandyopadhyayet al. (21) 发现了一些相同的links,无论如何,这里的module注释是手工完成的,然而我们的分析则完全是自动完成的,并产生一个更大的map。此外,我们的map通过显示那3个hubs是链接在一起的和为Rpd3L提供额外的links,来扩展到先前的观察数据上(our map extends on the previousobservations)。在图3,我们集中于map(P<1E-70)中的那三个最重要的links。图3A显示了Rpd3L和Set3复合物之间、Rpd3L和Swr1复合物之间的connections。Rpd3L和Set3复合物都是histonedeacetilases(组蛋白乙酰),并且他们之间的阴性GI在(20)中有记录。Rpd3L复合物被分裂成两个互不相交的modules,然而在我们的map中却是作为一个单个module被检测出来,包含该复合物所有的14个genes。图3B 显示了proteasomecomplex (39) 中的两个已经确认的次级单元之间的connections。这个例子展示了PPIs和GIs综合分析(jointanalysis)如何正确检测到核心功能单元(corefunctional subunits),甚至他们被很多PPIs链接在一起的时候。
酵母中DNA损伤应答网络的分析 Analysis of DNA damageresponse networks in yeast
本文先前描述的module map是通过分析已知的阴性GIs的整个集合来得到的。最近的研究已经超越使用静态分析(gone beyondstatic analysis)检测 对应DNA damage (21,31)的GIs网络中的变化。在这些研究中,GIs被用于测定未经治疗的细胞和DNA损伤造成的MMS(methyl methanesulfonate,甲磺酸甲酯)所带来的扰动。我们结合两种这样的数据集合(21,31)来检测‘DNAdamage-specific positive GIs’,比如,在以治疗的细胞发现的,但是在未治疗的细胞中未发现的差异阳性GIs(differential positive GIs)(详情看 ‘Materials and Methods’ section)。阴性GIs尤其发现(typically observed)于并行工作的genes之间,比如,包含在在两个互补(或互为备份)的复合物或pathways(compensatory complexes or pathways that backup each other)之间的基因,因此,其中一个受到损伤时另一个可以提供一定的减轻作用(buffered)。阳性GIs在来自相似complex或 pathway之间的genes的更有可能被检测到,这里的genes中,大部分表型效应(phenotypic effect)已经在每个single-knockout中找到。因此,DNA特定损伤阳性GIs(DNA damage-specific positive GIs)被希望用与于表示network针对MMS做出应答时的变化,来揭示pathway内或之间或complexes之间working in series 的DNA特定损伤相互作用。总计,1078 genes包含在两个研究中,他们之间带有2227 个DNA damage-specific positive GIs(附表S7)。共有6671PPIs在基因集合之内(within that gene set)。
我们将PPI网络看做H,将DNA damage-specific positive GIs 网络看做G,应用ModMap 算法。因为这些网络在先前的研究中都更小,所以我们设立一个module 的最小size为3.module的小size也能够影响针对links可得到的P-value值。在此,考虑到在改进阶段(improvement steps)该算法完成的所有统计学实验,如果在进行Bonferonnicorrection之后,一对modules的P-value<0.05,则这对modules被定义为linked。
生成的module map包含12 modules,其中共含有78 genes,在modules之间共有17 links。Module的size分布在3-15之间。附表S8-S10中有对此的详细描述。图4A中显示了一个使用GO terms来明显提升(significantly enriched with GO terms)的modules的 map。这个map的hub是一个enriched with DNA repair genes的module,链接到了6个涵盖了各式各样的functions的modules上。在图4B中,我们专注于DNA 修复相关模块(DNArepair-related module)和连接到他(hub)的三个modules上。DNA repair module包含了四个基因:RAD5,RAD18,HPR5,和UBC13.有趣的是,尽管UBC13已知和其他三个基因之间有物理相互作用(physically interact with。。。),阳性GIsquery依然在实验中一贯地保持稳定地(see ‘Materials and Methods’section)与其他三个基因相连,为这四个基因包含在一个共同的过程中提供更多的证据。RAD5,RAD18,和UBC13 三个基因已知包含在post-replication repair (41–43),HPR5 基因包含在checkpoint recovery (44,45)。
图中的DNA repair hubmodule链接到了一个与DNA damage应答相关的module上(amodule associated with response to DNA damage)。它包含5个基因:CTF4, ESC4, MMS1, MMS22 and Rt101。其中最后四个基因也是cullin-RING ubiquitin ligase complex (GO:0031461)的一部分。最后三个基因被用于引导形成复合物来稳定replication stress (46,47)过程中的复制体。通过它的GO注释,我们知道,CTF4基因和DNA修复和DNA复制启动有关(DNA repair and DNA replicationinitiation)。Links表示这个复合物可能与DNArepair module一起作用在损伤复制叉的复制(coping with damagedreplication forks)。有趣的是,这两个MMS基因都是最初在MMS敏感性测试中发现的,但是却在双链修复(double-stranded repair)中并不是需要的(47)。RAD52 module(RAD51, RAD52 and RAD59)与双链DNA damage repair 有关,并且被链接到DNA damage repairmodule 和 DNA damage response module中,这暗示着这些modules是一样地作用在DNAdamage来处理damaged replication forks 和double-stranded DNA breaks。第四个链接在一起的module,包含SuperKiller (SKI) complex (SKI2, SKI5 and SKI7)的三个基因。这些基因涉及到cytoplasmatic exosome 中的3-5 RNA degradation(49,50)。我们的分析显示这个复合物(SuperKiller (SKI)complex)可能也涉及到了DNA损伤应答。先前的研究表明,RNAdegradation cytoplasmatic genes 也可能在从他们的细胞质活动分开DNA损伤应答中起到了作用(play a role in DNA damage response separately from theircytoplasmatic activity)(51,52)。研究中表明的DNA degradation genes在DNA damage response起到的作用,包含DNA 稳定性和功能相关的端粒稳定性(DNAstability and telomere stability related functionality)(51),调节double-stranded breaks 中多蛋白复合物的组成(52),相关DNA损伤的特定mRNA 退化(specific mRNA degradation on DNA damage)(53)。因此,我们的发现和之前的研究结果是一致的,并且进一步强调了SKI complex 在DNA damageresponse 中的角色作用。
人类共表达和差分相关网络的分析 Analysis of humanco-expression and differential correlation networks
我们应用ModMap算法到NSCLC的个案对照基因表达数据(case-control geneexpression data)上,来揭示高关联度的gene modules之间的DC。这部分的贡献是两层的(two fold)。第一,在交叉验证实验中,我们揭示了gene modules之间的DC是可再生的。第二,我们分析被ModMap发现的gene modules之间的DC模式的map。
给定一个来自cases和controls 的基因共表达谱(gene expression profiles)数据集合,我们使用(24)中的方式来计算每一个基因对的两个评分:CC score,如果基因对在显型之间是一直相关的(consistently correlated across phenotypes),则CCscore是正值。DC score,如果cases和controls之间的correlation difference意外地(by chance)比期望值更高,则DC score为正值。这些scores都作为边的权重各自在network H和G中使用,在两个网络中一个module map可以学习得到。这些方法使用交叉验证的方式进行评估:给定一个基于一个profiles集合(the training set)和一个不相交的samples集合(the test set)建立的module map,predicted map的质量通过对比使用Wilcoxon rank-sum test 得到的links 和 non-links的DC ,其中零假设(nullhypothesis)是links和non-links之间的DC是没有区别的。这个测量方法是无参数的,而且反应所有的DC 变化。
我们使用2-fold-cross-validation来测试了几种算法的变型。Local improver方法得到的maps得到了较低的P-value值,但是却遭受了低覆盖范围(low coverage)。例如,对于MBC-DICER initiator,local improver实现了P-value值为4.43E-4,但是map覆盖了仅仅197 genes。相反地,当应用ModMap(MBC-DICER with the global improver)算法时,覆盖了1289个genes,P-value值为1.54E-10。附加文本中包含了更多的结果,测试了几种不同的global improver的参数,并且包含了Alzheimer’s disease(54)的测试,这些测试得出了相似的交叉验证结果。完整的肺癌数据在附表S11,Alzheimer’s disease数据在附表S12。合在一起,ModMap算法产生了更大的maps,并且他们在测试独立数据集合时,更加具有robust。
接下来,我们分析了通过在所有的NSCLC数据样本上运行ModMap算法得到的module map。它覆盖了1921 genes,包含在被405 links链接在一起的 76 modules(详情看附表S13,S14)。为了集中在modules之间correlation的strong change,我们比较了map中每一个link的DC,和通过200次重复计算相同size的随机的基因集合之间的DC,还有通过计算真实链接和像(24)中提到的最好的随机链接的fold-change。Link fold-change scores在附表S14中有给出。总计,150 links具有fold-change ≥ 1.5,最大的5个links甚至超过了2.3。这证明了linkedmodules之间的DC意外地比预计的要强的多。我们也使用pathwayenrichment analysis 和 microRNA enrichment analysis (seeSupplementary Table S15 for details)分析了最大的几个links的modules。其中一个links链接到了两个和免疫反应活化作用相关的modules(Oneof the links connected two modules related to immune response activation)。The linked modules are shown in Figure 5.在图5A中,我们观察了许多在control class中在modules((gene pairs with r > 0.4))之间高共表达边(high co-expressionedges)。Module 11 包含了B细胞受体信号通路基因(6 genes, P=3.1E-8)。Module 12 包含T细胞受体信号通路基因(4 genes, P=1.37E-4)(Module 12 is enriched with T-cell receptor signaling pathway genes(4genes, P=1.37E-4))。图5B显示了以上这10个genes的GeneMANIA 分析(7,55),这确认了他们之间是通过几种类型的相互作用(interactions)链接在一起的。图5C显示了在NSCLC class中几种相似modules之间的co-expression。在每一个modules内部,都有一个高水平的co-expression,但是modules之间的co-expression是彻底摧毁的,这显示了在NSCLC中不同的immune responses之间的co-regulation是丢失的。最后,module 11包含了很多microRNA34-a,b,c家族的目标(图5A中红点),它的成员是有注释的,因为NSCLC按照疾病数据库whose members are annotated ascausal to NSCLC according to the mir-2-disease database (56)。合在一起来讲,这些结果显示了我们的分析方法在不适用任何先验知识情况下来检测NSCLC-related functional modules的能力。
讨论DISCUSSION
在这篇文章中,我们展示了一个joint 分析方法来处理两个gene networks,它们各自表示了genes之间的一种类型的omic relation。这个方法将gene sets看做modules,并找到它们之间的复杂结构关系,并且最后综合成一个module map。在第一个network中,modules相当于相互作用基因集合,module map中的links相当于第二个network中的相互作用的modules(Modules correspond to interacting gene sets in the first network,and links in the module map correspond to interacting modules in the second)。Map 是同时基于两个网络来建造的,因此可以capture and reveal 在单独分析各个数据类型中所不能识别的structures。我们的新型算法重新获得了(recovered)在仿真数据中的植物map structure,甚至数据中当噪声水平相当高的时候。我们在三个生物应用中检测了我们的算法:(i) yeast PPIs and negative GIs, (ii) yeast PPIs and DNAdamage-specific positive GIs and (iii) DC analysis of human disease expressionprofiles。在所有的情况下,我们的maps的一些部分是被先前的生物知识所支持的,然而,剩下的部分则发现了新的structure或暗示了一些新的生物发现(reveal novel structureand suggest new biological findings)。Module mapparadigm 可以大体上用于基于相同节点的两个不同类型的networks。
关于yeast PPI andnegative GI data 分析,我们建立了一个大的map来描述复合物之间的上位关系(epistatic relations)。我们的发现不仅与之间的研究一致,显示了chromatinmodification-related complexes 之间相互关系的complex map,还提供了与其他功能相关的相互作用(interactions with otherfunctions),例如,蛋白修饰相关复合物(proteinmodification-related complexes)。yeast PPIs and DNAdamage-specific positive GIs的数据分析产生了一个稍小的map,其中包含了一个作为central hub 的DNA repair module。这个module的相互作用显示了(suggest)几个mechanisms在MMS响应中同时暴露出来,包括double strand repair, damagedreplication fork repair and exosome complex activity。在基于human NSCLC blood expression profiles数据建立的map中,modules 代表in cases and in healthy controls高度共表达(highly co-expressed)的gene sets,然而map links相当于在MSCLC patients中的co-expression network 内的特殊重新排线(specificrewiring)。特别地,我们发现two modules enriched with immuneactivation genes 在NSCLC patient的相互关系时有一个暴跌(manifesting a sharp drop),这暗示了在T细胞和B细胞富集模块(enriched modules)之间的协作减弱(diminished coordination)。
Module map的概念可以看做是聚类和双聚类的一种高层次的综合运用(higher levelcombination of clustering and biclustering)。每一个这样的问题都已经被广泛地研究,并且被成功地用于大量的单类型genomic and proteomic研究中(1,57-68)。通过演示两种不同类型数据之间的joint analysis,我们允许使用两个网络中各自的一些松弛的目标函数,来获取一个全局的更加具体的结构(genomic and proteomic)。因此,当一种类型数据的聚类或双聚类失败时,我们的新型分析方法依然可以获得结果。在聚类和双聚类中的一个难点是module size 必须足够大,使得可以获取高度显著的sets(highly significant sets)。就像我们的分析示范的那样,modulemap 算法的额外能量可以检测到比较小的精确的分组,这些是超出之前的方法的能力范畴的(identify relativelysmall precise groups that are beyond the detection ability of those priormethods)。
迄今为止,只有很少研究者联系到了modulemap problem,而且他们之间的大部分都关注于yeast PPI and GI networks的joint analysis。Ulitksy et al. (17)和 Bandyopadhyay et al. (69)改进了用于寻找map的聚类方法,簇内PPIs和GIs或链接的簇间GIs的边权重 的可能性高于一个给定背景分布的情况。Leisersonet al. (22,36)通过一个贪心增量方法,寻找GIs有权图中的局部最大的下降(local maximum cuts),产生一个链接在一起的modules对的集合。Kelley and Ideker (20)改进了一个基于图压缩的聚类方法,在此原始的GI图被压缩成为一个module map。因此,(22,36) 和 (20)都寻找链接gene modules的近似的双聚类算法。相反地,我们计算GIs 的最大 bicliques,通过直接考虑两种类型相互作用的数据来保证在初始解中获取到强烈密集相连的modules, 从而分析这些数据,并且使用我们的global improver来提升solutions 的质量。因为我们的算法是一般性的(generic),所以它不会像其他算法那样利用GI 数据中某些特定的概率学性质(22,36)。但是,我们已经显示了我们的算法在一些条件(several criteria)下的GI 数据上是优于那些利用特定概率学性质的算法和其他的现有算法的。另外,因为我们的算法并不局限于几种数据数据的类型,所以我们可以在分析中结合使用其他的许多种heterogeneous data sets (e.g. using all GIs of BioGRID)。
当我们仔细检查humanexpression profiles of disease patients and healthy controls, DC 分析被提议作为一个方式来检测gene modules,这些modules之间的相关水平在disease中被改变的(12,14,24,70)。我们先前改进了使用local 方式来检测module pairs的DICER方法(24)。在此,我们尝试通过寻找DC图中的最大bicliques和同时建立modules 的global map来超越这个方法。就像在此展示的这样,在大部分情况下,maplinks都是高度显著的(highly significant)。但是,我们也观察了一些情况:尽管module pair的DC是显著的,modules的绝对关联变化也可能是轻微的。一种可能的纠正方式(possible remedy)就是给map links的high absolute DC 更多的emphasis,以便来更好地观察DC信号。另一种可能的提高措施就是使用确定的启发式函数来计算bicliques(enumerate bicliques using established heuristics [e.g. (68)])。
ModMap算法性能的一个关键因素就是最优化的目标函数。在此,我们选择最大化modules内部权重和加上module links的权重和之后的总和,并基于概率模型来指定这些权重。关于无权网络,比如,PPI和GI yeast network,我们设定edge的权重为1,non-edge的权重为-1,从而提升强烈相连的modules和links(promoting strongly connected modules and links)。这个设定产生了很好的结果,并且显示了蛋白复合物之间的功能性相互作用。通过给图中的non-edges设定不同的权重,将来的分析可以提升稀疏的modules,因此能够更好地进行complete pathways之间的相互作用的检测。