Chapter 01無料公開

事前知識

Ryota Chijimatsu
Ryota Chijimatsu
2025.02.23に更新

【 グラフ 】

空間オミクス解析やwhole slide imageでの病理画像解析が盛んになり、その細胞や組織の分布、位置関係を解析するニーズが増えている。
細胞や組織のセグメンテーションを行った後は、その二次元空間(場合によっては三次元)での座標が得られる。基本的には近い座標距離にある者同士の相互作用に着目することが多いので、まずは距離に基づく隣接関係・接続関係をまとめていく。細胞を頂点(ノード)、細胞間の近い関係を辺(エッジ)としてグラフを築いておくと、「どことどこが近い位置関係にあるのか」、「組織上のどこに細胞が密集しているのか」、などを簡単に取り出せるようになる。

本書では二次元座標から作られるグラフとして、kNN、半径距離、ドロネー図を用いた方法が登場する。

ノード(頂点)とエッジ(辺)を用いてグラフを表す。

本書では細胞がノード、細胞と細胞の接続をエッジとしたグラフを考える。エッジには接続の有無で0/1の値を使うことが多いが、距離や接続確率のような値をエッジ重みとして使うこともできる。本書では基本は接続の有無を使用する。

隣接行列

どのノードとどのノードが繋がっているかを示す際に、ノード数*ノード数の行列で表現する。接続無しの組み合わせには0、接続有りの組み合わせには1を入れる。

https://www.geeksforgeeks.org/adjacency-matrix/

edge_list / edge_index

グラフが巨大になると(ノード数,ノード数)の巨大な行列に大量の0を保持する疎行列となる。(遠くのノードは繋がらないため、巨大になるほど0比率が増えていく。) データサイズを小さくするため、[接続元ノード番号, 接続先ノード番号]を一組として、(2,エッジ数)という形状で保持する。
NetworkXパッケージでは「edgelist」、PyTorch Geometric(以下PyG)では「edge_index」と称している。本記事ではedge_indexと記載する。

隣接行列 ⇔ edge_index

この2つは互いに簡単に変換可能である。

隣接行列の用意
import numpy as np

adj_matrix = np.array([
    [0,1,1,1,1,0],
    [1,0,0,1,0,0],
    [1,0,0,1,1,1],
    [1,1,1,0,0,1],
    [1,0,1,0,0,0],
    [0,0,1,1,0,0]
])

接続有り(値が1)のindexを取得すれば、edge_indexである。

隣接行列からedge_index作成
edge_index = np.argwhere(adj_matrix > 0)



edge_indexから隣接行列へは以下のように変換できる。

edge_indexから隣接行列作成
num_nodes = np.unique(edge_index).max() + 1
adj_matrix = np.zeros((num_nodes, num_nodes), dtype=int)
adj_matrix[edge_index[:,0], edge_index[:,1]] = 1

有向グラフ/無向グラフ

その名の通り、接続に向きが有るグラフを有向グラフ、向きが無いグラフを無向グラフと呼ぶ。上図は無向グラフで隣接行列は対称行列である。edge_indexでは[接続元ノード番号,接続先ノード番号]という順で扱うため、無向グラフでは同じノード番号の組み合わせでも、[0,1]と[1,0]という風に入れ替えた2つを記録しておく。

有向グラフの隣接行列を無向グラフ版に変換

有向グラフの隣接行列を無向グラフ版に変換するには以下のように、隣接行列に転置した隣接行列を足せばよい。

# 有向グラフの隣接行列
adj_matrix = np.array([
    [0,1,1,1,1,0],
    [1,0,0,1,0,0],
    [1,0,0,1,1,1],
    [0,0,0,0,0,1],
    [0,0,0,0,0,0],
    [0,0,0,0,0,0]
])

adj_matrix_undirected = adj_matrix + adj_matrix.T
adj_matrix_undirected[adj_matrix_undirected > 0] = 1

有向グラフのedge_indexを無向グラフ版に変換

有向グラフのedge_indexを無向グラフ版に変換するには以下のように、転送元・転送先を入れ替えて連結する。

edge_index_undirected = np.concatenate([
    edge_index,edge_index[:,[1,0]]
])

edge_index_undirected = np.unique(edge_index_undirected, axis=0)

次数 / 次数行列

1つのノードが持つ接続先の数(エッジ数)を次数と呼ぶ。ノードの接続先の情報を集計して計算する際に、ノードが沢山の接続を持つのか少ない接続を持つのかを補正することがある。
計算の上では次数行列として扱うことが多い。これはノード数*ノード数の行列の対角成分に次数、他が0の行列である。

edge_indexから次数行列へ

接続先が何個あるかは、edge_indexの接続元ノードの登場回数を調べればよい。np.unique()return_counts=Trueオプションが使える。

unique_nodes, degree = np.unique(edge_index[:,0],return_counts=True)

これを対角成分に配置した行列を作る。

