RNA-seq Count的数据分布
参考
- Differential expression analysis
- Gamma-Poisson Distribution
- https://hbctraining.github.io/DGE_workshop/lessons/01_DGE_setup_and_overview.html
- https://hbctraining.github.io/Intro-to-DGE/lessons/01b_DGE_setup_and_overview.html
- 2007, Moderated statistical tests for assessing differences in tag abundance
- 2020, ComBat-seq: batch effect adjustment for RNA-seq count data
- 2014, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
- 2010, Differential expression analysis for sequence count data, DeSeq
- https://arxiv.org/abs/1104.3889
RNA-seq Count分布的重要性
RNA-seq通常用于寻找两组不同处理样本的差异表达基因。其做法是,将NGS产生的reads回帖到参考基因组,落在基因区域的reads count代表基因表达高低的信号,接下来就是比较不同处理下基因表达的信号是否发生变化。
假设我们有两组处理的样本,乳腺癌和正常样本,并且关注一个感兴趣的基因叫做HER2。我们希望通过RNA-seq的分析,能够有信心的说它在乳腺癌中高表达。如何有信心得出这个结论是本文要讨论的内容。
RNA-seq产生基因表达的信号——reads count,本质上是一个随机数,我们不能简单的说HER2基因在乳腺癌中的reads count大,就说HER2基因在乳腺癌高表达。我们需要使用假设检验的想法,将我们期待的假设作为null hypothesis(例如这个例子中null hypothesis是,HER2基因在乳腺癌中高表达),当null hypothesis发生的概率小于一个极小的阈值(通常是0.01或0.05)时,我们就拒绝null hypothesis,并说明我们犯错的概率很小,可以忽略不计。
像t-test这类的工具并不是为count数据设计的,并且RNA-seq中一组的样本数量通常不是很大,我们也不能相信中心极限定量,而使用t-test。因此我们需要为我们的count数据建模,来准确的进行差异基因表达分析。我们建立广义线性模型,其中因变量就是Count,自变量就是处理条件(不同处理条件会影响Count的大小),我们可以通过数学方法解出自变量(处理条件)的系数,你可以想象在这例子中,我们处理条件就是乳腺癌和正常样本,在广义线性模型中我们就1表示乳腺癌组样本,0表示正常组样本,根据之前的假设,我们预期这个自变量的系数将是大于0的正数,因为只有这个系数大于0,乳腺癌组的count数会相对比较大。那么我们就可以建立null hypothesis,不同处理这个自变量是等于零,接下来,就可以通过一些test得出拒绝零假设的p value。也就是说通过test估计广义线性模型的自变量的系数是否为0可以进行差异基因表达分析,因此接下来的就是分析如何估计这个系数。
由于我们Count不是连续的并且其variance不是常数,因此我们不能使用最小平方法(least squares)对系数进行估计。这里我们需要使用极大似然估计(Maximum Likelihood Estimate,MLE)对自变量系数进行估计,其做法是:
- 第一步:写出观察到的随机变量的机率(这里也就是观察到的reads count的机率)
- 第二步: 将机率转化为(自变量系数)的函数
- 第三步: 计算哪些会使得我们观察到的Y发生的机率最大,这些的估计值就是MLE
为了方便理解,下面我使用logistic regression的极大似然估计来解释这一过程:
- 第一步:写出logistic regission的Y的机率
- 假设我们观察到一组Y的值为(1,0,1),其中,我们用符号表示为,那么其机率表示为(表示第i个随机变量的机率),其中
- 第二步: 将转化为的函数
- =>
- 将每一个观察的机率(Pr(Y_i=y_i)=)转化为,每一个的机率是独立的,那么总的机率,也就是likelihood是将所有观察到的随机比变量的几率相乘,如果使用向量表示,那么
- 第三步:使用
Fisher's scoring
、Newton–Raphson
、Gradient descent
等方法找出使最大的likelihood
极大似然(MLE)的方法对自变量系数估计的一个要点是,我们需要写出观察到的随机变量的机率,对于离散型随机变量来说,就是我们需要能写出它的Probability mass function,也就说我们需要知道随机变量的分布。对于RNA-seq的count数据来说,我们需要知道count数据的分布。
RNA-seq Count的分布是什么?
那么现在的问题就是,RNA-seq count数据的是什么样的分布?
首先我们需要清楚,RNA-seq的Count本质是一个随机数,为了能准确估计这个随机数来源于哪一个分布,我们需要有重复。目前来说,RNA-seq的重复一般是三个或三个以上的样本。
假设在三次重复实验中,我们的测序深度是100条reads。在第一个样本中,落在GeneA上的reads数为10条;在第二个样本中,落在GeneA上的reads数为11条;在第三个样本中,落在GeneA上的reads数为9条。也就是说,测序仪连续100测捕获体系中游离的基因片段,在第一个样本中,总过捕获到了10条GeneA的片段;在第二个样本中,总过捕获到了11条GeneA的片段;在第三个样本中,总过捕获到了9条GeneA的片段。回忆一下,二项分布的描述,在n次伯努利试验中,成功的次数X发生的机率是p,那么随机变量X来源于一个B(n,p)的二项分布。在上面的例子中,n就是我们的测序深度,X就是我们观察到的GeneA的reads count,p就是GeneA ,reads的占比。在上面的例子中,我们可以猜测GeneA的占比为,也就是说GeneA来源于一个的二项分布。
为了方便解释,我们使用表示,基因i(i=1,...,m)在样本j(j=1,...,n)count,其中表示gene的数量,n表示一个条件下样本重复的个数,进一步,可以使用公式表示为:
这代表,样本j下的基因1到m,来源于一个Multinomial的分布,
其中表示测序深度或文库大小;表示基因从1到m的比例,那么。
当时,基因从i到m来源于独立的 Poisson Distribution,表示为:
在解释上面公式之前,首先回忆下Poisson Distribution表示为,其中,也就说Poisson Distribution只有一个参数,并且随机变量的均值于方差相等。
上式中,表示基因1在样本j中的count,表示,测序深度乘以基因1的比例,也就是基因1 count的均值,就是Poisson Distribution中的。
在Multinomial Distribution中N是固定的,当为一个基因分配reads后,其他gene的文库就变小了,因此gene之间的count不是独立的。而当时,这一依赖就可以忽略不计。也就是说来源于一个独立的Poisson Distribution。
independent这一性质在写MLE的likelihood时可方便去写联合概率,也就是直接将每一个观察到的随机变量概率相乘即可。
然而当我们观察真实的RNA-seq数据时,我们绘制一个处理条件下,每个基因的均值与方差的散点图,如下图所示,我们会发现在高表达的基因中,方差大于均值。这与Poisson Distribution的性质不符,也就说使用Poisson Distribution不能很好的拟合我们的数据。
这里我们可以考虑hierarchical model,对于同一条件下的每一个Gene i,其均值服从一个Gamma Distribution。用公式表示为;同一处理条件的下的样本j基因i服从一个Poisson Distribution,用公式表示为,整体用公式表示如下:
或者也可以说服从Negative Binomial Distribution,用公式表示为。