基于贝叶斯框架的单细胞可变剪切算法

2020-08-03 本文已影响0人小潤澤

本推送同步上线于微信公众号“单细胞组学”

可变剪切事件

可变剪切指的是一个基因由于剪切方式的不同从而产生了不同的转录本。正常基因表达的前体RNA需要经过去内含子，外显子拼接等一系列步骤，最终形成mRNA。而同一基因的外显子拼接具有选择性，从而形成不同转录本。

上图中，左边的转录本具有1，2，3号外显子，而右边的转录本只有1，3号外显子，这就是可变剪切中的外显子跳跃。左边的拼接情况称为lnclusion isoform，右边的拼接情况称为Skipping isoform
那么，可变剪切具有多种方式，常见的分为以下7类：

可变剪切的研究意义

正是由于可变剪切的存在，使得同一个基因表达出不同的转录本，这也是导致蛋白质功能多样性的重要原因。因此，同一个基因由于转录本的不同，直接导致翻译出来的蛋白不同，其发挥的功能不同。间接导致生物体产生多种多样的性状，从而演化除了生物多样性。
另外一个方面，同一物种某些基因在雌雄间有着不同的可变剪切模式，在不易从外形区分性别的物种研究中，可以通过该基因雌雄间差异可变剪切来判断性别。

常用的可变剪切软件

对于不同的测序，有不同的鉴定可变剪切的软件。在bulk RNA-seq中，常用的鉴定可变剪切的软件有：MISO，rMAT，MAT等。
而在单细胞RNA-seq中，常用的鉴定可变剪切软件有BRIE，outrigger等。而目前大部分鉴定可变剪切的软件都基于贝叶斯模型框架（不同软件的细节处理可能有所不同）。那么我们就来看一下贝叶斯框架的大致原理。

贝叶斯因子

贝叶斯因子解释起来即为面对同一套观测值（数据），利用哪一个参数的概率更大。

贝叶斯框架

前面我们说到，目前检测可变剪切事件的软件（包括单细胞）大多都是基于贝叶斯框架下而制作的。
主要分为一下下几个步骤：

第一，软件先统计两个samples中测序的reads mapping到exon-exon junction处的counts数，统计结果包括lnclusion isoform和Skipping isoform两种情况的。接下来估算这两个samples中的exon inclusion levels（A图）。

第二，软件利用所有可变剪切的exon inclusion levels构造一个多变量的先验均匀分布（均匀分布的概率函数为正比例函数，与图B相符），这个信息主要衡量两个samples中可变剪切事件的总体相似性。（B图）

第三，软件计算每一个基因ψ1 - ψ2的值，并通过MCMC（拒绝抽样算法）对所有基因进行抽样，确定后验概率的最佳参数。由图中可以看到大部分的基因的ψ1 - ψ2都为0，所以大部分的基因都没有差异可变剪切。（C图）

第四，进行假设检验显著性分析，并确定阈值c（这个阈值c要满足仅有少部分基因服从 |ψ1 - ψ2| > c）。（D图）

接下来我们就一步一步具体看：

1. exon inclusion level

exon inclusion level 是衡量exon inclusion isoforms占总的 isoform（包括 inclusion isoforms 和 skipping isoforms）的比例
假设测序的read counts服从伯努利分布，那么 exon inclusion level 的极大似然估计可以表示为下图算式：

我们定义 :

I 是“ exon inclusion isoforms”的数量；

S 代表“ skipping isoforms” 的数量。

我们又定义：

UJC: 上游junction的数量；

DJC: 下游junction的数量；

SJC: skipping junction的数量。

2.计算贝叶斯后验概率

比较两个samples的可变剪切模式差异，我们分别定义sample 1 的exon inclusion level为ψ1，而sample 2 的exon inclusion level为ψ2。
于是假设原假设H0：| ψ1 - ψ2 | ≤ c ；备择假设H1：| ψ1 - ψ2 | ＞ c ，其中 c 是我们设定的阈值。
根据上一小节所说的，利用所有可变剪切的exon inclusion levels构造一个多变量的先验均匀分布，满足于下图：

其拟合模型图像为（图像上的点为不同的基因）：

其中，横坐标为ψ1，纵坐标为ψ2
那么根据阈值 c ，我们在这个多元的先验均匀分布中就可以划分出H0区域和H1区域，如下图：

