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）对自变量系数进行估计，其做法是：

• 第一步：写出观察到的随机变量Y_1,...,Y_n的机率（这里也就是观察到的reads count的机率）
• 第二步: 将机率转化为\beta(自变量系数)的函数
• 第三步: 计算哪些\beta会使得我们观察到的Y发生的机率最大，这些\beta的估计值就是MLE

为了方便理解，下面我使用logistic regression的极大似然估计来解释这一过程：

• 第一步：写出logistic regission的Y的机率
  • 假设我们观察到一组Y的值为(1,0,1)，其中Y_1=1,Y_2=0,Y_3=1，我们用符号表示为Y_i=y_i，那么其机率表示为Pr(Y_i=y_i)=p_i^{y_i}(1-p_i)^{1-y_i}（表示第i个随机变量的机率），其中y_i=0,1
• 第二步: 将p_i转化为\beta的函数
  • ln(\frac{Pr(y_i=1)}{1-Pr(y_i=1)}) = ln(\frac{p_i}{1-p_i})=\beta_0+\beta_1x_1+...+\beta_kx_k => p_i=\frac{exp(\beta_0+\beta_1x_1+...+\beta_kx_k)}{1+exp(\beta_0+\beta_1x_1+...+\beta_kx_k)}
  • 将每一个观察的机率(Pr(Y_i=y_i)=p_i^{y_i}(1-p_i)^{1-y_i})转化为Pr(Y_i=y_i)=\frac{exp(\beta_0+\beta_1x_1+...+\beta_kx_k)}{1+exp(\beta_0+\beta_1x_1+...+\beta_kx_k)}^{y_i}\times(1-\frac{exp(\beta_0+\beta_1x_1+...+\beta_kx_k)}{1+exp(\beta_0+\beta_1x_1+...+\beta_kx_k)})^{1-y_i},每一个Y_i的机率是独立的，那么总的机率，也就是likelihood是将所有观察到的随机比变量Y_i的几率相乘L(\beta)=\prod^{n}_{i=1}Pr(Y_i=y_i)=\prod_{n}^{i=1}p_i^{y_i}(1-p_i)^{1-y_i}=\prod_{n}^{i=1}\frac{exp(\beta_0+\beta_1x_1+...+\beta_kx_k)}{1+exp(\beta_0+\beta_1x_1+...+\beta_kx_k)}^{y_i}\times(1-\frac{exp(\beta_0+\beta_1x_1+...+\beta_kx_k)}{1+exp(\beta_0+\beta_1x_1+...+\beta_kx_k)})^{1-y_i},如果使用向量x_i\beta表示\beta_0+\beta_1x_1+...+\beta_kx_k，那么L(\beta)=\prod_{n}^{i=1}\frac{exp(x_i\beta)}{1+exp(x_i\beta)}^{y_i}\times(1-\frac{exp(x_i\beta)}{1+exp(x_i\beta)})^{1-y_i}
• 第三步：使用Fisher's scoring、Newton–Raphson、Gradient descent等方法找出使\beta最大的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的占比为\frac{10}{100}=\frac{1}{10}，也就是说GeneA来源于一个B(100,\frac{1}{10})的二项分布。

为了方便解释，我们使用X_{ij}表示，基因i(i=1,...,m)在样本j(j=1,...,n)count，其中表示gene的数量，n表示一个条件下样本重复的个数，进一步，可以使用公式表示为:

这代表，样本j下的基因1到m，来源于一个Multinomial的分布，
其中N_{j}= \sum_{m}^{i=1}X_{ij}表示测序深度或文库大小；(p_{1},..,p_{m})表示基因从1到m的比例，那么\sum_{m}^{i=1}p_{i}=1。

当N_{j} \rightarrow \infty, p_{1},...p_{m}\rightarrow 0时，基因从i到m来源于独立的 Poisson Distribution，表示为:

在解释上面公式之前，首先回忆下Poisson Distribution表示为X\sim Poi(\lambda),其中E(X)=Var(X)=\lambda，也就说Poisson Distribution只有一个参数\lambda，并且随机变量的均值于方差相等。

上式中，X_{1j}表示基因1在样本j中的count，N_j \times p_1表示，测序深度乘以基因1的比例，也就是基因1 count的均值，就是Poisson Distribution中的\lambda。

在Multinomial Distribution中N是固定的，当为一个基因分配reads后，其他gene的文库就变小了，因此gene之间的count不是独立的。而当N_{j} \rightarrow \infty, p_{1},...,p_{m}\rightarrow 0时，这一依赖就可以忽略不计。也就是说X_{ij}来源于一个独立的Poisson Distribution。

independent这一性质在写MLE的likelihood时可方便去写联合概率，也就是直接将每一个观察到的随机变量概率相乘即可。

然而当我们观察真实的RNA-seq数据时，我们绘制一个处理条件下，每个基因的均值与方差的散点图，如下图所示，我们会发现在高表达的基因中，方差大于均值。这与Poisson Distribution的性质不符，也就说使用Poisson Distribution不能很好的拟合我们的数据。

这里我们可以考虑hierarchical model，对于同一条件下的每一个Gene i，其均值服从一个Gamma Distribution。用公式表示为\lambda_{ij} \sim Gamma(\alpha_{i}, \beta_{i}) ;同一处理条件的下的样本j基因i服从一个Poisson Distribution，用公式表示为X_{ij}\sim Pio(\lambda_{ij})，整体用公式表示如下:

或者也可以说X_{ij}服从Negative Binomial Distribution，用公式表示为X_{ij} \sim NB(\mu,\phi)。