Deseq2的理论基础

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}\sim NB({\mu }_{ij},{\alpha }_i),i=1,\cdots ,n,j=1,\cdots ,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\sim 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},\cdots ,\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}\sim 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的理论基础的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。