空間トランスクリプトミクス解析：Monkeybreadチュートリアル解説

## はじめに：なぜ空間的な細胞解析が重要か

**従来の限界**

これまでの遺伝子発現解析では、主にバルクRNA-seqや単一細胞RNA-seq（scRNA-seq）が用いられてきました。しかしこれらの技術には重要な限界も存在します。

バルクRNA-seqでは組織をすりつぶして解析するため、組織全体の平均的な情報しか得られません。一方、scRNA-seqは個々の細胞の遺伝子発現を詳細に調べることができますが、細胞を組織から取り出す際に、その細胞が「組織のどこにいたか」という位置情報が失われてしまいます。

この空間情報の喪失は深刻な問題をもたらします。まず、細胞同士の実際の相互作用が分からなくなります。細胞は隣接する細胞とシグナル分子をやり取りして機能を調節していますが、位置情報がないと、どの細胞とどの細胞が実際に相互作用しているかを知ることができません。また、組織の構造と機能の関係も見えなくなります。例えば腫瘍では、中心部と辺縁部で細胞の性質が大きく異なることが知られていますが、空間情報がないとこのような場所による違いを捉えることができません。

さらに、細胞を取り巻く環境の影響を評価することも困難になります。細胞の振る舞いは周囲の環境に大きく左右されます。どのような細胞に囲まれているか、血管からどれくらい離れているかなど、これらの情報は空間的な文脈でのみ意味を持ちますが、従来の手法ではこうした情報が完全に失われてしまいます。

**空間解析で答えられる生物学的問い**

空間トランスクリプトミクス技術の登場により、組織切片上で遺伝子発現と位置情報を同時に取得できるようになりました。これにより従来は答えることが困難だった数多くの生物学的問いに答えられるようになります。

この解析から得られる知見は、基礎研究の深化だけでなく、臨床応用にも直結します。特にがん免疫療法の分野では、免疫チェックポイント阻害薬の効果が腫瘍微小環境の空間的構成に強く依存することが明らかになってきており、空間解析に基づいた患者層別化や新規治療標的の同定が進められています。

本記事では、空間トランスクリプトミクスデータを解析するためのPythonツールであるMonkeybreadを使って、これらの生物学的な問いにどのようにアプローチできるかを具体的に解説します。

### 空間トランスクリプトミクスの主な種類について

空間トランスクリプトミクスは、組織内のどこでどの遺伝子が働いているかを明らかにする技術ですが、そのアプローチの仕方によっていくつかの種類に大別できます。代表的には「画像ベースの手法」と「スポットベースの手法」があります。今回使用するMonkeybreadは**画像ベース**の手法であることに注意してください。

<details><summary>2種類の空間トランスクリプトームについての解説</summary>

**画像ベースの空間トランスクリプトミクス (Image-based Spatial Transcriptomics)**

**原理**

組織切片上で、特定のRNA分子を蛍光プローブなどで直接標識し、その光を顕微鏡で撮影（イメージング）することで、遺伝子の位置と種類を特定します。

**特徴**

この手法の最大の魅力は、細胞1つ1つ、あるいは細胞内部の分子まで見えるほどの高い空間解像度（単一細胞～サブセルラー解像度）を実現できることです。顕微鏡で直接観察するため、細胞の形状や大きさ、隣接する細胞との関係性も詳細に把握できます。

解析の過程では、画像情報から個々の細胞の境界を特定する「細胞セグメンテーション」という処理が重要な役割を果たします。これにより、検出された遺伝子情報を個々の細胞に正確に割り当てることが可能になります。

測定可能な遺伝子数については、初期の技術では数十～数百種類の限られた遺伝子（ターゲット遺伝子）を対象とすることが多かったものの、技術の進歩により現在では数千～数万種類の遺伝子を検出できる手法も登場しています。得られるデータには、個々の細胞ごとの遺伝子発現カウント、各細胞の正確なX, Y座標（場合によってはZ座標も含む）、そして細胞の面積、体積、バウンディングボックスなどの形態情報が含まれます。

**代表的な技術例**

現在市場で利用可能な主要な技術として、MERFISH (Multiplexed Error-Robust Fluorescence In Situ Hybridization) やその商用版であるVizgen社の**MERSCOPE**があります。また、NanoString社の**CosMx Spatial Molecular Imager**や10x Genomics社の**Xenium In Situ**も広く使用されています。その他にも、seqFISHやsmFISH (single molecule FISH) などの技術が研究現場で活用されています。

**どんな研究に向いているか**

この手法は、特定の細胞種が組織内のどこに微細に分布しているかを知りたい場合に特に有効です。また、細胞同士がどのように接触し、影響し合っているかを詳細に調べる研究、例えば免疫細胞とがん細胞の相互作用の解析などにも適しています。さらに、細胞の形状や状態と遺伝子発現の関係を解析したい場合にも威力を発揮します。

**課題・考慮点**

一方で、一度に解析できる組織の面積（視野）には限界があることを考慮する必要があります。ターゲット型の技術では、あらかじめ測定したい遺伝子を選定する必要があるため、未知の遺伝子の発見には向いていません。また、非常に高解像度なデータを扱うため、データ量が膨大になりがちで、解析には相応の計算資源が必要となります。

**スポットベースの空間トランスクリプトミクス (Spot-based Spatial Transcriptomics)**

**原理**

特殊なスライドガラス上に、位置情報を示すバーコードが付いた微小な領域（スポット）を多数配置します。このスライドに組織切片を乗せ、細胞を溶解させてmRNAを捕捉し、次世代シーケンサー（NGS）でバーコードと共にRNA配列を読むことで、どのスポット（位置）でどの遺伝子が発現していたかを網羅的に解析します。

**特徴**

この手法の最大の強みは、次世代シーケンサーを用いることで、原理的には細胞内で発現しているほぼ全ての遺伝子（トランスクリプトームワイド）を一度に測定できる**網羅的な遺伝子発現解析**が可能な点です。

スポットのサイズと配置はあらかじめ決まっており、直径は数マイクロメートルから100マイクロメートル程度です。解像度はこのスポットのサイズに依存し、小さなスポットでも複数の細胞を含むことが一般的です。そのため、厳密な単一細胞解像度ではないものの、組織領域ごとの遺伝子発現プロファイルを効率よく取得できます。得られるデータは、各スポットごとの網羅的な遺伝子発現カウントと、各スポットのX, Y座標情報から構成されます。

**代表的な技術例**

現在最も広く使用されている技術は、10x Genomics社の**Visium Spatial Gene Expression**です。この他にも、Slide-seqやDBiT-seqなどの技術が研究現場で開発・利用されています。

**どんな研究に向いているか**

この手法は、組織全体の遺伝子発現のパターンや傾向を大まかに把握したい場合に理想的です。特定の組織領域、例えばがん組織の中心部と周辺部で、どのような遺伝子が特徴的に発現しているかを網羅的に探索する研究に適しています。また、これまで知られていなかった遺伝子やパスウェイが、特定の組織構造や状態に関わっている可能性を探る探索的研究にも有用です。

**課題・考慮点**

