Deep-learning augmented RNA-seq analysis of transcript splicing

0. 简介

这篇文章是发表在Nature Methods上的使用深度学习来扩充可变剪切的分析。通讯作者是来自于UCLA的邢毅老师。这篇文章创新点在于将贝叶斯假设检验框架和深度学习结合在一起，首先使用贝叶斯假设检验来检验可变差异剪切，并用来作为训练模型的标签，然后根据深度学习预测的结果重新作为贝叶斯假设检验的先验概率从而估计后验概率；除此之外，作者针对于类别不平衡预测结果还做了rank-transformation。

文章链接
 文章代码

1. 摘要

可变剪切的RNA-Seq分析最大的局限在于非常依赖测序深度；
作者在这里提出DARTS, 通过整合先验的RNA-Seq证据以及深度学习预测结果来推断不同生物学样本中的差异可变剪切；
DARTS利用公共数据库中的大量的RNA-Seq数据通过深度学习提供可变剪切调节的knowledge base, 因此可以用来帮助研究人员即使在中等测序深度的情况下更好的刻画可变剪切。

2. 结果

2.1 DARTS方法简介

RNA-Seq测序使得分析转录组水平的可变剪切成为可能，RNA-Seq的快速积累也为分析可变剪切提供了前所未有的机遇。但是可变剪切的分析却受限于测序深度。
因此作者提出了DARTS，如下图所示：DARTS由两个成份组成，DNN（Deep neural network）基于外显子特异性的序列特征以及样本特异性的调控特征预测两个条件是不是差异可变剪切；BHT(Bayesian hypothesis testing)统计模型通过整合特定RNA-Seq数据的经验证据（由DNN预测得来）与可变剪切的先验概率推断差异可变剪切；
在训练的时候，首先使用无信息先验模型只基于RNA-Seq数据生成高置信度差异以及非差异的标签，然后再用于训练，然后预测的结果又用于DARTS BHT（info），用预测结果作为先验；
Figure 1B: 与其他方法不同的是，DARTS不仅使用了cis sequence还整合了两种情况下trans RNA-binding proteins(RBPs)的mRNA水平。这可以允许DARTS分析RBPs的表达如何影响剪切，作者使用了2926个序列特征以及1498个注释的RBPs的mRNA水平（trans RBP特征）；

20190327-Fig2.png

- Figure 1C: 作者使用了两个cell lines(K562和HepG2)对应的RBP-depletion RNA-seq，在两个cell lines中都有的196个 RBPs depleted by short-harpin RNA (shRNA), 对应着408对（knockdown-versus-control）;对于只在一个细胞系中出现的58组数据用于leave-out test;
Figure 1D: DARTS在leave-out数据集上的表现，与其他模型进行比较；

2.2 DARTS BHT框架的性能评估以及训练数据集对于DARTS DNN的影响

Figure 2A: DARTS BHT (info)与DARTS BHT (flat)的性能比较, 测序深度越浅，性能提升的越显著；
Figure 2B: DARTS在leave0out数据集上的预测效果；

2.3 DARTS在上皮-间质转化中的可变剪切分析

Figure 3A: 将DARTS DNN用于时序的RNA0Seq数据，底部的bar plot是无先验模型的结果，线图是DNN的预测结果；

3. 方法

3.1 DARTS贝叶斯假设检验框架

