single cell RNA-seq のdropout - 驚異のアニヲタ社会復帰の予備

MikuHatsune2017-06-14

読んだ。
MAGIC: A diffusion-based imputation method reveals gene-gene interactions in single-cell RNA-sequencing data
コードはPython で書かれている。

computational flowcytometry のDana Peer ラボ。RNA-seq のデータ行列が取るであろう高次元空間の形からデータ点(細胞)を復元しようという雰囲気。
single cell RNA-seq ではmRNA の読みの精度の悪さのせいで、全体の80-90% くらいは0 という結果が返ってくるdropout という現象に悩まされる。有効なデータが10% 程度しかなく、他が0 のときに元のデータがどのようなものだったかを推定するのが必要な作業だが、MAGIC(Markov Affinity-based Graph Imputation of Cells) という、マルコフ連鎖と近傍のデータ点からもともと読まれるはずだったmRNA のデータと細胞を再構成してデータを復元する方法があるらしい。
　

図はNat Methods. 2014 Jul; 11(7): 740–742.のFig.1a
　
n 遺伝子 m 細胞の行列 $D$ がある。このなかのほとんどはdropout のせいで0 であり、データを復元した $D_{imputed}$ がほしい。細胞で眺めれば遺伝子発現のベクトルになるので、このベクトルで表現型が分類でき、ベクトルが近いものはその近傍に細胞が固まるはずである(diffusion)。グラフ理論的にある表現型(ベクトルパターン)をもつ細胞の近傍k 個を考える。
$D$ からは簡単に距離行列 $D_{dist}$ が計算できる(適当にユークリッド)。そして、 $D_{dist}$ を適切なadaptive ガウスカーネル関数で変換するのだが、
$A=e^{-(\frac{D_{dist}}{\sigma})^2}$
というように変換してaffinity matrix を作る。このとき、 $\sigma$ はkernel width と呼ばれ、このアルゴリズムでは非常に重要である。 $\sigma$ 以下の $D_{dist}$ は速やかに0 affinity になるため、相当近い細胞同士でないとneighborhood が構成できない。
$A$ は対称行列でないので $M=A+A^T$ として強制的に対称行列化して、マルコフ過程を満たすために行和を1 になるようにする。つまり確率推移行列になるようにする。
ここで、 $t$ というべき乗パラメータを使って、 $D_{imputed}=M^tD$ とすると、目的だった $D_{imputed}$ が得られる。ここで $D_{imputed}$ は、dropout して80-90% が0 だったRNA-seq データが、もともと得られるはずだったRNA-seq リードカウントデータになっている、ということらしい。
( $M^t$ は行列積ではなく、ただ単に要素のべき乗である)
　
scImpute: Accurate And Robust Imputation For Single Cell RNA-Seq Data
MAGIC とは異なり、発現の多寡とdropout するかしないかを確率分布からモデル化して、それを説明するパラメータを推定する。
データの取得のされ方について以下のようなモデルを考える。
$f_{X_{i}}(x)=\lambda_i \textrm{Gamma}(x;a_i,b_i)+(1-\lambda_i)\textrm{Normal}(x;\mu_i,\sigma_i)$
$\lambda_i$ ：dropout rate
$a_i,b_i$ ：ガンマ分布のパラメータ。ガンマ分布はdropout をモデル化しており、遺伝子発現としては高くてばらつきの少ないものだが、これが0になるのはdropout によるもの、という感じ。
$\mu_i,\sigma_i$ ：正規分布のパラメータ。正規分布は通常の遺伝子発現をモデル化しており、低めから中程度で、ばらつきの大きな遺伝子発現をしており、これが0 になるということはまあ普通に考えて遺伝子発現がないから0 である、という感じ。
これをEM アルゴリズムでゴリゴリ解く。ベイズでやってもいいと思う。
実際に $i$ 遺伝子が $j$ 細胞でdropout する確率は
$d_{ij}=\frac{\lambda_i \textrm{Gamma}(x;a_i,b_i)}{\lambda_i \textrm{Gamma}(x;a_i,b_i)+(1-\lambda_i)\textrm{Normal}(x;\mu_i,\sigma_i)}$
となる(各パラメータはEM アルゴリズムで推定されたあとのハット付きである)
　
dropout 確率と、細胞遺伝子発現の隠れモデルが推定されたので、実際にdropout しているデータをimputation する。ある適当なdropout 確率閾値 $t$ により、細胞 $j$ についてimputation されないといけない遺伝子発現セット $A_j=\{i:d_{ij}\geq t\}$ と、dropout 確率閾値を下回っているため確実にデータとして信頼できる遺伝子発現セット $B_j=\{i:d_{ij}< t\}$ に分けて、 $B_j$ を用いてimputation する。これは $L_1$ 正則化で行い、
$\hat{\beta}^{(j)}=\underset{\beta^{(j)}}{\textrm{argmin}}\frac{1}{|B_j|}||w_{B_j}^T(X_{B_j,j}-X_{B_j,-j}\beta^{(j)})||_2^2+\theta||\beta^{(j)}||_1$
となる発現量を復元する。
dropout するかしないかというのは、例えば行動定量学で鳥の生態を観察するような、鳥は生きているが観察されないのか、鳥は死んでいるので観察されないのか、というのをモデル化するzero inflated model 的な感じがした。
　
Nat Methods. 2014 Jul; 11(7): 740–742.
もともとはこれがいっぱい引用されている。
ベイズ的にやるのだが、モデル式としては、dropout はポアソン分布、amplification は負の二項分布からサンプリングされるとしている。負の二項分布を使うと裾が伸びて分散が大きい値を表現できるから、これが高発現というかamplification を意図しているのだと思う。
あるsubpopulation で平均して $x$ の発現量がある遺伝子の事後確率を
$p_{S}=E[\Pi_{c\in B}p(x|r_c,\Omega_c)]$
とする。ただし、 $B$ は $S$ のブートストラップサンプル、 $r$ はRPM (RNA-seq で得られるリードカウントを適当に補正したようなもので、なんらかの定量値と思えばいい)、 $\Omega$ はエラー項であり、 $c$ の細胞について、 $p(x|r_c,\Omega_c)$ は事後確率であり、dropout を観測する確率 $p_d$ を用いて
$p(x|r_c,\Omega_c)=\p_d(x)p_{\textrm{Poisson}}(x)+(1-p_d(x))p_{\textrm{NB}}(x|r_c)$
と書く。
というようなことをやっているので、結局はdropout で0 になる、その0 になるのが通常発現パターンか、高発現パターンなのかを別のそれっぽい分布から出てきたようにモデル化するっぽい。