ContinuousVI

A Python library for analyzing single-cell RNA-seq data with continuous covariates using scVI.

ContinuousVI extends the popular single-cell Variational Inference (scVI) framework to incorporate one or more continuous factors (like pseudotime or aging metrics) while correcting for batch effects. It provides straightforward APIs for:

• Multiple model training (with different random seeds/initializations)
• Generating latent embeddings (e.g., UMAP, clustering)
• Regression against continuous covariates (linear, polynomial, spline)
• Sampling from the generative model for gene expression distributions

🧬 Key Features

1. Continuous Covariate Support: Include a single continuous factor (e.g., pseudotime) alongside batch/cell-type labels.
2. Multiple Model Training: Train N scVI models with identical hyperparameters but varying seeds, enabling robust downstream analyses.
3. Dimensionality Reduction & Clustering: Obtain latent embeddings, run UMAP or Leiden clustering, and easily visualize results.
4. Gene Expression Sampling: Sample expression parameters (px) from the learned generative models for posterior predictive analyses.
5. Regression Tools: Regress expression levels against the continuous covariate using OLS, polynomial, or spline models (including advanced multi-sampling approaches for uncertainty estimation).

📕 Installation

ContinuousVI will be published on PyPI. Once available, you can install it via:

pip install continuousvi

Or install directly from source:

git clone https://github.com/<your-org>/continuousvi.git
cd continuousvi
pip install .

Pip location: pip

🚀 Quick Usage Example

import scanpy as sc
from continuousvi import ContinuousVI

# Load AnnData
adata = sc.read_h5ad("my_data.h5ad")

# Initialize
vi_setup = ContinuousVI(
    adata=adata,
    batch_key="batch",
    label_key="cell_type",
    continuous_key="pseudotime"
)

# Train multiple models
trained_vi = vi_setup.train(n_train=5, n_latent=30)

# Calculate embeddings (UMAP, clustering)
trained_vi.calc_embeddings(resolution=0.5, n_neighbors=10, n_pcs=30)

# Perform a simple linear regression against the continuous covariate
df_regression = trained_vi.regression(mode="ols")
print(df_regression.head())

🧠 理論的背景：なぜContinuousVIは共変量を保存できるのか

ContinuousVIは、バッチ効果と生物学的共変量（年齢など）が相関している場合でも、共変量情報を98.8%以上保存しながらバッチ補正を行うことができます。これは以下の3つの理論的機構によって実現されています。

1. Extended Harmonyによる数学的分離

核心原理：直交化による共変量保存

Harmonyの論文で示されているように、バッチ効果と共変量を「数式的に」分離することができます。ContinuousVIのExtended Harmonyは、通常のHarmonyを拡張し、共変量（年齢など）に直交するバッチ成分のみを除去します。

実装の詳細 (continuousVI.py:727-791):

# 1. 基底空間の構築（保存したい共変量を含む）
phi_base = [1, age_normalized]  # 切片 + 年齢

# 2. バッチの直交化（residualization）
# batch_residual = batch - proj_age(batch)
# これにより、年齢と相関するバッチ成分を除去
batch_residual = batch_onehot - phi_base @ (phi_base^T @ phi_base)^{-1} @ phi_base^T @ batch_onehot

# 3. 直交化されたバッチ成分のみをHarmonyで除去
Z_corrected = Z - (batch_residual @ W_batch) * R_k

数学的意味:

• 生のバッチ効果: batch = age_component + technical_component
• Extended Harmony: 技術的成分のみを除去し、年齢成分は保持
• 標準的なscVI: 年齢とバッチが混在した成分を一括除去 → 年齢情報が94.6%喪失

ベンチマーク結果:

• 年齢-バッチ相関 ρ=0.866 の場合
• ContinuousVI: 年齢R²保持率 98.8%
• scVI: 年齢R²保持率 5.4%

重要な注意: なぜR²を主要指標とするのか

scVIは年齢AUROC保持率が98.9%と高いにもかかわらず、R²保持率は5.4%と壊滅的です。この違いは情報の粗視化によるものです：

• AUROC（二値分類）: 「若い群」vs「老い群」の区別だけを測定

  • scVI: 年齢情報の94.6%を失っても、残り5.4%の「おおまかな順序情報」で高いAUROCを達成
  • アナロジー: 10x10ピクセルに圧縮した画像でも「人 vs 犬」の区別は可能

