CLAMは、スライド画像からアノテーションせずに、特徴点を抽出し、クラス分類を行う弱教師あり学習です。
実行環境

構成	実行環境
OS	linux(Ubuntu18.04 LTS)
GPU	TITAN RTX （なくても実行できます。）
python	3.7.7

ライブラリ

パッケージ	バージョン
h5py	2.10.0
matplotlib	3.1.1
numpy	1.18.1
opencv-python	4.1.1.26
openslide-python	1.1.1
openslides	3.0
pandas	0.25.3
pillow	7.0.0
PyTorch	1.3.1
scikit-learn	0.22.1
scipy	1.4.1
tensorflow	1.14.0
tensorboardx	1.14.0
torchvision	0.1.8
smooth-topk	setup.py実行

（smooth-topkはCLAM-masterのsmooth-topk/setup.pyがあるので、ヒートマップの作成前までに実行すればいいです。）
smooth-topk: https://github.com/mahmoodlab/CLAM/blob/master/docs/INSTALLATION.md
ディレクトリ構成

CLAM
├── DATA_DIRECTORY/
│     ├── slide_1.svs
│     ├── slide_2.svs
│     └── ...
│
├── dataset_csv
│   ├── tumor_subtyping_dummy_clean.csv
│   └── tumor_vs_normal_dummy_clean.csv
├── heatmaps
│   ├── demo /
│   ├── configs
│   │   └── config_template.yaml
│   └── process_lists
│       └── heatmap_demo_dataset.csv
├── presets
│   ├── bwh_biopsy.csv
│   └── ...
├── splits
│   ├── task_1_tumor_vs_normal_75 /
│   └── ...
├── datasets /
├── models /
├── utils /
├── vis_utils /
├── wsi_core /
├── build_preset.py
├── create_heatmaps.py
├── create_patches.py
├── create_patches_fp.py
├── create_splits_seq.py
├── main.py
├── extract_features.py
├── extract_features_fp.py
└── eval.py

※ extract_features.py, create_patches.pyの語尾に"_fp"がついたものがあります。これは、パイプライン処理をしたもので、基本的に"_fp"がついているほうを使用します。

前提条件

CLAMで用いられる画像はバーチャルスライドとして使われるWSI(.svs, .ndpi, .tiff等)を想定しています。
今回は、.ndpiファイルで二値分類を行いました。

手順

完全自動実行
パラメータのカスタマイズ
特徴抽出
データセットの準備
学習
テスト・評価
ヒートマップの作成
感想
完全自動実行

　DATA_DIRECTORY内に入っているすべてのスライドに対して、デフォルトのパラメータで分割し、スライド領域内の組織を含むパッチを抽出し、つなぎ合わせて復元画像を作成します。

python create_patches_fp.py --source DATA_DIRECTORY --save_dir RESULTS_DIRECTORY --patch_size 256 --seg --patch --stitch

実行すると、以下のように出力されます。

source:  入力画像フォルダパス
patch_save_dir:  出力パッチフォルダパス
mask_save_dir:  出力マスクフォルダパス
stitch_save_dir:  出力スティッチフォルダパス
{'seg_params': ... } デフォルト設定

...

progress: 0.98, i/N    <- N枚のうちi番目
processing ファイル名
segmentation took 0.1509876251220703 seconds
patching took -1 seconds
stitching took -1 seconds
average segmentation time in s per slide: 2.1...
average patching time in s per slide: -0.8...
average stiching time in s per slide: -0.8...

復元された画像はRESULTS_DIRECTORYに以下のフォルダ構成で出力されます。masksでは組織を緑線で、空洞を青線で囲っています。stitchesはマスキング画像をもとに組織部分のみを抽出し、ほかの部分を黒くマスクしています。

RESULTS_DIRECTORY/
    ├── masks
    │    ├── slide_1.png
    │    ├── slide_2.png
    │    └── ...
    ├── patches
    │    ├── slide_1.h5
    │    ├── slide_2.h5
    │    └── ...
    ├── stitches
    │    ├── slide_1.png
    │    ├── slide_2.png
    │    └── ...
    └── process_list_autogen.csv

