Seurat、楽しんでますか？最近では rPCA が出たり、SCTransform が改良されたり、version5 ではオブジェクト構造自体がガラッと変わったり...などなど、進化が目まぐるしいですね。

今回は、そもそも SCTransform って何してるの？という疑問が生じたので、SCTransform 関数のデバッグ結果から、得られた理論を辿っていこうと思います。

使用パッケージ
Seurat: 5.3.0, sctransform: 0.4.2, glmGamPoi: 1.16.0

使用データ (Zheng Y et al., Cell Biosci. 12(1): 168 (2022))
GSM5293869: 32738 genes × 8447 cells

2025/11/20: より細かく、R で可能な限界まで debug し、追記しました。
追記箇所一覧
① カーネル密度推定の使われ方
② 負の二項回帰に使用される遺伝子の条件
③ その他、分かりにくい式を修正

前提として、負の二項モデルは 0 が豊富に含まれ、かつ平均よりも分散が大きいデータにおいて適しているとされています。scRNA-seqデータはこのモデルが適合します。しかし、以前よりパラメータの不安定性(過学習)が指摘されており、修正が必要です。そこで考案された方法が SCTransform です。

SCTでは主に４つの段階を踏みます。

① 各セルのシーケンス深度と遺伝子 i の発現カウントの間の関係を、負の二項回帰によりモデル化する。
② 各遺伝子の幾何平均発現と①で得られた各パラメータの関係を、それぞれカーネル回帰することで、パラメータの正則化(安定化)を行う。
③ 正則化パラメータを用いてピアソン残差を計算し、その分散から変動性の高い遺伝子を特定する。
④ 修正後のカウント値及びそのlog1p値を計算する。
⑤ 回帰パラメータを用いて線形回帰を行う。(オプション)

① 各セルのシーケンス深度と遺伝子 i の発現カウントの間の関係を、負の二項回帰によりモデル化します。モデルは遺伝子ごとに作成します。横軸には各細胞ごとのシーケンス深度(mj), 縦軸には各遺伝子ごとにカウント値(yij)を取り、回帰モデルを作成します。

ここで使用されるリンク関数、線形予測子、対数尤度関数は以下の通りです。

リンク関数と線形予測子
$$
log_{}(E(y_{ij})) = β_{i} + log_{}(m_j)
$$
対数尤度関数(負の二項分布の式をガンマ関数に変形)
$$
P(α_i, y_i) = \sum_{j=1}^{n} log \left \lbrace \frac{Γ(θ_i + y_{ij})}{Γ(θ_i)Γ(y_{ij} + 1)}
\left(\frac{θ_i}{θ_i + α_{ij}}\right)^{θ_i}
\left(\frac{α_{ij}}{θ_i + α_{ij}}\right)^{y_{ij}}
\right \rbrace
$$
$$
※　log_{}(α_{ij}) = {β_{i} + log_{}(m_j)}
$$
　　　　　　　　　　　　𝒚𝒊𝒋 : 細胞jにおける遺伝子iのカウント値
　　　　　　𝒎𝒋 : 細胞jにおけるカウント値の合計(シーケンス深度)
　　　　　　　　　　　𝜷𝒊, 𝜽𝒊 : 未知パラメータ(求めるパラメータ)
　　　　　　　　　　　　　　　　　　　　　　　Γ( ) : ガンマ関数

実際の実装は以下のようになっています。

①-1
モデル化に使用する細胞 (5000 が上限) を固定シード値でランダムに選びます。また、デフォルトパラメータであれば、5 細胞以上で検出がある遺伝子のみを以降で使用します。

①-2
モデル構築に使用する 2000 遺伝子を選びます。ただし、以下のような条件を設け、偏りが無いように選びます。