num_nodes = unique_nodes.max() + 1
degree_matrix = np.zeros((num_nodes, num_nodes))
degree_matrix[unique_nodes, unique_nodes] = degree



ただし、次数0のノードがある場合、edge_indexにはそのノード番号が記録されていない。edge_indexをnp.unique()すると次数があるノードしかカウントされないため、ノード番号がずれてしまう。

次数0のノードがある場合

次のコードは全てのノード番号の次数を集計する一例である。edge_indexに登場しないノードの次数に0が入るようにしている。

from collections import Counter

num_nodes = np.unique(edge_index).max()
degree_dict = { node_id:0 for node_id in range(num_nodes)}

degree_dict.update(Counter(edge_index[:,0]))

degree_matrix = np.diag(list(degree_dict.values()))

隣接行列から次数行列へ

隣接行列から行方向(列方向)で合計を計算すればそのノードが持つ次数となる。これを対角成分に配置した行列を作る。

degrees = np.sum(adj_matrix, axis=0)
degree_matrix = np.diag(degrees)


【 細胞コミュニティ特有のグラフ 】

生物学の領域で扱われるグラフには遺伝子制御ネットワーク、代謝反応ネットワークなどが有名である。

遺伝子制御ネットワーク 代謝反応ネットワーク
10.1007/s12031-022-02068-w 10.1186/1475-2859-8-27

遺伝子制御ネットワークでは極端に次数が多いハブ遺伝子がコミュニティを形成したグラフをよく見かける。ハブ遺伝子は多くの他の遺伝子に繋がっていることが特徴なので、その次数に着目することが手掛かりになる。化学反応をネットワークで示した代謝反応ネットワークでは、少ないエッジで繋がったノードが方向性を持って繋がっている。そのコミュニティ解析では、ネットワーク全体から解糖系やTCA回路のような一連の反応系を部分に分ける際に用いられることがある。
このようにグラフ全体から部分的なコミュニティを見つける解析には生物学的な意義を見つける上で一つのアプローチとなる。

二次元(三次元)空間での細胞間の距離に基づいてエッジを考えると、大抵は画像上の細胞全てが繋がった1つのグラフになる。(画像上に細胞がいない領域で分断されていなければ)

QuPathのドロネーグラフの例


kNNグラフやドロネーグラフなどでは細胞が持つエッジ数に偏りは少なく、エッジ数は細胞の密集度の指標にはならない。一方、細胞間距離に閾値を設けてグラフを作ればエッジが多い少ないは細胞が密集しているか否かで決まる。単純に考えると、これを空間的なコミュニティ・クラスターとして捉えることもできる。

QuPathのドロネーグラフの例 距離20 μm以内に接続を限定

エッジが多い少ないの議論を病理空間に適応しても、細胞が密集しているか疎かを議論するだけになってしまう。本書では発展的な方法として、エッジで繋がった局所コミュニティでどのように細胞の種類が混ざり合っているかを考える。同じ細胞種が混じった領域なのか、異なる細胞種が入り混じっているのか、もし入り混じっているならどんな細胞で構成されているのかを定量的に解析するためにコミュニティ探索に挑戦する。


【 解析環境 】

本書は以下の環境を使用している。

  • Windows
  • Anaconda
  • Python 3.11
  • GPU RTX3090 24GBメモリ
  • torch 2.2.0+cu118

Python

本書では、グラフニューラルネットワークを使用した解析も含んでおり、それに先立ってPython環境の構築、必要なライブラリが必要となる。
グラフニューラルネットワークには、PyTorchおよびPyTorch Geometricをメインで使用する。

こちらでPythonとPyTorchの導入を紹介しているので参考にしてほしい。
https://zenn.dev/rchiji/books/016f8c2a6c0ae0/viewer/84c4c0


PyTorch Geometric

PyTorchのTensorを使ってグラフニューラルネットワークを行うのにメジャーなライブラリ。

公式サイトからインストールコマンドが得られる。自身の環境に合ったインストールコマンドを使用すればよい。

https://pytorch-geometric.readthedocs.io/en/latest/install/installation.html

本書では次のコマンドでインストールした。

pip install torch_geometric

pip install pyg_lib torch_scatter torch_sparse torch_cluster torch_spline_conv -f https://data.pyg.org/whl/torch-2.4.0+cu118.html


GNNの訓練に必要なGPUメモリ

グラフニューラルネットワークはモデルが軽量なものが多く、GPUメモリの量は小さくても使用ができる。8GBあれば問題はない。もっと少なくても動作は可能。


R

細胞コミュニティの解析は、空間シングルセルトランスクリプトーム解析と共にニーズが高まっており、特にトランスクリプトーム解析にはRが主流であったことから、本書でもRでの解析が登場する。

導入は広く解説されているので割愛させて頂くが、RをCRANからダウンロードし、RStudioから使用するのが一般的。

CRAN: https://cran.r-project.org/
RStudio: https://posit.co/download/rstudio-desktop/