スポット内に複数の細胞が含まれるため、得られるデータは混合した細胞集団の平均的な発現情報となります。単一細胞レベルでの情報を得るためには、計算によるディコンボリューション（混合データから個々の成分を推定する手法）などの追加解析が必要になる場合があります。また、細胞の正確な境界や形態を直接捉えることはできないため、細胞の形態学的な解析には適していません。

**まとめと比較**

| 項目           | 画像ベースの手法                                   | スポットベースの手法                                   |
| :------------- | :------------------------------------------------- | :----------------------------------------------------- |
| **基本原理** | 蛍光標識したRNAを直接イメージング                     | バーコード付きスポットでRNAを捕捉しNGS解析                 |
| **空間解像度** | 高い (単一細胞～サブセルラー)                        | スポットサイズに依存 (通常、複数細胞/スポット)              |
| **遺伝子網羅性** | ターゲット型が多い (網羅的手法も発展中)                 | 高い (トランスクリプトームワイド)                         |
| **細胞形態情報** | 取得可能                                           | 限定的                                                 |
| **主な利点** | 微細構造、細胞間相互作用の詳細解析                     | 組織全体の網羅的遺伝子発現プロファイリング                  |
| **主な課題** | 視野、ターゲット遺伝子選定 (一部除く)、データ量        | 単一細胞解像度の限界 (追加解析で補完の試みあり)           |

</details>

## Monkeybreadライブラリについて

Monkeybreadは、空間トランスクリプトミクスのデータを解析するために開発されたPythonツールです。空間トランスクリプトミクスが生み出す大量かつ複雑なデータを扱いやすくし、細胞の配置パターンや細胞間の相互作用を効率的に見つけ出すことを可能にします。特に画像ベースの空間トランスクリプトミクスデータの解析に適しており、高解像度な空間情報を活かした詳細な解析を実現します。

### Monkeybreadの特徴

を使うことで特定の細胞が組織内のどこに集まっているかという細胞の局在を明らかにできます。これは単に細胞の位置を示すだけでなく、なぜその場所に集まるのか、どのような生物学的意味があるのかを考察する出発点となります。

さらに重要なのは、隣り合う細胞同士がどのように影響し合っているかという細胞間相互作用の解析です。細胞は孤立して存在するのではなく、周囲の細胞と常にコミュニケーションを取りながら機能しています。Monkeybreadは、このような細胞間の相互作用を空間的な文脈の中で理解することを可能にします。

また、特定の細胞の周囲でどのような遺伝子が特徴的に働いているかという微小環境の解析も可能です。同じ種類の細胞でも、置かれた環境によって異なる振る舞いをすることがあり、この微小環境の違いを理解することは、細胞の機能や運命を決定するメカニズムの解明につながります。

## 本記事で扱う解析の全体像

### 全体像

```
【データ準備】
データ読み込み
    ↓
【空間分布解析】  
細胞密度推定
    ↓
【微小環境解析】
ニッチの同定 → ニッチのクラスタリング
    ↓
【相互作用解析】
共在在パターン → リガンド-受容体ペアの同定
    ↓
【統合的理解】
免疫微小環境の全体像
```

### 解析の流れ

最初に行うのは、データの読み込みと前処理です。空間トランスクリプトミクスでは、各細胞の遺伝子発現データと位置座標が得られます。組織全体では数十万個の細胞が含まれることもあるため、解析したい領域にズームインしてデータをサブセット化することで、計算を効率化し、局所的な現象に焦点を当てます。

データの準備ができたら、まず細胞の空間分布を理解することから始めます。特定の細胞タイプ、例えばCD8陽性T細胞やB細胞がどこに多く存在するかを、細胞密度として可視化します。これにより、免疫細胞の浸潤パターンや腫瘍細胞の分布の偏りが明らかになります。

次に、個々の細胞ではなく、細胞の集団としての振る舞いに注目します。ニッチ分析では、各細胞の周囲にどのような細胞が存在するかに基づいて、組織内の微小環境を同定します。
最後に、細胞間の直接的な相互作用を解析します。どの細胞タイプが空間的に共在しているかを調べ、意味のある近接関係を同定します。さらに、リガンド受容体解析により、物理的に近接した細胞間で起こりうるシグナル伝達を推定します。

## 環境とコード
実行環境は以下の通りです。
~~~
OS - 5.4.0-67-generic #75-Ubuntu SMP Fri Feb 19 18:03:38 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux
CPU - Intel(R) Core(TM) i9-10920X CPU @ 3.50GHz　×23
メモリ - 131.6 GB
~~~
コードは以下です。

https://github.com/immunitastx/Monkeybread/blob/main/docs/notebooks/tutorial.ipynb

##  データ

どのようなデータを扱っているか確認するためにコードを確認してみましょう。
~~~python
# Use gsutil to download the data from Google Cloud where the public data is stored
subprocess.run(
    "mkdir HumanMelanomaPatient2",
    shell = True
)
subprocess.run(
    "gsutil -m cp -n gs://vz-ffpe-showcase/HumanMelanomaPatient2/cell_by_gene.csv \
    gs://vz-ffpe-showcase/HumanMelanomaPatient2/cell_metadata.csv \
    ./HumanMelanomaPatient2/",
    shell = True
)
# データのロード
adata = mb.util.load_merscope("./HumanMelanomaPatient2")
~~~

ダウンロードしたデータは以下の形になっていることがわかります。
~~~
HumanMelanomaPatient2/
├── cell_by_gene.csv
└── cell_metadata.csv
~~~
それぞれのファイルの中身を見てみます。

* `cell_by_gene.csv`
~~~
<bound method NDFrame.head of           
        cell  PDK4  CX3CL1  PGLYRP1  CD4  SNAI2  TNFRSF17  ICAM3  TBX21  \
0            0   0.0     0.0      0.0  0.0    0.0       0.0    0.0    0.0   
1            1   0.0     0.0      0.0  0.0    0.0       0.0    0.0    0.0   
2            2   0.0     0.0      0.0  0.0    0.0       0.0    0.0    0.0   
3            3   0.0     0.0      0.0  0.0    4.0       0.0    0.0    0.0   
4            4   0.0     0.0      0.0  0.0    0.0       0.0    0.0    0.0   
...        ...   ...     ...      ...  ...    ...       ...    ...    ...   
207864  207864   0.0     1.0      0.0  0.0    1.0       0.0    0.0    0.0   
207865  207865   2.0     1.0      0.0  0.0   16.0       0.0    0.0    0.0   
207866  207866   0.0     1.0      0.0  0.0    7.0       0.0    0.0    0.0   
207867  207867   0.0     0.0      0.0  0.0    3.0       0.0    0.0    0.0   
207868  207868   0.0     0.0      0.0  1.0   11.0       0.0    0.0    1.0   

        FAP  ...  Blank-29  Blank-38  Blank-46  Blank-47  CXCR1  NOS2  WNT5A  \
0       0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0   
1       0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0   
2       0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0   
3       0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0   
4       0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0   
...     ...  ...       ...       ...       ...       ...    ...   ...    ...   
207864  0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0   
207865  0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0   
207866  0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0   
207867  0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0   
207868  0.0  ...       0.0       0.0       0.0       0.0    0.0   0.0    0.0  
~~~


