可变剪切的可视化 ggsashimi.py

从github https://github.com/guigolab/ggsashimi 下载ggsashimi.py
安装对应的依赖包和环境，python和R以及里面的包

一个区间的可变剪切可视化

ggsashimi.py -b input_bams.tsv -c chr10:27040584-27048100 -g annotation.gtf -M 10 -C 2 -O 3 --shrink --alpha 0.25 --base-size=20 --ann-height=4 --height=3 --width=18 -P palette.txt

多个基因的可变剪切可视化

genes.list里面是tab分割的两列，第一列是基因ID,第2列是对应的基因的区间，格式是chr10:10000-200000

cat genes.list|while read gene region;do
ggsashimi.py -b input_bams.tsv -c ${region} -M 10 -C 3 -O 3 --shrink --alpha 0.25 --base-size=20 --ann-height=4 --height=3 --width=18 -P palette2.txt --fix-y-scale -o ${gene} -F pdf
done

-b 指定输入bam文件的信息
input_bams.tsv要求是三列，第一列是ID信息，第二列是bam文件的路径位置，第三列是分组信息。文件里可以放很多个bam的信息。
-c是要绘制的区间。格式是chr10:10000-200000
-g 基因的注释gtf文件,此文件必须包含exon信息，也可以不提供此文件,则不会绘制转录本
-M 画图时，最小覆盖的reads数量，默认是1
-C 绘图时的颜色分组
-O 列的索引
--shrink 将连接处缩小为一个系数，以获得更好的显示效果。该参数指定则收缩，不指定则不收缩。后面没有参数值
--alpha 密度直方图的透明度，默认是0.5
-base-size 基本字号
--ann-height 下面转录本的高度 英寸 默认是1.5
--height 每个bam样本的高度 英寸，默认是2
--width 宽度 英寸，默认是10
-P palette.txt 里面是一列不同的颜色，用于绘图的颜色。支持R的颜色名称和16进制的颜色名称
-F 输出的图片类型，支持svg,pdf,jpeg,png,tiff

image.png

二代可变剪切的鉴定

cash鉴定可变剪切位点

安装下载简单
https://sourceforge.net/projects/cash-program/files/2.2.1/cash_v2.2.1.zip/download
下载cash v2.2.1版本，已经是5年前的版本了，也是最新版。

unzip cash_v2.2.1.zip
cd cash_v2.2.1
java -jar -Xmx10g cash.jar --help

用法如下：

java -jar -Xmx10g /share/home/chaim/soft/cash/cash_v2.2.1/cash.jar --Case:Mutant Mutant_1.bam,Mutant_2.bam,Mutant_3.bam --Control:WT Normal_1.bam,Normal_2.bam,Normal_3.bam --GTF genome.gtf --Output samples

遇到

结果解析：
samples.ControlvsTreat.alldiff.statistics.txt 主要是统计分析结果
samples.MutationvsWildType.alldiff.txt 具体剪切信息文件
提取显著的外显子保留的可变剪切
awk '$11=="IR" && $10<0.05' samples.MutantvsWT.alldiff.txt >IR.list.cash
结果只有9个基因。

rmats鉴定可变剪切位点

安装不太方便
直接安装rmats V4.1.2,请参考https://blog.csdn.net/yaya_bioinfo/article/details/129047948
一堆依赖，一般人装不了。所以直接用conda.

conda create -n rmats
conda activate rmats
conda install -c conda-forge -c bioconda rmats=4.1.2

最新版是V4.1.2，所以指定版本是V4.1.2,如果不指定的话，默认是V4.0.2.

使用rmats分析可变剪切

#使用rmats分析可变剪切
#激活conda环境(shell脚本里，激活conda比较麻烦，需要先source)
condapath=`conda info | grep 'base environment'|cut -d : -f2|cut -d " " -f2`
source ${condapath}/etc/profile.d/conda.sh
conda deactivate
conda activate rmats
rmats.py --b1 WT.txt --b2 Mutant.txt --gtf genome.gtf --od AS --tmp tmp -t paired --readLength 150 --cstat 0.0001 --nthread 24

参数说明：
--b1 WT.txt文件内是逗号分割的 对照材料的bam文件名
--b2 Mutant.txt文件内是逗号分割的 处理材料的bam文件名
--gtf 基因组gtf文件
--od 输出文件夹
--tmp 缓存文件路径
-t 序列类型：单端single还是双端paired
--readLength 序列长度，二代是150bp,三代得看报告
--cstat cutoff的阈值，默认就是0.0001，即0.1%
--nthread 进程数量

rmats的输出结果比较复杂，每个文件的具体讲解可以看https://www.jianshu.com/p/d09b95a98c64
提取显著的RI类型的结果
awk '$20<0.05' AS/RI.MATS.JC.txt >RI.list.rmats
rmats输出的结果里面不是一个基因一行，而是一个可变剪切一行。所以行数会比较多。cash的输出结果是一行一个基因。但是rmats的输出结果还是比cash的多太多，rmats去除重复后鉴定到1950个RI的基因。

可变剪切的分类

image.png

SE 外显子跳跃Skipped exon
A5SS 可变5'截切位点
A3SS 可变3'截切位点
MXE mutually exclusive exons 同源互斥外显子
RI 外显子保留