from pathlib import Path

import scanpy as sc
import anndata as ad
import numpy as np


def load_and_qc(
    h5ad_path: Path,
    min_genes_per_cell: int = 500,
    max_genes_per_cell: int = 8000,
    max_pct_mito: float = 15.0,
    max_pct_ribo: float = 50.0,
) -> ad.AnnData:
    """K562データセットに対する標準QC。"""
    adata = sc.read_h5ad(h5ad_path)
    print(f"読み込み中: 細胞数 {adata.n_obs}, 遺伝子数 {adata.n_vars}")

    # ミトコンドリアおよびリボソーム遺伝子タグ
    adata.var["mt"] = adata.var_names.str.startswith("MT-")
    adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
    sc.pp.calculate_qc_metrics(adata, qc_vars=["mt", "ribo"], inplace=True)

    # フィルタリング
    sc.pp.filter_cells(adata, min_genes=min_genes_per_cell)
    sc.pp.filter_cells(adata, max_genes=max_genes_per_cell)
    adata = adata[adata.obs["pct_counts_mt"] < max_pct_mito, :].copy()
    adata = adata[adata.obs["pct_counts_ribo"] < max_pct_ribo, :].copy()

    # 正規化と対数変換
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)

    # HVG & スケーリング（Geneformerは生のカウントを使用するため、別のレイヤーを保持します）
    if "raw_counts" not in adata.layers:
        adata.layers["raw_counts"] = adata.X.copy()  # Geneformerの入力用

    print(f"QC後: 細胞数 {adata.n_obs}, 遺伝子数 {adata.n_vars}")
    return adata

import torch


class GeneformerEmbedder:
    """事前学習済みのGeneformer 95Mモデルのラッパー。

    公式のGeneformerリポジトリのEmbExtractorとTranscriptomeTokenizerを使用します。
    """

    MODEL_VARIANTS = {
        "gf-6L-30M-i2048": "Geneformer 30M small、コンテキスト2048",
        "gf-12L-95M-i2048": "Geneformer 95M standard、コンテキスト2048",
        "gf-12L-95M-i4096": "Geneformer 95M extended、コンテキスト4096",
    }

    def __init__(
        self,
        device: str = "cuda",
        model_variant: str = "gf-12L-95M-i2048",
    ):
        # Geneformerは独自のトークナイザーとBertForMaskedLM構造を持っています。
        # 実際の実装では、helicalまたは公式リポジトリの例を参照してください。
        from geneformer import TranscriptomeTokenizer, EmbExtractor
        self.tokenizer = TranscriptomeTokenizer(
            custom_attr_name_dict={"cell_type": "cell_type", "target_gene": "target_gene"},
            nproc=4,
        )
        self.extractor = EmbExtractor(
            model_type="Pretrained",
            num_classes=0,
            emb_mode="cell",  # 細胞埋め込み（CLSトークンに類似）
            filter_data=None,
            max_ncells=None,
            emb_layer=-1,
            forward_batch_size=8,
            nproc=4,
        )
        self.device = device
        self.model_variant = model_variant

    def tokenize_adata(self, adata: ad.AnnData, output_dir: Path) -> Path:
        """AnnDataをGeneformerトークンファイル（.dataset）に変換します。

        Ensembl遺伝子IDが必要です（adata.var["ensembl_id"）。
        """
        if "ensembl_id" not in adata.var.columns:
            # mygeneまたはpyensemblを使用して、シンボルをEnsembl IDにマッピングします。
            raise ValueError("adata.var['ensembl_id']が必要です。mygeneまたはpyensemblを使用してマッピングしてください。")

        output_dir.mkdir(parents=True, exist_ok=True)
        # anndataをloom形式に変換し、次にトークン化します（公式のGeneformerワークフロー）
        loom_path = output_dir / "cells.loom"
        adata.write_loom(loom_path)
        self.tokenizer.tokenize_data(
            data_directory=str(output_dir),
            output_directory=str(output_dir),
            output_prefix="tokenized",
            file_format="loom",
        )
        return output_dir / "tokenized.dataset"

    def extract_baseline_embeddings(self, tokenized_path: Path, output_dir: Path) -> np.ndarray:
        """ベースラインの細胞埋め込み（摂動前）を抽出します。"""
        embs = self.extractor.extract_embs(
            model_directory=self.model_variant,
            input_data_file=str(tokenized_path),
            output_directory=str(output_dir),
            output_prefix="baseline_embs",
        )
        return np.asarray(embs)

