Chapter 13無料公開

🐅1⃣ 【 UTAG - 隣接现胞の特城量も䜿っお现胞のクラスタリング - 】

Ryota Chijimatsu
Ryota Chijimatsu
2024.11.19に曎新

【 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」フォルダにアクセスしやすいような堎所に眮いおおく。

https://github.com/ElementoLab/utag/tree/main/utag

こちらに䟝存パッケヌゞ類の情報がある。深局孊習モデルでは無いのでPytorchなどは䞍芁。シングルセル解析に良く䜿甚されるScanPyず、空間オミクス解析ラむブラリのSquidPyが必芁ずなる。

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

https://scanpy.readthedocs.io/en/stable/
https://squidpy.readthedocs.io/en/stable/



以䞋の順で入れれば䜿甚できた。その他の関連パッケヌゞも同時にむンストヌルされる。

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でフルで䞊列実行を行う。