次元削減の手法としてVAE（変分オートエンコーダ）がありますが、マルチオミクス解析で使用できるライブラリを探していたところ、MOVE というライブラリに出会いましたので使い方を記載します。
MOVE
Git hub
公式ドキュメント
説明が少ないので試行錯誤しました。
https://move-dl.readthedocs.io/
MOVEを使ってVAEを行おう
利用データ
- 今回はGDC Data Portalから以下のデータをダウンロードしました
- RNA発現（以下、RNA）： Gene Expression Quantification / STARでフィルター、*_gene_counts.tsvを使用（517症例）
- DNAメチル化（以下、METH）：Methylation Beta Value / Illumina 450Kでフィルター、*level3betas.txtを使用（579症例）
- コピー数変異（以下、CNV）：Gene-level Copy Numberでフィルター、gene_level_copy_number*.tsvを使用（518症例）
前処理
- 以下の前処理を事前に行いました
モダリティ 解析用行列 主な加工 欠損補完 RNA CPM→log₂ 行列 標準化、分散0除外 各 gene の中央値 METH β→M 値 変換 同上 同上 CNV −2~2 正規化 同上 正常コピー数 (=2) - sample_id（TCGAの症例バーコード）は共通化し、444症例を解析で扱います
準備
- MOVEはコマンドでいろいろな処理を行っていきます
- 行いたい処理は
taskで指定します - taskの種類は以下
Available options in 'task':
analyze_latent
analyze_latent_schema
encode_data
identify_associations_bayes
identify_associations_bayes_schema
identify_associations_ks
identify_associations_ks_schema
identify_associations_ttest
identify_associations_ttest_schema
tcga_luad_latent # 自分で作ったtask（後で使用します）
tune_model_reconstruction
tune_model_reconstruction_schema
tune_model_stability
tune_model_stability_schema
- 今回は
encodeとanalyze_latentを使います
task: encode
- VAEのためにデータの前処理を行ってくれます
- 公式のconfig設定例
# DO NOT EDIT DEFAULTS
defaults:
- base_data
# FEEL FREE TO EDIT BELOW
raw_data_path: data/ # where raw data is stored
interim_data_path: interim_data/ # where intermediate files will be stored
results_path: results/ # where result files will be placed
sample_names:
random.small.ids # names/IDs of each sample, must appear in
# the other datasets
categorical_inputs: # a list of categorical datasets
- name: random.small.drugs
continuous_inputs: # a list of continuous datasets
- name: random.small.proteomics
log2: true #apply log2 before scaling
scale: true #scale data (z-score normalize)
- name: random.small.metagenomics
log2: true
scale: true
-
defaultsは絶対変えてはいけませんが、それ以降は自分で変えられます
-
今回自分で設定した
config/data/tcga-luad.yamlは以下
raw_data_path: data/
interim_data_path: interim/
results_path: results/
sample_names: com_sample # data/com_sample.txt を読む
continuous_inputs:
- name: rna_z_comsam_row
log2: false
scale: false
- name: meth_z_comsam_row
log2: false
scale: false
- name: cnv_z_comsam_row
log2: false
scale: false
continuous_names:
- rna_z_comsam_row
- meth_z_comsam_row
- cnv_z_comsam_row
# 記載しないとエラーになるので空リストで指定
categorical_inputs: []
categorical_names: []
-
log2やscaleも揃えてくれるらしいですが、今回は自分で行ってしまったので、
falseとしました -
com_sampleにはtxtファイルでindex情報（細胞）のみを記載（参考下図）
-
コマンドで実行
move-dl data=tcga_luad task=encode_data
[INFO - encode_data]: Beginning task: encode data
[INFO - encode_data]: Reading 'rna_z_comsam_row'
[INFO - encode_data]: Reading 'meth_z_comsam_row'
[INFO - encode_data]: Reading 'cnv_z_comsam_row'
- interimディレクトリにnumpy_array形式およびtxt形式でファイルが追加されました
task: analyze_latent
さて、VAEを行ってみます。
まずは何もハイパラの設定を行わず、defaultで学習を行います。
move-dl data=tcga_luad task=analyze_latent
[INFO - analyze_latent]: Beginning task: analyze latent space
[INFO - analyze_latent]: Generating visualizations
[INFO - analyze_latent]: Projecting into latent space
[INFO - analyze_latent]: Reconstructing
[INFO - analyze_latent]: Computing reconstruction metrics
[INFO - analyze_latent]: Computing feature importance
Killed
メモリが足らず、最後のfeature_importance（特定の特徴量を摂動させたときの潜在空間の変化を可視化。どの特徴が潜在空間に強く影響しているかを評価）が取得できませんでした。
実行すると勝手に以下のディレクトリ、ファイルが追加されます。
ターミナルの出力が味気ない代わりに、logファイルに過去の学習を含めて追記されていきます。
analyze_latent.logを見てみると、今回は潜在空間に1,000次元へと圧縮されたようです。
Model: VAE (135679 ⇄ 1000 ⇄ 150)
- los_curve
task: tcga_luad_latent（ハイパラを設定）
自分でハイパラメータを設定し学習を制御できるので行ってみます。
以下のようにconfig/task/tcga_luad_latentを設定してみます。
defaults:
- analyze_latent # パッケージ内のanalyze_latent を使う
- _self_
batch_size: 64
model:
num_hidden: [2048, 1024]
cuda: True
training_loop:
num_epochs: 200
early_stopping: True
patience: 10
コマンドで実行する際は、taskに先ほど設定したtcga_luad_latentを指定します。
move-dl data=tcga_luad task=tcga_luad_latent
[INFO - analyze_latent]: Beginning task: analyze latent space
[INFO - analyze_latent]: Generating visualizations
[INFO - analyze_latent]: Projecting into latent space
[INFO - analyze_latent]: Reconstructing
[INFO - analyze_latent]: Computing reconstruction metrics
[INFO - analyze_latent]: Computing feature importance
Killed
またkilledとなりました。。
モデル情報がログに追記され、正しく反映されていることを確認できました。
Model: VAE (135679 ⇄ 2048 ⇄ 1024 ⇄ 150)
- los_curve
defalutのハイパラで学習を行ったほうがlossが良いという結果になりました。
念のため、もう少し丁寧にモデルの中身を見てみます。
MOVEにmodelをloadする属性があるか探しましたが、なさそうなのでtorch.loadでモデルの中を見てみます。
import torch, re
ckpt = torch.load("results/latent_space/model.pt", map_location="cpu", weights_only=False)
# state-dict の key を確認
for k, v in ckpt.items():
# encoder の Linear 層だけ抽出
if re.search(r"encoder.*weight", k):
print(f"{k:40s} {tuple(v.shape)}")
結果
python check_model.py
encoderlayers.0.weight (2048, 135679)
encoderlayers.1.weight (1024, 2048)
encodernorms.0.weight (2048,)
encodernorms.1.weight (1024,)
なお、results_pathを変え忘れると、前の学習で作成されたモデルptを読み込もうとしてエラーが出てしまうので、resultディレクトリの名前を毎回変えておく必要があります。設定はtcga_luadのresults_pathで設定すれば簡単です。
補足：ハイパラチューニング
ドキュメントのtuningによると、experimentでconfigを指定すればHydraのmultirunを実施することが可能です。
その際は、設定をoverrideすればできそうです。
defaults:
- override /data: tcga_luad
- override /task: tune_model_reconstruction
taskには以下の2つを指定できます。
-
tune_model_reconstruction：異なるハイパーパラメータの組み合わせで学習したモデルの再構成精度を報告する -
tune_model_stability：異なるハイパーパラメータで学習されたモデルの潜在空間の安定性を報告する
実際には、以下のようにハイパラを複数設定します。
hydra:
mode: MULTIRUN
sweeper:
params:
task.batch_size: 10, 50
task.model.num_hidden: "[500],[1000]"
task.training_loop.num_epochs: 40, 60, 100
クラスタリング
実際に次元削減されたデータを用いて、クラスタリングまで行なってみます。
今回はdefaultのハイパラメータで設定した学習のほうがlossが下がったのでそちらを使用します。
データを読み込みadataに変換
latent = pd.read_csv("results/latent_space/latent_space.tsv", sep="\t", index_col=0)
adata_latent = sc.AnnData(latent, obs=pd.DataFrame(index=latent.index))
adata_latent
近傍グラフ → t-SNE → leidenクラスタリング
- t-SNEにしたのは、MOVEのdefaultが採用していたからです。シングルセル解析ではUMAPがよく使われているように感じます
- defaultでは
num_latent: 150でした
sc.pp.neighbors(adata_latent, n_neighbors=10, use_rep="X") # 150次元を10次元にする
sc.tl.tsne(adata_latent, use_rep="X", random_state=42)
sc.tl.leiden(adata_latent, resolution=1.0, key_added="latent_leiden")
結果の可視化
sc.pl.tsne(
adata_latent,
color="latent_leiden",
title="MOVE Latent Leiden",
)
右はPCAで次元削減を行い同条件でクラスタリングして可視化したものです。
VAEを行ったほうがクラスタがよく分けられているように見えます。
実際に数値でも確認します。
# --- Silhouette Score ---
sil_vae = silhouette_score(adata_latent.obsm["X_tsne"], labels_vae)
sil_pca = silhouette_score(adata_pca.obsm["X_tsne"], labels_pca)
print(f"Silhouette Score (VAE): {sil_vae:.3f}")
print(f"Silhouette Score (PCA): {sil_pca:.3f}")
# --- Calinski‐Harabasz ---
ch_vae = calinski_harabasz_score(adata_latent.obsm["X_tsne"], labels_vae)
ch_pca = calinski_harabasz_score(adata_pca.obsm["X_tsne"], labels_pca)
print(f"Calinski-Harabasz (VAE): {ch_vae:.1f}")
print(f"Calinski-Harabasz (PCA): {ch_pca:.1f}")
Silhouette Scoreは各クラスタの凝集度と分離度のバランスを評価した手法で、Calinski‐Harabaszはクラスタ内分散とクラスタ間分散の比を評価した手法です。どちらも大きいほどより優れていると解釈ができるので、今回はPCAよりVAEで次元削減しクラスタリングを行ったほうが良いという結論になりました。
以上です。読んでいただきありがとうございました。