这部分使用贝叶斯假设检验来确定差异可变剪切事件，DARTS的核心是先利用贝叶斯假设检验确定样本的标签（differential splicing or unchanged splicing）,然后使用深度学习建模预测，把深度学习预测结果作为先验信息，再结合RNA-Seq数据（实验证据）得到样本最终的标签。以最简单的情况（检测两个条件下exon-inclusion levels; PSI values）为例, 认为在特定PSI的情况下，外显子保留的read count符合二项分布，其中 $n$ 为外显子保留的read数目加上外显子跳跃read数目，概率 $p=f_i(\psi_{ij})\$ ，详情如下：
**首先声明: ** $i$ 表示第 $i$ 个外显子， $j\in{1,2}$ 分别代表差异splicing和unchanged splicing；
$I_{ij}|\psi_{ij}\sim Binomial(n=I_{ij}+S_{ij},p=f_i(\psi_{ij}))$
$\psi_{i1}=\mu_i$
$\psi_{i2}=\mu_i+\delta_i$
$\mu_i\sim Unif(0,1)$
$\delta_i\sim N(0,\tau^2)$
$I_{ij}$ : 外显子保留的read count；
$S_{ij}$ : 外显子跳跃的read count；
$\psi_{ij}$ : 外显子 $i$ 在sample group $j\in{1,2}$ 的外显子保留水平,也就是PSI值；
$f_i$ : 外显子 $i$ 的长度标准化函数，这个函数考虑了外显子保留和跳跃的有效长度；
$\mu_i$ : 外显子 $i$ 的baseline inclusion level；
$\delta_i$ : 两个条件下期望的外显子保留水平之差；
差异可变剪切分析的目标就是检验两个条件下是否有很高的概率使得 $\delta_i$ 大于指定的阈值 $c$ (例如：0.05)。对应如下公式：
$p(|\delta_i|>c|I_{ij},S_{ij})\approx1$
在贝叶斯统计中，可以通过使用参数 $\delta$ 的spike-and-slab先验来估计结果。spike-and-slab先验是双组分混合分布，其中“spike”描绘了模型参数 $\delta$ 被约束在零附近的概率，而“slab”分量描绘了模型参数 $\delta$ 的无约束分布。
在贝叶斯统计框架中，为了考虑到随机的技术干扰对PSI值的影响，作者利用了spike prior $H_0$ with a small variance $\tau=\tau_0$ ; slab prior $H_1$ with a small variance $\tau=\tau_1$ 。设置 $\tau_0=0.03$ , 对应着90%的密度约束在 $\delta\in[-0.05,0.05]$ , $\tau_1=0.3$ , 因此后验概率可以写做：
$p(H_1|I_{ij},S_{ij})=\frac1Zp(H_1)*p(I_{ij},S_{ij}|H_1)$
$p(I_{ij},S_{ij}|H_1)=\int_\tau\int_\mu p(I_{ij},S_{ij}|\mu_i,\delta_i)*p(\mu_i,\delta_i|H_1)d\mu_id\delta_i$
$p(H_0|I_{ij},S_{ij})=\frac1Zp(H_0)*p(I_{ij},S_{ij}|H_0)$
$p(I_{ij},S_{ij}|H_0)=\int_\tau\int_\mu p(I_{ij},S_{ij}|\mu_i,\delta_i)*p(\mu_i,\delta_i|H_0)d\mu_id\delta_i$
说明：
- $p(H_1)$ ：外显子 $i$ 差异剪切的先验概率，在模型训练之前这个采用无信息先验分布，训练模型之后用预测分值作为先验信息, 因此这里还是无信息先验；
- $p(H_0)=1-p(H_1)$ : 外显子 $i$ unchanged的先验概率；
- $p(I_{ij},S_{ij}|H_1)$ 和 $p(I_{ij},S_{ij}|H_0)$ 分别表示在differential splicing and unchanged splicing模型下的likelihoods；
- $Z$ ：标准常数；
由于是比较两个模型，因此又可以写成以下形式：
$\frac{p(H_1|I_{ij},S_{ij})}{p(H_0|I_{ij},S_{ij})}=\frac{p(H_1)}{p(H_0)}*\frac{p(I_{ij},S_{ij}|H_1)}{p(I_{ij},S_{ij}|H_0)}$
由于采用无信息先验，所以： $p(H_0)=p(H_1)=0.5$ ; 因此上式等价于似然比检验，这里作者称为DARTS BHT (flat); 当整合基于深度学习的预测结果作为先验信息时，作者称为DARTS BHT (info)。
最终，使用上述公式，我们可以得到参数 $\delta_i$ 的边际后验概率 $p(\delta_i|I_{ij},S_{ij})$ ，即两个模型的后验混合：
$p(\delta_i|I_{ij},S_{ij})=p(\delta_i|H_1,I_{ij},S_{ij})*p(H_1|I_{ij},S_{ij})+p(\delta_i|H_0,I_{ij},S_{ij})*p(H_0|I_{ij},S_{ij})$
最终就是推断 $p(|\delta_i|>c|I_{ij},S_{ij})$ , 作者在分析中，设置 $c=0.05$ , 即：外显子保留水平有5%的变化。认为 $p(|\delta_i|>0.05|I_{ij},S_{ij})>0.9$ 为显著差异剪切事件； $p(|\delta_i|>0.05|I_{ij},S_{ij})<0.1$ 为显著的unchanged splicing事件。