def in_silico_knockout(
    embedder: GeneformerEmbedder,
    tokenized_path: Path,
    target_gene_ensembl: str,
    output_dir: Path,
    max_ncells: int = 100,
) -> np.ndarray:
    """特定の遺伝子に対してin silicoノックアウトを実行し、摂動された細胞埋め込みを取得します。

    公式のGeneformer InSilicoPerturber API [3]を使用します。
    """
    from geneformer import InSilicoPerturber

    perturber = InSilicoPerturber(
        perturb_type="delete",  # delete: 遺伝子を削除、overexpress: 上に移動
        perturb_rank_shift=None,  # deleteの場合、None
        genes_to_perturb=[target_gene_ensembl],
        combos=0,  # 単一遺伝子
        anchor_gene=None,
        model_type="Pretrained",
        num_classes=0,
        emb_mode="cell",  # 細胞埋め込みの変化を観察します。
        cell_emb_style="mean_pool",
        filter_data=None,
        cell_states_to_model=None,
        max_ncells=max_ncells,  # 細胞のサブセット（計算量の削減のため）
        emb_layer=-1,
        forward_batch_size=8,
        nproc=4,
    )

    perturbed_embs = perturber.perturb(
        model_directory=embedder.model_variant,
        input_data_file=str(tokenized_path),
        output_directory=str(output_dir),
        output_prefix=f"perturbed_{target_gene_ensembl}",
    )
    return np.asarray(perturbed_embs)

def cosine_shift(baseline_embs: np.ndarray, perturbed_embs: np.ndarray) -> np.ndarray:
    """各細胞のコサインシフトを計算します。出力の形状は（N_cells、）です。"""
    baseline_norm = baseline_embs / (np.linalg.norm(baseline_embs, axis=1, keepdims=True) + 1e-8)
    perturbed_norm = perturbed_embs / (np.linalg.norm(perturbed_embs, axis=1, keepdims=True) + 1e-8)
    cos_sim = np.sum(baseline_norm * perturbed_norm, axis=1)
    return 1.0 - cos_sim  # 距離


def euclidean_shift(baseline_embs: np.ndarray, perturbed_embs: np.ndarray) -> np.ndarray:
    """ユークリッド距離（補助的なメトリック）。"""
    return np.linalg.norm(perturbed_embs - baseline_embs, axis=1)


def screening_by_shift(
    embedder: GeneformerEmbedder,
    tokenized_path: Path,
    baseline_embs: np.ndarray,
    target_genes: list[str],  # Ensembl IDのリスト
    output_dir: Path,
    max_ncells_per_gene: int = 100,
) -> dict[str, dict]:
    """複数の遺伝子に対してin silico KOを連続的に実行し、シフト統計を計算します。"""
    results = {}
    for i, gene in enumerate(target_genes):
        if i % 20 == 0:
            print(f"[{i}/{len(target_genes)}] {gene}")
        try:
            perturbed = in_silico_knockout(
                embedder, tokenized_path, gene, output_dir,
                max_ncells=max_ncells_per_gene,
            )
            cos = cosine_shift(baseline_embs[:len(perturbed)], perturbed)
            euc = euclidean_shift(baseline_embs[:len(perturbed)], perturbed)
            results[gene] = {
                "mean_cosine_shift": float(np.mean(cos)),
                "median_cosine_shift": float(np.median(cos)),
                "std_cosine_shift": float(np.std(cos)),
                "mean_euclidean_shift": float(np.mean(euc)),
                "n_cells": int(len(perturbed)),
            }
        except Exception as e:
            results[gene] = {"error": str(e)}

    return dict(sorted(
        results.items(),
        key=lambda x: -(x[1].get("mean_cosine_shift", 0.0) if "error" not in x[1] else 0.0),
    ))

import umap
import matplotlib.pyplot as plt


