Deseq的理论基础

原文：Differential Expression Analysis for Sequence Count Data by Anders and Huber 2010

写篇博文简单总结一下这个模型的统计方法。首先，Deseq的目标是给定基因，检验不同组的read counts是否存在显著差别。如果reads是互相独立的，那么read counts服从二项分布，可以由Poisson分布近似（当the probability of read足够小，且样本数足够大时），所以Poisson分布曾被广泛用于count data analysis。但Poisson分布是单参数分布，它的参数除了决定reads的频率外，同样决定了其均值和方差，所以让read counts服从Poisson的假设非常严格，并且均值和方差相等的假设也会造成overdispersion problem，为了避免Poisson分布可能带来的问题，学界采用了负二项分布作为替换。原文也继承了这样的思想，用负二项分布对read counts建模，下面我们开始讨论Deseq的细节。

模型

Number of reads in sample jjj that are assigned to gene iii记为KijK_{ij}Kij，假设
Kij∼NB(μij,σij2)K_{ij} \sim NB(\mu_{ij},\sigma^2_{ij})Kij∼NB(μij,σij2)

其中μij\mu_{ij}μij代表均值，σij2\sigma_{ij}^2σij2代表方差。

通常我们用试验成功的次数rrr与成功概率ppp作为负二项分布的参数，mass function为
f(k)=(k+r−1r−1)pr(1−p)kf(k)=\left( \begin{matrix} k+r-1 \\ r-1 \end{matrix} \right)p^r(1-p)^kf(k)=(k+r−1r−1)pr(1−p)k将参数用均值和方差替换：p=μijσij2,r=μij2σij2−μijp=\frac{\mu_{ij}}{\sigma_{ij}^2},r = \frac{\mu_{ij}^2}{\sigma_{ij}^2-\mu_{ij}}p=σij2μij,r=σij2−μijμij2可以得到read counts的mass function f(kij)=(kij+μij2σij2−μij−1μij2σij2−μij−1)(μijσij2)μij2σij2−μij(1−μijσij2)kijf(k_{ij})=\left( \begin{matrix} k_{ij}+ \frac{\mu_{ij}^2}{\sigma_{ij}^2-\mu_{ij}}-1 \\ \frac{\mu_{ij}^2}{\sigma_{ij}^2-\mu_{ij}}-1 \end{matrix} \right)\left(\frac{\mu_{ij}}{\sigma_{ij}^2} \right)^{\frac{\mu_{ij}^2}{\sigma_{ij}^2-\mu_{ij}}}\left(1-\frac{\mu_{ij}}{\sigma_{ij}^2} \right)^{k_{ij}}f(kij)=⎝⎛kij+σij2−μijμij2−1σij2−μijμij2−1⎠⎞(σij2μij)σij2−μijμij2(1−σij2μij)kij
假设μij=qi,ρ(j)⋅sj\mu_{ij}=q_{i,\rho(j)} \cdot s_jμij=qi,ρ(j)⋅sj，其中sjs_jsj是size factor，它代表样本jjj对基因的覆盖程度，qi,ρ(j)q_{i,\rho(j)}qi,ρ(j)中下标ρ(j)\rho(j)ρ(j)代表sample jjj的experiment condition，所以qi,ρ(j)q_{i,\rho(j)}qi,ρ(j)是condition-dependent value of gene iii，代表expression strength，它与true concentration of fragments of gene iii under condition ρ(j)\rho(j)ρ(j)的期望成正比；
假设σij2=μij⏟shot noise+sj2vi,ρ(j)⏟raw variance\sigma^2_{ij}=\underbrace{\mu_{ij}}_{\text{shot\ noise}}+\underbrace{s_j^2v_{i,\rho(j)}}_{\text{raw\ variance}}σij2=shot noiseμij+raw variancesj2vi,ρ(j) 其中vi,ρ(j)v_{i,\rho(j)}vi,ρ(j)是关于qi,ρ(j)q_{i,\rho(j)}qi,ρ(j)的平滑函数，这个假设的作用是pool the data，因为在估计单个gene的方差时，样本数很可能不足以保证估计量的质量，所以加这个假设使得在做单个gene的方差估计时可以borrow information from entire dataset。