那么当两个samples中某基因的exon inclusion levels满足这个先验均匀分布时，那么该基因是没有发生差异可变剪切的。

所以，所得的这个多元的先验均匀分布，满足以下关系：

而ψ1和ψ2的相关性满足均匀分布uniform (0,1)

对于P(D | ψ1 - ψ2)分布，软件认为满足伯努利分布，其中伯努利分布的概率为ψ。
可以这样理解，对于某个基因A来说，软件经过计算会得到一个ψ值，假设比对到A基因的reads有100条，那么每一个reads就有两种可能性，一种是exon inclusion isoforms，那么另一种是skipping isoforms，概率为ψ。正好满足伯努利实验，即做100次重复实验，有两种可能性，其中一种的概率为ψ。因此可以看成是已知概率为ψ的条件下，每个reads为exon inclusion isoforms或者是skipping isoforms的伯努利概率分布。
那么Ii1 | ψi1表示已知某基因的ψi1值的条件下，其中某一条reads为exon inclusion isoforms的似然值；Ii2 | ψi2表示已知某基因的ψi2值的条件下，其中某一条reads为exon inclusion isoforms的似然值。

根据伯努利分布的线性关系，我们可以得到似然函数P(D | ψ1)和P(D | ψ2)的分布。

因此，先验概率表示为 P(ψ)，这是具有普遍性的，那么，如果我任意给一套Data（简称：D），这里的
Data, D 是可变剪切事件的集合：

集合D
集合D包括 I 的count数和 S 的count数

对于某一个基因，我们的目的是通过该基因的先验分布和似然函数计算 P(ψ | D) 的期望（均值），以及对应的ψ的均值：

并求得 P(ψ1-ψ2 | D)均值 ，以此判断是否发生可变剪切。

上图的常数c为高维积分，很难计算，因此通常采用MCMC的方法估算后验概率的期望。

MCMC—MH算法（插播）：

MCMC(MH)

由公式：

我们可以得到简单的一个正比关系：

由于这种正比关系的存在，后验分布的形状与先验分布和似然函数的乘积所构成的分布函数形状相似。

而 P( ψ )分布可以从所有基因的 ψ 值来拟合

于是ψ 服从P(D | ψ)的分布，在P(D | ψ)P(ψ)采样，前后两个采样状态的参数ψ一定满足下式：

{注：这里的ψ(new)和ψ(current)是同一分布下的不同参数(或不同自变量)}

那么，求后验概率分布P(ψ | D)的关键是求出参数ψ，根据前面说的正比关系，我们可以在该分布中：

取样，从而模拟 P(ψ | D) 的分布。并且M-H算法的判别是否取样的条件为：

因此经过MCMC多次迭代后有：

那么，取迭代中的最优值作为后验概率分布参数 ψ 的均值，并且该ψ值作为后验概率分布P(ψ | D)中的ψ值（唯一值）。根据MCMC总的采样情况并做拟合，就可以得到后验概率分布的期望 P(ψ | D) 。
P(ψ | D)这个后验分布表示，已知某个基因 I 的count值和 S 的count值，而它的 ψ为某一个值的概率。那么对于某个基因来说，我们选后验抽样分布的 均值ψ 作为该基因的 ψ值，从而得到该基因的后验概率分布的期望 P(ψ | gene)

以DNMT3B为例，抽样结果为：

其中左边是用sashimi plot可视化的结果；右边蓝色曲线表示先验的分布，直方图表示后验抽样分布。
对于其中某一个基因（DNMT3B），已知该基因 I 的count值和 S 的count值，其 ψ 的后验抽样分布。那么，根据后验抽样分布，我们可以得到 ψ 的均值，即图中红色竖线对应的 ψ 值。所以我们就确定了DNMT3B的后验概率分布的期望 P(ψ | DNMT3)

回到下图：

根据贝叶斯理论，我们得到一套数据D，并且给定阈值 c，根据后验概率分布的期望 P(ψ1-ψ2 | D) 就可以判断某个基因是否存在差异可变剪切。即某个基因在已知 I 的count值和 S 的count值的条件下，| ψ1-ψ2 | ≤ c的概率

对于某个基因，这里的 D 特指某一个基因（或者某一基因的可变剪切isoform的情况），比较下面的大小：