この `cell_by_gene.csv` ファイルは、その名の通り、**各細胞（cell）が、どの遺伝子（gene）をどれだけ発現しているか**を示しています。

* **行 (row)**: 一つ一つの細胞を表しています。この例では、0番から207868番までの細胞IDが振られています。
* **列 (column)**: `cell`列以外は、測定対象となった遺伝子の名前（例: `PDK4`, `CX3CL1`など）や、実験の品質管理などに使われる可能性のあるブランクコントロール（例: `Blank-29`）を表しています。
* **値**: 各細胞における各遺伝子の発現量（カウント数）を示しています。例えば、3番の細胞では遺伝子`SNAI2`が`4.0`という値で検出されています。

これは**シングルセルRNAシーケンス解析（scRNA-seq）でよく見られる「カウントマトリックス（count matrix）」や「遺伝子発現マトリックス」に相当するもの**です。scRNA-seqでは、組織から細胞を一つ一つ分離して遺伝子発現を測定しますが、空間トランスクリプトミクスでは、この情報に加えて、**各細胞が組織内のどこに存在したか**という位置情報（`cell_metadata.csv`に含まれる）が重要になります。

つまり、`cell_by_gene.csv` は「どの細胞で、どの遺伝子が、どれくらい働いているか」という基本的な情報を提供するファイルであり、これと細胞の位置情報を組み合わせることで、組織構造と遺伝子機能の関連性を探ることが可能になります。

* `cell_metadata.csv`：各細胞の空間情報と属性データ

```
         Unnamed: 0   fov        volume     center_x      center_y        min_x  \
0            183932     0   1130.140916   152.274028   6249.049955   142.748725
1            183934     0    138.886608     3.877718   6243.144839     2.109583
2            183935     0    904.852155   202.958029   6247.245333   195.628544
3            183936     0    392.052611   134.961771   6257.559419   129.688256
4            183938     0    179.407528   133.132464   6254.298536   130.274355
...             ...   ...           ...          ...           ...          ...
207864        34219  1180    251.003761  7600.293902     16.622712  7597.543673
207865        34220  1180    321.750910  7717.710980     35.380081  7714.540597
207866        34221  1180    540.939982  7581.991590     55.118039  7578.434254
207867        34222  1180    128.577054  7564.624416    111.726298  7561.956910
207868        34225  1180    163.863254  7672.158707    151.154438  7669.524286

               max_x        min_y        max_y
0         161.799330  6239.801266  6258.298645
1           5.645852  6240.431482  6245.858195
2         210.287514  6242.523300  6251.967366
3         140.235286  6251.266833  6263.852005
4         135.990574  6251.846677  6256.750394
...              ...          ...          ...
207864  7603.044131    13.144266    20.101158
207865  7720.881363    31.421764    39.338398
207866  7585.548927    49.298724    60.937353
207867  7567.291923   109.716782   113.735813
207868  7674.793127   148.432741   153.876135
```

この `cell_metadata.csv` ファイルは、各細胞に関する追加情報、特に空間的な位置情報を格納しています。これにより、`cell_by_gene.csv` の遺伝子発現データと合わせて、組織内での細胞の振る舞いを理解することができます。

各列の意味は以下の通りと推測されます。

* **`Unnamed: 0`**: 元のデータフレームのインデックス、あるいは細胞を識別するための一意のID。`cell_by_gene.csv` の細胞IDと対応していると考えられます。
* **`fov` (Field of View)**: 視野。顕微鏡などで組織を観察・測定する際に、どの視野領域でその細胞が検出されたかを示す番号やID。大きな組織サンプルを複数の視野に分けて測定する場合に重要になります。
* **`volume`**: 細胞の体積。画像解析などから推定された細胞の大きさに関する情報です。
* **`center_x`, `center_y`**: 細胞の中心座標（X座標、Y座標）。組織切片上での細胞の2次元的な位置を示します。これが空間情報を代表する最も重要なデータの一つです。
* **`min_x`, `max_x`, `min_y`, `max_y`**: 細胞のバウンディングボックス（境界ボックス）の座標。細胞が占める領域のX座標の最小値・最大値と、Y座標の最小値・最大値を示します。これにより、細胞の形状や広がり、他の細胞との接触具合などをより詳細に把握できます。

この `cell_metadata.csv` に含まれる空間座標（特に `center_x`, `center_y`）と、`cell_by_gene.csv` の遺伝子発現情報を組み合わせることで、初めて空間トランスクリプトミクス解析が可能になります。「どの細胞が、どこにあって、どんな遺伝子を発現しているか」という情報が入っているファイルです。

##  細胞密度の推定と可視化

カーネル密度推定（KDE）を用いることで、特定細胞タイプの空間的な密度分布を連続的な分布として描画できます。この手法は、個々の細胞の位置情報から、組織全体における細胞密度の勾配を滑らかに推定します。

この解析の生物学的意義は、免疫細胞や腫瘍細胞など、特定細胞集団の集積や空白地帯の存在を客観的に明らかにできることにあります。例えば、T細胞が高密度に集積している領域は活発な免疫応答が起きている可能性を示唆し、逆に免疫細胞が存在しない領域は免疫抑制的な環境である可能性があります。このような密度分布パターンの理解は、細胞間相互作用の予測や微小環境の機能的異質性の発見につながります。

まずは今回登場する細胞について以下の表でみてみましょう。

| 英語表記     | 日本語表記・簡単な説明                                                                 |
|--------------|------------------------------------------------------------------------------------|
| B lineage    | B細胞系（抗体を産生し体液性免疫を担うリンパ球）                                     |
| CD4 T        | CD4陽性T細胞（免疫応答の司令塔として他の免疫細胞を活性化するT細胞）                |
| CD8 T        | CD8陽性T細胞（ウイルス感染細胞やがん細胞を直接攻撃・破壊するT細胞）                |
| Treg         | 制御性T細胞（免疫応答を抑制し恒常性を維持するT細胞サブセット）                      |
| NK           | NK細胞（感染細胞や腫瘍細胞を非特異的に攻撃する自然免疫の細胞）                     |
| pDC          | 形質細胞様樹状細胞（ウイルス感染時に大量のI型インターフェロンを産生する樹状細胞）   |
| cDC          | 古典的樹状細胞（抗原提示によりT細胞を活性化し獲得免疫を誘導する樹状細胞）           |
| mregDC       | 調節性樹状細胞（免疫応答を調節しバランス維持に関与する樹状細胞）                    |
| mono/mac     | 単球／マクロファージ（異物を貪食・消化し抗原提示や炎症制御も担う免疫細胞）          |
| granulocyte  | 顆粒球（主に細菌や寄生虫などの異物を貪食・排除する白血球の一群）                    |
| stromal      | 間質細胞（組織や臓器の構造を支え免疫細胞の活動を助ける支持細胞）                    |
| malignant    | 悪性細胞（異常増殖し周囲組織に浸潤・転移するがん細胞）                              |
| endothelial  | 内皮細胞（血管やリンパ管の内側を覆い物質交換やバリア機能を担う細胞）                |
| other        | その他（上記以外の細胞種）                                                           |

