Chapter 13無料公開

🐅1️⃣ 【 UTAG - 隣接細胞の特徴量も使って細胞のクラスタリング - 】

2024.11.19に更新

このチャプターの目次

【 UTAG 】
【 UTAG Tutorial 】

【 UTAG 】

原著（Nature Method, 2022）: https://www.nature.com/articles/s41592-022-01657-2
GitHub: https://github.com/ElementoLab/utag
チュートリアル：https://github.com/ElementoLab/utag/blob/main/documentation/UTAG Tutorial.ipynb
クラスタリング後の解析例（非常に参考になる）: https://github.com/ElementoLab/utag/blob/main/documentation/IMC Healthy Lung.ipynb

既にアノテーション付けされた細胞のコミュニティを見つけるのではなく、入力特徴量と細胞位置グラフで新たなクラスター（アノテーション）を見つける。マーカー発現や細胞形態といった特徴量行列を隣接行列とかけて、1つの細胞に周囲の細胞の情報を加える。その結果を基にクラスタリングへ。クラスタリング手法は既報のLeidenやPARCなどを使用したとのこと。

ScanPyを内部で使う。Rでもできる。

原著 Fig.1

Python環境

GitHubにはPython 3.9でtestしたとある。本書では3.12で行ってみたが次のエラーが出たので言われた通りPython 3.9を使用した。

AttributeError: module 'pathlib' has no attribute '_windows_flavour'
Python version 3.12ではエラーが出るようだ。
https://github.com/devopshq/artifactory/issues/430

仮想環境作り

version 3.9のpython環境を用意する。

conda create --name utag python=3.9

 JupyterLabを使うならインストールしておく。
conda install -c conda-forge jupyterlab
conda install conda-forge::jupyterlab-variableinspector
pip install --upgrade ipywidgets

UTAGのインストール

方法１） `pip install`で一発で導入

次のコマンドによりGitHub内のutagパッケージおよび関連パッケージがインストールされる。

pip install git+https://github.com/ElementoLab/utag.git@main

方法２）パッケージフォルダをダウンロード + 関連パッケージのインストール

GitHub内の「utag」フォルダがメインのパッケージなので、これをローカルにダウンロードして使用してもよい。（git cloneしておくとチュートリアルのipynbファイルが入った「documentation」フォルダも用意できる。）
pythonで解析する際に「utag」フォルダにアクセスしやすいような場所に置いておく。

こちらに依存パッケージ類の情報がある。深層学習モデルでは無いのでPytorchなどは不要。シングルセル解析に良く使用されるScanPyと、空間オミクス解析ライブラリのSquidPyが必要となる。

https://github.com/ElementoLab/utag/blob/main/setup.cfg

以下の順で入れれば使用できた。（その他の関連パッケージも同時にインストールされる。）

pip install scanpy
pip install squidpy
pip install parmap
pip install setuptools-scm
pip install parc

デモデータ

Zenodoにデータがdepositされている。
https://zenodo.org/records/6376767

GitHubの「documentation」フォルダに幾つかのチュートリアルのipynbファイルがある。各がん種ごとの解析例があり、Zenodoで配布されているScanPyのanndataを事前にダウンロードして進めるよう想定されている。
「UTAG Tutorial.ipynb」ではノートブックの中でZenodoからデモデータをダウンロードするので事前にダウンロードしておく必要はない。

【 UTAG Tutorial 】

GitHubの「documentation」フォルダにある「UTAG Tutorial.ipynb」ファイルを流してみる。

0. パッケージのimport

git cloneでパッケージを取得した場合は、importしやすいように「utag」フォルダにパスを通しておく。（pip installでutsgを取得した場合はパスが通った場所にインストールされている。）

import sys
sys.path.append(r"C:\Users\...\utag")

from utag import utag
import scanpy as sc
import squidpy as sq
import anndata
import numpy as np
import matplotlib.pyplot as plt

1. デモデータの読み込み

ZenodeのリンクからScanPyのanndataをダウンロードする。
※ ScanPyのAnnDataの扱いについては詳細には説明しない。

adata = sc.read(
    'data/healthy_lung_adata.h5ad',
    backup_url='https://zenodo.org/record/6376767/files/healthy_lung_adata.h5ad?download=1'
)

71946細胞×33マーカーの発現マトリクスを含むデータである。

データの中身確認

X属性の発現マトリクスを確認してみる。
発現マトリクスは平均が0でSDが1なので標準化されたデータであることがわかる。値の最小値・最大値は-3・3にclipされているようだ。

発現マトリクスを確認