まず、遺伝子ごとに分散と平均を計算し、分散 > 平均であるような遺伝子を過分散であるとします。次に、過分散遺伝子における以下の値を計算し、カーネル密度推定によって分布密度を推定します。
$$
log_{10} \left (e^{\frac{1}{N_j}\sum_{j=1}^{N_j}log(y_{ij} + 1)}-1 \right)
$$
　　　　　　　　　　　　　　　　　　　　　　　Nj: 細胞数(今回は 5000)

log 内部は遺伝子ごとの幾何平均発現値で、補正として -1 を加えることで、実数全体の値を取ることが出来ます。

最後に、推定密度の逆数を各遺伝子の抽出確率とし、密度的に偏りが無いように 2000 遺伝子を取り出します。以降、5000 細胞 × 2000 遺伝子がモデル構築に用いられます。

①-3
分散パラメータ, および回帰係数 β の初期値を求めます。初期値は以下のように計算されます。
$$
dispersion_{i, init} = {\frac{σ^2_i + \frac{μ_i}{\frac{1}{N_j}\sum_{i=1}^{N_j}m_{j}}}{μ_i^2}}
$$
　　　　　　　　　　　　　　　　　　　　　μi: 遺伝子ごとの算術平均値
　　　　　　　　　　　　　　　　　　　　　　　σi^2: 遺伝子ごとの分散
$$
β_{i, init} = log \left (\frac{1}{N_j}{\sum_{j=1}^{N_j}{\frac{y_{ij}}{m_j}}} \right )
$$
dispersion の式は、以下に示す負の二項分布の性質より導くことができますが、実際の初期値には、μ に対する補正としてシーケンス深度の平均値が入っています。また、値は 0 以上として定義されるようです。
$$
σ^2 = μ+\frac{μ^2}{θ} ⇔ \frac{1}{θ} = \frac{σ^2-μ}{μ^2}
$$
また、β に関しては log 正規化した発現値を使用しています。これも先ほど同様、シーケンス深度の影響を緩和する役割があると思われます。

①-4
β の推定を行います。推定には "_glmGamPoi_fitBeta_one_group" 関数が使用されており、R では debug 出来ませんでしたが、勾配降下のような方法で、繰り返し推定を行うようです。一般に、解を導くことが困難な場合には、このような方法が取られるとのことです。

①-5
推定した β を用い、dispersion の推定を行います。推定には _glmGamPoi_estimate_overdispersions_fast 関数が使用されており、先ほど同様 debug 出来ませんでしたが、繰り返し推定を行うという方法は変わらないようです。また、先ほど β を先に推定したのは、μ に β が関係しているためだと思われます。

①-6
dispersion の収縮を行います。低発現遺伝子の dispersion が非常に大きくなってしまうことがあるため、緩和措置として行うようです。

①-7
最後に、β の再推定を行います。なお、θ は dispersion の逆数として、最終的に格納されます。

②各遺伝子の幾何平均発現と①で得られた各パラメータの関係を、それぞれカーネル回帰することで、パラメータの正則化(安定化)を行います。

回帰における分散パラメータとしては、以下が設定されています。また、幾何平均発現の計算方法も再掲します。

幾何平均発現(対数値)
$$
m(y_i) = log_{10}(e^{\frac{1}{N_j}\sum_{j=1}^{N_j}log(y_{ij} + 1)}-1)
$$

分散パラメータ
$$
v(y_i) = log_{10} \left \lbrace 1+\frac{e^{\frac{1}{N_j}\sum_{j=1}^{N_j}log(y_{ij} + 1)}-1}{θ_i} \right \rbrace
$$
負の二項分布の性質より、分散パラメータは以下のように変形できます。
$$
log_{10} \left(1+\frac{μ}{θ}\right)
= log_{10}\left(\frac{μ+\frac{μ^2}{θ}}{μ}\right)
= log_{10} \left(\frac{σ^2}{μ}\right)
= log_{10} \left(CV^2 \times μ\right)
$$
(この部分はあまりよく理解できていないのですが、平均当たりの分散とすることで、分散を安定化しているのでしょうか？もしくは、ポアソン分布からどの程度過分散かをみているのでしょうか？)