参数估计
考虑data {kij}i=1,j=1i=n,j=m\{k_{ij}\}_{i=1,j=1}^{i=n,j=m}{kij}i=1,j=1i=n,j=m，上述模型中需要估计的参数及其估计量如下

size factor {sj}j=1m\{s_j\}_{j=1}^m{sj}j=1m，估计量为
s^j=medianikij(∏v=1mkiv)1m\hat s_j = \text{median}_i \frac{k_{ij}}{(\prod_{v=1}^m k_{iv})^{\frac{1}{m}}}s^j=mediani(∏v=1mkiv)m1kij这是一个没用统计方法的，人为构造出来的估计，首先要明确sjs_jsj代表样本jjj的数据能覆盖多少基因片段，因此给定基因iii，用kij(∏v=1mkiv)1m\frac{k_{ij}}{(\prod_{v=1}^m k_{iv})^{\frac{1}{m}}}(∏v=1mkiv)m1kij代表样本jjj与所有样本在此基因上的coverage的比值，再取这一系列{kij(∏v=1mkiv)1m}\{\frac{k_{ij}}{(\prod_{v=1}^m k_{iv})^{\frac{1}{m}}}\}{(∏v=1mkiv)m1kij}的中位数作为样本jjj对所有基因片段coverage的平均水平，用median是因为count data的variation可能比较大，为了让估计量更稳健而选择中位数；
expression strength {qi,ρ}i=1n\{q_{i, \rho}\}_{i=1}^n{qi,ρ}i=1n given ρ\rhoρ，估计量为q^i,ρ=1mρ∑j:ρ(j)=ρkijs^j\hat q_{i,\rho}=\frac{1}{m_{\rho}}\sum_{j:\rho(j)=\rho} \frac{k_{ij}}{\hat s_j}q^i,ρ=mρ1j:ρ(j)=ρ∑s^jkij其中mρm_{\rho}mρ代表number of replicates of condition ρ\rhoρ，这个估计量就是对基因iii的所有用size factor做了rescaling的count的平均值
smooth function vi,ρ(qi,ρ(j)):R+→R+v_{i,\rho}(q_{i,\rho(j)}):\mathbb R^+ \to \mathbb R^+vi,ρ(qi,ρ(j)):R+→R+，它的估计量为v^ρ=wρ(q^i,ρ)−zi,ρ\hat v_{\rho}=w_{\rho}(\hat q_{i,\rho})-z_{i,\rho}v^ρ=wρ(q^i,ρ)−zi,ρ 其中wρw_{\rho}wρ是sample variance on the common scale，zi,ρz_{i,\rho}zi,ρ是一个bias correction termwρ(q^i,ρ)=1mρ−1∑j:ρ(j)=ρ(kijs^j−q^i,ρ)2,zi,ρ=q^i,ρmρ∑j:ρ(j)=ρ1s^jw_{\rho}(\hat q_{i,\rho})=\frac{1}{m_{\rho}-1}\sum_{j:\rho(j)=\rho} \left(\frac{k_{ij}}{\hat s_j} -\hat q_{i,\rho}\right)^2,z_{i,\rho}=\frac{\hat q_{i,\rho}}{m_{\rho}}\sum_{j:\rho(j)=\rho} \frac{1}{\hat s_j}wρ(q^i,ρ)=mρ−11j:ρ(j)=ρ∑(s^jkij−q^i,ρ)2,zi,ρ=mρq^i,ρj:ρ(j)=ρ∑s^j1