var属性を見てみると33のマーカーの内、最後の5つは細胞の形状に基づくパラメーターであった。（これらも標準化されている。）

var属性

obs属性から1つのAnnDataのメタデータが確認できた。各細胞の由来サンプルやXY座標、ROI IDなどの情報が得られる。

obs属性

その他のobsm、obsp、varm、uns属性にはデータ形状が不確定のものを入れることが多い。これらは辞書型のようにadata.属性[キー]のようにアクセスする。

このキーの名前は自由に設定できるが、ScanPyの解析ワークフローに従うと、キーの名前はScanPyが用意したものになる。ScanPyのPlotなどでは特定のキーを基にデータを探すので従っておいた方が無難。

UMAP

obsm属性のX_umapキーにUMAP情報も入っていた。ScanPyのscanpy.pl.umap()コマンドではこのキーのUMAP値を参照して次元削減データが視覚化できる。ドットの色をサンプルで分けてみると、4サンプルの細胞を統合処理などした後に見える。

sc.pl.umap(adata, color="sample")

細胞のアノテーションも既に付いている。

sc.pl.umap(adata, color=["cell type"])

Spatial plot

scanpy.pl.spatial()コマンドは各細胞をXY座標に投影できる。obsm属性のspatialキーを自動で参照する。

adata.obsm["spatial"]に細胞のXY座標がある

デモデータのAnnDataは26個のROIをひとまとめにしているので、その内一つを取り出して可視化してみる。

# ユニークなROI ID一覧取得
roi_ids = np.unique(adata.obs["roi"])

# 1つ目のROIのAnnData取得
tmp_adata = adata[ adata.obs["roi"] == roi_ids[0]]
# plot
sc.pl.spatial(tmp_adata, color="cell type", spot_size=10)

2. グラフ作りの距離閾値を事前確認

グラフ作り自体はutag()機能の内部で行われるのだが、どこまで遠い細胞を接続有りとみなすかの閾値決めは事前にやっておく。

グラフ作りはSquidPyパッケージのsq.gr.spatial_neighbors()機能で行われる。kNNグラフ、Radiusグラフ、ドロネーグラフに対応しており、その内部では本書でも紹介したようなsklearn.neighbors.NearestNeighborsやscipy.spatial.Delaunayが使われている。

https://squidpy.readthedocs.io/en/stable/api/squidpy.gr.spatial_neighbors.html

UTAGではRadiusグラフを使用する。

まずは1症例取り出してグラフを可視化する。

# 1つのROIの細胞だけのAnnData取り出し
tmp_adata = adata[adata.obs["roi"] == roi_ids[1]]

fig, ax = plt.subplots(1, 1, figsize = (5,5), dpi = 100)

sc.pl.spatial(tmp_adata, color = 'cell type', spot_size = 7, 
              ax = ax, frameon = False,
              edges = False)

グラフを作る前の細胞位置を確認

sq.gr.spatial_neighbors()機能は入力に渡したAnnDataを直接更新する。返り値は無い。

max_dist = 15

sq.gr.spatial_neighbors(tmp_adata, radius=max_dist,
                        coord_type="generic")

sc.pl.spatial()機能で可視化する。

fig, ax = plt.subplots(1, 1, figsize = (5,5), dpi = 100)

sc.pl.spatial(tmp_adata, color = 'cell type', spot_size = 7, ax = ax, 
              title = 'Spatial Connectivity', frameon = False, 
              edges = True, neighbors_key = 'spatial_neighbors',
              edges_width = 0.5, edges_color = 'black')

距離閾値15では隣り合う細胞の接続にとどまっており、隣の隣までは繋がっていない距離間。

距離閾値15のRadiusグラフ

ちなみに距離閾値30も見てみる。

tmp_adata = adata[adata.obs["roi"] == roi_ids[1]]

max_dist = 30

sq.gr.spatial_neighbors(tmp_adata, radius=max_dist)

fig, ax = plt.subplots(1, 1, figsize = (5,5), dpi = 100)

sc.pl.spatial(tmp_adata, color = 'cell type', spot_size = 7, ax = ax, 
              title = 'Spatial Connectivity', frameon = False, 
              edges = True, neighbors_key = 'spatial_neighbors',
              edges_width = 0.5, edges_color = 'black')

隣の隣まで繋がることで接続数が増大する。ただ実際は距離閾値30が30 μmだとすると細胞間相互作用が無いとは言えないので幾つかの閾値で試してみると良いだろう。

距離閾値30のRadiusグラフ