• R²（連続予測）: 22歳 vs 25歳 vs 28歳...の細かい年齢値を測定

  • scVI: 連続的な年齢勾配が失われ、R²が5.4%に低下
  • アナロジー: 10x10ピクセル画像では顔のしわやそばかすは見えない

下流解析（遺伝子発現予測、疑似時間解析、連続的な共変量回帰など）では連続的な情報が必要なため、R²が真の保存率を示す指標となります。

実験データの構造:

Batch 0: age ~ Uniform(20, 35)  ← 若い群
Batch 1: age ~ Uniform(35, 50)  ← 老い群

scVI補正後:
✅ 「Batch 0由来（若）」vs「Batch 1由来（老）」の痕跡は残る → AUROC 98.9%
❌ しかし年齢の連続的な勾配は消失 → R² 5.4%

結論: AUROCが高い ≠ 共変量情報が保存されている

2. 潜在空間とグラフベース補正

核心原理：バッチに依存しない大域的な重要構造の獲得

VAEの潜在空間は、バッチ効果に関連する・しないに関わらず、「全体において重要な空間」を取得します。その後、グラフベースのHarmony補正を適用することで、VAEにバッチ構造を学習させます。

2段階プロセス:

Stage 1: VAEによる大域的構造の学習 (continuousVI.py:442-520)

# Pearson residualsに変換（分散安定化）
X_residuals = PearsonResidualTransform.fit_transform(X_counts)

# VAEで潜在空間へエンコード
Z = VAE.encode(X_residuals)  # Z: (n_cells, n_latent)

# 潜在空間Zは以下を捉える:
# - 細胞タイプ構造（最も重要な生物学的変動）
# - 年齢シグナル（共変量）
# - バッチ効果（技術的変動）
# これらがバッチ依存性とは独立に混在

Stage 2: グラフベースHarmonyによるバッチ構造の学習 (continuousVI.py:545-567)

# k-meansクラスタリングで潜在空間を分割
clusters = KMeans(n_clusters=10).fit_predict(Z)

# ソフトアサインメント（グラフ構造）
R = softmax(-dist(Z, cluster_centers) / sigma)  # (n_clusters, n_cells)

# クラスタごとにバッチ効果を推定して除去
for k in range(n_clusters):
    W_batch = solve_ridge_regression(Z[cluster_k], batch_residual[cluster_k])
    Z_corrected -= (batch_residual @ W_batch) * R_k

なぜこれが機能するのか:

1. 大域的構造の保存: VAEは全データの構造を学習するため、バッチ固有の局所的なアーティファクトに過適合しない
2. 局所的補正: グラフベース補正は潜在空間の局所的な領域（クラスタ）ごとにバッチ効果を推定
3. 直交化: Extended Harmonyにより、年齢と相関しないバッチ成分のみを除去
4. ソフトアサインメント: R_kによる重み付けで、クラスタ境界での不連続性を回避

対比: scVIの問題点

• scVIは p(X | z, batch) を学習し、E[X | z, batch=reference] を返す
• 年齢とバッチが相関している場合、zには g(age + batch) の混合表現が入る
• バッチを参照値に固定すると、年齢-バッチ相関が壊れ、年齢情報が失われる

3. Inverse Transformによるバッチ効果漏出の防止

核心原理：バッチ中立的なライブラリサイズによる逆変換

Pearson residual変換の逆変換（inverse_transform）において、バッチ中立的なライブラリサイズを使用することで、カウント空間への変換時にバッチ効果が再導入されることを防ぎます。

完全なパイプライン (continuousVI.py:442-613):

# [Forward] カウント → Pearson residuals
avg_lib_size = library_sizes.mean()
mu = (library_sizes / avg_lib_size) * gene_means
variance = mu + mu^2 / theta
X_residuals = (X - mu) / sqrt(variance)  # 正規分布に近似

# [VAE + Harmony] 潜在空間で補正
Z = VAE.encode(X_residuals)
Z_corrected = ExtendedHarmony(Z, age, batch)
X_residuals_corrected = VAE.decode(Z_corrected)

