PythonのnetworkXを使ってpathway解析をする

ある遺伝子(機能)セットと代謝ネットワークデータを用いて、上流の化合物から下流の化合物に行くpathwayがあるか調べたい。
PythonのnetworkXを使って最短経路を出すドキュメントはたくさんあったが、エッジの有無によってFlow解析についてはあまりなかったのでまとめた。

Flow解析をするには、edgeに"capacity"を設定する必要がある。
0→1→(2,3)→4 という枝分かれを1つ持つ単純なネットワークを、以下のように設定する。

>>> import networkx as nx
>>> Graph = nx.DiGraph()
>>> Graph.add_nodes_from([0,1,2,3,4])
>>> Graph.add_edges_from([[0,1], [1,2], [1,3], [2,4], [3,4]])
>>> Graph[0][1]["capacity"] = 1.0
>>> Graph[1][2]["capacity"] = 1.0
>>> Graph[1][3]["capacity"] = 1.0
>>> Graph[2][4]["capacity"] = 1.0
>>> Graph[3][4]["capacity"] = 1.0

edgeの"capacity"を1.0に設定することは、その遺伝子が含まれることを意味する。
今、すべてのedgeの"capacity"が1.0なので、

>>> nx.maximum_flow(Graph, 0, 4)
(1.0, {0: {1: 1.0}, 1: {2: 1.0, 3: 0}, 2: {4: 1.0}, 3: {4: 0}, 4: {}})

化合物0から化合物4までの経路は1.0であり、すなわち最下流の物質を生成するpathが存在することを意味する。
(複数存在する場合、もっとも短い経路で、かつ若いインデックスを返すことに注意)

ここで、枝分かれのうちの一つの"capacity"を0.0に設定する(遺伝子を欠損させる)。

>>> Graph[2][4]["capacity"] = 0.0
>>> nx.maximum_flow(Graph, 0, 4)
(1.0, {0: {1: 1.0}, 1: {2: 0, 3: 1.0}, 2: {4: 0}, 3: {4: 1.0}, 4: {}})

もう片方の枝分かれが生きているので、化合物0から化合物4までの経路は依然として1.0である。

もう一つの枝分かれを欠損させると、

>>> Graph[3][4]["capacity"] = 0.0
>>> nx.maximum_flow(Graph, 0, 4)
(0, {0: {1: 0}, 1: {2: 0, 3: 0}, 2: {4: 0}, 3: {4: 0}, 4: {}})

化合物0から化合物4までの経路は0になった。

実際にはKEGGなどからPathwayを取ってきて、目的のゲノムなどの機能組成を構築し、networkXに埋め込む必要があるがここでは割愛する。

参考