最后证明一下v^ρ\hat v_{\rho}v^ρ是一个无偏估计：因为s^j\hat s_js^j是构造的，与我们熟悉的参数估计方法无关，所以为了简化问题，我们假设它就是s^j≈sj\hat s_j \approx s_js^j≈sj。因为E[Kij]=μij=qi,ρ(j)⋅sjE[K_{ij}]=\mu_{ij}=q_{i,\rho(j)} \cdot s_jE[Kij]=μij=qi,ρ(j)⋅sj，所以
E[q^i,ρ]=1mρ∑j:ρ(j)=ρEKijs^j=1mρ∑j:ρ(j)=ρqi,ρ(j)⋅sjs^j≈qi,ρE[\hat q_{i,\rho}]=\frac{1}{m_{\rho}}\sum_{j:\rho(j)=\rho} \frac{EK_{ij}}{\hat s_j}=\frac{1}{m_{\rho}}\sum_{j:\rho(j)=\rho} \frac{q_{i,\rho(j)} \cdot s_j}{ \hat s_j} \approx q_{i,\rho}E[q^i,ρ]=mρ1j:ρ(j)=ρ∑s^jEKij=mρ1j:ρ(j)=ρ∑s^jqi,ρ(j)⋅sj≈qi,ρ

因此，当我们认可s^j=sj\hat s_j= s_js^j=sj的假设时，可以将q^i,ρ\hat q_{i,\rho}q^i,ρ视为qi,ρq_{i,\rho}qi,ρ的无偏估计。根据
(mρ−1)wρ(q^i,ρ)=∑j:ρ(j)=ρ(kijs^j−q^i,ρ)2=∑j:ρ(j)=ρ(kijs^j−1mρ∑j:ρ(j)=ρkijs^j)2=(1−1mρ)∑j:ρ(j)=ρkij2s^j2−1mρ∑j:ρ(j)=ρl:ρ(l)=ρj≠lkijkils^js^l\begin{aligned}(m_{\rho}-1)w_{\rho}(\hat q_{i,\rho})& =\sum_{j:\rho(j)=\rho} \left(\frac{k_{ij}}{\hat s_j} -\hat q_{i,\rho}\right)^2 \\ &=\sum_{j:\rho(j)=\rho} \left(\frac{k_{ij}}{\hat s_j} -\frac{1}{m_{\rho}}\sum_{j:\rho(j)=\rho} \frac{k_{ij}}{\hat s_j}\right)^2 \\ & = \left(1-\frac{1}{m_{\rho}} \right)\sum_{j:\rho(j)=\rho} \frac{k_{ij}^2}{\hat s_j^2}-\frac{1}{m_{\rho}}\sum_{\begin{matrix}j:\rho(j)=\rho \\ l:\rho(l)=\rho \\ j \ne l \end{matrix}} \frac{k_{ij}k_{il}}{\hat s_j \hat s_l}\end{aligned}(mρ−1)wρ(q^i,ρ)=j:ρ(j)=ρ∑(s^jkij−q^i,ρ)2=j:ρ(j)=ρ∑⎝⎛s^jkij−mρ1j:ρ(j)=ρ∑s^jkij⎠⎞2=(1−mρ1)j:ρ(j)=ρ∑s^j2kij2−mρ1j:ρ(j)=ρl:ρ(l)=ρj=l∑s^js^lkijkil