たくさんの細胞がありますが、チュートリアルではCD4 T, CD8 T, TregをまとめてT細胞とし、B細胞との位置関係について解析しています。
この解析の結果は以下です。

![スクリーンショット 2025-05-14 233149.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/2e257382-37cd-4be1-8395-e610208a2f4d.png)
B細胞の方が分散して存在していることがみてとれます。B細胞は液性免疫であり、病変部に対する場所特異性はあまりないことが予想されます。一方でCD4細胞もCD8細胞の様に直接細胞を攻撃することは基本的になく、病変部に対する場所特異性はあまりないことが考えられるためにCD4細胞の方がより特異的であるのは面白い結果だと思います。

一方CD8細胞は細胞性免疫で直接攻撃するという機能があり、病変の中あるいは近傍に存在していること考えられます。したがってB細胞とCD8細胞、CD4細胞とCD8細胞を比較するといずれの場合でもCD8細胞がより特異的に集まっていることが予想されます。
結果は以下です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/1f91969d-d64e-4b15-bca3-2bbae52b7b87.png)

これを見るとCD8細胞が最も分散しているようにみてとれます。予想と逆の結果となり面白いです。


##  ニッチ解析（細胞近傍プロファイルのクラスタリング）

### ニッチ計算
ここでは各細胞の近傍（周囲一定距離内）にどの細胞タイプがどれだけ存在するかをカウントし、それをクラスタリングして各クラスターを「ニッチ（niche）」と定義します。  
組織内では、特定の細胞タイプが集まって機能的なユニットを形成します。ニッチ解析により、空間的・機能的な細胞コミュニティの実体を推定できると考えられます。

チュートリアルでは免疫細胞について解析しています。免疫細胞がどのような微小環境（ニッチ）に存在するかについて調べているということです。

免疫細胞に対するニッチを計算していきます。
~~~python
# 免疫細胞のみを対象にフィルタリング
immune_mask = ~adata.obs['niche_subtypes'].isin(['malignant/other', 'endothelial', 'stromal'])
# 近傍プロファイルの計算（細胞ごとの隣接細胞タイプのカウント）
adata_neighbors = mb.calc.cellular_niches(
    adata,
    cell_type_key='niche_subtypes',
    radius=75,  # 近傍の定義範囲（腫瘍微小環境の典型的な相互作用距離）
    normalize_counts=True,   # 総細胞数で正規化（相対化）
    standard_scale=True,     # z-score化（細胞タイプ間のカウント差を平準化）
    clip_min=-5, clip_max=5, # 外れ値のクリッピング
    mask=immune_mask,        # 免疫細胞のみを解析対象
    n_neighbors=100,         # クラスタリング時の近傍数
    resolution=0.25,         # Leidenアルゴリズムの解像度
    min_niche_size=300,      # 最小ニッチサイズ
    key_added='niche'
)
~~~

つぎに免疫細胞の周囲に存在するB細胞のプロットをみてみましょう。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/026b27e2-ba1b-4c0d-aefd-be68e9f381b4.png)

X軸は各免疫細胞の周囲に存在するB lineage細胞の正規化/Zスコア値でY軸はその値を持つ免疫細胞の数となります。
例えば一番左の棒について考えるとB細胞についてのnormalized neighbor countsが-0.7くらいの免疫細胞が8500個くらいいることがわかります。

<details><summary>x軸のNormarized neighbor countsがマイナスとはどういうことか</summary>

まずはNormarized neighbor countsの定義から見てみましょう。
**Neighborhood counts（近傍カウント）：**
ある免疫細胞の周囲（半径75μm）に存在する特定の細胞タイプ（例：B細胞）の数をカウントした値。
ある免疫細胞を中心に、半径75μmの円内に存在する**すべての細胞**を「近傍細胞」と定義します。
この円内には、B細胞・マクロファージ・悪性細胞が含まれるとします。

この時点での「B細胞のNeighborhood count」は、単純に円内に存在するB細胞の数です。例えば、B細胞が3個あれば、カウントは「3」です。この値は必ず0以上の整数になります。

では、なぜX軸の「Normalized neighbor counts」はマイナスになるのでしょうか？これは、チュートリアルのコード中で指定された以下の2つの処理が関係しています。

1.  **`normalize_counts=True` （正規化）**
    まず、各免疫細胞について計算された「特定の細胞タイプ（例: B細胞）のNeighborhood count」を、その免疫細胞の周囲（半径75μm以内）に存在する「**全細胞タイプ**の総数」で割ります。

    * **計算例:**
        * ある免疫細胞Aの周囲75μm以内に、B細胞が3個、マクロファージが5個、悪性細胞が2個いたとします。
        * 近傍の全細胞数は 3 + 5 + 2 = 10個 です。
        * B細胞の正規化カウント = B細胞の数 / 近傍の全細胞数 = 3 / 10 = 0.3
        * マクロファージの正規化カウント = 5 / 10 = 0.5
        * 悪性細胞の正規化カウント = 2 / 10 = 0.2

    この正規化により、近傍の細胞密度が高い場所と低い場所の影響をある程度補正し、特定の細胞タイプが「どれくらいの割合で存在するか」という相対的な指標になります。この段階では、値は通常0から1の間（または0以上）になります。

2.  **`standard_scale=True` （標準化、Z-score化）**
    次に、上記で計算された「正規化されたカウント値」を、**解析対象となる全ての免疫細胞にわたって**標準化（Z-score化）します。Z-scoreは、以下の式で計算されます。

    $$Z = \frac{(個々の細胞の正規化カウント) - (全免疫細胞における正規化カウントの平均値)}{(全免疫細胞における正規化カウントの標準偏差)}$$


    * **Z-scoreがマイナスになる場合:**
        ある免疫細胞の周囲におけるB細胞の「正規化カウント」が、**「全免疫細胞の周囲におけるB細胞の正規化カウントの平均値」よりも低い**場合、Z-scoreはマイナスになります。
        つまり、その免疫細胞の周りには、**平均的な状況と比較してB細胞が相対的に少ない**ことを意味します。

    * **Z-scoreがプラスになる場合:**
        その免疫細胞の周囲におけるB細胞の「正規化カウント」が、平均よりも高い場合、Z-scoreはプラスになります。これは、その免疫細胞の周りに**平均よりも相対的に多くのB細胞が存在する**ことを意味します。

    * **Z-scoreが0に近い場合:**
        その免疫細胞の周囲におけるB細胞の「正規化カウント」が、全体の平均値に近いことを意味します。

したがって、プロットのX軸が「-0.7」ということは、その免疫細胞の周囲75μmには、**B細胞が平均よりも明らかに少ない環境である**ということを示しています。「B細胞があまりいない場所を好む（あるいは結果的にそのような場所にいる）免疫細胞群」と解釈できます。

この標準化処理により、異なる細胞タイプ間（例えば、元々数の多い細胞タイプと少ない細胞タイプ）でのカウント数の絶対的な差の影響を抑え、各細胞タイプがそれぞれの平均的な存在量と比べて、特定の免疫細胞の周囲で「相対的に濃縮しているか」「相対的に枯渇しているか」を比較しやすくしています。ニッチを定義する上で、このような相対的な偏りを捉えることが重要になるため、標準化が行われています。
</details>

