前言
igraph 提供了许多图论算法，例如我们熟知的最短路径、最大流最小割、最小生成树等。
我们针对下面的图，介绍一些简单的图论相关的函数
g <- graph_from_data_frame(d = sub_path, vertices = node_attr)
V(g)$color <- if_else(V(g)$type == "other", "#80b1d3", "#fb8072")
V(g)$size <- (V(g)$mut_class + 1) * 6
E(g)$weight <- abs(rnorm(n = ecount(g)))
E(g)$arrow.size <- .2
E(g)$edge.color <- "gray80"
graph_attr(g, "layout") <- layout_with_lgl
plot(g)
1. 距离和路径
计算平均路径长度，即每对基因之间的最短距离的均值
> mean_distance(g)
[1] 2.625
当图是无向图时的平均路径长度
> mean_distance(g, directed = FALSE)
[1] 2.571429
计算所有位点之间的最短路径，默认会考虑边的权重
> distances(g, weights = NA)
EGFR SOS1 SOS2 GRB2 EGF MAP2K1 MAP2K2 ARAF RAF1 BRAF HRAS KRAS NRAS MAPK1 MAPK3
EGFR 0 2 2 1 1 5 5 4 4 4 3 3 3 6 6
SOS1 2 0 2 1 3 3 3 2 2 2 1 1 1 4 4
SOS2 2 2 0 1 3 3 3 2 2 2 1 1 1 4 4
GRB2 1 1 1 0 2 4 4 3 3 3 2 2 2 5 5
EGF 1 3 3 2 0 6 6 5 5 5 4 4 4 7 7
MAP2K1 5 3 3 4 6 0 2 1 1 1 2 2 2 1 1
MAP2K2 5 3 3 4 6 2 0 1 1 1 2 2 2 1 1
ARAF 4 2 2 3 5 1 1 0 2 2 1 1 1 2 2
RAF1 4 2 2 3 5 1 1 2 0 2 1 1 1 2 2
BRAF 4 2 2 3 5 1 1 2 2 0 1 1 1 2 2
HRAS 3 1 1 2 4 2 2 1 1 1 0 2 2 3 3
KRAS 3 1 1 2 4 2 2 1 1 1 2 0 2 3 3
NRAS 3 1 1 2 4 2 2 1 1 1 2 2 0 3 3
MAPK1 6 4 4 5 7 1 1 2 2 2 3 3 3 0 2
MAPK3 6 4 4 5 7 1 1 2 2 2 3 3 3 2 0
设置 weights = NA，表示不考虑边的权重
使用 distances，我们可以计算我们感兴趣的基因与其他基因之间的距离，例如，我们想知道 EGFR 与其他基因之间的距离
dist_egfr <- distances(g, v = V(g)[1], to = V(g), weights = NA)
# 设置颜色板
pal <- colorRampPalette(c("#fb8072", "#80b1d3"))
col <- pal(max(dist_egfr) + 1)
plot(g, vertex.color = col[dist_egfr+1], vertex.label = dist_egfr,
edge.arrow.size = .4, vertex.label.color = "white")
我们也可以只计算某两个基因之间的最短路径，例如，计算 EGFR 到 MAPK1 之间的最短距离
e2m_path <- shortest_paths(g, from = V(g)[name == "EGFR"],
to = V(g)[name == "MAPK1"],
output = "both")
# 设置边的颜色
edge_col <- rep("grey80", ecount(g))
edge_col[unlist(e2m_path$epath)] <- "#80b1d3"
# 设置边的宽度
edge_width <- rep(2, ecount(g))
edge_width[unlist(e2m_path$epath)] <- 4
# 设置节点的颜色
vcol <- rep("gray40", vcount(g))
vcol[unlist(e2m_path$vpath)] <- "#fb8072"
plot(g, vertex.color=vcol, edge.color=edge_col,
edge.width=edge_width, edge.arrow.size = .6,
vertex.label.cex = 1
)
我们可以获取与一个基因相关的所有边，例如
> incident(g, V(g)[name=="KRAS"], mode="in")
+ 2/29 edges from d01cfa6 (vertex names):
[1] SOS1->KRAS SOS2->KRAS
> incident(g, V(g)[name=="KRAS"], mode="out")
+ 3/29 edges from d01cfa6 (vertex names):
[1] KRAS->ARAF KRAS->RAF1 KRAS->BRAF
也可以一次性获取多个基因的边
> incident_edges(g, V(g)[1:2], mode = "all")
$EGFR
+ 2/29 edges from d01cfa6 (vertex names):
[1] EGFR->GRB2 EGF ->EGFR
$SOS1
+ 4/29 edges from d01cfa6 (vertex names):
[1] SOS1->HRAS SOS1->KRAS SOS1->NRAS GRB2->SOS1
类似地，我们可以获取与一个基因直接互作的基因
> neighbors(g, V(g)[name == "KRAS"], mode = "all")
+ 5/15 vertices, named, from d01cfa6:
[1] SOS1 SOS2 ARAF RAF1 BRAF
获取调控 KRAS 基因的基因
> neighbors(g, V(g)[name == "KRAS"], mode = "in")
+ 2/15 vertices, named, from d01cfa6:
[1] SOS1 SOS2
获取被 KRAS 调控的基因
> neighbors(g, V(g)[name == "KRAS"], mode = "out")
+ 3/15 vertices, named, from d01cfa6:
[1] ARAF RAF1 BRAF
如果我们想获取多个感兴趣的基因的直接互作基因，可以使用 adjacent_vertices()
> adjacent_vertices(g, v = V(g)[c(1, 5)], mode = "all")
$EGFR
+ 2/15 vertices, named, from d01cfa6:
[1] GRB2 EGF
$EGF
+ 1/15 vertex, named, from d01cfa6:
[1] EGFR
上面两个函数获取的都是直接互作的基因，即一步互作基因。如果想要获取多步互作基因，可以使用 ego() 函数，并设置 order 参数。
> ego(g, order = 1, nodes = V(g)["KRAS"])
[[1]]
+ 6/15 vertices, named, from d01cfa6:
[1] KRAS SOS1 SOS2 ARAF RAF1 BRAF
> ego(g, order = 2, nodes = V(g)["KRAS"])
[[1]]
+ 11/15 vertices, named, from d01cfa6:
[1] KRAS SOS1 SOS2 ARAF RAF1 BRAF GRB2 HRAS NRAS MAP2K1 MAP2K2
order = 2 或包含 order = 0 和 order = 1 的结果
# 获取 KRAS 的所有边
inc.edge <- incident(g, V(g)[name=="KRAS"], mode="all")
# 设置这些边的颜色
ecol <- rep("grey80", ecount(g))
ecol[inc.edge] <- "#33a02c"
# 获取直接互作基因
nei.node <- neighbors(g, V(g)[name == "KRAS"], mode = "all")
# 设置这些基因的颜色
vcol <- rep("gray40", vcount(g))
vcol[nei.node] <- "#fb8072"
vcol[V(g)[name == "KRAS"]] <- "#80b1d3"
plot(g, vertex.color=vcol, edge.color=ecol)
2. 连通性
判断一个图是不是连通的、是强连通图还是弱连通图，可以使用 is_connected
> is_connected(g)
[1] TRUE
> is_connected(g, mode = "weak")
[1] TRUE
> is_connected(g, mode = "strong")
[1] FALSE
可以使用 count_components 计算连通分量的数目
> count_components(g)
[1] 1
因为我们的图是连通的，所以只有一个连通分量
获取图中的连通分量信息
> components(g)
$membership
EGFR SOS1 SOS2 GRB2 EGF MAP2K1 MAP2K2 ARAF RAF1 BRAF HRAS KRAS NRAS MAPK1 MAPK3
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
$csize
[1] 15
$no
[1] 1
是不是有向非循环图（DAG）
> is_dag(g)
[1] TRUE
对于如下图
g1 <- sample_gnp(20, 1/20)
plot(g1)
> count_components(g1)
[1] 12
> is_connected(g1)
[1] FALSE
将图分解为连通分量列表
> decompose(g1)
[[1]]
IGRAPH dd688c4 U--- 3 2 -- Erdos renyi (gnp) graph
+ attr: name (g/c), type (g/c), loops (g/l), p (g/n)
+ edges from dd688c4:
[1] 1--2 1--3
...
[[12]]
IGRAPH 7f3046a U--- 1 0 -- Erdos renyi (gnp) graph
+ attr: name (g/c), type (g/c), loops (g/l), p (g/n)
+ edges from 7f3046a:
根据连通分量进行分组
igraph::groups(components(g1))
$`1`
[1] 1 8 9
...
$`4`
[1] 4 6 12 13 17 19
...
$`12`
[1] 20
找到所有与给定顶点连通的点
> subcomponent(g1, 1)
+ 3/20 vertices, from 360a2eb:
[1] 1 8 9
> subcomponent(g1, 4)
+ 6/20 vertices, from 360a2eb:
[1] 4 13 17 6 12 19
3. 社区检测
有一些算法用于检测一些群组，即在群组内密集连接的节点，而群组之间的连接较少。
- 基于
Newman-Girvan
通过不断地删除网络中的边来检测网络中的社区。在最终剩余的网络中的连通分量也就是社区
计算图的社区结构
ceb <- cluster_edge_betweenness(g)
绘制社区结构树状图
dendPlot(ceb, mode="hclust")
在图中展示社区结构
plot(ceb, g)
我们可以来看下社区检测对象
> class(ceb)
[1] "communities"
获取社区结构的数目
> length(ceb)
[1] 3
查看每个社区结构包含的成员
> membership(ceb)
EGFR SOS1 SOS2 GRB2 EGF MAP2K1 MAP2K2 ARAF RAF1 BRAF HRAS KRAS NRAS MAPK1 MAPK3
1 2 2 1 1 3 2 2 2 2 2 2 2 3 3
分区的高度模块化反映了社区内的紧密连接和跨社区的稀疏连接。即基因之间的紧密互作关系。
- 基于贪婪算法
这个函数试图通过直接优化模块化得分来寻找密集子图
cfg <- cluster_fast_greedy(as.undirected(g))
plot(cfg, as.undirected(g))
- 基于随机游走
基于随机游走的思想，即短的随机游走往往停留在同一个社区
cw <- cluster_walktrap(g)
plot(cw, g)
还有其他的一些社区检测函数，例如 cluster_leading_eigen、cluster_label_prop、cluster_louvain 等，就不再一一介绍了
我们也可以自定义设置每个社区的绘制，例如
V(g)$community <- cfg$membership
col <- c("#fb8072", "#80b1d3", "#fdb462", "#b3de69")
plot(g, vertex.color = col[V(g)$community])
4. 节点的度
degree 可以计算节点的度，其中 mode 参数的值可以是：
-
in：入度 -
out：出度 -
all或total：所有的度
plot(g, vertex.size = degree(g, mode = "all") * 3)
查看节点度的分布
deg.dist <- degree_distribution(g, cumulative=T, mode="all")
plot(
x = 0:max(deg),
y = 1 - deg.dist,
pch = 19, cex = 1.2,
col = "orange", xlab = "Degree",
ylab = "Cumulative Frequency"
)
5. 网页分析算法
PageRank 算法是 Google 使用的网页排名算法。其基本假设是，重要的页面往往会被很多其他页面引用，即指向该节点的节点越多，则该节点越重要。
我们可以将该算法应用于基因互作网络的分析，越多的基因来调控的某些基因可能行使的功能越重要
pg <- page_rank(g)$vector
plot(g, vertex.size=pg*400)
Hyperlink-Induced Topic Search(HITS) 算法也是一个网页排序算法，与 PageRank 不同的是，它将网页分为：
-
Hub：如门户网站，会指向很多其他的网址。 -
Authority：用户实际希望访问的网站，如淘宝、京东等
一个页面的 Authority 值和 Hub 值的计算公式为
$$
\begin{align}
auth(p) = \sum_{q \in p_{to}}{hub(q)}\
hub(p) = \sum_{q \in p_{from}}{auth(q)}
\end{align}
$$
hs <- hub_score(g, weights=NA)$vector
as <- authority_score(g, weights=NA)$vector
par(mfrow=c(1,2))
plot(g, vertex.size=hs*50, main="Hubs")
plot(g, vertex.size=as*30, main="Authorities")
6. 图形的集合操作
我们可以使用 %du% 操作符合并两张图
g1 <- make_star(10, mode="undirected")
V(g1)$name <- letters[1:10]
g2 <- make_ring(10)
V(g2)$name <- letters[6:15]
plot(g1 %du% g2)
使用 %u% 取两张图的并集
plot(g1 %u% g2)
可以看到，由于两张图之间有交叠的节点，最后会合并为一张图
使用 %m% 取两张图的差集
plot(g1 %m% g2)
使用 %d% 取两张图的交集
par(mfcol = c(1, 2))
net1 <- graph_from_literal(D-A:B:F:G, A-C-F-A, B-E-G-B, A-B, F-G,
H-F:G, H-I-J)
net2 <- graph_from_literal(D-A:F:Y, B-A-X-F-H-Z, F-Y)
plot(net1 %du% net2, main = "disjoint union")
plot(net1 %s% net2, main = "intersection")