def visualize_perturbation_trajectory(
    baseline_embs: np.ndarray,
    perturbed_embs: np.ndarray,
    gene_symbol: str,
    output_path: str = "perturbation_umap.png",
    n_arrows: int = 20,
) -> None:
    """UMAPを使用して、ノックアウト後の細胞状態の変化を可視化します。"""
    combined = np.vstack([baseline_embs, perturbed_embs])
    reducer = umap.UMAP(n_components=2, n_neighbors=15, min_dist=0.1, random_state=42)
    embedding_2d = reducer.fit_transform(combined)

    n_cells = len(baseline_embs)
    baseline_2d = embedding_2d[:n_cells]
    perturbed_2d = embedding_2d[n_cells:]

    fig, ax = plt.subplots(figsize=(12, 8))
    ax.scatter(baseline_2d[:, 0], baseline_2d[:, 1],
               color="lightgray", alpha=0.5, s=20, label="ベースライン（WT）")
    ax.scatter(perturbed_2d[:, 0], perturbed_2d[:, 1],
               color="crimson", alpha=0.7, s=30, label=f"ノックアウト({gene_symbol})")

    # サンプル細胞の移動を示す矢印
    n_arrows_actual = min(n_arrows, n_cells)
    idx = np.random.choice(n_cells, n_arrows_actual, replace=False)
    for i in idx:
        ax.annotate("", xy=perturbed_2d[i], xytext=baseline_2d[i],
                    arrowprops=dict(arrowstyle="->", color="black", alpha=0.4, lw=0.8))

    ax.set_title(f"細胞状態の変化：{gene_symbol}のin silicoノックアウト\n(細胞数={n_cells}, 矢印数={n_arrows_actual})")
    ax.set_xlabel("UMAP1")
    ax.set_ylabel("UMAP2")
    ax.legend()
    plt.tight_layout()
    plt.savefig(output_path, dpi=150)
    plt.close()

import pandas as pd
from scipy.stats import spearmanr


def load_replogle_ground_truth(
    replogle_h5ad_path: Path,
    baseline_ctrl_label: str = "non-targeting",
) -> pd.DataFrame:
    """Replogle 2022 Perturb-seqからの遺伝子ごとのKO効果を読み込みます。

    GEO GSE168191からK562必須ライブラリをダウンロードし、処理します。
    各ターゲット遺伝子の細胞状態が、ノントargetingコントロールとどれだけ異なるかを定量化します。
    """
    adata = sc.read_h5ad(replogle_h5ad_path)
    # 'target_gene'列が必要です
    if "target_gene" not in adata.obs.columns:
        raise ValueError("adata.obs['target_gene']が必要です")

    # コントロール細胞の埋め込み（例：PCA）
    if "X_pca" not in adata.obsm:
        sc.pp.pca(adata, n_comps=50)

    ctrl_mask = adata.obs["target_gene"] == baseline_ctrl_label
    ctrl_centroid = adata.obsm["X_pca"][ctrl_mask].mean(axis=0)

    # 各ターゲット遺伝子のセントロイドとコントロールのセントロイドの間の距離
    effects = []
    for target in adata.obs["target_gene"].unique():
        if target == baseline_ctrl_label:
            continue
        target_mask = adata.obs["target_gene"] == target
        target_centroid = adata.obsm["X_pca"][target_mask].mean(axis=0)
        # ユークリッド距離とコサイン距離
        eucl_dist = float(np.linalg.norm(target_centroid - ctrl_centroid))
        cos_sim = float(np.dot(target_centroid, ctrl_centroid) /
                        (np.linalg.norm(target_centroid) * np.linalg.norm(ctrl_centroid) + 1e-8))
        effects.append({
            "target_gene": target,
            "measured_shift": eucl_dist,
            "measured_cosine": 1.0 - cos_sim,
            "n_cells": int(target_mask.sum()),
        })
    return pd.DataFrame(effects)


def compare_with_wetlab(
    in_silico_results: dict[str, dict],
    wet_lab_df: pd.DataFrame,
    top_k: int = 20,
    gene_symbol_mapping: dict[str, str] | None = None,  # Ensemblからシンボルのマッピング
) -> dict:
    """in silico予測と実実験の測定値を比較します。"""
    is_records = []
    for ensembl, stats in in_silico_results.items():
        if "error" in stats:
            continue
        symbol = gene_symbol_mapping.get(ensembl, ensembl) if gene_symbol_mapping else ensembl
        is_records.append({
            "target_gene": symbol,
            "silico_shift": stats["mean_cosine_shift"],
        })
    is_df = pd.DataFrame(is_records)
    merged = is_df.merge(wet_lab_df, on="target_gene", how="inner")
    print(f"共通の遺伝子: {len(merged)}")

    if len(merged) < 3:
        return {"error": "共通の遺伝子が少なすぎます"}

    # スピアマン相関（シフトの大きさ）
    rho, pval = spearmanr(merged["silico_shift"], merged["measured_cosine"])

    # 上位Kにおける再現率
    is_topk = set(merged.nlargest(top_k, "silico_shift")["target_gene"])
    wl_topk = set(merged.nlargest(top_k, "measured_cosine")["target_gene"])
    recall_at_k = len(is_topk & wl_topk) / min(top_k, len(merged))

    # 上位Kにおける適合率
    precision_at_k = len(is_topk & wl_topk) / len(is_topk)

    return {
        "n_common_genes": len(merged),
        "spearman_rho": float(rho),
        "spearman_p": float(pval),
        f"recall_at_top{top_k}": float(recall_at_k),
        f"precision_at_top{top_k}": float(precision_at_k),
        "is_topk_genes": list(is_topk),
        "wl_topk_genes": list(wl_topk),
        "overlap_genes": list(is_topk & wl_topk),
    }