因为
EKij2=μij2+σij2=sj2qi,ρ2+sjqi,ρ+sj2vρEKijKil=sjslqi,ρ2EK_{ij}^2=\mu_{ij}^2+\sigma_{ij}^2=s_j^2q_{i,\rho}^2+s_jq_{i,\rho}+s_j^2v_{\rho} \\ EK_{ij}K_{il}=s_js_lq_{_{i,\rho}}^2EKij2=μij2+σij2=sj2qi,ρ2+sjqi,ρ+sj2vρEKijKil=sjslqi,ρ2
计算期望可得，(mρ−1)Ewρ(q^i,ρ)=(1−1mρ)∑j:ρ(j)=ρEKij2s^j2−1mρ∑j:ρ(j)=ρl:ρ(l)=ρj≠lEKijKils^js^l=vρ+qi,ρmρ∑j:ρ(j)=ρ1sj⏟=E[zi,ρ]=E[q^i,ρmρ∑j:ρ(j)=ρ1s^j]\begin{aligned}(m_{\rho}-1)Ew_{\rho}(\hat q_{i,\rho})& =\left(1-\frac{1}{m_{\rho}} \right)\sum_{j:\rho(j)=\rho} \frac{EK_{ij}^2}{\hat s_j^2}-\frac{1}{m_{\rho}}\sum_{\begin{matrix}j:\rho(j)=\rho \\ l:\rho(l)=\rho \\ j \ne l \end{matrix}} \frac{EK_{ij}K_{il}}{\hat s_j \hat s_l} \\ & = v_{\rho}+\underbrace{\frac{q_{i,\rho}}{m_{\rho}}\sum_{j:\rho(j)=\rho}\frac{1}{s_j}}_{=E[z_{i,\rho}]=E[\frac{\hat q_{i,\rho}}{m_{\rho}}\sum_{j:\rho(j)=\rho} \frac{1}{\hat s_j}]}\end{aligned}(mρ−1)Ewρ(q^i,ρ)=(1−mρ1)j:ρ(j)=ρ∑s^j2EKij2−mρ1j:ρ(j)=ρl:ρ(l)=ρj=l∑s^js^lEKijKil=vρ+=E[zi,ρ]=E[mρq^i,ρ∑j:ρ(j)=ρs^j1]mρqi,ρj:ρ(j)=ρ∑sj1

假设检验
考虑ρ∈{A,B}\rho \in \{A,B\}ρ∈{A,B}，replicates分别为mA,mBm_A,m_BmA,mB，我们想检验
H0:qi,A=qi,BHa:qi,A≠qi,BH_0:q_{i,A}=q_{i,B} \\ H_a:q_{i,A} \ne q_{i,B}H0:qi,A=qi,BHa:qi,A=qi,B

为此，我们引入检验统计量
KiA=∑j:ρ(j)=AKij,KiB=∑j:ρ(j)=BKij,KiS=KiA+KiBK_{iA}=\sum_{j:\rho(j)=A}K_{ij},K_{iB}=\sum_{j:\rho(j)=B}K_{ij},K_{iS}=K_{iA}+K_{iB}KiA=j:ρ(j)=A∑Kij,KiB=j:ρ(j)=B∑Kij,KiS=KiA+KiB

给定observation (kiA,kiB)(k_{iA},k_{iB})(kiA,kiB)，用pip_ipi表示这个检验的p值，下面介绍p值的计算方法。

假设H0H_0H0为真，记p(a,b)=P(KiA=a,KiB=b)=P(KiA=a)P(KiB=b)p(a,b)=P(K_{iA}=a,K_{iB}=b)=P(K_{iA}=a)P(K_{iB}=b)p(a,b)=P(KiA=a,KiB=b)=P(KiA=a)P(KiB=b)，KiAK_{iA}KiA与KiBK_{iB}KiB均服从负二项分布，分布的参数估计方法如下：

Pooled mean estimation: q^iO=∑j:ρ(j)∈{A,B}kijs^j\hat q_{iO}=\sum_{j:\rho(j) \in \{A,B\}} \frac{k_{ij}}{\hat s_j}q^iO=j:ρ(j)∈{A,B}∑s^jkij
Estimation of mean: μ^iA=∑j:ρ(j)=As^jq^iO,μ^iB=∑j:ρ(j)=Bs^jq^iO\hat \mu_{iA} = \sum_{j :\rho(j)=A} \hat s_j\hat q_{iO},\hat \mu_{iB} = \sum_{j :\rho(j)=B} \hat s_j\hat q_{iO}μ^iA=j:ρ(j)=A∑s^jq^iO,μ^iB=j:ρ(j)=B∑s^jq^iO
Estimation of variance: σ^iA2=∑j:ρ(j)=As^jq^iO+s^j2vA(q^iO),σ^iB2=∑j:ρ(j)=Bs^jq^iO+s^j2vB(q^iO)\hat \sigma_{iA}^2 = \sum_{j :\rho(j)=A} \hat s_j\hat q_{iO}+\hat s_j^2 v_{A}(\hat q_{iO}) ,\hat \sigma_{iB}^2 = \sum_{j :\rho(j)=B} \hat s_j\hat q_{iO}+\hat s_j^2 v_{B}(\hat q_{iO}) σ^iA2=j:ρ(j)=A∑s^jq^iO+s^j2vA(q^iO),σ^iB2=j:ρ(j)=B∑s^jq^iO+s^j2vB(q^iO)