UTAGのグラフは自己ループ有りのグラフが計算に使用される。。sq.gr.spatial_neighbors()のset_diag=Trueとする。距離閾値の確認の時には見づらくなるので自己ループ無しでよい。

tmp_adata = adata[adata.obs["roi"] == roi_ids[1]]

max_dist = 15

sq.gr.spatial_neighbors(tmp_adata, radius=max_dist, set_diag=True)

fig, ax = plt.subplots(1, 1, figsize = (5,5), dpi = 100)

sc.pl.spatial(tmp_adata, color = 'cell type', spot_size = 7, ax = ax, 
              title = 'Spatial Connectivity', frameon = False, 
              edges = True, neighbors_key = 'spatial_neighbors',
              edges_width = 0.5, edges_color = 'black')

自己ループありのグラフ

余談だが、neighbors_key = 'spatial_neighbors'が無いと、"neighbors"が採用され、これはUMAPの計算で作られた近傍細胞となる。

3. 必要なデータのみに削減

UTAGに必要なのはマーカー発現データと座標のみである。細胞アノテーション情報も使用しない。
AnnDataから不要な情報を除いておく。

# obs属性からROI IDと細胞のXY座標を取得
obs = adata.obs[['roi', 'X_centroid', 'Y_centroid']]
X = adata.X
var = adata.var

minimal_adata = anndata.AnnData(
    X = np.array(X).astype(np.float64),
    obs = obs,
    var = var
)

obsm属性は辞書型のように扱える。keyをspatial、valueをXY座標のnumpy.arrayとして保存しておく。

minimal_adata.obsm['spatial'] = np.array(minimal_adata.obs[['Y_centroid', 'X_centroid']])

さらに不要な情報を削除しておく。var属性はマーカー名だけの情報になる。

del minimal_adata.var['mean']
del minimal_adata.var['std']
del minimal_adata.obs['Y_centroid']
del minimal_adata.obs['X_centroid']

var属性マーカー名だけになった

minimal_adata

obs属性：ROI IDのみ

X属性：発現マトリクスは元と同じ

4. UTAGの実行

utag()機能で、グラフ作り、接続細胞間の特徴量の伝播、次元削減、クラスタリングまで一気に行われる。チュートリアル通り進めてみる。（引数の解説は後述する。）

1サンプル処理

1サンプル（1領域）のみで行うには、目的サンプルのみのAnnDataを用意し、utag()機能に渡す。
デモデータでは数秒で終了した。

single_slide = minimal_adata[minimal_adata.obs['roi'] == minimal_adata.obs['roi'].unique()[2]].copy()

# Run UTAG on provided data
utag_single_results = utag(
    single_slide,
    slide_key=None,
    max_dist=15,
    normalization_mode='l1_norm',
    apply_clustering=True,
    clustering_method = 'leiden', 
    resolutions = [0.3]
)

返り値のオブジェクトもAnnDataで、各属性に主成分分析、接続情報（neighbors、distances、connectivities）、クラスタリング結果（obs['UTAG Label_leiden_0.3']）などが追加されている。

sc.pl.spatial()でクラスタリング結果をXY座標上に投影して確認できる。

sc.pl.spatial(utag_single_results, 
              color = 'UTAG Label_leiden_0.3',
              spot_size = 10, 
              title = '',
              frameon = False # 枠無し
             )

複数サンプル同時処理

同時に比較解析したいサンプルは1つのAnnDataに含んだままでUTAGを実行することで、サンプル間（領域間）で統一されたクラスター番号が得られるので、サンプル（領域）によりどのクラスターが多い・少ないを述べることができる。

slide_key="roi"オプションでサンプル（領域）を分ける情報を持ったキー名を指定する。

utag_results = utag(
    minimal_adata,
    slide_key="roi",
    max_dist=15,
    normalization_mode='l1_norm',
    apply_clustering=True,
    clustering_method = 'leiden', 
    resolutions = [0.3]
)

デモデータでは3分ほどで終了した。AnnDataの各属性に主成分分析、接続情報（neighbors、distances、connectivities）、クラスタリング結果（obs['UTAG Label_leiden_0.3']）などが追加されている。

チュートリアル通りのplot機能でクラスタリング結果をXY座標上で確認する。

