Deseq2的理論基礎

原文：Moderated estimation of fold change and dispersion for RNA-seq data with Deseq2 by Love, Anders and Huber 2014

這是對Deseq的延申，簡單總結一下這個模型的統計方法。

---

模型
Number of reads in sample $j$ that are assigned to gene $i$記為$K_{ij}$，假設
$$\begin{gathered} K_{ij}～NB({\mu }_{ij},{\alpha }_i),i=1,?,n,j=1,?,m \\ {\mu }_{ij}=s_j{q}_{ij},\log q_{ij}=\sum _{r=1}^p{x}_{jr}{\beta }_{ir} \end{gathered}$$

其中$s,q$的含義與Deseq中$s,q$的含義相同，$[x_{jr}]$為design matrix，$[{\beta }_{ir}]$是系數矩陣，${\alpha }_i$是dispersion parameter，
$$Var[K_{ij}]={\mu }_{ij}+{\alpha }_i{\mu }_{ij}^2$$

${\alpha }_i$越接近0，$K_{ij}$的方差越接近均值，$s_j$作為size factor，用與Deseq中一樣的方法確定$$s_j={\text{median}}_i\frac{k_{ij}}{({\prod }_{v=1}^m{k}_{iv}{)}^{\frac{1}{m}}}$$

Inference on Dispersion
假設dispersion的先驗為$\log {\alpha }_i～N(\log {\alpha }_{tr}({\bar{\mu }}_i),{\sigma }_d^2)$，${\bar{\mu }}_i=\frac{1}{m}{\sum }_j\frac{K_{ij}}{s_j}$，${\alpha }_{tr}(\bar{\mu })=\frac{a_1}{\bar{\mu }}+{\alpha }_0$，dispersion估計分為三步：

估計gene-wise dispersion ${\alpha }_i^{gw}$, 用MLE估計，${\alpha }_i^{gw}=\underset{\alpha }{\mathrm{argmax}}{l}_{CR}(\alpha )$，其中$l_{CR}$代表用了Cox-Reid Adjustment的對數似然，$$\begin{gathered} {\alpha }_i^{gw}=\underset{\alpha }{\mathrm{argmax}}\sum _{j=1}^m\log f_{NB}(K_{ij};{\mu }_{ij},{\alpha }_i)?\underset{\text{cox-Reid?Bias?Adjustment}}{\frac{1}{2}\log \det (X^{T}WX)} \\ W=\text{diag}\left(\frac{1}{\frac{1}{{\mu }_{i1}}+{\alpha }_i},?,\frac{1}{\frac{1}{{\mu }_{im}}+{\alpha }_i}\right) \end{gathered}$$

擬合dispersion trend ${\alpha }_{tr}$: gamma-family GLM of ${\alpha }_i^{gw}$ on ${\bar{\mu }}_i$ to get estimations of $a_1$ and ${\alpha }_0$.

結合似然與trend prior得到dipersion的MAP估計，$$\begin{gathered} {\alpha }_i^{MAP}=\underset{\alpha }{\mathrm{argmax}}[l_{CR}(\alpha )+\underset{\text{log-Normal?prior}}{{\Lambda }_i(\alpha )}] \\ {\Lambda }_i(\alpha )=\frac{?(\log \alpha ?\log {\alpha }_{tr}({\bar{\mu }}_i){)}^2}{2{\sigma }_d^2} \\ {\sigma }_d^2=\max (0.25,s_{lr}^2?{\psi }_1(\frac{m?p}{2})),s_{lr}={\text{mad}}_i(\log {\alpha }_i^{gw}?\log {\alpha }_{tr}({\bar{\mu }}_i)) \end{gathered}$$ 其中${\psi }_1$是trigamma function，mad表示median absolute deviation，$s_{lr}$為standard logrithm residual，如果$\log {\alpha }_i^{gw}>\log {\alpha }_{tr}({\bar{\mu }}_i)+2s_{lr}$，則認為基因$i$是一個dispersion outlier。

Fold change (系數${\beta }_{ir}$代表fold change)
假設系數的先驗為${\beta }_{ir}～N(0,{\sigma }_r^2)$，用empirical method確定$${\sigma }_r=\frac{Q_{|{\beta }_r|}(1?p)}{Q_N(1?p/2)}$$ 原文默認值$p=0.05$，$Q_N(1?p/2)$代表標準正態分布的$1?p/2$上分位點，$Q_{|{\beta }_r|}$代表$\{{\hat{\beta }}_{ir}^{MLE}\}$的$1?p$ empirical quantile，其中${\hat{\beta }}_{ir}^{MLE}$可以由最開始的模型用IRLS得到。系數的MAP為
$${\beta }_i=\underset{\beta }{\mathrm{argmax}}\left[\sum _{j=1}^m\log f_{NB}(K_{ij};{\mu }_{ij},{\alpha }_i)+\Lambda (\beta )\right]$$

其中
$${\mu }_{ij}=s_j{e}^{{\sum }_{r=1}^p{x}_{jr}{\beta }_{ir}},\Lambda (\beta )=?\sum _{r=1}^p\frac{{\beta }_{ir}^2}{2{\sigma }_r^2}$$

使用IRLS計算，迭代方程為
$$\begin{gathered} {\beta }_i\leftarrow (X^{T}WX+\lambda I{)}^{?1}{X}^{T}Wz \\ {\lambda }_r=\frac{1}{{\sigma }_r^2},z_j=\log \frac{{\mu }_{ij}}{s_j}+\frac{K_{ij}?{\mu }_{ij}}{{\mu }_{ij}} \end{gathered}$$

從迭代方程可以看出，與標準的IRLS不同，這里的迭代方程盡管也有WLS的形式，但由于系數有一個正態先驗，所以$(X^{T}WX+\lambda I{)}^{?1}$繼承了ridge regression的特點，因此最后得到的估計量與標準IRLS估計相比會有fractional shrinkage。

總結

以上是生活随笔為你收集整理的Deseq2的理论基础的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。