したがってこのプロットは以下のように解釈できます。
**一部の免疫細胞（Normalized neighbor countsが５）の周りにB細胞がよく集まっている。**

### 近傍プロファイル（ニッチ）のクラスタリング
ここで重要なのは、必ずしも細胞の種類によってクラスター化されるわけではないことです。むしろ、近傍の構成がどれだけ類似しているかによってクラスター化されます。
例えば同じ「マクロファージ」（niche_subtypesがマクロファージ）でも、あるマクロファージはがん細胞に囲まれた環境に、別のマクロファージはB細胞に囲まれた環境にいるかもしれません。この場合、これらのマクロファージはUMAP上では離れた位置にプロットされる可能性があります。
逆に、異なる種類の免疫細胞（例：T細胞とNK細胞）でも、もし両者が同じような「がん細胞が多い」環境にいれば、UMAP上で比較的近い位置にプロットされるかもしれません。

ニッチごとに色分けした組織のプロットは以下です。
横軸（spatial1）と縦軸（spatial2）は組織上の位置を表しており、1つ1つの点が個々の細胞の位置を示しています。
グレーで表示されている領域は、悪性細胞（がん細胞）やその他の細胞が占める部分です。一方、0から9の番号は、それぞれ異なる種類の免疫細胞ニッチを表しています。
これらがどのような種類のニッチかは下で同定します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/bf4f2a2f-791b-4b76-be39-1d9322217a32.png)

次に各ニッチにおける細胞タイプのエンリッチメントのプロットを見ていきます。
**エンリッチメントとは**
エンリッチメント（濃縮/集積）は、特定のニッチにおいて、ある細胞タイプが組織全体での平均的な割合と比べて多いか少ないかを評価する指標です。
ある特定のニッチ（例：ニッチA）を取り出したとき、そのニッチAを構成している細胞たちの「細胞タイプ」の内訳を見ます。
そして、その内訳が、組織全体での細胞タイプの平均的な割合と比べて、特定の細胞タイプだけが際立って多くなっている場合、その細胞タイプは「ニッチAにエンリッチ（濃縮/集積）している」と言います。


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/26c34cb1-f705-4184-bf71-bd4a13242c3f.png)

例えばB細胞はニッチクラスター４にエンリッチしていることがわかります。
ニッチクラスターがどの細胞種に対応しているかを調べるためにチュートリアルではこのままDEG解析を行っています。
DEG解析とは異なるグループの細胞間で、遺伝子の発現量が大きく違う遺伝子（＝DEG）を見つける解析です。
グループ間の違いが、どの遺伝子の働きの違いによって生まれているのかを知ることができます。
例えば、「がん細胞の近く（ニッチA）」と「遠く（ニッチB）」にいる免疫細胞でDEG解析をすれば、がん細胞が免疫細胞の遺伝子発現をどう変えるかについての知見を得ることができます。

## 細胞種間の共在性の探索

組織の中で、異なる種類の細胞が「どこで、どのように一緒に集まっているか（あるいは互いに避けているか）」という空間的な配置パターンを調べることが共在性解析です。細胞はランダムに散らばっているわけではなく、特定の細胞同士が特定の場所でグループを形成しています。この共在のパターンを、細胞の位置情報と種類から明らかにします。

細胞種間の共在性を調べることで、組織内の機能的な構造が見えてきます。特定の細胞が集まって特別な役割を果たす微小環境やニッチが組織内には存在しており、共在性解析によってこれらの機能単位がどこにあり、どのような細胞から構成されているかが明らかになります。

また、病態の理解においても重要な知見が得られます。がんや炎症などの病変部では、健康な組織とは異なる細胞の集まり方をしていることがあり、この違いが病気の原因や進行の程度を反映している可能性があります。さらに、特定の細胞の空間的な配置パターンが治療効果と相関する場合、新しい治療法や創薬ターゲットの発見につながることも期待されます。

### 共局在のプロット

チュートリアルでは樹状細胞とCD4細胞の共局在についてプロットしています。この組み合わせを選択したのには、免疫応答の開始における両細胞の重要な役割があるからです。

**免疫応答における物理的接触の必要性**

樹状細胞とCD4 T細胞の相互作用は、効果的な免疫応答の開始に不可欠ですが、この相互作用には必ず物理的な接触が必要です。その理由を段階的に見ていきましょう。

まず、抗原提示の段階では、樹状細胞が取り込んだ抗原の断片をMHCクラスII分子に乗せて細胞表面に提示します。CD4 T細胞は、自身のT細胞受容体（TCR）を使ってこの抗原断片を認識しますが、この認識は分子レベルでの直接的な結合を必要とするため、両細胞が物理的に接触していなければ成立しません。

しかし、TCRによる抗原認識だけでは不十分です。CD4 T細胞の完全な活性化には、樹状細胞表面の共刺激分子（CD80、CD86など）とCD4 T細胞表面のCD28分子の結合という「第二のシグナル」が必要です。これも細胞膜上の分子同士の直接的な相互作用であり、細胞間の密接な接触なしには起こりません。

さらに、樹状細胞は活性化の際にサイトカイン（IL-12、IL-4、IL-6、TGF-βなど）を産生します。これらの分泌性因子は理論的には遠くまで拡散可能ですが、実際には産生源の近くで最も高濃度となり、効率的に作用します。このサイトカイン環境がCD4 T細胞の分化方向（Th1、Th2、Th17、Tregなど）を決定するため、両細胞の空間的近接性は機能的な意味を持ちます。

**共局在解析の生物学的意義**
このように、樹状細胞とCD4 T細胞の相互作用は複数のレベルで物理的接触を必要とするため、両細胞の共局在パターンを解析することで、組織内で実際に免疫応答が開始されている場所を特定できます。共局在が観察される領域は、活発な抗原提示とT細胞活性化が起きている「免疫応答のホットスポット」である可能性が高く、逆に共局在が見られない領域では、免疫応答が抑制されているか、異なるメカニズムが働いている可能性があります。

プロットを見てみましょう。
なおこのプロットは近接している(<20 μm)樹状細胞とT細胞をプロットしたものです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/ca4d59c9-f67f-4850-8378-5c3dbd5c2fbe.png)

空間的な配置を見ると、両細胞種は組織全体にランダムに散在しているのではなく、特定のパターンを示しています。特に注目すべきは、組織の辺縁部に沿って両細胞が分布していることです。このことはこの領域で活発な抗原提示とT細胞の活性化が起きている可能性を示唆しています。2つ前のプロットでこの両細胞がぬけている真ん中の部分が癌であることが明らかになっています。したがって免疫細胞が腫瘍の辺縁部に集積している典型的なパターンであることがわかります。

### 共局在に関する統計的検定

この解析では共局在がたまたま起こったものではなく意味のあるものであることを見ていきます。