K562_ESSENTIAL_GENES_EXAMPLE = [
    # 例。実際には、DepMap、MAGeCK、または必須遺伝子データベースから選択します。
    "ENSG00000141510",  # TP53
    "ENSG00000012048",  # BRCA1
    "ENSG00000139618",  # BRCA2
    "ENSG00000186092",  # MYC（例）
    # ... 200遺伝子
]


def full_perturbation_pipeline(
    k562_h5ad_path: Path,
    essential_gene_ensembls: list[str],
    replogle_ref_path: Path | None,
    output_dir: Path,
    device: str = "cuda",
) -> dict:
    """K562必須遺伝子のin silicoスクリーニングを実行し、Perturb-seqで検証します。"""
    output_dir.mkdir(parents=True, exist_ok=True)

    print("[1/6] K562データのQC")
    adata = load_and_qc(k562_h5ad_path)

    print("[2/6] Geneformerをロード")
    embedder = GeneformerEmbedder(device=device)

    print("[3/6] トークン化とベースライン埋め込み")
    tokenized = embedder.tokenize_adata(adata, output_dir / "tokens")
    baseline = embedder.extract_baseline_embeddings(tokenized, output_dir / "baseline")
    print(f"      ベースラインの形状={baseline.shape}")

    print(f"[4/6] {len(essential_gene_ensembls)}個の遺伝子のin silico KO")
    silico_results = screening_by_shift(
        embedder, tokenized, baseline, essential_gene_ensembls,
        output_dir / "perturbations", max_ncells_per_gene=50,
    )
    import json
    with open(output_dir / "silico_shifts.json", "w") as f:
        json.dump(silico_results, f, indent=2, ensure_ascii=False)

    print("[5/6] 上位遺伝子のUMAP軌跡を可視化")
    top_genes = [g for g in list(silico_results.keys())[:5] if "error" not in silico_results[g]]
    for gene in top_genes:
        perturbed_dir = output_dir / "perturbations" / f"perturbed_{gene}"
        # ... 埋め込み結果をembedderまたはキャッシュから再読み込みします。
        pass  # 実際には、計算量を削減するために、キャッシュを使用して再計算を避けるようにします。

    print("[6/6] Perturb-seqの測定値と比較")
    validation = {}
    if replogle_ref_path:
        wet_lab_df = load_replogle_ground_truth(replogle_ref_path)
        validation = compare_with_wetlab(silico_results, wet_lab_df, top_k=20)
        with open(output_dir / "validation.json", "w") as f:
            json.dump(validation, f, indent=2, ensure_ascii=False)
        print(f"      スピアマンρ={validation.get('spearman_rho', 0):.3f}、再現率@20={validation.get('recall_at_top20', 0):.2f}")

    return {"silico_ranking": silico_results, "validation": validation}

## 性能、コスト、および既知の失敗事例

### 性能の指標（公開ベンチマークを使用）

| モデル | ベンチマーク（Perturb-seq リコール@20） | スピアマンの ρ | ソース |
|------|:----------------------------:|:----------:|------|
| ランダムベースライン | K562 エッセンシャル | 0.10 | ベースライン |
| GENIE3（相関ベース） | K562 エッセンシャル | 0.25 | 従来の方法 |
| scGPT in silico KO | K562 エッセンシャル | 0.45〜0.55 | Cui et al., Nat Methods 2024 [2] |
| Geneformer in silico KO | Replogle 2022 サブセット | 0.50〜0.65 | Theodoris et al., Nature 2023 [3] |
| scFoundation | マルチティッシュ | 0.55〜0.70 | Hao et al., Nat Methods 2024 [4] |
| UCE | クロススペシャス | 0.50〜0.60 | Rosen et al. 2024 [5] |
| アンサンブル（scGPT + Geneformer） | K562 | ～0.70（推定） | コミュニティベンチマーク |