如果计算出来的| ψ1-ψ2 | ≤ c那么就没有发生差异可变剪切；如果计算出来的| ψ1-ψ2 | > c那么就发生差异可变剪切.
我们依据这个就可以判断两个samples中，这个是否有差异可变剪切事件了。

3.如何确定阈值 c

比方说软件对5万个基因进行可变剪切计算，那么对应就会有5万个 ψ1-ψ2 ，软件通过做 ψ1-ψ2 的频率分布直方图，来确定阈值 c；首先，软件认为大部分基因是不存在差异可变剪切的，因此存在差异可变剪切的只有一小部分，因此在频率分布直方图里那一小部分对应的ψ1-ψ2就为阈值c

图中阴影部分即为那一小部分的差异可变剪切基因，虚线处即代表阈值 c

主体思想

我们说P(gene | ψ)的ψ只能看作是总体 ψ 的一个样本，并不一定是最优的ψ。并且总体的ψ分布还不知道，因此用MCMC去估计总体的ψ值(最优ψ)，从而求得P(ψ | gene)，这个思想可以类比二元函数来理解固定其中一个变量讨论另一个变量来理解，当给定ψ的条件下讨论gene的分布，随着ψ不同，那么gene的分布也就不同，因此我们要求解的是最优的ψ值，故用MCMC去估计

单细胞可变剪切软件

BRIE就是根据贝叶斯模型进行统计推断的一款单细胞可变剪切分析软件。主要是判断两个cell cluster间的差异可变剪切事件。常用的软件有BRIE，其基本用法：

#1. Generate exon-skipping events from gene annotation.
briekit-event -a gencode.vM17.annotation.gtf.gz -o
$DATA_DIR/AS_events

#2. Filter splicing events.
briekit-event-filter -a AS_events/SE.gff3.gz --anno_ref
gencode.vM17.annotation.gtf.gz \
-r GRCm38.p6.genome.fa

#3.Extract the sequence features.
briekit-factor -a AS_events/SE.filtered.gff3.gz -r
GRCm38.p6.genome.fa \
-c mm10.60way.phastCons.bw -o mouse_features.csv -p 10 --
bigWigSummary ./bigWigSummary

#4.Quantify the annotated exon�skipping events on the aligned reads files
brie -a SE.filtered.gff3.gz -s cell_n.sorted.bam -f mouse_-
features.csv.gz -o $OUT_DIR -p 10

#5. All cells with comma separation, and for setting shared weights for individual cells it is just by using -w
brie -a SE.filtered.gtf -f mouse_features.csv.gz \
-s cell1.sorted.bam,cell2.sorted.bam,cell3.sorted.bam -o
$OUT_DIR/cell_1to4 -p 10

brie -a SE.filtered.gtf -f mouse_features.csv.gz \
-s cell1.sorted.bam -w $OUT_DIR/cell_1to4/weights.tsv -o $OUT_DIR/cell1_fixedW -p 10

#6.Detect differential splicing between any pair of cells (or cell clusters)
$fileList=cell1/samples.csv.gz,cell2/samples.csv.gz,cell3/
samples.csv.gz,cell4/samples.csv.gz

brie-diff -i $fileList -o cell1_cell4.diff.tsv

参考文献：
[1].Huang Y , Sanguinetti G . BRIE: transcriptome-wide splicing quantification in single cells[J]. Genome Biology, 2017, 18(1).
[2].Katz Y , Wang E T , Airoldi E M , et al. Analysis and design of RNA sequencing experiments for identifying isoform regulation[J]. Nature Methods, 2010, 7(12):1009.
[3].Shihao S , Won P J , Jian H , et al. MATS: a Bayesian framework for flexible detection of differential alternative splicing from RNA-Seq data[J]. Nuclc Acids Research, 2012(8):e61.
[4].Yuanhua H , Guido S . Statistical modeling of isoform splicing dynamics from RNA-seq time series data[J]. Bioinformatics(19):2965.
[5].Trapnell C , Williams B A , Pertea G , et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[J]. Nature Biotechnology, 2010, 28(5):511-515.
[6].Bayesian Inference: Gibbs Sampling
[7].Bayesian Inference: Metropolis-Hastings Sampling
[8].Computational Methods for Single-Cell Data Analysis
[9].Deep-learning augmented RNA-seq analysis of transcript splicing
[10]. MCMC与贝叶斯估计