最後に、下記のパラメータに対し m(yi) を個別に回帰することで、各パラメータをそれぞれ正則化します。

$$
β_{i}, v(y_i)
$$
実際の実装は以下のようになっています。

②-1
正則化の前に、ポアソン分布に従うような遺伝子を取り除きます。条件は以下の通りです。

・σ^2 < μ (過分散でない)
・算術平均が非常に小さい(< 0.001)
・以下の式で計算される θ が、モデルで推定された θ よりも比率として非常に小さい(< 0.001)　※ この時、θ は実際には非常に大きく、μ^2/θ ≒ 0 であると考え、ポアソン分布であるとみなします。
$$
θ = \frac{μ^2}{σ^2-μ}
$$

また、推定された β が、m(yi) に対して外れ値である場合も、取り除かれます。

②-2
保持された過分散遺伝子について、βi, v(yi) に m(yi) をカーネル回帰します。最終的に、①-1 で得られた全ての遺伝子について各パラメータを得ることが出来ますが、ポアソン分布に従う遺伝子の dispersion は 0 とされます。

③ 正則化パラメータを用いてピアソン残差を計算し、その分散から変動性の高い遺伝子を特定します。平均、分散は上記同様、負の二項分布に基づきます。
$$
Z_{ij} = \frac{y_{ij}-E(y_{ij})}{\sqrt{E(y_{ij}) + \frac{E(y_{ij})^2}{θ_i}}}
$$

また、外れ値のクリッピングも行います。N は細胞数です。
$$
Z_{ij}^{*} =
\begin{cases}
\sqrt{N} & \text{if } Z_{ij} > \sqrt{N} \
,-\sqrt{N} & \text{if } Z_{ij} < -\sqrt{N} \
, Z_{ij} & \text{otherwise}
\end{cases}
$$
発現がシーケンス深度のみによる場合、ピアソン残差は正規性を持つと思われます。しかし、今回は生物学的多様性があるため、そのようにはならないはずです。よって、この Zij の分散が大きいものから順に、変動性が高いとみなしていると考えられます。