パラメータのカスタマイズ

　csvにデータセットに適用するパラメータ名と値を記載し、コマンド実行時に引数としてパスを指定することで、
実行時に指定したパラメータで処理させることができます。
(スライドに応じて像の境界や穴部分を検知させるために有効です。)

コマンド実行の際に特定の引数をつけ、csvファイルパスを指定することで利用可能ですが、
以下の２つの方法があります。

１．コマンド実行時に--preset引数を付け、csvファイルパスを指定する。

　　➡テンプレートとしてpresets/bwh_biopsy.csvを用いて値を変更します。
　　　各行に、各スライド毎のパラメータが記載してあります。

２．コマンド実行時に--process_list引数を付け、csvファイルパスを指定する。

　　➡初回実行後に出力されるRESULTS_DIRECTORY/process_list_autogen.csvを
　　　ベースに値を変更します。１．と異なる点は主にprocess列があることです。
(process列については後述。)
　　　
どちらのcsvも基本的には同じですが、２．にはスライドを処理するかどうかを指定するprocess列があります。

　データセットに適用したパラメータのテンプレートを定義します。テンプレートはpresets/bwh_biopsy.csvをもとにRESULTS_DIRECTORY/process_list_autogen.csvを変更していきます。（process_list_autogen.csvはパラメータ調整をする前に複製して別名で編集していったほうがいいです。）

python create_patches_fp.py --source DATA_DIRECTORY --save_dir RESULTS_DIRECTORY --patch_size 256 --seg --process_list process_list_edited.csv

今回は、process_list_autogen.csvをprocess_list_edited.csvに変更して、パラメータ調整を行いました。
パラメータの種類には以下のようなものがあります。

パラメータ	デフォルト	説明
seg_level	-1	WSIを分割する際のダウンサンプリングレベル
sthresh	8	セグメンテーション閾値。より高い閾値を用いることで、前景の検出が少なくなり、背景の検出が多くなる。
mthresh	7	平均フィルターサイズ(自然数、奇数)
use_otsu	False	大津の二値化フィルタ
close	4	後の初期閾値化(二値化)の前に行う追加のモルフォロジークロージング処理(自然数 or -1)
a_t	100	組織の領域フィルターの閾値(自然数, 考慮すべき検知される前景の輪郭の最小サイズ) レベル0での512x512の参照パッチサイズを基準とする。
a_h	16	穴用の領域フィルターの閾値 (自然数, 避けるべき前景輪郭の検知される最小の穴/空洞のサイズ)
max_n_holes	10	検知された前景輪郭ごとに考慮すべき穴の最大数(自然数)
vis_level	-1	セグメンテーション結果を視覚化する際のダウンサンプリングレベル -1: 1/64にするダウンサンプリングに最も近いWSIのダウンサンプリングを用いる
line_thickness	250	レベル0でのセグメンテーション結果を描画する線の太さ
use_padding	True	スライドの境界を埋めるかどうか
contour_fn	four_pt	パッチが前景か背景か判断する輪郭チェック機能

モルフォロジークロージング処理：物体の外郭(境界線)を得るための処理です。

process_list_autogen.csvのprocess列について：
　process列のセルを0にすることで、その行のスライドの処理をスキップすることができ、1にすることでその行のスライドは処理されます。
　調整したい画像が書かれた行のprocessセルを1にして、パラメータの変更→create_patches_fp.py実行→mask画像の確認... を繰り返して組織がちゃんと囲われているように調整していきます。

　すべての画像でセグメンテーションの結果に満足できたら、process列をすべて1にして保存します。保存先は、presetsフォルダにしました。

python create_patches_fp.py --source DATA_DIRECTORY --save_dir RESULTS_DIRECTORY --patch_size 256 --seg --process_list CSV_FILE_NAME --patch --stitch

特徴抽出

CUDA_VISIBLE_DEVICES=0,1 python extract_features_fp.py --data_h5_dir DIR_TO_COORDS --data_slide_dir DATA_DIRECTORY --csv_path CSV_FILE_NAME --feat_dir FEATURES_DIRECTORY --batch_size 512 --slide_ext .svs