**帰無仮説によるランダムシャッフルとの比較**
細胞タイプ間の共局在が偶然によるものか、生物学的に意味のあるパターンかを判定するため、統計的検定を行います。検定では「group1の特定の細胞タイプとgroup2の特定の細胞タイプの間には、特別な空間的な共局在（集まりやすさや避け合い）はない」という帰無仮説を設定します。
この帰無仮説のもとで期待される最短距離の分布（帰無分布）を生成するため、以下の手順を踏みます。group1の細胞タイプの位置は固定したまま、それ以外の細胞の位置情報をランダムにシャッフルし、その都度group1の各細胞から最も近いgroup2の細胞までの距離を計算します。この操作をn_perms回（ここでは100回）繰り返すことで、ランダムな配置での距離分布を得ることができます。
**p値による有意性の評価**
実際に観測された最短距離の分布と、設定した距離の閾値（threshold）を用いてp値を計算します。このp値は、group1の細胞から閾値以内の距離にgroup2の細胞が存在する頻度が、実際の観測値と同等またはそれ以上になる確率を表しています。言い換えれば、観察されたパターンがランダムな細胞配置でも生じうる確率を示しています。
p値が低い場合（慣例的に0.05未満）、観測された近接性は偶然では説明しにくく、2つの細胞タイプ間に生物学的に意味のある空間的な引き寄せ合い（共局在）が存在する可能性が高いと判断されます。逆にp値が高い場合は、観察されたパターンはランダムな配置でも十分起こりうることを意味し、特別な相互作用は存在しない可能性が示唆されます。

コードは以下です。
~~~python
# Compute distance from each myeloid cell type to the nearest T cell subtype
ct1_to_ct2_to_dists = mb.calc.shortest_distances_pairwise(
    adata_test, 
    'Cell Subtype', 
    MOMAC_SUBTYPES, 
    T_CELL_SUBTYPES
)

# Calculate significance of these shortest distances
ct1_to_ct2_to_perms, ct1_to_ct2_to_pval = mb.stat.shortest_distances_pairwise(
    adata_test, 
    'Cell Subtype', 
    MOMAC_SUBTYPES, 
    T_CELL_SUBTYPES,
    n_perms=100,
    observed=ct1_to_ct2_to_dists, 
    threshold=20
)
~~~

プロットは以下です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/3cfbc8cd-8179-4f2b-8fac-a702c972081c.png)

このヒートマップは、単球/マクロファージおよび樹状細胞の各サブタイプ（縦軸）と、T細胞の各サブタイプ（横軸）との間の空間的共局在の統計的有意性を示しています。各マスの数値はp値を表し、色の濃淡で視覚的に有意性の程度を把握できるようになっています。
* 有意な共局在パターン
黄色で表示されている低いp値（0.01～0.05）は、統計的に有意な共局在を示しています。特に注目すべきは、CXCL9を高発現する単球/マクロファージとCD4 T細胞（Tfh）の組み合わせでp値が0.01となっており、これらの細胞が偶然を超えて空間的に近接していることを示しています。同様に、mregDC（調節性樹状細胞）もCD4 T（Tfh）と有意な共局在を示しています（p=0.01）。
これらの共局在は生物学的に理にかなっています。CXCL9はT細胞を引き寄せるケモカインであり、これを産生する単球/マクロファージの近くにT細胞が集まることは、ケモカイン勾配に沿った細胞遊走を反映していると考えられます。また、Tfh（濾胞性ヘルパーT細胞）は抗体産生応答において重要な役割を果たすため、抗原提示細胞である樹状細胞との共局在も機能的な意味を持ちます。
* 共局在が見られないパターン
一方、濃い紫色で示される高いp値（1.00）は、ランダムな配置と区別がつかないことを意味します。CD163やSPP1を高発現する単球/マクロファージは、ほとんどのT細胞サブタイプとの間で有意な共局在を示していません。これらのマーカーは一般的にM2型（抗炎症性）マクロファージと関連しており、T細胞との直接的な相互作用よりも組織修復や免疫抑制的な機能を担っている可能性があります。

##  細胞間のリガンド-受容体相互作用解析

この統計検定の目的は、組織内の特定の細胞タイプ間（例：腫瘍細胞とT細胞）で観察されるリガンド・受容体（LR）ペアの発現に基づく相互作用が、単なる偶然によるものなのか、それとも生物学的に意味のある空間的な関連性に基づいているのかを統計的に明らかにすることです。
ここでは単球／マクロファージとT細胞間のエンリッチメントのP値を計算します。
~~~python
# 受容体のペア
lr_pairs = []
for g1, g2 in lrs:
    if g2 in ['CXCL9', 'CXCL10']:
        lr_pairs.append((g2, g1))
    else:
        lr_pairs.append((g1, g2))
        
# For each pDC, find the neighboring Tregs（注目したい細胞タイプペアを選ぶ）
cell_to_neighbors = mb.calc.cell_neighbors(
    adata,
    groupby='Cell Type',
    group1=['mono/mac'],
    group2=['CD8 T', 'CD4 T', 'Treg'],
    radius=20
) 

# 共発現スコアの計算
lr_pair_to_score = mb.calc.ligand_receptor_score(
    adata,
    cell_to_neighbors,
    lr_pairs=lr_pairs
)

# P値の計算
res = mb.stat.ligand_receptor_score(
    adata,
    cell_to_neighbors,
    actual_scores=lr_pair_to_score
)
~~~

結果のプロットは以下になります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/43903036-cc3b-48d9-9fbe-d6a629688964.png)

**X軸: 相互作用スコア**
詳しくは「統計検定の解説」のトグルの中で説明します。簡単には特定のリガンドと受容体のペアがでどれだけ強く共発現しているかの強度を数値化したものと考えればよいです。
X軸の右側にあるほど、相互作用スコアが高く、そのL-Rペアによる潜在的な相互作用が強い（またはより多くの細胞ペアで観察される）ことを意味します。

**Y軸: -log10(pval)（p値の対数変換値）**
これは、各L-Rペアの相互作用スコアが統計的にどれだけ有意であるかを示す指標です。
-log10変換であるため、Y軸の上側にあるほど、p値が小さく、そのL-Rペアの相互作用は統計的に有意であることを意味します。

**各点（データポイント）**
プロット上のそれぞれの点は、個別のリガンド・受容体ペアを表しています。点に付いているラベルがそのL-Rペアの名前です（例: CXCL9/CCR4）。
点のX座標がそのL-Rペアの「相互作用スコア」を、Y座標が「-log10(p値)」を示します。

**水平なグレーの線 (y ≈ 1.3)**
この線は、多くの場合、統計的な有意性の一般的な基準である **p値 = 0.05** のラインを示しています（-log10(0.05) ≈ 1.301）。
この線よりも上にある点は、p値が0.05未満であり、統計的に有意な相互作用を持つL-Rペアであると一般的に解釈されます。

**解釈**
プロットの右上に位置するL-Rペアは、相互作用スコアが高く（X軸が大きい）、かつ統計的有意性も高い（Y軸が大きい）ことを意味します。これらは、生物学的に重要で、かつ信頼性の高い相互作用の候補として最も注目すべきペアです。
例えば、図の最も右上に位置するCXCL9/CCR4は、相互作用スコアが約0.5と最も高く、-log10(pval)も約2.0（p値 ≈ 0.01）と非常に有意性が高いです。