④ 各細胞におけるカウント値の log 中央値、正則化したパラメータ、およびピアソン残差を用いて、修正後のカウント値及びそのlog1p値を計算します。ピアソン残差はクリップ前のものです。
$$
log(E(y_i')) = β_{i} + log(median(m_1,m_2, ..., m_j))
$$
$$
y_{ij, correct} = E(y_i') \times Z_{ij} \times \sqrt{E(y_{ij}) + \frac{E(y_{ij})^2}{θ_i}}
$$
この correct の計算は見たことが無いのですが、変形すると以下のようになります。直観的には、遺伝子ごとの発現中央値に、その細胞固有の生物学的変動をかけたものという感じでしょうか？
$$
y_{ij, correct} = E(y_i') \times (y_{ij}-E(y_{ij}))
$$

なお、このままだと負の値や少数を持ちますが、負の値は 0 に、少数は丸められます。

⑤ 回帰パラメータを用いて線形回帰を行います。(オプション)

オプションとして、ミトコンドリア由来遺伝子の割合や、細胞周期差の影響を取り除くことが出来ます。遺伝子ごとに、ピアソン残差 Zij に対して、回帰したい変数を説明変数とした線形回帰を行います。なお、ピアソン残差は √(N/30) を堺にクリッピングされた後、使用されます。最終的な scale.data には、変動性の高い 3000 遺伝子のみが含まれます。

以上となります。いかがでしたでしょうか。負の二項分布や回帰分析の知識があれば理解できると思ったのですが、correct カウントの計算や分散パラメータの計算理論が曖昧なままです。今後も継続した理解に努めていこうと思います。次回は PCA, KNN あたりでも見てみようかな？

参考文献(ゼロから理解するにはこれくらいが必要かも...?)
リンク切れしてたらすみません。
Hafemeister, et al. Genome Biol 20, 296 (2019)
David W. Scott. Make Your Publications Visible. Papers, No. 2004,16 (2004)
Sheather, et al. Journal of the Royal Statistical Society. Series B(Methodological).Volume 53, issue 3, 683-690(1991)
https://www.rdocumentation.org/packages/Seurat/versions/4.0.1/topics/SCTransform
https://jp.mathworks.com/help/stats/residuals.html
https://atmarkit.itmedia.co.jp/ait/articles/2105/24/news019.html
https://masamunetogetoge.com/norm-dist
https://ai-trend.jp/basic-study/normal-distribution/nature/#i
https://ocw.u-tokyo.ac.jp/lecture_files/engin_05/2/notes/ja/ishikawa2.pdf
https://masamunetogetoge.com/fisher
https://masamunetogetoge.com/kernel-linearspeace
https://www.math.kyoto-u.ac.jp/~ichiro/lectures/11st.pdf
https://ja.wikipedia.org/wiki/負の二項分布
 https://qiita.com/wsuzume/items/09a59036c8944fd563ff
https://hoxo-m.hatenablog.com/entry/20151012/p1
https://ja.wikipedia.org/wiki/ガンマ関数
 https://www.hulinks.co.jp/software/math_develop/mj/mj_1308
https://qiita.com/tsnry7913/items/e56e08aa23d141175a04
https://qiita.com/tjmnmn/items/3aed6fb85f75446f74ca
https://www.agr.nagoya-u.ac.jp/~seitai/document/R2009/t090614glmIntro.pdf
http://bin.t.u-tokyo.ac.jp/summercamp2017/file/2-7.pdf
https://support.minitab.com/ja-jp/minitab/18/help-and-how-to/modeling-statistics/regression/supporting-topics/logistic-regression/link-function/
https://www.ism.ac.jp/editsec/toukei/pdf/61-2-271.pdf
https://www.bananarian.net/entry/rnbinomMLE
https://cpp-learning.com/glm/#i-8
https://stats.biopapyrus.jp/bayesian-statistics/estimation/negative-binomial-regression.html
https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/ksmooth
https://runebook.dev/ja/docs/r/library/stats/html/density
http://ibis.t.u-tokyo.ac.jp/suzuki/lecture/2015/dataanalysis/L9.pdf
https://www.iwanttobeacat.com/entry/2019/02/24/011111#fn-ce7d4500
https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/bandwidth
https://ja.wikipedia.org/wiki/カーネル密度推定
 https://data-science.gr.jp/implementation/ida_r_kde.html
https://github.com/satijalab/seurat/blob/13b615c27eeeac85e5c928aa752197ac224339b9/R/preprocessing.R
https://github.com/satijalab/seurat/issues/2414
https://github.com/satijalab/seurat/blob/4e868fcde49dc0a3df47f94f5fb54a421bfdf7bc/vignettes/sctransform_vignette.Rmd
https://github.com/ChristophH/sctransform/blob/master/R/vst.R
https://www.rdocumentation.org/packages/sctransform/versions/0.3.2/topics/row_gmean
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/61.html
http://joemphilips.com/post/r_model_matrix/
http://r-statistics.co/Linear-Regression.html
https://rdrr.io/cran/Seurat/src/R/preprocessing.R#sym-RegressOutMatrix
https://rdrr.io/cran/Seurat/src/R/preprocessing.R#sym-RegressOutMatrix
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/70.html
https://decafish.blog.ss-blog.jp/2013-12-11-1
https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/qr
https://bellcurve.jp/statistics/course/9704.html
https://github.com/satijalab/seurat/issues/1679
https://www.rdocumentation.org/packages/ggplot2/versions/1.0.0/topics/cut_number
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/13.html
https://qiita.com/kazutan/items/5406f0b40f8921f5a75e
https://www.rdocumentation.org/packages/Seurat/versions/4.0.4/topics/CellCycleScoring
https://github.com/satijalab/seurat/issues/728
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/37.html