本推送同步上線于微信公眾號“單細胞組學”

可變剪切事件

可變剪切指的是一個基因由于剪切方式的不同從而產生了不同的轉錄本。正?；虮磉_的前體RNA需要經過去內含子，外顯子拼接等一系列步驟，最終形成mRNA。而同一基因的外顯子拼接具有選擇性，從而形成不同轉錄本。

可變剪切的研究意義

正是由于可變剪切的存在，使得同一個基因表達出不同的轉錄本，這也是導致蛋白質功能多樣性的重要原因。因此，同一個基因由于轉錄本的不同，直接導致翻譯出來的蛋白不同，其發(fā)揮的功能不同。間接導致生物體產生多種多樣的性狀，從而演化除了生物多樣性。
另外一個方面，同一物種某些基因在雌雄間有著不同的可變剪切模式，在不易從外形區(qū)分性別的物種研究中，可以通過該基因雌雄間差異可變剪切來判斷性別。

常用的可變剪切軟件

對于不同的測序，有不同的鑒定可變剪切的軟件。在bulk RNA-seq中，常用的鑒定可變剪切的軟件有：MISO，rMAT，MAT等。
而在單細胞RNA-seq中，常用的鑒定可變剪切軟件有BRIE，outrigger等。而目前大部分鑒定可變剪切的軟件都基于貝葉斯模型框架（不同軟件的細節(jié)處理可能有所不同）。那么我們就來看一下貝葉斯框架的大致原理。

貝葉斯因子

貝葉斯因子解釋起來即為面對同一套觀測值（數據），利用哪一個參數的概率更大。

貝葉斯框架

前面我們說到，目前檢測可變剪切事件的軟件（包括單細胞）大多都是基于貝葉斯框架下而制作的。
主要分為一下下幾個步驟：

第一，軟件先統計兩個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 的極大似然估計可以表示為下圖算式：

我們定義 :

1. I 是“ exon inclusion isoforms”的數量；
2. S 代表“ skipping isoforms” 的數量。

我們又定義：

1. UJC: 上游junction的數量；
2. DJC: 下游junction的數量；
3. 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 是我們設定的閾值。
根據上一小節(jié)所說的，利用所有可變剪切的exon inclusion levels構造一個多變量的先驗均勻分布，滿足于下圖：

所以，所得的這個多元的先驗均勻分布，滿足以下關系：

對于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(ψ)，這是具有普遍性的，那么，如果我任意給一套Data（簡稱：D），這里的
Data, D 是可變剪切事件的集合：

集合D

由于 P(D | ψ1) 和 P(D | ψ2) 無法很好的計算 P(D | ψ1-ψ2) ，并且 P(D | ψ1-ψ2) 也不好直接的判斷某基因是否發(fā)生差異可變剪切。而后驗概率分布P(ψ1-ψ2 | D)可以較好判斷某基因是否發(fā)生差異可變剪切。所以依據貝葉斯理論，需要求出后驗概率分布 P(ψ1-ψ2 | D)。

對于某一個基因，我們的目的是通過該基因的先驗分布和似然函數計算 P(ψ | D) 的期望（均值），以及對應的ψ的均值：

接下來，為了便于計算和判斷，我們需要計算后驗概率P(ψ1 | D)和P(ψ2 | D)的期望，才能計算出P(ψ1-ψ2 | D)進行判斷。所以接下來的問題是如何計算P(ψ1 | D)和P(ψ2 | D)。
那么我要計算后驗概率的期望，即P(ψ1 | D)和P(ψ2 | D)的期望才能夠判斷這套數據是否有可變剪切事件發(fā)生。那么在正常情況下，我們求解P(ψ1 | D)和P(ψ2 | D)的期望是非常困難的：

MCMC—MH算法（插播）：

MCMC(MH)

由公式：

而 P( ψ )分布可以從所有基因的 ψ 值來擬合

于是ψ 服從P(D | ψ)的分布，在P(D | ψ)P(ψ)采樣，前后兩個采樣狀態(tài)的參數ψ一定滿足下式：

那么，求后驗概率分布P(ψ | D)的關鍵是求出參數ψ，根據前面說的正比關系，我們可以在該分布中：

以DNMT3B為例，抽樣結果為：

因此通過上述計算出 P(ψ1 | D)和P(ψ2 | D) 這兩個后驗概率分布的期望，而通過這兩個概率分布以及概率分布運算性質，在相同條件下可以求出P(ψ1-ψ2 | D)分布的期望。因此，給定相應的閾值 c ，我們依據后驗概率的似然值就可以計算當某個基因在已知 I 的count值和 S 的count值的條件下，| ψ1-ψ2 | ≤ c的概率，以此判斷某個基因是否存在差異可變剪切

回到下圖：

對于某個基因，這里的 D 特指某一個基因（或者某一基因的可變剪切isoform的情況），比較下面的大?。?br>

3.如何確定閾值 c

比方說軟件對5萬個基因進行可變剪切計算，那么對應就會有5萬個 ψ1-ψ2 ，軟件通過做 ψ1-ψ2 的頻率分布直方圖，來確定閾值 c；首先，軟件認為大部分基因是不存在差異可變剪切的，因此存在差異可變剪切的只有一小部分，因此在頻率分布直方圖里那一小部分對應的ψ1-ψ2就為閾值c

圖中陰影部分即為那一小部分的差異可變剪切基因，虛線處即代表閾值 c

主體思想

我們說P(gene | ψ)的ψ只能看作是總體 ψ 的一個樣本，并不一定是最優(yōu)的ψ。并且總體的ψ分布還不知道，因此用MCMC去估計總體的ψ值(最優(yōu)ψ)，從而求得P(ψ | gene)，這個思想可以類比二元函數來理解固定其中一個變量討論另一個變量來理解，當給定ψ的條件下討論gene的分布，隨著ψ不同，那么gene的分布也就不同，因此我們要求解的是最優(yōu)的ψ值，故用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與貝葉斯估計