<details><summary>統計検定の解説</summary>
この統計検定の解説では公式チュートリアルに参照するように書いてあった以下の図（のパネルB）を参照します。下のリンクが元論文です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/b1ccc0a2-f222-426b-9cc7-5d003e3a48de.png)

https://www.biorxiv.org/content/10.1101/2021.11.03.467020v3.full

### 統計検定のワークフロー

パネルに示されたリガンド・受容体相互作用の統計的有意性を評価するための5つのステップを解説します

**Step 1: Build Delaunay Network（ドロネーネットワークの構築）**
ドロネー三角形分割アルゴリズムを用いて、細胞の位置座標から細胞間を結ぶエッジ（辺）のネットワークを構築します。
ドロネーネットワークは細胞の位置情報に基づいて、どの細胞がどの細胞と物理的に「隣り合っている」可能性が高いか、つまり直接的な細胞間相互作用の候補となるかを、解析者の主観を排して客観的に定義するために用いられます。例えば、「半径〇µm以内の細胞」という基準では半径の大きさを決める恣意性が残りますが、ドロネーネットワークはデータ自身の構造から近接関係を導き出します。

**Step 2: Extract Linkages（リンケージの抽出）**

構築されたネットワーク（実際には距離ベースの近傍関係）から、解析対象となる特定の細胞タイプペアを抽出します。例えば、単球/マクロファージとT細胞サブタイプの組み合わせなど、生物学的に意味のある相互作用が期待される細胞ペアに焦点を絞ります。チュートリアルでは、`group1=['mono/mac']`と`group2=['CD8 T', 'CD4 T', 'Treg']`として、これらの細胞タイプ間の相互作用に注目しています

チュートリアルでは以下が対応します。
~~~python
# For each pDC, find the neighboring Tregs（注目したい細胞タイプペアを選ぶ）
cell_to_neighbors = mb.calc.cell_neighbors(
    adata,
    groupby='Cell Type',
    group1=['mono/mac'],
    group2=['CD8 T', 'CD4 T', 'Treg'],
    radius=20
) 
~~~
このようにチュートリアルではgroup1, group2として単球、T細胞が選ばれこの間の接続に注目することがわかります。しかしここでパラメーターに`radius=20`があるのはなぜでしょうか。ドロネーネットワークの良さを潰してしまうことになります。ソースコードを見てみましょう。
~~~python
"""
Calculate a ligand-receptor co-expression score between neighboring cells.
"""

import itertools
from typing import Dict, List, Set, Tuple, Union, Optional
from omnipath.interactions import import_intercell_network

import numpy as np
from anndata import AnnData


def ligand_receptor_score(
    adata: AnnData,
    cell_to_neighbors: Dict[str, Set[str]],
    lr_pairs: Optional[Union[Tuple[str, str], List[Tuple[str, str]]]] = None,
) -> Dict[Tuple[str, str], float]:
    """
    Calculates an average co-expression score of a 
    ligand-receptor pair between neighboring cells..

    Statistical test is as described in :cite:p:`He2021.11.03.467020` (See 
    Figure 4).

    Calculates scores for ligand-receptor pairs in neighboring cells.

    Parameters
    ----------
    adata
        Annotated data matrix.
    cell_to_neighbors
        A mapping of cells to their neighbors, as calculated by 
        :func:`monkeybread.calc.cell_neighbors`.
    lr_pairs
        One or multiple tuples corresponding to (ligand, receptor). If `None` then
        ligand/receptor pairs will be downloaded from Omnipath via the 
        :func:`omnipath.interactions.import_intercell_network`.

    Returns
    -------
    A mapping from ligand-receptor tuple pairs to pair scores.
    """
    # Convert to list if only one pair provided
    if lr_pairs is None:
        print("Using ligand/receptor pairs from omnipath...")
        interactions = import_intercell_network(
            transmitter_params={"categories": "ligand"},
            receiver_params={"categories": "receptor"}
        )
        genes_in_data = set(adata.var_names)
        lr_pairs = [
            (s, t)
            for s, t in zip(
                interactions['genesymbol_intercell_source'], 
                interactions['genesymbol_intercell_target']
            )
            if s in genes_in_data and t in genes_in_data
        ]
    elif isinstance(lr_pairs, tuple):
        lr_pairs = [lr_pairs]
    
    # Convert contacts dictionary to list of cell_id tuples
    neighbor_pairs = np.array(
        list(
            itertools.chain.from_iterable(
                itertools.product([cell], neighbors) 
                for cell, neighbors in cell_to_neighbors.items()
            )
        )
    )

    # Define score calculation function
    calculate_score = (
        lambda l, r: np.sum(np.sqrt(
            adata[neighbor_pairs.T[0], l].X.toarray() * adata[neighbor_pairs.T[1], r].X.toarray()
        ))
        / neighbor_pairs.shape[0]
    )

    lr_pair_to_score = {(ligand, receptor): calculate_score(ligand, receptor) for ligand, receptor in lr_pairs}

    return lr_pair_to_score


def ligand_receptor_score_per_niche(
    adata: AnnData,
    niche_to_cell_to_neighbors: Dict[str, Set[str]],
    lr_pairs: Optional[Union[Tuple[str, str], List[Tuple[str, str]]]] = None,
) -> Dict[str, Dict[Tuple[str, str], float]]:
    """
    Calculates an average co-expression score of a
    ligand-receptor pair between neighboring cells within each cellular niche
    calculated by :func:`monkeybread.calc.cellular_niches`. Statistical test is as described 
    in :cite:p:`He2021.11.03.467020` (See Figure 4). This function is a wrapper around 
    :func:`monkeybread.calc.ligand_receptor_score` and calls this function separately for 
    each niche.

    Parameters
    ----------
    adata
        Annotated data matrix.
    cell_to_neighbors
        A mapping of cells to their neighbors, as calculated by
        :func:`monkeybread.calc.cell_neighbors`.
    lr_pairs
        One or multiple tuples corresponding to (ligand, receptor). If `None` then
        ligand/receptor pairs will be downloaded from Omnipath via the
        :func:`omnipath.interactions.import_intercell_network`.

    Returns
    -------
    A dictionary mapping each niche to a sub-dictionary mapping each ligand-receptor 
    pair to its co-expression score.
    """

    niche_to_lr_pair_to_score = {}
    for niche, cell_to_neighbors in niche_to_cell_to_neighbors.items():
        
        lr_pair_to_score = ligand_receptor_score(
            adata,
            cell_to_neighbors,
            lr_pairs,
        )
        niche_to_lr_pair_to_score[niche] = lr_pair_to_score
    return niche_to_lr_pair_to_score

~~~
ドロネーネットワークはどこにも実装されていません。確かに公式ーチュートリアルには統計検定に関しては参照することと書いてあったのでこの部分はまだ？実装されていないということでしょう。