确定了分布的参数后可以得到p(a,b)p(a,b)p(a,b)。以KiA,KiBK_{iA},K_{iB}KiA,KiB为检验统计量，cut-off point为
{(a,b):a+b=kiS,p(a,b)=p(kiA,kiB)}\{(a,b):a+b=k_{iS},p(a,b)=p(k_{iA},k_{iB})\}{(a,b):a+b=kiS,p(a,b)=p(kiA,kiB)}
pi=∑a+b=kiS,p(a,b)≤p(kiA,kiB)p(a,b)∑a+b=kiSp(a,b)p_i=\frac{\sum_{a+b=k_{iS},p(a,b) \le p(k_{iA},k_{iB})}p(a,b)}{\sum_{a+b=k_{iS}}p(a,b)}pi=∑a+b=kiSp(a,b)∑a+b=kiS,p(a,b)≤p(kiA,kiB)p(a,b)

Deseq的理论基础相关推荐

Deseq2的理论基础
Deseq2的理论基础原文:Moderated estimation of fold change and dispersion for RNA-seq data with Deseq2 by Lo ...
【转载】图像缩放与插值理论基础
图像的缩放图像经过缩放后有可能在原图中招不到对应的像素点,这需要用图像插值来解决. 1.理论基础假设图像的X轴方向缩放比例是Kx,Y轴方向的缩放比是Ky,则缩放后输出图像的点(x' , y')对应 ...
机器学习理论基础到底有多可靠？
机器学习领域近年的发展非常迅速,然而我们对机器学习理论的理解还很有限,有些模型的实验效果甚至超出了我们对基础理论的理解. 目前,领域内越来越多的研究者开始重视和反思这个问题.近日,一位名为 Aidan ...
信息理论基础周炯槃常迥
通信的数学原理作者乡农 https://wenku.baidu.com/view/080af641ddccda38376baf4e.html 信息理论是应用十分广泛的基础学科.它随着通讯技术的 ...
人大魏哲巍：图神经网络的理论基础
[专栏:研究思路]近年来,由于图结构数据的强大表现力,用机器学习方法分析图的研究越来越受到重视.图神经网络是一类基于深度学习的处理图结构数据的方法,在众多领域展现出了卓越的性能,因此已成为一种广泛应用 ...
魏哲巍：图神经网络的理论基础 | 青源 Talk 第 7 期
活动议程日期:11月12日(周五) 时间主题 13:30-13:35 开场简介张峰中国人民大学副教授,青源会会员 13:35-14:20 主题:图神经网络的理论基础魏哲巍中国人民大学教授, ...
汪昭然：构建“元宇宙”和理论基础，让深度强化学习从虚拟走进现实
作者 | 陈彩娴深度强化学习的故事,可以追溯到2015年: 当时,位于英国伦敦的一家小公司 DeepMind 在<Nature>上发表了一篇文章"Human-level con ...
从2019 AI顶会最佳论文，看深度学习的理论基础
2020-01-27 13:15:38 如果能有一种理论告诉我们什么样的模型架构.运算方式能最好地表示某种数据,什么样的损失函数.迭代方式能最高效地学习到某种能力,什么样的设置又使这种能力能处理各种意 ...
科普丨可拓学，诞生于中国的人工智能的理论基础之一
[作者简介] 蔡文,1942 年生,广东澄海人,广东工业大学研究员, 可拓学创立者,国家级有突出贡献的专家杨春燕,1964 年生,山东泰安人,研究员,中国人工智能学会可拓工程专业委员会副主任可 ...

Deseq的理论基础

Deseq的理论基础

Deseq的理论基础相关推荐

最新文章

热门文章