歡迎關(guān)注”生信修煉手冊”!

MISO是一款經(jīng)典的可變剪切分析工具磅废，和rmats類似锅棕，該軟件也支持對可變剪切事件進(jìn)行定量和差異分析拙泽，網(wǎng)址如下

https://miso.readthedocs.io/en/fastmiso/index.html#

這個軟件支持exon和transcript兩種水平的可變剪切分析，在rmats的文章中裸燎，我們也提到了rmats是從exon水平給出的可變剪切結(jié)果顾瞻，因為二代測序讀長短的特點，無法有效得到轉(zhuǎn)錄本全長德绿，從exon水平得到的結(jié)果更加的準(zhǔn)確荷荤，而且陽性結(jié)果更容易通過RT-PCR驗證出來，但是無法詳細(xì)的探究某個基因不同isoform之間的變化移稳；transcript水平直接給出不同isoform間的定量和差異蕴纳，能有效的探究基因不同isofrm的變化情況，但是結(jié)果準(zhǔn)確性較差个粱。

該軟件是一個python包古毛，直接通過pip就可以安裝，分析的pipeline如下

1. 對參考基因組的GFF文件建索引

對于transcript水平的分析而言，只需要提供轉(zhuǎn)錄本的GFF文件稻薇，可以從Ensembl等數(shù)據(jù)庫下載參考基因組的gtf文件嫂冻，然后自己轉(zhuǎn)換成GFF3格式；對于exon水平而言塞椎，需要提供已知的可變剪切事件的GFF格式文件桨仿，示意如下

chr1 ?SE ? ? ?gene ? ?4772649 4775821 . ? ? ? - ? ? ? . ? ? ? ID=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-;Name=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-
chr1 ?SE ? ? ?mRNA ? ?4772649 4775821 . ? ? ? - ? ? ? . ? ? ? ID=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.A;Parent=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-
chr1 ?SE ? ? ?mRNA ? ?4772649 4775821 . ? ? ? - ? ? ? . ? ? ? ID=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.B;Parent=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-
chr1 ?SE ? ? ?exon ? ?4775654 4775821 . ? ? ? - ? ? ? . ? ? ? ID=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.A.up;Parent=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.A
chr1 ?SE ? ? ?exon ? ?4774032 4774186 . ? ? ? - ? ? ? . ? ? ? ID=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.A.se;Parent=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.A
chr1 ?SE ? ? ?exon ? ?4772649 4772814 . ? ? ? - ? ? ? . ? ? ? ID=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.A.dn;Parent=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.A
chr1 ?SE ? ? ?exon ? ?4775654 4775821 . ? ? ? - ? ? ? . ? ? ? ID=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.B.up;Parent=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.B
chr1 ?SE ? ? ?exon ? ?4772649 4772814 . ? ? ? - ? ? ? . ? ? ? ID=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.B.dn;Parent=chr1:4775654:4775821:-@chr1:4774032:4774186:-@chr1:4772649:4772814:-.B

第二列表示可變剪切的類型，以外顯子跳躍為例案狠，ID的格式如下

chr1:4775654:4775821:-@chr1:4774032:4774186:@chr1:4772649:4772814

包含了用@符號隔開的3個外顯子服傍，中間的exon的跳過的外顯子，第一個為上游的外顯子骂铁，第二個為下游的外顯子吹零，對應(yīng)如下示意圖中的3個exon

transcript水平的GFF文件從數(shù)據(jù)庫中下載即可，而exon水平的GFF文件是需要自己先識別可變剪切的不同isoform,然后整理得到的从铲，對于人和小鼠等常見物種瘪校，官網(wǎng)提供了exon水平的GFF文件，鏈接如下

https://miso.readthedocs.io/en/fastmiso/annotation.html

準(zhǔn)備好GFF文件之后名段，就可以建立索引了，命令如下

index_gff --index ensGene.gff3 index_db

index_db為索引保存的目錄泣懊。

2. 運(yùn)行miso

運(yùn)行miso需要第一步建好的索引以及樣本對應(yīng)的bam文件伸辟，該bam文件必須是經(jīng)過排序處理的，而且有對應(yīng)的bai索引馍刮，對于雙端數(shù)據(jù)信夫，用法如下

miso --run
index_db \
algin.sorted.bam \ ?
--output-dir out_dir \
--read-len 150 \
--paired-end 250 15 \
--settings-filename miso_settings.txt

read-len是reads的平均長度，paired-end代表插入片段長度的平均值和方差卡啰，miso_settings.txt是配置文件静稻，內(nèi)容如下