**Step 3a: リンケージの方向付け**
リガンド-受容体相互作用には明確な方向性があります。リガンドを産生する細胞から受容体を発現する細胞へのシグナル伝達という一方向の流れです。この方向を決めます。
チュートリアルでは以下が対応します。
~~~python
# 受容体のペア
lr_pairs = []
for g1, g2 in lrs:
    if g2 in ['CXCL9', 'CXCL10']:
        lr_pairs.append((g2, g1))
    else:
        lr_pairs.append((g1, g2))
~~~
これはLRペアのリストを準備している段階であり、個々の細胞間リンケージに方向性を与える処理そのものではない点に注意が必要です。しかし、方向性はスコア計算時にリガンド側細胞と受容体側細胞を区別することで暗黙的に考慮されます。
**Step 3b & 4b: LR相互作用スコアの計算**
**計算の詳細**:
1.  対象とする各「方向付けされたリンケージ $i$」（例：リガンド供給細胞 $cell_L$ から受容体発現細胞 $cell_R$ へ）を考えます。
2.  このリンケージ $i$ について、リガンドの発現量 $L_{cell_L}$ と受容体の発現量 $R_{cell_R}$ を用いて、個別のリンケージスコア $s_i = \sqrt{L_{cell_L} \times R_{cell_R}}$ を計算します。この平方根を取る形式は、リガンドと受容体の両方がバランス良く発現している場合にスコアが高くなりやすく、片方だけ極端に高くてもスコアの上昇が抑えられる効果があります。発現量が0の細胞が関与する場合、そのリンケージのスコアは0になります。
3.  全ての対象リンケージ（総数 $N_{linkages}$）について計算された個別のスコア $s_i$ を合計し、リンケージの総数で割ることで、特定のLRペアに関する平均的な相互作用強度（LR Interaction Score）を算出します
 $$ \text{LR Interaction Score} = \frac{\sum_{i=1}^{N_{\text{linkages}}} s_i}{N_{\text{linkages}}} $$

この計算により実際の組織サンプルにおいて、注目している細胞タイプ間で、特定のリガンドと受容体を介した潜在的なコミュニケーションが平均してどれくらいの強度で起こっているかを定量的に示します。4bでは帰無仮説（”観測されたLRペアの相互作用は偶然であり、細胞の空間配置やLRの発現パターンには特別な関係性はない”という仮説）でのランダムな状況でスコアが計算されます。

**Step 4a: リンケージのランダム化**
帰無分布を作成するために、元のデータの構造を一部ランダムにシャッフルします。具体的には、細胞のラベル（種類）をランダムに入れ替えたり、リンケージをランダムに再接続したり、遺伝子発現値をシャッフルするなどの方法が考えられます。このランダム化されたデータに対して、ステップ4bと同様にLR相互作用スコアを計算します。この操作を多数回（例: 1000回、10000回など）繰り返すことで、偶然によって生じうるスコアの分布（帰無分布）が得られます。

**Step 5: Calculate p-value（p値の計算）**
最後に、ステップ3bで計算された「実際に観測されたLR相互作用スコア」が、ステップ4aと4bで作成した「ランダムな場合に期待されるスコアの分布（帰無分布）」の中で、どれくらい稀な値なのかを評価します。
具体的には、観測スコア以上に極端なスコアが、帰無分布の中で出現する割合を計算し、これを**p値**とします。図では、観測スコア（Interaction Score）が帰無分布の右側の非常に端に位置しており、これはp値が小さいことを示唆しています。
p値が事前に設定した有意水準（例えば0.05）よりも小さい場合、「帰無仮説は棄却され、観測されたLR相互作用は偶然では説明できず、統計的に有意である」と結論づけられます。つまり、そのLRペアを介した細胞間相互作用が、その組織の特定の場所で実際に起こっている可能性が高いと考えられます。

</details>
共発現スコアのプロットは以下です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3566009/eb7c838f-b306-40b4-bed2-48ddbf221b8c.png)

このプロットは、各ニッチにおけるCXCL9（ケモカイン）とCXCR3（その受容体）の共発現スコアを示しています。

X軸はニッチ番号でY軸は上のトグル内で説明したリガンド/受容体スコアです。

最も高いスコアを示すのはニッチ9（約0.42）で、次いでニッチ4（約0.38）となっています。一方、最も低いのはニッチ10（約0.10）で、ニッチ6（約0.21）も比較的低い値を示しています。この4倍以上の差は、ニッチによってCXCL9/CXCR3を介した細胞間コミュニケーションの活発さが大きく異なることを示しています。

CXCL9は炎症性ケモカインの一つで、主に活性化されたマクロファージや樹状細胞から産生されます。その受容体であるCXCR3は、エフェクターT細胞、特にTh1細胞やCD8陽性T細胞に高発現しています。このリガンド-受容体ペアは、炎症部位へのT細胞の遊走と集積において中心的な役割を果たします。

高スコアを示すニッチ9やニッチ4は、活発な抗腫瘍免疫応答が起きている領域である可能性が高いと考えられます。これらのニッチでは、CXCL9を産生する抗原提示細胞とCXCR3を発現するエフェクターT細胞が空間的に近接し、かつ両分子が高発現していることを示しています。

対照的に、低スコアのニッチ10（悪性細胞主体）やニッチ6は、免疫学的に比較的静穏な領域、あるいは免疫抑制的な環境である可能性があります。特にニッチ10の低スコアは、腫瘍細胞周辺でのT細胞遊走シグナルが限定的であることを示唆しています。

## まとめ：解析の意義

ここで紹介したMonkeybreadを用いた解析により、メラノーマ組織における免疫微小環境の複雑な空間構造が明らかになりました。従来の単一細胞解析では失われていた空間情報を保持したまま解析することができました。

特に重要な発見は、組織内に機能的に異なる複数のニッチが存在し、それぞれで特異的な細胞間相互作用が展開されていることです。免疫活性化ニッチではCD80/CD28共刺激シグナルとCXCL9/CXCR3ケモカインシグナルが協調的に働き、効果的な抗腫瘍免疫応答を形成している一方、免疫抑制的ニッチではこれらのシグナルが著しく低下していることが定量的に示されました。

Monkeybreadの強みは、観察された空間パターンが偶然ではなく生物学的に意味があることを統計的に検証できる点にあります。樹状細胞とT細胞の共局在解析では、物理的近接性が統計的に有意であることが示され、実際に抗原提示が起きている場所を客観的に同定できました。また、リガンド-受容体解析により、空間的近接性と分子発現の両方を考慮した真の相互作用ネットワークを描き出すことができました。

免疫チェックポイント阻害薬などの免疫療法において、治療効果を予測するバイオマーカーの開発は喫緊の課題です。Monkeybreadによる解析は、単なる免疫細胞の数や種類だけでなく、それらの空間的配置と相互作用パターンが治療応答性を規定する可能性を示しています。

例えば、CXCL9/CXCR3スコアが高いニッチが優勢な患者では既存の免疫療法が有効である可能性が高く、逆に免疫抑制的ニッチが優勢な患者では、まず免疫細胞の浸潤を促進する前処置が必要かもしれません。このような空間的視点に基づいた治療戦略の個別化が可能になります。

(この記事は研究室インターンで取り組みました：https://kojima-r.github.io/kojima/)