3.2 DARTS预测差异可变剪切的深度神经网络模型(DARTS DNN model for predicting differential alternative splicing)

DARTS的一个核心就是利用DNN检测两个条件下一个特定的外显子是不是差异剪切，这里作者利用了两个层面的特征，第一个是序列层面的特征，第二个是在两个条件下样本特异性的反式RBP表达水平的特征，因此预测问题可以表述为以下公式：
$p(Y_{ik}=1)=F(Y_{ik};E_i,G_k)$
其中 $Y_{ik}$ 就是外显子 $i$ 在比较 $k$ 中的标签； $E_i$ 是一个向量，包含2936个进化保守型特征、2973个剪切点的长度、2971个调控motif的组成以及1748个RNA二级结构特征（ skipped exons, alternative 5′ splice sites, alternative 3′ splice sites, and retained introns）；而 $G_k$ 是2996维的向量代表1498个RBPs在两个条件下的标准化基因表达量；DNN的预测结果可以背进一步整合为贝叶斯假设检验框架中的先验信息；
模型由4个全连接隐层分别包含1200，500，300和200个神经元均使用ReLu激活函数；RMSprop优化，minibatch设置为1000，加入Dropout防止过拟合。

3.3 RNA-seq数据处理以及DARTS DNN模型的训练

使用rMATS处理bam文件检测可变剪切；
Kallisto定基因表达量，使用Gencode生成index；
使用DARTS BHT(flat)生成训练标签，后验概率 $p(|\Delta\psi|>0.05)>0.9$ 作为正样本；后验概率 $p(|\Delta\psi|>0.05)<0.1$ 作为负样本；
10%的样本用于测试数据，剩下90%数据用于五倍交叉验证；
同时也收集了只在一个细胞系中存在的ENCODE RBP-knockdown实验作为leave-out datasets来评估模型的精度；

3.4 DARTS信息先验的Rank-transformation

由于unchanged splicing要比differential splicing大几个数量级，因此机器学习算法会偏向于类别多的样本，作者为了解决这个问题，使用了非监督的rank-tranformation重新scale DARTS DNN的预测分值来获得DARTS BHT框架的先验信息；
具体来说，作者将所有DARTS DNN的预测分值拟合成的双峰（两个component）的高斯混合模型，然后得出两个混合高斯峰的均值和方差，以及属于特定成分的每个DARTS DNN得分的后验概率 $\lambda$ ；
每个高斯成分的均值和方差分别设置为: $\mu_0$ , $\mu_1$ , $\sigma_0$ , $\sigma_1$ ，然后每个原始预测分值通过rank-transformed成新的高斯成分然后通过权重参数 $\lambda$ 加权平均；
最后，为了保持有限的先验概率，把transformed DARTS DNN score scale到 $[\alpha, 1-\alpha]$ , 其中 $\alpha\in[0,0.5)$ ；
在实际操作中，用到的参数设置为： $\mu_0=0.05$ , $\mu_1=0.95$ , $\sigma_0=\sigma_1=0.1$ , $\alpha=0.05$ 。