単細胞RNA-seqデータ解析における主要課題の克服

近年、Geneformerはトランスフォーマー構造を応用したゲノムスケールの遺伝子発現モデルとして注目を集めています。しかし、生データから生物学的知見を得るまでのプロセスには依然として多くの技術的障壁があります。以下のような課題に直面している研究者は少なくありません：

• loomやH5AD形式の生データをモデル入力に変換する手間
• 分類タスクにおけるハイパーパラメータ調整の試行錯誤
• 特定遺伝子の機能喪失効果を定量的に評価する手段の欠如
• 高次元埋め込み表現の解釈困難性
• 複数フェノタイプを同時に予測する統合モデルの構築難易度

本稿では、これらの課題を解決するGeneformerエコシステム内の5つのコアコンポーネントについて、実践的なコード例とともに紹介します。

1. TranscriptomeTokenizer：自動化されたトークン化エンジン

単細胞発現マトリクスは通常3万以上の遺伝子を含みますが、Geneformerは最大4096遺伝子までしか処理できません。この次元削減と標準化を一括で行うのがTranscriptomeTokenizerです。

from geneformer import TranscriptomeTokenizer

tokenizer = TranscriptomeTokenizer(
    custom_attr_name_dict={"cell_type": "cell_type", "tissue": "organ"},
    nproc=16,
    model_version="V2"
)

tokenizer.tokenize_data(
    data_directory="input/scrna",
    output_directory="output/tokens",
    output_prefix="heart_cells",
    file_format="h5ad"
)

内部では以下の処理が行われます：

1. ID正規化：Ensembl IDへの統一変換（ensembl_mapping_dict_gc104M.pkl使用）
2. 順位変換：各遺伝子の発現量を中央値ベースの相対順位に置き換え
3. シーケンス整形：上位4096遺伝子を選択し、必要に応じてパディングまたは切り捨て

並列処理により、百万セル規模のデータセットも数時間以内に処理可能です。

2. Classifier：標準化された分類パイプライン

心筋症など疾患状態の予測モデル構築を簡素化するフレームワークです。事前学習済みモデルの一部をフリーズすることで、小規模データでも過学習を抑えて高精度な分類が可能になります。

from geneformer import Classifier

clf = Classifier(
    classifier="cell",
    cell_state_dict={
        "state_key": "disease",
        "states": ["healthy", "dcm", "hcm"]
    },
    filter_data={"cell_type": ["Cardiomyocyte"]},
    training_args={
        "num_train_epochs": 5,
        "learning_rate": 5e-4,
        "per_device_train_batch_size": 16,
        "weight_decay": 0.2
    },
    freeze_layers=2,
    model_version="V2",
    nproc=16
)

clf.prepare_data(
    input_data_file="output/tokens/heart_cells.arrow",
    output_directory="output/classifier",
    output_prefix="cardio_disease"
)

results = clf.validate(
    model_directory="models/Geneformer-V2-104M",
    prepared_input_data_file="output/classifier/cardio_disease_labeled_train.dataset",
    id_class_dict_file="output/classifier/cardio_disease_id_class_dict.pkl",
    output_directory="output/classifier",
    output_prefix="eval"
)

評価結果はROC曲線や混同行列として可視化され、AUCスコアはV2モデルで0.95以上に達することがあります。

3. InSilicoPerturber：仮想的遺伝子ノックアウトシミュレーション

特定遺伝子の発現を計算機上で抑制または増強し、それが細胞状態に与える影響を定量化します。これにより、候補治療標的のスクリーニングが可能になります。

from geneformer import InSilicoPerturber, EmbExtractor

# 基準となる埋め込みの取得
emb_extractor = EmbExtractor(
    model_type="CellClassifier",
    num_classes=3,
    filter_data={"cell_type": ["Cardiomyocyte"]},
    max_ncells=500,
    model_version="V2"
)

baseline_embs = emb_extractor.get_state_embs(
    cell_states_to_model={
        "state_key": "disease",
        "start_state": "dcm",
        "goal_state": "healthy"
    },
    model_directory="output/classifier/fine_tuned",
    input_data_file="output/tokens/heart_cells.dataset",
    output_directory="output/embs",
    output_prefix="dcm_baseline"
)