GPUを使わない場合は、"CUDA_VISIBLE_DEVICES=0,1"を削除してください。
create_patches_fp.py実行後に出力された.h5ファイルがDIR_TO_COORDSフォルダ下に保存されていることを想定しており、実行するとバッチサイズ512で各スライドの各組織パッチから1024次元の特徴を抽出し、以下のフォルダ構造を作成します。

FEATURES_DIRECTORY/
    ├── h5_files
    │      ├── slide_1.h5
    │      ├── slide_2.h5
    │      └── ...
    └── pt_files
          ├── slide_1.pt
          ├── slide_2.pt
          └── ...

各.h5ファイルには、抽出された特徴配列およびパッチの座標が含まれます。
(学習を高速化するために、各スライドのパッチの特徴だけが含まれる.ptファイルも作成されます。)

データセットの準備

学習のための準備をします。

データの整理（データセットが複数ある場合のみ）
データセットオブジェクト作成のためのコードの編集
データセットの分割（train, validation, test用）

データの整理は複数のデータセットがある場合に下のように一つのフォルダにまとめてください。

DATA_ROOT_DIR/
    ├── DATASET_1_DATA_DIR/
    │    ├── h5_files
    │    │      ├── slide_1.h5
    │    │      ├── slide_2.h5
    │    │      └── ...
    │    └── pt_files
    │           ├── slide_1.pt
    │           ├── slide_2.pt
    │           └── ...
    ├── DATASET_2_DATA_DIR/
    │    ├── h5_files
    │    │         ├── slide_a.h5
    │    │         ├── slide_b.h5
    │    │         └── ...
    │    └── pt_files
    │            ├── slide_a.pt
    │            ├── slide_b.pt
    │            └── ...
    └── ...

　各データセットはDATA_ROOT_DIRのサブフォルダに入っており、抽出された.ptファイルがpt_filesに保存されていることを想定しています。
　DATA_ROOT_DIR内に保存されている学習用画像ファイルをまとめたテーブルを作成します。作成するテーブルは、case_id, slide_id, とラベル情報が入った3つの列で作成し、.csvで保存してください。(例：tumor_vs_normal_dummy_clean.csv)
slide_idは.ptファイルと対応するため、画像ファイルから拡張子を削除した文字列にすれば問題ありません。

もし、自分で分類するモデルを作成する場合は、main.pyとeval.pyに以下のようなコードの追加を行います。

(例："TASK NAME"という分類クラスを作る場合)
if args.task == 'TASK NAME':
    args.n_classes=2
    dataset = Generic_MIL_Dataset(csv_path = 'dataset_csv/DATASET.csv',
                            data_dir= os.path.join(args.data_root_dir, 'tumor_vs_normal_feat_resnet'),
                            shuffle = False, 
                            seed = args.seed, 
                            print_info = True,
                            label_dict = {'LABEL_1':0, 'LABEL_2':1},
                            label_col = 'label',
                            ignore=[])

二値分類なら、TASK NAME, DATASET.CSV,label_dictの変更でいけると思います。多クラス分類の場合には、data_dirなどの変更が必要になる？と思います。

その後、TASK NAMEを引数として渡すため、--taskのリストに追加します。

parser.add_argument('--task', type=str, choices=[
    'task_1_tumor_vs_normal',
    'task_2_tumor_subtyping',
    'TASK NAME' <--- 追加
])

データセットの分割します。アルゴリズムのパフォーマンスを評価するため、学習/検証/テスト用に分割し、複数のfoldに分けていきます。分割は以下のコマンドを実行します。

python create_splits_seq.py --task task_1_tumor_vs_normal --seed 1 --label_frac 0.75 --k 10

