- 著者: Christoph Hafemeister, Rahul Satija
- Corresponding author: Christoph Hafemeister (New York Genome Center, New York, NY, USA); Rahul Satija (New York Genome Center, New York, NY, USA)
- 雑誌: Genome Biology
- 発行年: 2019
- Epub日: 2019-12-23
- Article種別: Original Article
- PMID: 31870423
背景
single-cell RNA-seq (scRNA-seq; 単一細胞RNAシーケンシング) 技術の進歩により、個々の細胞レベルでの遺伝子発現プロファイリングが可能となったが、細胞ごとのシーケンス深度 (検出される総分子数) の極めて大きな変動が、下流解析における重大な技術的ノイズとして立ちはだかっている。同一の細胞型を対象とした実験であっても、細胞溶解効率、逆転写効率、およびシーケンス時の確率的な分子サンプリングのばらつきに起因し、観測される UMI (Unique Molecular Identifier; 固有分子識別子) カウント数には桁違いの差が生じる。この技術的ノイズと真の生物学的異質性を分離するための正規化手法の開発は、単一細胞解析における最重要課題の一つである。
従来の主流なアプローチは大きく二系統に分類される。第一の系統は、細胞ごとにサイズファクターを計算して一律にスケーリングを行い、log(x+1) 変換を施す手法であり、Satija et al. NatBiotechnol 2015 や Butler et al. NatBiotechnol 2018、Wolf et al. GenomeBiol 2018 などで広く採用されてきた。しかし、これらのサイズファクターに基づく手法は「すべての遺伝子に対して共通のスケーリング定数を適用できる」という仮定を置いており、低発現遺伝子においてこの仮定が破綻することや、擬似カウント (pseudocount) の加算と対数変換というアドホックな操作が下流解析の結果を歪めることが指摘されていた。第二の系統は、負の二項 (NB; negative binomial) 分布や、ゼロ膨張負の二項 (ZINB; zero-inflated negative binomial) 分布などの確率モデルを用いてカウントデータを直接モデリングする手法 (ZINB-WaVE や DCA (deep count autoencoder) など) である。しかし、これらの非制約的な確率モデルは、遺伝子ごとに独立してパラメータを推定しようとするため、データが極めて疎 (sparse) である scRNA-seq においては過適合 (overfitting) を引き起こしやすく、安定したパラメータ推定が困難であるという課題が残されていた。
このように、従来のサイズファクターベースの手法や非制約的確率モデルには、技術的ノイズを生物学的シグナルから完全に分離できないという限界が存在した。実際に、log 変換を施したデータであっても、遺伝子発現量とシーケンス深度の間に系統的な相関が残存することが定量的に示されており、より原理的な正規化・分散安定化フレームワークが決定的に不足していた。この問題は、バルク RNA-seq データの正規化手法である Love et al. GenomeBiol 2014 や Robinson et al. Bioinformatics 2010 を scRNA-seq データに直接適用することを困難にする要因でもあり、技術的交絡を排除しつつ生物学的異質性を正確に保持する新しいアプローチの確立が強く求められていた。特に、低発現遺伝子と高発現遺伝子で異なるスケーリングが必要であるにもかかわらず、一律のサイズファクターを適用する従来法では対応できないという技術的ギャップが未解明のまま放置されており、単一細胞特有の平均と分散の依存関係を適切に制御する数理モデルが不足しているという課題が残されていた。
目的
本研究の目的は、UMI カウントデータに基づく scRNA-seq 解析において、(i) 正規化後の遺伝子発現量が細胞のシーケンス深度と完全に無相関となり、(ii) 正規化された遺伝子の分散が遺伝子発現量やシーケンス深度に依存せず、真の生物学的異質性のみを反映する、という二つの厳密な条件を満たす新しい正規化・分散安定化手法「sctransform」を開発することである。さらに、擬似カウントの加算や log 変換といった従来のヒューリスティックなステップを完全に排除し、10x Chromium、Drop-seq、inDrop、CEL-Seq2 (single-cell RNA-seq protocol 2) などの多様な UMI プラットフォームに汎用的に適用可能なオープンソースの R パッケージとして実装・公開し、Seurat などの標準的ワークフローに統合することを目指す。
結果
log 正規化におけるシーケンス深度依存性の残存: 従来のサイズファクターに基づく log(x+1) 正規化 (log 正規化) を施した後であっても、遺伝子発現量と細胞の総 UMI 数との間に有意な系統的相関が残存することを定量的に明らかにした (Fig 1)。特に高発現遺伝子群においてこの傾向が顕著であり、n=33148 cells を含む PBMC データセットでは、高発現遺伝子群の平均発現量とシーケンス深度の Pearson 相関係数が log 正規化後も 0.3 を超えるケースが観察された (Fig 1d)。さらに、細胞を総 UMI 数に応じて5つのグループに分割して分散の寄与度を解析したところ、log 正規化データでは、総 UMI 数が少ない細胞グループにおいて高発現遺伝子の分散が不当に高く見積もられ、細胞間での分散の均一性が保たれていないことが示された (Fig 1e)。
非制約的 NB モデルの過適合と正則化による克服: 遺伝子ごとに独立して非制約的な NB 回帰モデルを当てはめた場合、特に低発現から中発現の遺伝子において、分散パラメータ θ や傾き β1 の推定値が極めて不安定になり、過適合が生じることをブートストラップ法 (n=13 replicates) およびフィッシャー情報行列を用いた信頼区間解析により実証した (Fig 2)。これに対し、遺伝子平均発現量に対するカーネル平滑化を適用してパラメータを正則化することで、推定値の標準偏差が平均して約 10-fold 減少することを確認した (Fig 2b)。この正則化パラメータを用いることで、低発現遺伝子におけるノイズの過剰フィッティングが効果的に抑制された。
Pearson 残差によるシーケンス深度の完全な排除と分散安定化: sctransform によって算出された Pearson 残差は、遺伝子発現量の全範囲にわたって細胞の総 UMI 数との相関がほぼ 0 にまで低減され、技術的交絡が完全に排除された (Fig 3a, c)。また、細胞のシーケンス深度が異なるグループ間でも、遺伝子の残差分散は極めて一貫していた (Fig 3b)。生物学的異質性を持たないコントロールデータセット (ヒト脳 RNA テクニカルレプリケート) においては、大半の遺伝子の残差分散が理論値通り平均 0、分散 1 の分布に収束した (Fig 4a, b)。一方、実際の PBMC データセットにおいては、PPBP や GNLY などの既知の細胞型特異的マーカー遺伝子のみが、生物学的変動を反映して 1 を大きく超える高い残差分散を示した (Fig 4c)。
下流解析 (次元削減・クラスタリング) の精度向上: 従来の log 正規化では、MALAT1 などの高発現ハウスキーピング遺伝子が技術的要因によって高度可変遺伝子 (HVG) の上位に誤って選出されていたが、sctransform では真の生物学的可変遺伝子のみが正確に選択された (Fig 4c, Fig 5)。PCA および UMAP による次元削減において、log 正規化データでは同一の細胞クラスター内部で細胞がシーケンス深度の勾配に沿って整列する技術的バイアスが生じていたが、sctransform の Pearson 残差を用いた解析ではこのバイアスが著しく消失した (Fig 6a)。これにより、CD4 記憶 T 細胞とナイーブ T 細胞、あるいは CD8 エフェクターサブセットといった、微細な細胞亜集団の分離能が大幅に向上した。
差次的発現 (DE) 解析における偽陽性の抑制と堅牢性: 生物学的に同一である CD14+ 単球 (n=5551 cells) をランダムに2群に分け、一方の群の UMI カウントを 50% にダウンサンプリングした擬似的な技術的バイアス実験を実施した (Fig 6c)。log 正規化データに対して t 検定を適用した場合、不適切な正規化に起因して 2,000 個以上の遺伝子が誤って差次的発現遺伝子 (DEG; differentially expressed gene) として検出された (FDR < 0.01、p<0.001) (Fig 6d)。これに対し、sctransform の Pearson 残差を用いた DE 解析では、偽陽性として検出された遺伝子はわずか 11 個にとどまり、シーケンス深度の差による偽陽性を極めて強力に抑制できることが示された。さらに、CD14+ 単球と CD16+ 単球 (n=1475 cells) の比較において、CD16+ 細胞のシーケンス深度を 20% にダウンサンプリングした場合でも、sctransform を用いた DE 解析の結果はダウンサンプリング前と極めて高い一致率 (Pearson 相関 > 0.9) を維持し、発現変動の検出において log2FC 1.8 以上の主要なマーカー遺伝子が安定して同定された (Fig 6e)。
考察/結論
先行研究との違い: 本研究で提案された sctransform は、従来のサイズファクターに基づく一律なスケーリングや、非制約的な確率モデルによる個別遺伝子フィッティングと異なり、細胞のシーケンス深度を共変量とする正則化負の二項回帰モデルを導入している。これにより、log 変換や擬似カウントの加算といった従来のヒューリスティックな前処理が抱えていた「高発現遺伝子における技術的相関の残存」や「低発現遺伝子における分散の歪み」という限界を根本的に克服した。
新規性: 本研究は、scRNA-seq の UMI カウントデータにおいて、遺伝子ごとの負の二項モデルが過適合を起こす問題を本研究で初めて定量的かつ体系的に示し、類似する発現量を持つ遺伝子間で情報を共有するカーネル回帰によるパラメータ正則化を新規に導入した。このアプローチにより、データ駆動型で極めて安定した分散安定化変換 (VST) を実現し、Pearson 残差を正規化発現量として下流解析に直接利用する道を拓いた。
臨床応用: sctransform は、シーケンス深度やクオリティのばらつきが非常に大きい臨床検体由来の scRNA-seq データ (例: がん微小環境における腫瘍浸潤リンパ球や骨髄系細胞の解析) において、技術的ノイズに惑わされることなく、真の疾患関連細胞状態や治療抵抗性クローンを正確に同定することを可能にする。この堅牢な正規化は、バイオマーカー探索や単一細胞レベルでの薬効評価といった臨床現場におけるトランスレーショナル研究の信頼性を担保する基盤技術となり、臨床応用への展開が期待される。
残された課題: 今後の検討課題として、(i) Pearson 残差は標準化された統計量であるため、絶対的な発現コピー数の比較や直感的な解釈が求められるシナリオでは従来の log 変換値の方が好まれる場合があること、(ii) UMI カウントを前提としたモデルであるため、フルレングススマートシーケンス (Smart-seq2 など) のリードカウントデータには直接適用できないこと、(iii) データセット内の大部分の遺伝子が有意な生物学的変動を示さないという暗黙の仮定を置いているため、極めて均一な細胞集団や逆に極端に不均一な系では正則化曲線が影響を受ける可能性があることが挙げられる。しかし、後続の Seurat v3/v4 への統合や、計算速度を劇的に改善した glmGamPoi バックエンドの導入などにより、本手法は scRNA-seq 解析における事実上のデファクトスタンダードとして確立されている。
方法
sctransform は、各遺伝子について、細胞あたりの総 UMI 数 (シーケンス深度) の対数変換値 log10(depth) を独立変数とする正則化負の二項回帰 (regularized negative binomial regression) モデルを構築する。具体的には、遺伝子ごとに NB 一般化線形モデル (GLM; generalized linear model) を当てはめ、切片 (intercept; β0)、傾き (slope; β1)、および分散パラメータ (dispersion; θ) の3つのパラメータを推定する。
低発現遺伝子や細胞数が限られるデータセットにおいて、非制約的な GLM フィッティングが過適合を起こす問題を解決するため、カーネル回帰 (kernel regression) による正則化を導入する。全遺伝子のパラメータ推定値を遺伝子の平均発現量に対して平滑化 (smoothing) することで、類似する発現量を持つ遺伝子間で情報をプールし、安定したパラメータ曲線 (regularization curve) を得る。この平滑化には、データ駆動型のバンド幅選択法 (bandwidth selection) を用いた正規カーネルを適用する。
最終的な正規化値として、正則化されたモデルから算出される Pearson 残差 (Pearson residuals) を採用する。Pearson 残差は、観測カウントからモデルの期待値を差し引き、モデル由来の標準偏差で除算することで得られ、自動的に分散安定化変換 (VST; variance-stabilizing transformation) として機能する。極端な外れ値の影響を排除するため、残差は細胞総数 N の平方根である √N にクリップされる。
本手法の性能評価およびベンチマークには、33,148個のヒト末梢血単核細胞 (PBMC; 10x Genomics v2) データ、ヒト脳コントロール RNA データセット、Drop-seq による網膜双極細胞、CEL-Seq2、inDrop などの多様なデータセットを使用した。また、技術的検証として、HEK293T (Human Embryonic Kidney 293T) 細胞株から得られたデータセットも活用した。比較対象として、従来の log 正規化 (Seurat の LogNormalize) および scran を用いた。評価指標には、遺伝子発現量と総 UMI 数の相関、分散の安定性、高度可変遺伝子 (HVG; highly variable gene) 選択の正確性、主成分分析 (PCA; principal component analysis) および UMAP (Uniform Manifold Approximation and Projection) による次元削減、Louvain クラスタリングの再現性、および差次的発現 (DE; differential expression) 解析における偽陽性率と検出感度を設定した。
統計的解析およびパラメータ評価においては、独立した遺伝子ごとのフィッティング精度を評価するために Pearson correlation (ピアソン相関係数) や Spearman correlation (スピアマン順位相関係数) を用いた。計算速度を最適化するため、初期 GLM フィッティングステップでは全遺伝子ではなく、遺伝子平均発現量の密度推定に基づいてサンプリングされた 2,000 遺伝子のみを使用し、Poisson 誤差分布を仮定して β 係数を推定した後、最尤法を用いて NB の θ パラメータを推定する高速化アルゴリズムを実装した。