### 学習者が再現するための推定コスト

- APIコスト：0（ローカルGPUを使用する場合）。
- 十分な性能を持つGPUがない学習者は、クラウドのオンデマンド時間単位の価格を参照する必要があります（Geneformer 95Mは24GBのVRAMで実行できます）。
- 200個の遺伝子をスクリーニングするには、約30〜60分かかります。
- データダウンロード：Geneformerのサイズは4GB、Replogle Perturb-seqサブセットは5〜10GBです。

### 5つの既知の失敗事例（コミュニティおよび論文から収集）

1. **in silico KOとウェットラボの結果の乖離（特定の細胞タイプ）**
   症状：in silico KO予測は、事前学習データに存在しない細胞タイプに対して意味のない値を生成します（例：特定の脳ニューロンのサブタイプ）。
   原因：基盤モデルの事前学習データのカバー範囲の限界。
   対策：（a）ターゲットの細胞タイプを含むデータでモデルをファインチューニングする、（b）複数の基盤モデルをアンサンブルする（Geneformer + scGPT）、（c）信頼度に基づいて予測をフィルタリングする（埋め込みシフトの大きさが非常に小さい場合）、（d）K562は事前学習データに含まれているため、比較的安定しています。
   ソース：Kedzierska et al. "Assessing the limits of zero-shot foundation models in single-cell biology." bioRxiv 2023 [6]。

2. **Gene Ensembl IDマッピングエラー**
   症状：遺伝子シンボル（TP53）をEnsembl ID（ENSG00000141510）にマッピングできない→in silico KOが誤った遺伝子に適用される。
   原因：遺伝子シンボルにはエイリアスがあり、Geneformerは特定のEnsemblバージョン（例：GRCh38.p13）でトレーニングされています。
   対策：（a）`pyensembl`または`mygene.info`を使用して、正確にマッピングし、バージョンを指定する、（b）Ensemblバージョンを固定する、（c）マッピングに失敗した遺伝子をログに記録してスキップする、（d）更新されたGeneformer再トレーニングバージョンを確認する。
   ソース：Geneformer GitHub Issues [7]。

3. **バッチ効果がin silicoシフトに影響を与える**
   症状：異なるバッチの細胞をまとめて分析すると、バッチ効果がシフトとして解釈され、誤った遺伝子効果の解釈につながる。
   原因：事前QCと正規化によって、バッチ効果が完全に除去されない。
   対策：（a）バッチ効果の補正（Harmony、scVI、Scanorama）、（b）in silico KOを単一のバッチ内でのみ実行する、（c）コントロール遺伝子（中性コントロール）のシフトをベースラインノイズとして使用して、信号対雑音比を計算する、（d）シフト値の再現性を複数のバッチで検証する。
   ソース：Peidli S et al. "scPerturb: harmonized single-cell perturbation data." Nat Methods 2024 [8]。

4. **ランクドロッププロトコルと実際のKOのギャップ**
   症状：ランクドロップシミュレーションは、実際のgRNA-Cas9 KOの結果と異なる。
   原因：ランクドロップは単に発現ランキングをシフトさせるだけであり、実際のKOはタンパク質、複合体、下流のカスケード、および時間的ダイナミクスを反映する。
   対策：（a）シフトが大きい遺伝子のみを強いシグナルとして解釈する、（b）Geneformerの複数の摂動モード（削除、過剰発現、ノックダウン）を試す、（c）結果を仮説生成のみに使用し、ウェットラボでの検証を必要とする、（d）Perturb-seqベンチマークとの整合性を定量的に評価する。
   ソース：Theodoris et al. Nature 2023 discussion [3]; Kedzierska 2023 [6]。

5. **大規模なデータセットのダウンロードと保存にかかる負担**
   症状：Replogle 2022 Perturb-seqデータセット全体は数十GB、K562エッセンシャルサブセットも5〜10GBである。
   原因：シングルセルデータの生データは大きく、スパース行列であっても大きい。
   対策：（a）サブセットを使用する（エッセンシャル遺伝子200サブセットのみ）、（b）GEOから必要なサンプルのみをダウンロードする、（c）圧縮して保存する（h5ad + zstd）、（d）学術的なクラウドリソースを利用する（例：CZI Chan Zuckerberg Initiative BioHub）。
   ソース：GEO GSE168191 Replogle 2022 [1]。

