本文档包含关于如何使用SANTA包中包含的功能的说明。SANTA建立在Ripley的k函数基础上,这是空间统计学中一个成熟的方法,通过将其应用于网络来测量平面上点的聚类。使用这种方法,Knet
而且Knode
函数被定义。Knet
检测网络上命中的聚类,以测量网络和表型之间的关联强度。Knode
根据与命中点的距离对基因进行排序,为后续分析提供了一种自然的优先排序顶点的方法。
SANTA解决了两个互补的目标:它可以(i)使用策展基因集功能性地注释实验衍生的网络,以及(ii)使用策展网络注释实验衍生的基因集。
SANTA包还包含可用于测量顶点对之间距离的函数(GraphDiffusion
而且GraphMFPT
)和一个功能,可用于可视化Knet
函数(情节。Knet
).igraph}包用于处理网络。
这个小插图,以及SANTA包中提供的数据,旨在允许用户重现论文中显示的结果:Cornish, A.和Markowetz, F. (2014)量化分子网络的功能内容。《科学公共图书馆·计算生物学》,10:9,e1003808。为了减少运行小插图所需的时间,许多案例研究都设置为减少试验、排列或基因集的数量。它指出了案例研究与论文分析的不同之处和方式。文中进一步解释了SANTA包所使用的方法。
关联负罪感(GBA)原理指出,相互作用的基因更有可能共享生物功能。例如,在合成基因阵列(SGA)中相互作用的两个基因比不相互作用的两个基因更有可能编码涉及同一途径或复合物的产物。Knet
使用这一原则来衡量网络和表型之间的关联强度。通过测量基因集在生物网络中的聚类,就有可能确定最能解释该基因集的机制和过程的网络。这个原则也被Knode
优先考虑后续研究的基因。如果一组基因被确定与特定的细胞功能相关,那么GBA原则表明,与这组基因相互作用的基因也可能参与该功能。
Knet
函数的Knet
Function通过量化网络上高权重顶点的聚类来衡量网络和表型之间的关联强度。高权重顶点代表那些基因是一个基因集的一部分或与表型联系最紧密的基因。GBA原则意味着聚类越强,网络就越能解释产生和维持表型的机制。
该函数基于Ripley的k统计量,这是一种空间统计工具,用于测量平面上点的聚类。雷普利k统计量的工作原理是计算每个点周围半径圈内包含的点的数量。增加\(s\)的值,并测量结果曲线下的面积(AUK)。如果点聚集在一起,那么它们将更快地被圆圈吞没,函数将更快地增加,AUK将更大。通过将平面上的欧几里得距离与网络上的距离测量值交换,这个统计值可以适用于网络。我们称这个函数为Knet
函数。
不同于Ripley的k统计量Knet
函数可以度量顶点集合的聚类或顶点权重的连续分布。如果使用生物网络,那么高权重将表明某个基因与正在研究的表型密切相关。在某些情况下,可能不是所有顶点的权重都可用。缺少权重的顶点可以包含在网络中。但是,应赋予它们NA}的权重,以确保它们的分布不影响显著性分析。
聚类的重要性通过排列网络中的顶点权重来量化。缺少权重的顶点被排除在外。下面的面积Knet
函数的观测集,然后比较的面积下Knet
函数,用z检验来计算p值。产生的p值表示顶点权重聚类至少与观察到的一样强的概率,假设顶点权重在网络中随机分布。
的Knet
函数是:
\ [K ^{净}(s) = \压裂{2}{酒吧(\ p {} n) ^ {2}} \ sum_{我}p_{我}\ sum_ {j} (p_ {j} - {p} \酒吧)~ \ textbf{我}(d ^ {g} (i, j) \ leq s) \]
其中\(p_i\)是在顶点\(i\)上观察到的表型,\(\textbf{i}(d^{g}(i,j)\leq s)\)是一个恒等函数,如果顶点\(i\)与\(j\)之间的距离小于\(s\),则等于\(1\),否则等于0,\(\bar{p} = \frac{1}{n}\sum_{i=1}^n p_i\)。
Knode
函数取内和Knet
方程可以通过它们与高权重顶点的距离来对顶点进行排序。我们称这个函数为Knode
函数。这一功能为后续研究提供了一种自然的方法来确定基因的优先级。
的Knode
函数是:
\ [K_{我}^{节点}(s) = \压裂{2}{p} \ sum_ {j} (p_ {j} - {p} \酒吧)~ \ textbf{我}(d ^ {g} (i, j) \ leq s) \]
SANTA包含三种测量顶点对之间距离的方法。每种方法都包含了网络的不同方面。
Knet
成功识别模拟网络上的聚类在本节中,我们将使用模拟网络来演示Knet
-function成功地量化了网络命中的聚类。
为了演示如何Knet
函数可以用来度量网络上的点击的聚类,我们将模拟一个网络,将不同聚类强度的点击的聚类应用到网络上,并使用Knet
函数来量化聚类的强度。
我们将使用偏好依附的Barabasi-Albert模型来模拟网络,因为已经观察到这种方法产生的网络具有与许多现实世界网络相似的属性。由500个顶点组成的网络既大到足以包含不同强度的命中集群,又小到足以快速运行大量试验。在Barabasi-Albert模型中,顶点是逐步添加的,并连接到先前添加的顶点。这是在基于每个顶点程度的概率分布下完成的。如果\(幂=1\),则概率分布与顶点度之间是线性关系。我们将使用这个值,因为它创建了一个具有多个高级枢纽的网络,就像许多生物网络一样。每个新顶点都连接到以前添加的\(m\)个顶点。为了创建一个中等密度的网络,我们将设置\(m=1\)。
#生成模拟网络要求(SANTA)
##加载所需的包:SANTA
##加载所需的包:图
## ##附件:“igraph”
以下对象从'package:stats'中屏蔽:## ##分解,频谱
下面的对象从'package:base'屏蔽:## ## union
要求(题图)set.seed(1) #用于重现性g <- barabasi。游戏(n=500,幂=1,m=1,定向=F)
通过首先随机选择一个种子顶点来模拟一个命中群集。然后,网络中的所有其他顶点根据它们到种子顶点的距离进行排序。距离可以使用前面描述的任何距离度量(最短路径,扩散核和MFPT)来测量。在本例中,我们将使用最短路径度量,因为它是运行最快的度量。但是,使用其他方法运行分析将产生相同的结果。
#测量g dist.method中顶点对之间的距离<- "short。D <- DistGraph(g, dist.method=dist。方法详细= F)
的Knet
函数计算\(g\)中每对顶点之间的距离。但是,为了减少运行函数所需的时间,可以在函数外部计算这些距离。的Knet
函数不能在原始连续距离上运行。因此,有必要将顶点对放入离散的容器中。这是使用BinGraph}函数完成的。
#将距离放入离散的bin B <- BinGraph(D, verbose=F)
接下来,我们将生成命中的聚类,并使用Knet
函数。最接近种子顶点的\(s\)顶点形成样本集。从这个集合中,随机选择5个顶点形成命中集。与网络中的顶点数量相比,命中集的大小需要较小。增加\(s\)的值会降低聚类的强度。我们将测试5个不同的\(s\)值,以模拟不同的聚类强度。等于网络中顶点总数的\(s\)值将生成一个包含从整个网络中随机选择的顶点的命中集。对于\(s\)的每一个值,我们将进行10次试验。在相关论文中,进行了1000次试验。这里我们也将使用100种排列,而不是1000种,以减少所需的时间。
集群。Size <- 5 s.to.use <- c(10, 20, 50, 100, 500) n.trials <- 10 pvalues <- array(0, dim=c(n。审判,长度(s.to.use)), dimnames=list(NULL, as.character(s.to.use))) #为s的每个值运行审判for (s in s.to.use) {for (i in 1:n.trials){#生成命中集种子。顶点<- sample(vcount(g), 1) #标识种子样本。set <- order(D[种子。顶点,])[1:s]命中。Set <- sample(sample。Set, cluster.size) #测量关联强度g <- Set .vertex。attribute(g, name="hits", value=as.numeric(1:vcount(g) %in% hit.set)) pvalues[i, as.character(s)] <- Knet(g, nperm=100, dist.method=dist. numeric)方法,顶点。attr="hits", B=B, verbose=F)$pval}}
图1为模拟研究结果。的Knet
函数正确地识别使用更小的\(s\)值创建的命中集,使聚类更加显著。
Knet
对密实度接下来,我们将展示Knet
将其与一个更简单的度量进行比较。
Glaab等人将紧凑度分数定义为网络中命中对之间的平均距离。通过排列属于命中集的顶点,并将观察到的紧凑性得分与排列集的紧凑性得分进行比较,可以产生一个经验p值。这很像Knet
。但是,与Knet
~时,自适应紧凑度评分不能应用于顶点权重的连续分布,并且在测量关联的显著性时不考虑顶点度。这个玩具示例演示了这一点。
首先,我们创建网络。这个网络中的顶点有一个程度范围,允许我们演示之间的差异Knet
而且密实度
。我们将应用2组顶点权重,每一组都包含两个相互作用的命中。
#创建网络n.nodes <- 12 edges <- c(1,2,1,3,1,4,2,3,4,3,4,1,5,5,6,2,7,7,8,4,9,9,10,3,11,11,12) weights1 <- weights2 <- rep(0, n.nodes) weights1[c(1,2)] <- 1 weights2[c(5,6)] <- 1 g <- graph.empty(n。节点,directed=F) g <- add.edges(g, edges) g <- set.vertex。属性(g, "weights1", value=weights1) g <- set.vertex。属性(g, "weights2", value=weights2)
图2显示了网络上每组命中点的位置。第1组和第2组中击球的距离是相同的。然而,直观地,我们期望集合2中命中的聚类更显著,因为命中的交互伙伴更少,因此不太可能被观察到偶然交互。
现在我们将应用Knet
而且密实度
以每组命中值来衡量聚类显著性。我们将运行带有100个排列的函数。增加排列的数量可以提高函数的准确性。然而,这也增加了运行时间。使用了1000种排列来创建论文中显示的结果。将使用最短路径距离度量。然而,其他方法也会产生同样的结果。
# set 1 Knet(g, nperm=100,顶点。attr =“weights1 verbose = F) pval美元
## [1] 0.3986885
紧凑性(g, nperm=100,顶点。attr =“weights1 verbose = F) pval美元
## [1] 0.08767958
# set 2 Knet(g, nperm=100,顶点。attr =“weights2 verbose = F) pval美元
## [1] 0.002557689
紧凑性(g, nperm=100,顶点。attr =“weights2 verbose = F) pval美元
## [1] 0.09822864
的密实度
函数在度量命中聚类的重要性时不考虑顶点度。因此,它为每个命中集输出一个相似的p值。Knet
合并顶点的全局分布,因此将集合2比集合1更明显地识别为集群。因此,最好使用它Knet
而不是密实度
当顶点的程度在网络中变化时。
Knet
以证明GI剖面中的相关性可以产生功能更丰富的网络在本节中,我们将演示如何Knet
函数可以用来确定构建交互网络的最有效方式。
有证据表明,蛋白质功能与遗传相互作用(GI)之间的整体相似性比与个体相互作用更密切相关。为了测试这一点,我们将使用Knet
函数来量化GO术语和两个酿酒酵母网络之间的关联强度:一个是原始相互作用的网络,一个是由相互作用剖面中的相关性创建的网络。这两个都将使用Costanzo等人产生的数据来创建。如果GO术语与网络的关联更强,则这表明网络在功能上提供的信息更多。
首先,我们将加载两个预处理网络,形式为igraph
对象的圣诞老人
包中。圣诞老人
要求所有的网络形式都是igraph
对象。用于创建这些网络的数据可从http://drygin.ccbr.utoronto.ca/~costanzo2009/使用的方法在相关的SANTA论文中进行了描述。
Knet
在量化高权重顶点的聚类时,也能够合并不同的边缘距离。边缘距离可以代表不同的生物学特性,例如两个基因产物之间物理相互作用的强度,或者在我们的两个网络的情况下,遗传相互作用的强度。边缘距离越小,相互作用越强或越显著。在创建网络时,我们通过取分数的\(-log10\)将相关系数转换为边缘距离。这些距离被存储在edge属性下距离
。
#加载金图对象数据(g.costanzo.raw)数据(g.costanzo.cor)网络<- list(raw=g.costanzo. cor)Raw, cor=g.c stanzo.cor) network.names <- names(网络)网络。基因<- V(网络$raw)$name
##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…
#跨网络的基因相同
基因和GO术语之间的关联将从基因本体论(GO)项目中获得org.Sc.sgd.db
包中。为了减少运行此分析所需的时间,我们将只使用5个GO术语。然而,像相关论文中所做的那样,使用更多的GO术语来运行这种分析很容易。
#从org.Sc.sgd.db包库中获取GO术语关联
##加载所需的包:AnnotationDbi
##加载所需的包:stats4
##加载所需的包:BiocGenerics
##加载所需的包:并行
## ##附加包:“BiocGenerics”
以下对象将从'package:parallel'中屏蔽:## ## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, ## clusterExport, clusterMap, parApply, parCapply, parApply, ## parapplylb, parRapply, parSapply, parSapplyLB
以下对象从'package:igraph'中屏蔽:## ## normalize, union
以下对象从'package:stats'中屏蔽:## ## IQR, mad, xtabs
以下对象从'package:base'中屏蔽:## ## Filter, Find, Map, Position, Reduce, anyduplication, append, ## as.data.frame, cbind, colnames, do。调用,复制,eval, ## evalq, get, grep, grepl, intersect, is。Unsorted, lapply, ## length, mapply, match, mget, order, paste, pmax, pmax.int, ## pmin, pmin.int, rank, rbind, rownames, sapply, setdiff, sort, ## table, tapply, union, unique, unsplit, which, which。Max, ## which.min
##加载所需的包:Biobase
##欢迎访问Bioconductor ## ##小插图包含介绍性材料;查看## 'browseVignettes()'。要引用Bioconductor,请参见##“citation(“Biobase”)”,以及软件包的“citation(“pkgname”)”。
##加载所需的包:IRanges
##加载所需的包:S4Vectors
## ##附加包:“S4Vectors”
下面的对象从'package:igraph'中屏蔽:## ## compare
以下对象从'package:base'中屏蔽:## ## colMeans, colsum, expand。grid, rowMeans, rowsum
## ##附件:“IRanges”
下面的对象从'package:igraph'中屏蔽:## ## simplify
# #
xx <- as.list(org.Sc.sgdGO2ALLORFS) go。terms <- c("GO:0000082", "GO:0003682", "GO:0007265", "GO:0040008", "GO:0090329") #将GO术语应用到网络for (name in network.names) {for (GO .names)在go.terms中的术语){网络[[name]] <- set.vertex. net。属性(networks[[name]], name=go.)项,价值= as.numeric(网络。% xx[[go.term]])))}}
##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…
##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它… ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly... ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly...
我们现在将应用Knet
函数到每个网络上的每个GO项,以量化关联的强度。如果一个基因与GO术语相关,它的权重将为1。如果它没有关联,它将被赋予0的权重。我们现在将这些分数存储在两个图的顶点属性下抽样
。由于需要时间,这段代码的结果被注释。
# results <- list() # for (network.names中的名称){# results[[name]] <- Knet(networks[[name]]], nperm=1000, # vertex.attr=go. name . #而言,优势。attr="distance", verbose=F) # results[[name]] <- sapply(results[[name]], # function(res) res$pval) #}
此代码可用于绘制结果。每个点代表一个GO项。p值的\(-log_{10}\)已被绘制出来。
# p.values <- array(unlist(results), dim=c(length(go.terms), # length(network.names)), dimnames=list(go. terms)。# network.names)) # p.values。Ml10 <- -log10(p.values) #轴。范围<- c(0, max(p.values.ml10)) # plot(p.values. ml10)Ml10 [, "raw"], p.values。Ml10 [, "cor"], asp=1, # xlim=轴。范围,ylim =轴。range, bty="l", # xlab="原始GI网络p值的-log10 ", # ylab="相关网络p值的-log10 ", # main="") # abline(0,1, col="red")
图3比较了每个GO术语与每个网络的关联强度。通过增加测试的GO术语的数量,可以证明在GI配置文件中使用相关性而不是原始交互可以产生功能上信息更丰富的网络。
Knet
研究酵母相互作用网络在mms处理下的功能重组在本节中,我们将演示如何Knet
-function可以通过比较2个网络之间基因集的关联强度来研究细胞的功能重布线。
Bandyopadhyay等人在正常的实验室条件下绘制了酵母的遗传相互作用(GI)网络,以及暴露于dna破坏剂甲基甲烷磺酸(MMS)的酵母。为了研究这些功能网络之间的差异,我们将使用Knet
-功能,以量化一组参与响应DNA损伤的基因的关联强度。虽然我们不能从这个例子中了解到很多东西,但也可以测试参与其他过程的基因集,以确定参与DNA损伤反应的其他过程。
现在我们将加载两个交互网络。这些网络是使用GI配置文件中的相关性创建的,如前一节所述。这是因为,如前所述,在GI轮廓中使用相关性作为边可以得到功能上信息更丰富的网络。
#加载金图对象data(g.b idyopadhyay .treated) data(g.b idyopadhyay .treated) networks <- list(g.b idyopadhyay .treated)治疗,治疗= g.bandyopadhyay。network.names <- names(网络)
我们将开展基因本体论(GO)项目,以识别那些参与响应DNA损伤的基因。通过改变下面的代码,可以测量不同功能涉及的基因集的关联强度。在相关论文中,我们测量了194个GO术语与每个网络的关联强度。
#获取GO术语关联库(org.Sc.sgd.db) xx <- as.list(org.Sc.sgdGO2ALLORFS) #更改为使用关联的替代GO术语。genes <- xx[["GO:0006974"]] associations <- sapply(网络,函数(g) as.numeric(V(g)$name %in% associated.genes), simplify=F)
##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它… ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly...
- sapply(network.names, function(name))attribute(networks[[name]], " rds ", value=associations[[name]]), simplify=F)
##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它… ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly...
我们现在将应用Knet
函数到每个网络上的GO项,以测量关联的强度。排列的次数越多,p值的准确性就越高。
我们也会运行Knet
利用最短路径法计算顶点对之间的距离。然而,正如相关论文所示,扩散
而且mfpt
距离测量也会产生同样的结果。可以通过替换dist.method
用其中一种替代距离测量方法进行论证。通过运行该分析与额外的GO项,稳健性Knet
跨越不同的距离可以看到测量。
# results <- sapply(networks, function(g) Knet(g, nperm=1000, # dist.method="最短。路径”,顶点。attr="rdds", # edge.attr="distance"), simplify=F) # p.values <- sapply(results, function(res) res$pval) # p.values
因为测试的GO项是对DNA损伤刺激的反应
,这并不奇怪,基因集更显著地与处理网络相关(由较低的p值所示)。暴露于dna损伤剂的酵母激活或上调了参与对该剂反应的途径,从而增加了dna损伤反应相关基因之间的遗传相互作用的强度。
我们现在将观察到的和排列的形象化Knet
GO项在dna损伤和未损伤网络上的曲线和AUKs。
# plot(结果$未处理)
观察到的Knet
-函数曲线(红线)和排列Knet
函数曲线分位数(黄色区域)为对DNA损伤刺激的反应
dna受损(左上)和未受损(左下)的GI网络上的基因集。观察曲线相对于排列分位数的位置表明,基因集与两个网络强烈相关。(右)观察到的Knet
-function AUK(红线)和排列Knet
-在dna受损(右上)和未受损(右下)的GI网络上的功能AUKs(灰色直方图)。观察到的和排列的AUKs之间的差异表明,该基因集与两个网络都有很强的关联,但与dna受损网络的关联最强。
Knet
研究酵母相互作用网络在紫外线损伤下的功能重组在本节中,我们将进行与前一节类似的分析。然而,在本节中,我们将使用暴露于高剂量紫外线下的酿酒酵母的地理位置图和未处理的酿酒酵母的地理位置图所创建的网络。由于本研究中映射的相互作用数量减少,因此不可能从GI中的相关性构建网络。因此,网络是由原始交互构建的。
# laod图表对象数据(g.srivas.high)数据(g.srivas.不平)网络<- list(high=g.srivas. high)高,治疗= g.srivas。network.names <- names(网络)
我们将再次只测试一个GO项(GO:0000725,重组修复)。在相关论文中,对更多的GO术语进行了测试。
#获取GO术语关联库(org.Sc.sgd.db) xx <- as.list(org.Sc.sgdGO2ALLORFS) associated。genes <- xx[["GO:0000725"]] associations <- sapply(网络,函数(g) as.numeric(V(g)$name %in% associated.genes), simplify=F)
##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它… ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly...
- sapply(network.names, function(name))attribute(networks[[name]], "dsbr", value=associations[[name]]), simplify=F)
##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它…##此图表是由旧的(er) igraph版本创建的。##在它上面调用upgrade_graph()来使用当前的igraph版本##现在我们在飞行中转换它… ## This graph was created by an old(er) igraph version. ## Call upgrade_graph() on it to use with the current igraph version ## For now we convert it on the fly...
的Knet
函数将再次用于比较GO项与2个网络的关联强度。此代码被注释以减少运行小插图所需的时间。
# p.values <- sapply(网络,函数(g) # Knet(g, nperm=1000, dist.method="最短。路径”,#顶点。Attr ="dsbr", edge.attr="distance")$pval) # p.values
正如预期的那样,GO术语与UV处理的网络关联更强。这可能是由于暴露在破坏dna的紫外线辐射下,酵母内部发生了功能性重组。
Knet
找出最能提供癌症细胞系信息的网络在本节中,我们将应用顶点权重的连续分布,描述基因参与癌细胞系维持的证据强度,以确定最能解释细胞系维持机制的网络。
生物网络可以使用广泛的数据源构建,包括共表达数据、物理相互作用数据或遗传相互作用数据。还可以通过组合数据源来创建网络。HumanNet网络由21个数据源组合而成。完整网络是使用物理交互的精选数据创建的。
在下一节中,我们将把来自癌细胞系的预处理RNAi数据作为顶点权重应用于这两个网络。GENE-E程序已用于将每个shRNA的重要性分数转换为基因的p值。这个p值包含在加载的数据集中。我们还将这些p值转换为可以应用于网络的顶点权重。通过取p值的\(-log_{10}\),我们确保最重要的基因被赋予了最高的顶点权重。
#导入和转换RNAi数据data(RNAi .cheung)张<- -log10(rnai. Cheung) rnai. Cheung [!is.finite(rnai. Cheung)] <- max(rnai. Cheung [is.finite(rnai. Cheung)]) rnai. Cheung [rnai. Cheung]张< 0]<- 0
接下来,我们将加载两个网络的交互数据并创建一个igraph
对象。对于完整的网络,没有给出边权值。HumanNet网络中的每个交互都与一个对数似能分数(LLS)相关联,描述给定所使用的数据,功能交互发生的概率。为了降低网络的密度,与一个\(LLS < 2\)的相互作用已经被删除。其余边的距离是通过将权重除以最大权重并取这个分数的\(-log_{10}\)来产生的。我们已经从每个网络中删除了所有不属于最大群集的基因。
#导入并创建完整的网络数据(edgelist.完整)g.完整<- graph.edgelist(as.matrix(edgelist.完整),directed=FALSE) #导入数据并创建HumanNet网络数据(edgelist.humannet) g. HumanNet <- graph.edgelist(as.matrix(edgelist.humannet)[,1:2], directed=FALSE) g. HumanNet <- set.edge.attribute(g. HumanNet)Humannet, "distance", value=edgelist.humannet$distance)网络<- list(完整=g.完整,Humannet =g.humannet)
接下来,我们将对网络应用顶点权重。并不是所有的网络基因都有顶点权重。
network.names <- names(网络)网络。基因<- sapply(网络,get.vertex。attribute, name="name", simplify=F) rnai.cheung.genes <- rownames(rnai.cheung) cancer <- colnames(rnai.cheung) for (cancer in cancer) {for (name in network.names) {vertex. genes <- rownames(rnai.cheung)weights <-rep(NA, vcount(networks[[name]])) names(vertex.weights) <- network. weights)基因[[名字]]常见。基因<- rnai.cheung. Genes [rnai.cheung. Genes]基因%在%网络。基因[[名字]]]vertex.weights[常见。<- rnai.cheung[普通的。基因,癌症]网络[[name]] <- set.vertex。属性(networks[[name]], cancer, value=vertex.weights)}}
我们现在将应用Knet
-函数到这些分数。由于计算距离矩阵所需的时间,这段代码已被注释掉。但是,如果希望的话,可以很容易地运行它。
# knet。res <- sapply(networks, Knet, nperm=100, # dist.method="最短。路径”,顶点。Attr =癌症,#边缘。attr="distance", simplify=F) #p。值<- sapply(knet。Res,函数(i) sapply(i, #函数(j) j$pval)
来自结肠癌和卵巢癌细胞系的RNAi数据与这两个网络密切相关。然而,它们与HumanNet网络的联系更紧密,这表明HumanNet网络更好地解释了维持这些细胞系所涉及的细胞机制。
Knode
和BioNet使用模拟数据在下一节中,我们将比较Knode
-function to BioNet。生物网络采用节点信息覆盖分子网络,识别高分子网络。的Knode
函数标识靠近大量高分节点的节点,因此也可用于标识高分子网络。然而,由于Knode
-function单独考虑每个节点,它可以用于跨多个集群识别高分节点。
我们将模拟一个包含多个高分节点集群的网络。然后我们将应用Knode
功能和BioNet,以比较性能。在相关论文中,使用2,3和4个聚类进行了1000次试验,重复了该分析。在下面的代码中,默认情况下只运行1个trial。由于需要时间,代码被注释。
BioNet工具可作为r包使用。
#库(生命网络)
由于前面解释的原因,我们将使用偏好依恋的Barabasi-Albert模型来创建网络。我们将使用SpreadHits
函数将3个10命中的集群应用到网络。在相关论文中,也对2和4个聚类进行了分析。的SpreadHits
函数首先为每个集群标识一个种子节点,它们之间的距离最小。我们将使用12的距离截断,以确保每个命中集群位于网络的不同部分。的shortest.paths
方法将用于在传播命中时测量距离。方法也可以产生相同的结果扩散
而且mfpt
距离的措施。命中群集对参数的依赖程度有多大λ
,我们将其设置为10,以确保命中有很强的聚类。
# #需要的参数# n.nodes <- 1000 # n.hits <- 10 # clusters <- 3 # # #创建网络并在3个集群中传播hits # g <- barabasi.game(n=n. hits)# g <- SpreadHits(g, h=n, m=1, direct =FALSE)命中,集群=集群,距离。Cutoff =12, # lambda=10, dist.method="最短。路径”,verbose = FALSE)
现在我们将模拟p值的分布,并将它们作为顶点属性应用到网络中。非命中点的p值将从均匀分布中选择,如果这些基因与所研究的表型之间没有关联,则可以预期。命中的p值将从以0为中心的截断正态分布中选择。
# #模拟p值#库(msm) # hits <- where (get.vertex.)Attribute (g, "hits") == 1) # p.values <- runif(vcount(g)) # names(p.values) <- as.character(1:vcount(g)) # p.values[as.character(hits)] <- rtnorm(n。Hits * clusters, mean=0, # sd=10e-6, lower=0, upper=1)
我们现在将BioNet应用于网络和模拟的p值,以尝试识别丰富的子网络。BioNet工具的完整解释在BioNet R-package小插图中给出。模块
包含标识为位于丰富子网络中的每个顶点。
# #应用BioNet到网络和p值# bum <- fitBumModel(p。values, plot=F) # scores <- scoreNodes(network=g, fb=bum, fdr=0.1) # module <- runFastHeinz(g, scores)
我们现在将应用Knode
对网络的功能。首先,我们将p值转换为顶点权重。这将通过取p值的\(-log_{10}\)来完成,以确保具有最低p值的顶点具有最高的分数。Knode
返回网络中每个顶点的分数排名列表。
# #应用Knode到网络# g <- set.vertex。属性(g, name="pheno", value=-log10(p.values)) # knode。results <- Knode(g, dist.method="diffusion", #顶点。attr = "把" verbose = FALSE)
总共有30次命中网络。现在,我们将比较BioNet识别的丰富子网络中包含的这30个命中数,以及由Knode
函数。
# # BioNet识别的点击数# sum(hits %in% as.numeric(V(module)$name)) # # # Knode排名前30的点击数# sum(hits %in% as.numeric(names(knod .results))[1:(n. number)]点击*集群)]))
如前所述,虽然BioNet倾向于识别仅包含在单个集群中的命中,但事实是Knode
单独考虑每个顶点允许它识别包含在多个簇中的命中。图~\ref{fig:bionet_simulation}显示了仅从网络中的命中创建的子网。在第一个子网络和第二个子网络中,由BioNet和Knode识别的命中结果是彩色的。由BioNet成功识别的顶点在左侧网络中显示为蓝色,由BioNet成功识别的顶点为Knode
~在右边的网络上是绿色的。
# #创建子网# g.bionet <- g.knode <-诱导。Subgraph (g, hits) # color。Bionet <- color。Knode <- rep("grey", vcount(g.b onet)) # color。生命网络[hits %in% as.numeric(V(module)$name)] <- "blue" # color.knode[hits %in% as.numeric(names(knode.results)[1:(n.hits * clusters)])] <- "green" # g.bionet <- set.vertex.attribute(g.bionet, "color", value=color.bionet) # g.knode <- set.vertex.attribute(g.knode, "color", value=color.knode) # # # plot subnetworks # par(mfrow=c(1,2)) # plot(g.bionet) # plot(g.knode)
Knode
和BioNet使用真实数据上一节比较了Knode
BioNet使用模拟数据;在本节中,我们将使用真实数据进行比较。这些数据的分析最初是在BioNet R-package的小插图中进行的,Beisser等人在论文中进行了讨论。这使得比较两种方法更加容易。
在本分析中,我们将使用两种弥漫大b细胞淋巴瘤亚型(DLBCL, ABC和GCB亚型)的表达数据,来自生存分析的p值和来自HPRD数据库的网络。我们将使用这些数据来识别DLBCL涉及的交互模块。由于所需时间,此代码被注释。
# #加载所需的包#库(SANTA) #库(BioNet) #库(DLBCL) #数据(exprLym) #数据(dataLym) #数据(interactome)
为了同时使用表达式和生存数据,有必要将其合并为单个p值。首先,差异分析用于识别两种肿瘤类型之间表达显著差异的基因(结果包含在DLBCL包中)。然后,我们将使用顺序统计数据将差分表达式p值与生存p值结合起来。
# #提取entrez id #库(stringr) # # #聚合p-values # pvals <- cbind(dataLym$t.;pval, dataLym$s.pval) # pval <- aggrPvals(pval, order=2, plot=F) # names(pval) <- dataLym$label
HPRD网络中包含的许多基因,我们没有它们的表达或存活数据。因此,我们将移除这些基因。
的圣诞老人
包需要网络igraph
对象。因此,我们将网络从agraphNEL
对象。
# #派生淋巴芯片特定网络#网络<- subNetwork(featureNames(exprLym), interactome) #网络<- largestComp(网络)#只使用最大的组件#网络<- igraph.from.graphNEL(网络,name=T, weight=T) #网络<- simplify(网络)
我们现在将使用BioNet来识别网络中的丰富子网。
# fb <- fitBumModel(pval, plot=F) # scores <- scoreNodes(network, fb, fdr=0.001) # module <- runFastHeinz(network, scores) # extract。entrez < -函数(x) str_extract (str_extract (x, #”[(][0 - 9]+[0 - 9]+“),#生命网络。- extract.entrez(V(模块)$name)
如前所述,Knode
函数要求将p值转换为顶点权重。因此,我们将获得聚合p值的\(-log_{10}\),并在运行Knode
函数。
# #将p值转换为顶点权重权重<- -log10(pval)[get.顶点。属性(network, # "name")] # network <- set.vertex. set. name . name . name . name . name . name . name。attribute(network, name="pheno", # value=vertex.weights) # # #在淋巴芯片特定的网络上运行Knode, # #转换聚合p值# Knode。results <- Knode(network, dist.method="diffusion", #顶点。attr="pheno", verbose=F) # knode。基因<- extract.entrez(names(knod .results)[1:vcount(模块)])
由于BioNet和Knode
,识别出不同的富集子网络。Beisser等人注意到,当他们将BioNet应用于该数据时,他们确定了参与凋亡过程调节的基因模块的显著富集(GO:0042981)。经申请Knode
,鉴定出凋亡调控中富集的模块(\(p<1*10^{-7}\))。
# data(go.entrez) # sum(go. entrez)在% bionet.genes中输入# sum(开始。在% know .genes中输入
本文件使用以下方法制作:
sessionInfo ()
## R版本3.3.1(2016-06-21)##平台:x86_64-pc-linux-gnu(64位)##运行在:Ubuntu 16.04.1 LTS ## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_US。UTF-8 LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基础包:## [1]parallel stats4 stats graphics grDevices utils datasets ##[8]方法基础## ##其他附加包:## [1]org.Sc.sgd.db_3.4.0 AnnotationDbi_1.36.0 IRanges_2.8.0 ## [4] S4Vectors_0.12.0 Biobase_2.34.0 BiocGenerics_0.20.0 ## [7] SANTA_2.12.0 igraph_1.0.1 formatr1.4 ## [10] knitr_1.14 ## ##通过命名空间加载(且未附加):## [5] magrittr_1. 1.5 RSQLite_1.0.0 evaluate_0.10 stringi_1.1.2 ## [9] Matrix_1.2-7.1 tools_3.3.1 string_1 .1.0
[1] S. Bandyopadhyay, M. Mehta, D. Kuo等,“DNA损伤对基因网络的反应”。:科学330(2010),第1385-1390页。
[2] D. Beisser, G. W. Klau, T. Dandekar等,“BioNet:用于生物网络功能分析的R-Package”。:生物信息学26.8(2010),第1129-30页。ISSN: 1367 - 4811。DOI: 10.1093 /生物信息学/ btq089。
[3] H. W.张,G. S. Cowley, B. a. Weir等,“对癌细胞系遗传脆弱性的系统调查揭示了卵巢癌的谱系特异性依赖性”。:PNAS108.30(2011),第12372-7页。ISSN: 1091 - 6490。DOI: 10.1073 / pnas.1109363108。
M. Costanzo, a . Baryshnikova, J. Bellay等,“细胞的遗传景观”。:科学327.5964(2010),第425-31页。ISSN: 1095 - 9203。DOI: 10.1126 / science.1180823。
C.盖坦,X.盖恩。空间统计与建模。纽约:施普林格,2010。
[6] I. Lee, U. M. Blom, P. I. Wang等,“通过基于网络的全基因组关联数据增强来优先考虑候选疾病基因”。:基因组研究21(2011),第1109-21页。DOI: 10.1101 / gr.118992.110.7。
S. Orchard, M. Ammari, B. Aranda等,“MIntAct项目-完整作为11个分子相互作用数据库的公共管理平台”。:核酸研究42.1 (2014), pp. D358-63。ISSN: 1362 - 4962。DOI: 10.1093 / nar / gkt1115。
[8] S. Peri, J. D. Navarro, R. Amanchy等,“人类蛋白质参考数据库的开发,作为接近人类系统生物学的初始平台”。:基因组研究13.10(2003),第2363-71页。ISSN: 1088 - 9051。DOI: 10.1101 / gr.1680803。
R. Srivas, T. costele, A. Carvunis等,“紫外线诱导的遗传网络将RSC复合体与核苷酸切除修复联系起来,并显示出剂量依赖的重布线。”:细胞的报道5.6(2013),第1714-24页。ISSN: 2211 - 1247。DOI: 10.1016 / j.celrep.2013.11.035。
S.怀特和P.史密斯。估计网络中相对重要性的算法。:第九届ACM SIGKDD知识发现和数据挖掘国际会议论文集(2003),第266-75页。