# ノックアウトシミュレーション
perturber = InSilicoPerturber(
    perturb_type="delete",
    genes_to_perturb="top_100_variable",
    model_type="CellClassifier",
    num_classes=3,
    cell_states_to_model={
        "state_key": "disease",
        "start_state": "dcm",
        "goal_state": "healthy"
    },
    state_embs_dict=baseline_embs,
    max_ncells=1000,
    model_version="V2"
)

perturber.perturb_data(
    model_directory="output/classifier/fine_tuned",
    input_data_file="output/tokens/heart_cells.dataset",
    output_directory="output/perturb",
    output_prefix="knockout_sim"
)

出力される各遺伝子の「遷移促進スコア」は、DCM→正常状態への移行をどれだけ助けたかを示します。これをボルケーノプロットで可視化することで、候補遺伝子を絞り込めます。

4. EmbExtractor：高次元特徴の抽出と可視化

Transformerの中間層から得られる埋め込みベクトルは、細胞の状態を包括的に表現しています。これらの特徴量を次元削減して解釈可能な空間に投影できます。

from geneformer import EmbExtractor
import umap
import matplotlib.pyplot as plt

extractor = EmbExtractor(
    model_type="Pretrained",
    emb_mode="cls",
    emb_layer=-1,
    max_ncells=2000,
    model_version="V2"
)

embeddings = extractor.extract_embs(
    model_directory="models/Geneformer-V2-104M",
    input_data_file="output/tokens/heart_cells.dataset",
    output_directory="output/embeddings",
    output_prefix="cls_layer_last"
)

# UMAPによる可視化
reducer = umap.UMAP(n_neighbors=15, min_dist=0.1)
reduced = reducer.fit_transform(embeddings)

plt.figure(figsize=(10, 8))
scatter = plt.scatter(reduced[:, 0], reduced[:, 1], c=labels, cmap='Spectral', s=6)
plt.colorbar(scatter)
plt.title("UMAP of Cardiomyocyte States")
plt.savefig("umap_cardio.png", dpi=300, bbox_inches='tight')

この手法により、細胞集団の連続性やサブクラスター構造を視覚的に確認できます。

5. MTLClassifier：多タスク学習による知識共有

単一のモデルで複数の生物学的属性（例：細胞種、疾患ステータス、年齢群）を同時に予測することで、汎化性能とデータ効率を向上させます。

from geneformer import MTLClassifier

mtl_clf = MTLClassifier(
    task_columns=["cell_subtype", "disease_status", "age_group"],
    pretrained_path="models/Geneformer-V2-104M",
    model_save_path="output/mtl_model",
    batch_size=16,
    epochs=10,
    use_task_weights=True,
    nproc=16
)

mtl_clf.train()
performance = mtl_clf.evaluate()
mtl_clf.save_results("output/mtl_results")

タスクごとの損失重みを動的に調整することで、バランスの取れた学習が実現されます。特に稀なクラスの予測精度が向上する傾向が観察されています。

ツール間の連携と最適なワークフロー

これら5つのモジュールは独立して使用できるだけでなく、以下のように連鎖的に活用できます：

1. 前処理：TranscriptomeTokenizerでトークン化
2. 分類：Classifierで疾患状態を予測
3. 特徴抽出：EmbExtractorで潜在表現を取得
4. 仮想実験：InSilicoPerturberで候補遺伝子を同定
5. 統合予測：MTLClassifierで複数フェノタイプを関連付ける

このようなパイプラインにより、生データから創薬ターゲット候補の提案までを一貫して実施できます。

推奨設定とトラブルシューティング

大規模データセットでの安定した実行のためのガイドライン：

• メモリ不足：バッチサイズを8以下に下げ、gradient_accumulation_stepsを使用
• 学習不安定：学習率ウォームアップ（10%エポック）と高い重み減衰（0.2～0.3）を適用
• 過学習：早期停止（early stopping）、ドロップアウト率の増加、データオーグメンテーション

また、GPUリソースに応じた推定所要時間は次の通りです：

今後の展開として、空間的転写産物データとの統合や、注意力機構に基づく解釈可能性の強化が期待されています。