KEGGgraph : 根据kgml 文件从pathway中重构出基因互作网络

本文主要是介绍KEGGgraph : 根据kgml 文件从pathway中重构出基因互作网络，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

KEGGgraph 包可以解析kgml 文件，从中得到不同对象之间的网络结构，并在此基础上进一步挖掘其中的信息。

KEGGgraph 包提供了以下3种基本功能：

将kgml 文件转换为基因互作网络

用法示例

# 读取hsa00020xml 文件
> mapkG <- parseKGML2Graph("hsa00020.xml",expandGenes=TRUE, genesOnly = TRUE)
> mapkG
A graphNEL graph with directed edges
Number of Nodes = 30
Number of Edges = 101
>  nodes(mapkG)
[1] "hsa:1738"  "hsa:4967"  "hsa:55753" "hsa:1743"  "hsa:8801"  "hsa:8802"
[7] "hsa:8803"  "hsa:3417"  "hsa:3418"  "hsa:3419"  "hsa:3420"  "hsa:3421"
[13] "hsa:47"    "hsa:2271"  "hsa:48"    "hsa:50"    "hsa:1431"  "hsa:4190"
[19] "hsa:4191"  "hsa:5091"  "hsa:5160"  "hsa:5161"  "hsa:5162"  "hsa:1737"
[25] "hsa:5105"  "hsa:5106"  "hsa:6389"  "hsa:6390"  "hsa:6391"  "hsa:6392"
>  edges(mapkG)
$`hsa:1738`
[1] "hsa:4967"  "hsa:55753" "hsa:5160"  "hsa:5161"  "hsa:5162"  "hsa:1737"
$`hsa:4967`
[1] "hsa:3419" "hsa:3420" "hsa:3421" "hsa:3417" "hsa:3418"
$`hsa:55753`
[1] "hsa:3419" "hsa:3420" "hsa:3421" "hsa:3417" "hsa:3418"

在 parseKGML2Graph 中，有两个参数，expandGenes 和 genesOnly。

expandGenes 控制是否将基因进行展开，在pathway 中，会有1个KO 对应多个gene的情况，比如下面这种

    <entry id="32" name="hsa:8801 hsa:8802 hsa:8803" type="gene" reaction="rn:R00405"link="">http://www.kegg.jp/dbget-bin/www_bget?hsa:8801+hsa:8802+hsa:8803">;<graphics name="SUCLG2, G-SCS, GBETA, GTPSCS..." fgcolor="#000000" bgcolor="#BFFFBF"type="rectangle" x="260" y="574" width="46" height="17"/></entry>

expandGenes = TRUE 表示将基因展开，每个基因作为一个节点。

genesOnly 参数控制是否将其他类型的entry （比如compound等类型）展现在network 中，默认值为 TRUE，所以最终得到的network 中节点全部是基因。

通过parseKGML2Graph 这一步我们就可以从一张pathway 中得到基因产物（蛋白）的互作网络，还需要注意一点，整个网络是一个有向图，因为基因产物之间的互作关系是由方向性的。

对network 进行可视化

由于自带的可视化不够美观，我们把nodes和edges 写入文件，用cytoscape 进行可视化，用法示例

mapkNodes <- nodes(mapkG)
mapkEdges <- edges(mapkG)
mapkEdges <- mapkEdges[sapply(mapkEdges, length) > 0]
res <- lapply(1:length(mapkEdges), function(t){name <- names(mapkEdges)[t]len  <- length(mapkEdges[[t]])do.call(rbind, lapply(1:len, function(n){c(name, mapkEdges[[t]][n])}))
})
result <- data.frame(do.call(rbind, res))
write.table(result,  "edges.txt", sep = "\t", row.names = F, col.names = F, quote = F)
write.table(mapkNodes, "nodes.txt", sep = "\t", row.names = F, col.names = F, quote = F)

导入cytoscape 画出来的图

查询节点的`degree`信息

对于一个netwrok 而言，每个节点的degree 信息是我们最常用的信息，示例

> mapkGoutdegrees <- sapply(edges(mapkG), length)
> mapkGindegrees <- sapply(inEdges(mapkG), length)
> degrees <- data.frame(indegrees = mapkGindegrees, outdegrees = mapkGoutdegrees)
> head(degrees)indegrees outdegrees
hsa:1738          1          6
hsa:4967          2          5
hsa:55753         2          5
hsa:1743          3          3
hsa:8801          4          1
hsa:8802          4          1

由于是有向图，所以有入度 indegrees 和出度 outdegrees 的概念。

除了以上基础功能外，还可以借助其他的R包进一步挖掘信息，比如在整个基因互作网络，哪个基因是最关键的。

示例：

> library(RBGL)
> mapkG <- parseKGML2Graph("hsa00020.xml",expandGenes=TRUE, genesOnly = TRUE)
>  bcc <- brandes.betweenness.centrality(mapkG)
> rbccs <- bcc$relative.betweenness.centrality.vertices[1L,]
> toprbccs <- sort(rbccs,decreasing=TRUE)[1:4]
> toprbccshsa:1743   hsa:2271   hsa:1738     hsa:47
0.21597893 0.16177167 0.14965648 0.09880362

对于network 而言，我们一般认为degree 越大的点在这个网络中越重要，所以需要看节点的degree 信息。除了这种基本的认识外，还有很多成熟的算法，从network 中挖掘关键节点。 RBGL 包提供了Brandes 的算法，用来衡量节点在网络中的重要性，上面的结果中，toprbccs 就是我们筛选出的4个比较重要的节点。