def visualize_rois(adata, n_rois = 3, roi_key = 'roi', color_key = 'cell type', dpi = 300):
    # increase dpi for a more crisp visualization
    fig, axes = plt.subplots(1, n_rois, figsize = (n_rois * 3 + 1,3), dpi = 100)

    for i, roi in enumerate(adata.obs[roi_key].unique()[:n_rois]):
        a = adata[adata.obs[roi_key] == roi].copy()

        if i == n_rois - 1:
            sc.pl.spatial(a, color = color_key, spot_size = 10, ax = axes[i], title = '', frameon = False, show = False)
        else:
            sc.pl.spatial(a, color = color_key, spot_size = 10, ax = axes[i], title = '', frameon = False, legend_loc = 'none', show = False)
    plt.tight_layout()
    plt.show()

# increase dpi for a more crisp visualization
visualize_rois(utag_results, color_key = 'UTAG Label_leiden_0.3', dpi = 100)

左2つの領域は似たクラスターで構成されており、右の領域は区画を分けたクラスターが幾つも存在することが分かる。

`utag()`機能の引数

1つの機能で複数の処理を行うため、引数が多い。

adata: AnnData
解析に用いるAnnData。obsm["spatial"]に細胞の座標がある前提。
channels_to_use: Optional[Sequence[str]]
解析に用いるマーカー（変数）名。何も指定しなければ全ての変数を使用する。
slide_key: {str, None}
utag()機能では複数のサンプル・領域を扱う際に複数のAnnDataを渡すのではなく、1つのAnnDataに統合してから渡す。サンプル・領域を分けるのに必要なobs属性のキー名を指定する。デフォルトでslide_key="Slide"になっている。1サンプル・領域のAnnDataであればNoneとすればよい。
save_key: str
クラスタリング結果をAnnDataに保存する際のキーの名前。"{save_key}_{method}_{resolution}"という名前になる。デフォルトは"UTAG Label"。
filter_by_variance: bool
分散で変数をフィルタリングするかいなか。デフォルトはFalse。Trueにすると全変数のSDの中央値より大きなSDを持つ変数だけにフィルタリングされる。（なので入力特徴量が半分ぐらいになる。）
max_dist: float
グラフ作りの際に接続ありとする最大距離。obsm["spatial"]の座標の単位がpxなのかμmなのか、どんな解像度（μm/px）なのかで閾値は大きく変わるので自身のデータに合わせて指定する必要がある。
デフォルトは20。これはimaging mass cytometryで1 μm/pxの解像度のデータで良い結果が得られた閾値とのこと。
normalization_mode: str
隣接行列の正規化方法。デフォルトはnormalization_mode="l1_norm"。sklearn.preprocessing.normalize()が実行される。
これ以外の正規化方法は用意されていないみたい。
keep_spatial_connectivity: bool
グラフ作りの結果得られた隣接行列や距離行列を残すか否か。slide_key=Noneの場合は、必ず残る。
obsp["spatial_connectivities"]やobsp["spatial_distances"]として記録される。
pca_kwargs: Dict[str, Any]
接続細胞間の特徴量の伝播した後の情報を使って主成分分析を行う。その際に内部で使うscanpy.pp.pca()機能に渡す引数の辞書。デフォルトでpca_kwargs={"n_comps":10}。主成分10個になっているので、入力特徴量の次元数が多いときは増やしても良いかも。
apply_umap: bool
UMAPを行うか否か。Trueの場合、scanpy.tl.umap()を行う。デフォルトでFalse。
umap_kwargs: Dict[str, Any]
UMAPを行う場合にscanpy.tl.umap()機能に渡す引数の辞書。
apply_clustering: bool
クラスタリングを行うか否か。FalseだとPCAもUMAPも行われない。デフォルトはTrue。
clustering_method: Sequence[str]
クラスタリングに使用する方法。Leidenアルゴリズム（scanpy.tl.leiden()）、PARC（PARCパッケージ）、KMeans（sklearn.cluster.KMeans）が使用できる。デフォルトはclustering_method = ["leiden", "parc", "kmeans"]で、3種類のアルゴリズムでクラスタリングされる。
resolutions: Sequence[float]
クラスタリング解像度に関するパラメーター。LeidenとPARCには対応する引数に渡される。KMeansではint(np.ceil(resolution * 10))の値が渡される。デフォルトでresolutions=[0.05, 0.1, 0.3, 1.0]で、複数のクラスタリング解像度で実施される。
leiden_kwargs: dict[str, Any]
Leiden法でクラスタリングする場合にscanpy.tl.leiden()に渡す追加の引数の辞書。
parc_kwargs: dict[str, Any]
PARC法でクラスタリングする際にPARC()インスタンス化に渡す追加の引数の辞書。
parallel: bool
接続細胞間の特徴量の伝播の処理を並列実行するかどうか。デフォルトでTrue。
processes: int
並列実行数。デフォルトはprocess=-1でフルで並列実行を行う。