## 拡張アイデア

- **過剰発現シミュレーション：** ランクを上げる代わりにランクを下げることで、過剰発現をシミュレートする（例：転写因子、腫瘍抑制遺伝子）。
- **多遺伝子組み合わせ：** 2つの遺伝子を同時にノックアウトする（合成致死性候補をスクリーニングする）。
- **薬物応答予測：** 薬物標的KOの細胞シフトを実際の薬物治療データと比較する→薬物再利用（セクション07および08へのリンク）。
- **クロススペシャス転送：** マウスデータで学習したパターンをヒトデータに適用する（UCEを使用）。
- **ライブscRNA-seq統合：** リアルタイムの実験データをストリームする→in silico予測をリアルタイムデータと比較する。
- **HIFネットワークとシグナル伝達経路の詳細な予測：** 特定のシグナル伝達経路のKO組み合わせ。

## 次のセクション

- セクション04 `crispr-guide-scoring`：in silicoでスクリーニングされた遺伝子のgRNAを設計する（K562ウェットラボへのエントリー）。
- セクション07 `drug-target-gnn`：遺伝子KO予測を利用して、薬物標的を優先する。
- セクション13 `protein-design-multimodal`：KO遺伝子を置き換える人工アクチベータータンパク質を設計する。
- セクション14 `bio-mcp-agent`：in silico摂動をMCPツールとして公開する→「K562で大きなシフトを示す遺伝子を抽出する」を自律的に実行する。

## 参考文献

1.  Replogle JM, Saunders RA, Pogson AN, et al. "Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq." Cell 2022. `https://www.cell.com/cell/fulltext/S0092-8674(22)00597-9` · GEO GSE168191
2.  Cui H, Wang C, Maan H, et al. "scGPT: toward building a foundation model for single-cell multi-omics using generative AI." Nature Methods 2024. `https://www.nature.com/articles/s41592-024-02201-0`
3.  Theodoris CV, Xiao L, Chopra A, et al. "Transfer learning enables predictions in network biology (Geneformer)." Nature 2023. `https://www.nature.com/articles/s41586-023-06139-9`
4.  Hao M, Gong J, Zeng X, et al. "Large-scale foundation model on single-cell transcriptomics (scFoundation)." Nature Methods 2024. `https://www.nature.com/articles/s41592-024-02305-7`
5.  Rosen Y, Roohani Y, Agarwal A, et al. "Universal Cell Embeddings: A Foundation Model for Cell Biology (UCE)." bioRxiv 2024. `https://www.biorxiv.org/content/10.1101/2023.11.28.568918v2`
6.  Kedzierska KZ, Crawford L, Amini AP, Lu AX. "Assessing the limits of zero-shot foundation models in single-cell biology." bioRxiv 2023. `https://www.biorxiv.org/content/10.1101/2023.10.16.561085`
7.  Geneformer GitHub Issues: `https://huggingface.co/ctheodoris/Geneformer/discussions`
8.  Peidli S, Green TD, Shen C, et al. "scPerturb: harmonized single-cell perturbation data." Nature Methods 2024. `https://www.nature.com/articles/s41592-023-02144-y`
9.  Geneformer HuggingFace: `https://huggingface.co/ctheodoris/Geneformer`
10. scGPT GitHub: `https://github.com/bowang-lab/scGPT`
11. CELLxGENE (CZI): `https://cellxgene.cziscience.com/`
12. scanpy: `https://scanpy.readthedocs.io/`
13. anndata: `https://anndata.readthedocs.io/`
14. Human Cell Atlas: `https://www.humancellatlas.org/`
15. UMAP-learn: `https://umap-learn.readthedocs.io/`
16. helical SDK (Geneformer/scGPT/UCE 統合): `https://github.com/helicalAI/helical`
17. DepMap (essential gene DB): `https://depmap.org/`
18. MAGeCK (CRISPR screen 分析): `https://sourceforge.net/p/mageck/wiki/Home/`
19. Harmony batch correction: `https://github.com/immunogenomics/harmony`
20. scVI (batch correction + generative model): `https://scvi-tools.org/`