rMATS可视化-rmats2sashimiplot

之前的对rMATS软件做笔记的时候，提到过rmats2sashimiplot这款专门用来对rMATS分析结果做可视化的软件，最近尝试用了下，操作比较简洁，但是还是留了点问题没解决

rmats2sashimiplot，网址：rmats2sashimiplot

下载

git clone https://github.com/Xinglab/rmats2sashimiplot

or

https://pypi.org/project/rmats2sashimiplot/2.0.2/#files

安装

python setup.py install

ubuntu 16.04 报错：error while loading shared libraries: libgsl.so.0: cannot open shared object file: No such file or directory

这是因为ubuntu 16.04不支持libgsl0ldbl，而是替换为libgsl2，所以先安装libgsl2，然后在将其软链接/复制到/usr/lib目录下

sudo apt-get install libgsl-dev
dpkg -L libgsl2
sudo cp /usr/lib/x86_64-linux-gnu/libgsl.so.19 /usr/lib/libgsl.so.0

解决方法参照自：https://jingyan.baidu.com/article/63acb44a13cc7d61fdc17e6d.html

如果直接使用rMATS从fastq输入文件到结果输出文件，那么rMATS是不保留中间的bam文件的

后来看了下rmats.py中的代码，发现只要将shutil.rmtree(args.tmp)注释掉就会保留STAR比对的bam结果文件了，但是其是放在服务器上/tmp目录中，也不太方便

而rmats2sashimiplot可视化则需要bam文件作为输入，所以需要我们先用STAR比对得到bam文件再用rMATS做差异可变剪切分析，如下所示：

用STAR对各个样本做比对生成bam文件，比对参数参考rMATS软件调用STAR时所用参数，其实就是比默认比对参数多了：

--chimSegmentMin 2 --outFilterMismatchNmax 3 --outSAMstrandField intronMotif --alignEndsType EndToEnd --alignSJDBoverhangMin 6 --alignIntronMax 299999

STAR --readFilesIn A1_1.clean.fq A1_2.clean.fq --chimSegmentMin 2 --outFilterMismatchNmax 3 --alignEndsType EndToEnd --runThreadN 12 --outSAMstrandField intronMotif --outSAMtype BAM SortedByCoordinate --alignSJDBoverhangMin 6 --alignIntronMax 299999 --genomeDir ~/index/STAR/Ensembl/human/GRCh38.92/ --sjdbGTFfile ~/annotation/Ensembl/human/Homo_sapiens.GRCh38.92.chr.gtf --outFileNamePrefix A1

rMATS做差异可变剪切分析（例如A.txt包含了A组bam的文件名）

python rmats.py --b1 A.txt --b2 B.txt --gtf Homo_sapiens.GRCh38.92.chr.gtf --od ./test -t paired --readLength 150 --cstat 0.0001 --nthread 10

rmats2sashimiplot作图

rmats2sashimiplot --b1 A1.bam,A2.bam,A3.bam --b2 B1.bam,B2.bam,B3.bam -t SE -e ./SE.MATS.JC_top20.txt --l1 A --l2 B --exon_s 1 --intron_s 1 -o SE_plot

rmats2sashimiplot的注意事项：

可以预先用samtools对bam文件建索引，不然rmats2sashimiplot也会先建索引（比较慢）
rMATS的结果文件中，如SE.MATS.JC.txt中染色体都是默认加上chr的，而bam文件中的染色体根据基因组注释文件来源不同，有些的没chr，而是直接以数字表示，这时需要预先将两者保持一致
rmats2sashimiplot默认是出PDF图片的，如果需要出png，则需要自己修改下源码（由于我也不懂python作图，所以用了最无脑的方式就是将所有跟出图相关的脚本中pdf都替换为png。。。）

rMATS可视化的图片如下：

rmats2sashimiplot

但是有个问题我还是没想清楚，按照rmats2sashimiplot作者的意思，这图中的线上数字代表着 junction-spanning reads，那么这跟rMATS的分析结果中的IC_SAMPLE_1和SC_SAMPLE_1一致吗？

比如这张图中第一个样本A-1的Exon Inclusion Isoform的reads数分别为176和159，exon skipping isoform的reads为2

而rMATS分析结果中IC_SAMPLE_1的A-1样本数值为259，SC_SAMPLE_1为2；

我的理解图中Exon Inclusion Isoform上的两个数值加和应该等于259（有些图片确实相等），但事实好多图片都与我所想不一致，有时相差的还很奇怪，想不明白。。。

之前在用rMATS时，对其结果的文件的表头并不是理解的很清楚，最近查了下资料（主要是rMATS的原理文献，反正也没完全看懂），再整理下：

IncLevel1可被认作是exon inclusion level（ψ），是Exon Inclusion Isoform在总（Exon Inclusion Isoform + Exon Skipping Isoform）所占比例

For a skipped exon, the exon inclusion level (denoted as percent spliced in or ψ) is calculated by the number of reads uniquely mapped to the exon inclusion isoform or the exon skipping isoform

因此计算公式如下：

ψ = (I/lI)/(S/lS)

I 是指mapping到exon inclusion isoform上的reads数，对应结果表格中的IC_SAMPLE_1
S 是指mapping到exon skipping isoform上的reads数，对应结果表格中的SC_SAMPLE_1
lI 是指effective length of the exon inclusion isoform，对应结果表格中的IncFormLen，主要用于对I标准化，用于计算ψ
lS 是指effective length of the exon skipping isoform，对应结果表格中的SkipFormLen，主要用于对S标准化，用于计算ψ

对于effective length的计算公式则是：

For a segment whose length is l, and the length of the read is r, the effective length of the segment is defined by the number of unique read intervals in this region, which is l − r + 1

这里的计算可以参照https://www.biostars.org/p/256949/，我还是没搞清楚。。。

IncLevelDifference则是指两组样本IncLevel的差异，如果一组内多个样本，那么则是各自组的均值之间差值

rMATS采用的是Likelihood-ratio test来比较两组样本可变剪切是否有差异

If the maximum-likelihood estimations (MLEs) of ψi1,ψi2 have a difference smaller than or equal to the user-defined cutoff (i.e., |ψi1 −ψi2|≤ c), we set the P value to be 1

所以跟我们参数--cstat设置有关（默认是0.0001），基于marginal distribution（在零假设前提下，近似χ2 distribution），备择假设是|ψi1 −ψi2|> c，因此粗略的理解就是ψi1和ψi2差异比c越大P值就越小

虽然上面作图可视化的问题还是没搞清楚，但是也记录下了，之后如果搞清楚了再回来修改。。