如何使用MISO進行可變剪切的分析

如何使用MISO進行可變剪切的分析，很多新手對此不是很清楚，為了幫助大家解決這個難題，下面小編將為大家詳細講解，有這方面需求的人可以來學習下，希望你能有所收獲。

MISO是一款經典的可變剪切分析工具，和rmats類似，該軟件也支持對可變剪切事件進行定量和差異分析。

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

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

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

對于transcript水平的分析而言，只需要提供轉錄本的GFF文件，可以從Ensembl等數據庫下載參考基因組的gtf文件，然后自己轉換成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的跳過的外顯子，第一個為上游的外顯子，第二個為下游的外顯子，對應如下示意圖中的3個exon

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

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

準備好GFF文件之后，就可以建立索引了，命令如下

index_gff --index ensGene.gff3 index_db

index_db為索引保存的目錄。

2. 運行miso

運行miso需要第一步建好的索引以及樣本對應的bam文件，該bam文件必須是經過排序處理的，而且有對應的bai索引，對于雙端數據，用法如下

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是配置文件，內容如下

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

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

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

3. 樣本間的差異分析

進行樣本間差異分析的代碼如下

compare_miso --compare-samples control case/ comparisons/

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

4. 對結果進行過濾

用法如下

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是配置文件，其中設置了樣本的bam文件和可變剪切的輸出結果，示例如下

[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]

最終會產生如下所示的結果

這種圖稱之為sashimi plot , 是一種專用于可變剪切可視化的圖表，上述示意圖表示的是一個外顯子跳躍事件在不同樣本中的表達情況，左下方是GFF文件中共的exon結構，左上方是每個樣本中比對上exon的reads的可視化，采用了RPKM表示，不同剪切方式用曲線鏈接，曲線上標記的是比對上該區域的reads數目，不同分組的樣本用不同顏色表示，右側的圖片是樣本中對應的可變剪切的表達量值。

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

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

看完上述內容是否對您有幫助呢？如果還想對相關知識有進一步的了解或閱讀更多相關文章，請關注億速云行業資訊頻道，感謝您對億速云的支持。