前言
对跨物种的RNA-seq进行标准化和差异分析已知是一个问题,而目前对此类问题的相关研究还比较少,有用RPKM进行各物种之间标准化的,也有基于count文件利用DESeq2的标准化方法对各物种进行标准化的,而今天介绍的方法来自于文章:《A statistical normalization method and differential expression analysis for RNA-seq data between different species》
跨物种基因
首先,对于同一物种相同基因的比较,由于其基因长度和功能都一样,因此可以直接比较;而对于跨物种的RNA-seq比较来说,一般选取直系同源的基因来比较
跨物种RNA-seq基本模型
首先,作者定义基本模型如下:
其中:
- E(Xgkt) 代表物种 t 中文库 k 基因 g 基因观测到的count值的期望(均值;上式右边可以看作为对 μgk 求均值的过程,因此式子左边用 E(Xgkt) 表示,另外一层意思参照下面的泊松分布模型);
- Xgkt 代表物种 t 中文库 k 基因 g 观测到的count
- μgkt 代表物种 t 中文库 k 基因 g 的真实表达水平;
- Lgkt 代表物种 t 中文库 k 基因 g 的基因长度;
- St 代表
Nt 代表物种 t 中,文库 k 的所有基因count数总和;
上面的模型建立了物种 t 中文库 k 基因 g 的真实表达水平与观测值之间的关系,有助于下一步的标准化及差异分析
跨物种RNA-seq标准化及差异分析
首先对于两个物种的直系同源基因的比较,我们有如下假设:
那么H0对应该基因没有差异表达,H1对应该基因发生了差异表达;之前我们说 Xgkt 代表物种 t 中文库 k 基因 g 观测到的count,那么事实上对于其中两个物种的直系同源基因,我们需要对每一个基因的真实表达值(count值)假设一个分布,方便后续的假设检验,作者这里利用的是泊松分布
因此定义泊松分布的参数:
这里的泊松分布模型可以理解为对 Xgkt 做多次测量,最终对 Xgkt 做的一个频率分布(横坐标为 Xgkt ,纵坐标为频率)服从泊松分布
基于上面的模型,我们可以对H0做恒等变换,所以我们的假设问题就转变成为了:
所以满足H0的直系同源基因,我们认为是没有差异的;否则就是有差异的
跨物种RNA-seq差异分析p值计算
又由于当我们获得实际测的数据后,对于两个物种来说,满足:
即 Xgk1 + Xgk2 等于一个定值,因此我们可以引入伯努利实验的思想,构建二项分布:
- Xgk1 代表物种 t 中文库 k 基因 g 观测到的count
- xgk1 代表从 1—ngk 的count数
其中:
则p值计算如下:
所谓p值,本质上就是比较括号内前一项比后一项大的概率
Xgkt 表示实际观测到的在物种 1 中文库 k 基因 g 观测到的count数
xgk1 代表从 1—ngk 的count数;即计算 xgk1 取从 1—ngk 的count数时,(1)比(2)大的概率,即为p值
最后附上该文章的R包链接:SCBN
网友评论