# [Inverse] Pearson residuals → カウント (ここが重要!)
library_sizes_neutral = compute_batch_neutral_libsize(X, batch)
mu_neutral = (library_sizes_neutral / avg_lib_size) * gene_means
variance_neutral = mu_neutral + mu_neutral^2 / theta
X_corrected = X_residuals_corrected * sqrt(variance_neutral) + mu_neutral

バッチ中立的ライブラリサイズの計算 (continuousVI.py:793-844):

def compute_batch_neutral_libsize(X, batch_ids):
    """
    バッチ効果の再導入を防ぐためのバッチ中立的ライブラリサイズ
    """
    s = log(library_sizes)  # 対数ライブラリサイズ

    # バッチone-hot行列
    B = one_hot(batch_ids)  # (n_cells, n_batches)

    # バッチ成分を直交射影で除去
    # s* = s - P_B·s, where P_B = B(B^T B)^{-1}B^T
    P_B = B @ inv(B.T @ B) @ B.T
    s_star = s - P_B @ s

    # 絶対スケールを保存（平均を復元）
    s_star = s_star + (s.mean() - s_star.mean())

    library_sizes_neutral = exp(s_star)
    return library_sizes_neutral

なぜこれが重要なのか:

1. ライブラリサイズとバッチの相関除去:

   • バッチ1: 平均ライブラリサイズ 10,000 UMI
   • バッチ2: 平均ライブラリサイズ 20,000 UMI (2倍の技術的差)
   • library_sizes_neutralはこの差を除去し、細胞間の生物学的差のみを保持

2. Inverse transformでの一貫性:

   • 補正されたresiduals X_residuals_correctedには、もはやバッチ効果が含まれない
   • もし元のlibrary_sizesを使うと、カウント空間への変換時にバッチ効果が再導入される
   • library_sizes_neutralを使うことで、バッチ効果の「漏出」を防ぐ

3. 数学的な整合性:

   Forward:  X_batch1 (high lib) → residuals (normalized)
             X_batch2 (low lib)  → residuals (normalized)
   
   VAE + Harmony: residuals → corrected_residuals (batch-free)
   
   Inverse (BAD):  corrected_residuals + lib_batch1 → X_corrected_batch1 (batch effect returns!)
   Inverse (GOOD): corrected_residuals + lib_neutral → X_corrected (batch-free)
   

実験的検証 (benchmark結果):

• バッチ中立的ライブラリサイズを使用: 年齢R²保持率 98.8%
• 元のライブラリサイズを使用（仮説的）: バッチ効果が再導入され、年齢情報が部分的に失われる

まとめ：3つの機構の相乗効果

ContinuousVIの成功は、以下の3つの機構が相乗的に作用することによります:

1. Extended Harmony: 数学的直交化により、年齢と相関しないバッチ成分のみを除去
2. VAE + グラフ補正: 大域的構造を保ちながら、局所的にバッチ効果を学習・除去
3. Batch-neutral inverse transform: カウント空間への変換時に、バッチ効果の再導入を防ぐ

これらにより、年齢とバッチが強く相関している場合（ρ=0.866）でも、年齢情報を98.8%保持しながらバッチ補正を実現できます。

ベンチマーク詳細: /benchmark/SCIENTIFIC_VALIDITY.md を参照

🛠️ Developer Guide

🔧 Environment Setup with uv

If you use uv (a command-line tool for managing Python environments), you can set up a development environment as follows:

# Clone the repository
git clone https://github.com/<your-org>/continuousvi.git
cd continuousvi

# Create and activate a new uv environment (example name: 'contvi-dev')
uv new env contvi-dev
uv activate contvi-dev

# Install an editable version of ContinuousVI along with dev requirements
pip install -e .[dev]

Note: The [dev] extra (or similar) could include testing and linting dependencies if specified in setup.cfg or pyproject.toml.

📁 Project Structure

• ContinuousVI: Sets up your AnnData object and trains multiple scVI models.
• TrainedContinuousVI: Manages trained models, provides methods for embeddings, regression, and sampling.
• Utility Methods: Perform regression (linear, polynomial, spline), advanced regression with multi-sampling, and more.

🪄 Contributing

1. Fork the repository and create your feature branch from main.
2. Make your changes, ensuring that new code is tested and documented.
3. Create a Pull Request, describing your changes and the reason behind them.

📝 License

ContinuousVI is licensed under the MIT License (or the license relevant to your project). Please see the LICENSE file for details.