背景と目的

STRING はPPI (Protein-Protein Interaction) ネットワークの公開データベースで、生命情報学の研究でよく使われます。例えば、何らかのトランスクリプトームデータの解析をして数十個から数百個程度の発現変動遺伝子 (Differentially Expressed Genes, DEG) のリストが得られたとき、エンリッチメント解析 (DAVID, Enrichr, Metascape, clusterProfilerなど) と合わせてSTRINGのネットワーク見ることで、遺伝子リストの大まかな特徴を捉えられる場合があります。

しかし、STRINGのネットワーク図は、遺伝子数が多いとPCのモニタの一画面分に収まらないほど巨大になります。その場合、拡大して描画範囲を移動しながら眺めなくてはならなくなります。また、発表スライドや論文に載せる図としても使いづらいです。

そこで、下記の投稿で作った関数を利用して、STRINGの出力結果をコンパクトに再描画するPython関数を作ってみました。遺伝子数が100くらいまでなら、図を拡大せずそのまま全体を眺められると思います。

networkxの頂点を一様に配置する関数

手法

STRINGでPPIネットワークを取得する手順

遺伝子リストを用意 (遺伝子名が縦一列に並んだcsvファイル)
STRING (https://string-db.org/) にアクセス
左メニューから "Multiple Proteins" を選択
"Browse..." ボタンより遺伝子リストのファイルをアップロード
Organismsを適宜選択 (ヒトならHomo Sapiens、マウスならMus Musculus)
"SEARCH" ボタンを押す
次ページ以降は "CONTINUE" ボタン連打でOK
ネットワーク図が表示されたら、下にスクロールして、ネットワークの下に横一列で並んでいるタブの中から "Exports" を選択
ダウンロード可能なファイル形式の一覧が出てくるので、2つあるTSV (Tab-Separated Values) 形式の上の方の "download" を押す

ちなみに、STRINGのネットワークの頂点配置は、公式ページの説明によるとバネモデルを使用しているとのことです。しかし、networkxのバネモデルより遥かに見やすいので、おそらく内部で色々と工夫しているのだろうと思います。少なくとも頂点同士が重ならないような何らかの処理をしているはずです。

一様配置関数

以下に一様配置関数のソースコードを載せます。頂点同士が互いに反発し合うようにして、全体的に一様に散らばるようにする関数です。初期配置はバネモデルで決めます。各ステップで各頂点の最近傍点を求め、その反対方向に少し動かす、ということを繰り返しています。描画領域の端にぶつかったらそれ以上いかないように制限しています。

import numpy as np
import networkx as nx
from scipy.cluster.hierarchy import distance

def uniform_layout(G, alpha=0.1, n_iter=None, seed=None, circle=False,
                   **kwargs):
  pos = nx.spring_layout(G, seed=seed, **kwargs)
  X = np.array(list(pos.values()))
  if n_iter == None:
    n_iter = 10 * len(G)
  for _ in range(n_iter):
    D = distance.squareform(distance.pdist(X))
    np.fill_diagonal(D, None)
    X += alpha * (X - X[np.nanargmin(D, axis=0)])
    if circle:
      X = (X.T / np.linalg.norm(X,axis=1).clip(1,None)).T
    else:
      X = X.clip(-1,1)
  return dict(zip(pos.keys(), X))

STRINGの出力結果の再描画関数

以下に再描画関数のソースコードを載せます。STRINGの出力結果のTSVファイルをpandasのデータフレームとして読み込んで、networkxのGraphオブジェクトに変換し、最大連結成分だけ残して他は捨ててしまって、前述の一様配置関数にかけています。頂点の色は何でも構いませんが、ひとまず次数を割り当てています。ただし、頂点の描画順は必ずしも名前順とは限らないため、頂点の順と次数の順を明示的に指定しています。

細かい調整の説明は省きます。頂点の大きさ、文字の大きさ、頂点の色の範囲、余白の大きさなど、あまり深く考えず適当に調節しています。適宜修正して使って下さい。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import networkx as nx
sns.set_context('talk', font_scale=0.8)

def plot_string(file_name):

  # load a tsv file obtained from STRING (https://string-db.org/).
  df = pd.read_table(file_name, usecols=[0,1])

  # convert it to a 'Graph' object of networkx
  df.columns = ['source','target']
  G = nx.from_pandas_edgelist(df)

  # select the largest connected component
  lcc_node_set = max(nx.connected_components(G), key=len)
  G = G.subgraph(lcc_node_set)

  # calculate a uniform layout
  pos = uniform_layout(G, seed=12345, circle=True)

  # calculate node colors based on degrees
  degrees = G.degree()
  nodes = G.nodes()
  node_color = [degrees[x] for x in nodes]
  d_max = max(node_color)

  # plot
  n = len(G)
  fig, ax = plt.subplots(figsize=(8,8))
  nx.draw_networkx(G,
                   pos=pos,
                   ax=ax,
                   node_size = np.min([10**5 / n, 5000]),
                   font_size = np.min([120 / n**0.5, 20]),
                   nodelist = nodes,
                   node_color = node_color,
                   cmap = 'Blues',
                   vmin = -1 * d_max,
                   vmax = 2 * d_max,
                   edge_color = '0.7',
                   width=1)
  if n > 50:
    ax.margins(0.1,0.05)
  elif n > 30:
    ax.margins(0.15,0.1)
  else:
    ax.margins(0.2)
  ax.axis('off')
  fig.tight_layout()
  fig.show()
  fig.savefig('tmp.png')

if __name__ == '__main__':
  file_name = sys.argv[1]
  plot_string(file_name)

STRINGの出力結果の再描画関数その2 (2024/04/29追記)

昨日のバージョンでは隣接しない頂点同士が近くに配置される問題があったため、ネットワーク上の近さを色で表すバージョンも作ってみました。以下にソースコードを載せます。Communicabilityの逆数を階層的クラスタリングにかけ、その順にhuslの色を割り当てています。次数は頂点の大きさに対応させました。

技術的な工夫は、networkxのCommunicabilityの関数があまりにも遅かったので書き直したことと (単に隣接行列を固有値分解して固有値を指数関数にかけて元に戻すだけです)、huslがmatplotlibには無くてseabornにしかないため、ListedColormapを使って渡したことの2点です。

ついでに文字の大きさ、枝の色、余白の大きさも微調整しました。また、長い名前は改行するようにしました。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
import networkx as nx
import re
from scipy.cluster.hierarchy import linkage, leaves_list
sns.set_context('talk', font_scale=0.8)

def label_wrap(sr):
  out_list = []
  for x in sr:
    if len(x) >= 6:
      i = int(np.ceil(len(x)/2))
      m = re.match('[^\d]+', x)
      if m:
        i2 = m.end()
        if np.abs(i-i2) <= 2:
          i = i2
      x = x[:i] + '\n' + x[i:]
    out_list.append(x)
  return pd.Series(out_list)

def plot_string2(file_name):

  # load a tsv file obtained from STRING (https://string-db.org/).
  df = pd.read_table(file_name, usecols=[0,1])

  # convert it to a 'Graph' object of networkx
  df.columns = ['source','target']
  df['source'] = label_wrap(df.source)
  df['target'] = label_wrap(df.target)
  G = nx.from_pandas_edgelist(df)

  # select the largest connected component
  lcc_node_set = max(nx.connected_components(G), key=len)
  G = G.subgraph(lcc_node_set)

  # calculate a uniform layout
  pos = uniform_layout(G, seed=12345, circle=True)

  # calculate node colors based on communicability
  # NOTE: nx.communicability() is too slow and not recommended.
  A = nx.adjacency_matrix(G).toarray()
  w_arr, V = np.linalg.eigh(A)
  C = V.dot(np.diag(np.exp(w_arr))).dot(V.T)
  D = 1 / C[np.triu_indices_from(C,1)]
  Z = linkage(D, method='ward', metric='euclidean')
  node_arr = np.asarray(G.nodes())
  sr = pd.Series(np.arange(len(G)), index=node_arr[leaves_list(Z)])
  node_color = sr.loc[node_arr]
  cmap = sns.color_palette('husl', len(G))
  cmap = np.array(cmap) ** 0.5  # brighter
  cmap = mpl.colors.ListedColormap(cmap)

  # calculate node size based on degrees
  scale = 0.2
  degrees = G.degree()
  deg_arr = np.array([degrees[x] for x in G.nodes()])
  deg_arr = 2 * scale * (deg_arr - deg_arr.min()) / \
    (deg_arr.max() - deg_arr.min()) + 1 - scale

  # plot
  n = len(G)
  fig, ax = plt.subplots(figsize=(8,8))
  nx.draw_networkx(G,
                   pos=pos,
                   ax=ax,
                   node_size = np.min([150000 / n, 3000]) * deg_arr,
                   font_size = np.min([140 / n**0.5, 20]),
                   node_color = node_color,
                   cmap = cmap,
                   edge_color = '0.6',
                   width=1)
  ax.margins(np.clip(300 / n**2, 0.05, 0.2),
             np.clip(0.05 + 200 / n**2, 0.05, 0.2))
  ax.axis('off')
  fig.tight_layout()
  fig.show()
  fig.savefig('tmp.png')

if __name__ == '__main__':
  file_name = sys.argv[1]
  plot_string2(file_name)

結果

公開可能な例として、DAVIDのデモリスト1と2の結果を示します。

デモリスト1

ヒトの149個の遺伝子リストです。遺伝子記号は145個です。これをSTRINGにかけた結果が以下の図です。既定の設定よりフォントサイズを大きくし、遺伝子名の位置も各頂点の右上から中央に変えてみました。決して見づらい図ではありません。ただし、孤立遺伝子や2, 3個だけで繋がっている遺伝子が多いため、それらを削除すれば、最大連結成分をもっと大きく表示出来るだろうと期待されます。

デモリスト1のネットワークを再描画した結果が以下の図です。図全体の大きさに対する文字の大きさが格段に大きくなりました。ハブやコミュニティ構造も読み取れます。欠点としては、近くにあってもネットワーク上で隣接していない (直接つながっていない) 場合が多いことなどが挙げられます。なお、枝の太さは一定としましたが、STRINGの出力ファイルに各枝のスコアが含まれるため、元図と同様に枝の太さをスコアに応じて変えることは技術的には可能です。

[2024/04/29追記] デモリスト1の再描画関数その2の結果が以下の図です。中心部のコミュニティが比較的近い色 (ピンク色) で表されています。図上の距離が近くても色が大きく違うものは、実際のネットワーク上では離れている頂点です。ただし、使える色の数が限られているため、色が似ているからといって必ずしもネットワーク上で近いとは限りません。あくまで色は補助的なもので、基本的には枝でネットワーク構造を判断して下さい。

デモリスト2

デモリスト2はヒトの390個の遺伝子リストです。遺伝子記号は372個です。デモリスト1と比べて遺伝子数が2倍以上に増えています。これをSTRINGにかけた結果が以下の図です。明らかに中央部に1つの密に繋がったコミュニティがあります。おそらくここに、この遺伝子リストの中で中心的な役割を果たす遺伝子が集積しているだろうと考えられます。しかし、中心部は頂点と枝が密集していて遺伝子名が読みづらいです。

デモリスト2のネットワークを再描画したのが以下の図です。文字の大きさは小さいものの、頂点同士が重なっていないので元図よりずっと読みやすくなっています。

[2024/04/29追記] デモリスト2の再描画関数その2の結果が以下の図です。中心部のコミュニティが、枝、頂点の大きさ、頂点の色の3種類で強調されています。

まとめ

STRINGのPPIネットワークをPythonでコンパクトに再描画する関数を作りました。発現変動遺伝子などの遺伝子リストの意味を理解する上で便利なツールになると思います。特に、ある程度主要な遺伝子の名前と機能を知っている中級者以上にとっては、エンリッチメント解析よりも分かりやすいかもしれません。何故なら、STRINGのネットワークは既知の知見に基づいており、有名な遺伝子ほどネットワーク上で繋がりやすいためです。GO (Gene Ontology) の注釈を見ただけではピンと来なかったものの、STRINGでネットワークの中心部にいる遺伝子名を見たら、実際の生体内で何が起きていたかある程度予想出来た、という経験が私自身もあります。その意味で、ネットワークの構造が多少見づらくなっても、頂点の名前がはっきり読める表示法は大事だろうと思います。

最大連結成分以外を除いていますが、上記コードの該当箇所2行をコメントアウトすれば全遺伝子を一様配置することも可能です (2024/04/29追記: Communicabilityは非連結グラフだと0が現れ、その逆数を使っているのでゼロ割が発生します。適宜修正して下さい)。