label_frac 0.75では(全体-検証データ数-テストデータ数)×0.75が学習データ数になります。
このコマンドの場合には、10個のfoldを作成し、学習データに(該当ラベルのデータ全数 - 検証データ数 - テストデータ数)の75％を用いることを意味します。
検証データ数とテストデータ数はデフォルトで該当ラベルのデータ全数×0.1ですが、
それぞれコマンド実行時に--val_frac引数と--test_frac引数を指定することで変更できます。
実行するとsplitsフォルダ下に'TASK NAME'_'label_frac'のフォルダが作成され、3種類k個のsplitデータファイルが保存されます。

splits
   └── TASK NAME_'label_frac'
         ├── splits_0.csv
         ├── splits_0_bool.csv
         ├── splits_0_descriptor.csv
         ├── ...
         └── splits_k_descriptor.csv

学習（二値分類の場合）

以下のコマンドを引く数を変更します。--label_fracの閾値や--exp_codeのフォルダ名の後ろの数字は直前で使用した--label_fracの値と同じにしてください。もしわからなくなったら、CLAM/splits内に保存されているので確認しましょう。
引数が指定できたら、実行します。
（GPUを使用しない場合には、"CUDA_VISIBLE_DEVICES=0"を削除してください。）
出力はresults/--exp_codeにされていくことに注意してください。

CUDA_VISIBLE_DEVICES=0 python main.py --drop_out --early_stopping --lr 2e-4 --k 10 --label_frac 0.75 --exp_code task_1_tumor_vs_normal_CLAM_75 --weighted_sample --bag_loss ce --inst_loss svm --task task_1_tumor_vs_normal --model_type clam_sb --log_data --data_root_dir DATA_ROOT_DIR

テスト・評価

学習モデルのパフォーマンスをテスト・評価するコマンドを実行します。
(引数はこれまでの同様適宜変更して下さい。)

CUDA_VISIBLE_DEVICES=0 python eval.py --drop_out --k 10 --models_exp_code task_1_tumor_vs_normal_CLAM_75_s1 --save_exp_code task_1_tumor_vs_normal_CLAM_75_s1_cv --task task_1_tumor_vs_normal --model_type clam_sb --results_dir results --data_root_dir DATA_ROOT_DIR

結果は、eval_results/'--save_exp_code_'に保存されているので、その下のsummary.csvを確認します。このファイルにはtest_auc, test_accが各foldごとにのっています。もし、評価が悪いようなら、スライドを見直すかパラメータ調整の調整を行ってください。

その結果を確認し、問題なさそうなら、heatmap.pyの作成に移ります。

ヒートマップの作成

　ヒートマップの作成には、設定ファイルとして.yamlファイルを作成する必要があります。（/heatmaps/configs/config_template.yaml参照）
process_listはパラメータのカスタマイズで使ったprocess_list_edited.csvを複製してslide_id, labelの列を残して、heatmaps/process_list下に保存しました。
作成した設定ファイルは/heatmaps/configsに保存してください。

config.yamlが作成できたら、以下のコマンドを実行します。

CUDA_VISIBLE_DEVICES=0,1 python create_heatmaps.py --config config.yaml

途中結果はheatmaps/heatmap_raw_results
最終結果はheatmaps/heatmap_production_results
に保存されます。

感想

私たちのチームでは、windowsとLinuxで試していきましたが、OSによって動作が少し変わるみたいで「"Not found"パスが見つかりません。」が散見されました。もし、実行するさいにはフルパスで書いたほうがいいかもしれないですね。

わからない点・誤認がありましたら、質問してもらえると助かります。
作成者：坂口

参考文献URL

CLAM github

Lu, M.Y., Williamson, D.F.K., Chen, T.Y. et al. Data-efficient and weakly supervised computational pathology on whole-slide images. Nat Biomed Eng 5, 555–570 (2021). https://doi.org/10.1038/s41551-020-00682-w

@article{lu2021data,
  title={Data-efficient and weakly supervised computational pathology on whole-slide images},
  author={Lu, Ming Y and Williamson, Drew FK and Chen, Tiffany Y and Chen, Richard J and Barbieri, Matteo and Mahmood, Faisal},
  journal={Nature Biomedical Engineering},
  volume={5},
  number={6},
  pages={555--570},
  year={2021},
  publisher={Nature Publishing Group}
}