[data]
filter_results = True
min_event_reads = 20
strand = fr-unstranded
[sampler]
burn_in = 500
lag = 10
num_iters = 5000
num_processors = 4

配置文件中的參數(shù)很多，就不一一解釋了匈辱，每個參數(shù)的意義請參考官方文檔振湾。
通過上述方式得到的結(jié)果可以直接用于后續(xù)的差異分析，但是這個結(jié)果不利于我們查看亡脸，所以官方提供了匯總程序押搪，用法如下

summarize_miso \
--summarize-samples \
raw_out/ \
summary_out1

3. 樣本間的差異分析

進(jìn)行樣本間差異分析的代碼如下

compare_miso --compare-samples control case/ comparisons/

在輸出目錄，會生成一個后綴為bf的文件浅碾。

4. 對結(jié)果進(jìn)行過濾

用法如下

filter_events \
--filter ?case_vs_control.miso_bf \
--num-inc 1 \
--num-exc 1 \
--num-sum-inc-exc 10 \
--delta-psi 0.20 \
--bayes-factor 10 \
--output-dir filter_dir

5. 可視化

用法如下

sashimi_plot \
--plot-event "chr1:7778:7924:-@chr1:7096:7605:-@chr1:6717:6918:-" \
index_db/ \
sashimi_plot_settings.txt ?\
--output-dir out_dir

sashimi_plot_settings.txt是配置文件大州，其中設(shè)置了樣本的bam文件和可變剪切的輸出結(jié)果，示例如下

[data]
# directory where BAM files are
bam_prefix = ./test-data/bam-data/
# directory where MISO output is
miso_prefix = ./test-data/miso-data/

bam_files = [
 ? ?"heartWT1.sorted.bam",
 ? ?"heartWT2.sorted.bam",
 ? ?"heartKOa.sorted.bam",
 ? ?"heartKOb.sorted.bam"]

miso_files = [
 ? ?"heartWT1",
 ? ?"heartWT2",
 ? ?"heartKOa",
 ? ?"heartKOb"]

[plotting]
# Dimensions of figure to be plotted (in inches)
fig_width = 7
fig_height = 5
# Factor to scale down introns and exons by
intron_scale = 30
exon_scale = 4
# Whether to use a log scale or not when plotting
logged = False
font_size = 6

# Max y-axis
ymax = 150

# Whether to plot posterior distributions inferred by MISO
show_posteriors = True

# Whether to show posterior distributions as bar summaries
bar_posteriors = False

# Whether to plot the number of reads in each junction
number_junctions = True

resolution = .5
posterior_bins = 40
gene_posterior_ratio = 5

# List of colors for read denisites of each sample
colors = [
 ? ?"#CC0011",
 ? ?"#CC0011",
 ? ?"#FF8800",
 ? ?"#FF8800"]

# Number of mapped reads in each sample
# (Used to normalize the read density for RPKM calculation)
coverages = [
 ? ?6830944,
 ? ?14039751,
 ? ?4449737,
 ? ?6720151]

# Bar color for Bayes factor distribution
# plots (--plot-bf-dist)
# Paint them blue
bar_color = "b"

# Bayes factors thresholds to use for --plot-bf-dist
bf_thresholds = [0, 1, 2, 5, 10, 20]

最終會產(chǎn)生如下所示的結(jié)果

這種圖稱之為sashimi plot , 是一種專用于可變剪切可視化的圖表垂谢，上述示意圖表示的是一個外顯子跳躍事件在不同樣本中的表達(dá)情況厦画，左下方是GFF文件中共的exon結(jié)構(gòu)，左上方是每個樣本中比對上exon的reads的可視化滥朱，采用了RPKM表示根暑，不同剪切方式用曲線鏈接力试，曲線上標(biāo)記的是比對上該區(qū)域的reads數(shù)目，不同分組的樣本用不同顏色表示购裙，右側(cè)的圖片是樣本中對應(yīng)的可變剪切的表達(dá)量值懂版。

從這種圖中，可以直觀的看到兩組樣本間的可變剪切表達(dá)有無差異躏率，上圖中heartWT組中的表達(dá)量高于heartKO組躯畴。

實際分析時，由于需要手動整理可變剪切isofrom對應(yīng)的gff文件薇芝，所以使用的難度較大蓬抄，但是其提供的可視化功能是非常值得借鑒的。

·end·

—如果喜歡夯到，快分享給你的朋友們吧—

掃描關(guān)注微信號嚷缭，更多精彩內(nèi)容等著你！