1. 创建小环境
rMATS是常用的选择性剪切（alternative splicing）分析软件。目前最新版是rMATS 4.1.2，已经可以通过conda安装，步骤如下：
1.1 创建rmats小环境
conda create -n rmats -y
1.2 安装rMATS 4.1.2
conda install -c bioconda rmats
进入到rMATS软件目录
cd ~/miniconda3/envs/rmats/rMATS
1.3 测试rMATS
image.png
查看版本
image.png
2. 数据质控
2.1 质控软件安装
这里主要安装质控软件
conda install -y fastqc trim-galore multiqc
2.2 fastqc质控
使用fastqc对原始数据进行质量评估
# 激活conda环境
conda activate rmats
# 上传数据到自己的文件夹~/project/AS/data/rawdata
cd ~/project/AS/data
# 使用FastQC软件对单个fastq文件进行质量评估，结果输出到qc/文件夹下
# 目录改成自己的目录，否则会报错：permission deny.
qcdir=~/project/AS/data/qc
fqdir=~/project/AS/data/rawdata
# 多个数据质控
fastqc -t 20 -o $qcdir $fqdir/*.fq.gz
# 使用MultiQc整合FastQC结果
cd ~/project/AS/data/qc
multiqc *.zip
3. 数据过滤
3.1 trim_galore过滤
# 新建文件夹
cd ~/project/AS/data
mkdir -p cleandata/trim_galore
#多个样本
ls ../../rawdata/*.gz | awk -F'/' '{print$4}' | awk -F '_' '{print$1}'| uniq >sample.ID.txt
# 多个样本 vim trim_galore.sh，以下为sh的内容
# 目录改成自己的目录，否则会报错：permission deny.
rawdata=~/project/AS/data/rawdata
cleandata=~/project/AS/data/cleandata/trim_galore
cat ~/project/AS/data/cleandata/trim_galore/sample.ID.txt | while read id
do
trim_galore --phred33 -q 20 --length 36 --stringency 3 --fastqc --paired --max_n 3 -o ${cleandata} ${rawdata}/${id}_1.fq.gz ${rawdata}/${id}_2.fq.gz
done
# 提交任务到后台
nohup sh trim_galore.sh >trim_galore.log &
4. 数据比对
4.1 参考基因组
## 参考基因组准备:注意参考基因组版本信息
# 下载，Ensembl：http://asia.ensembl.org/index.html
# ftp://ftp.ensembl.org/pub/release-105/fasta/Mus_musculus/dna/
# 进入到参考基因组目录
cd ~/database/genome/Ensembl/Mus_musculus/GRCm39_release105
# 下载基因组序列
wget -c http://ftp.ensembl.org/pub/release-105/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz
# 下载基因组注释文件
wget -c ftp://ftp.ensembl.org/pub/release-105/gtf/Mus_musculus/Mus_musculus.GRCm39.105.gtf.gz
4.2 数据比对
4.2.1 建立索引
先进入参考基因目录
cd ~/database/genome/Ensembl/Mus_musculus/GRCm39_release105/
mkdir -p index/STAR
再建立索引
STAR \
--runMode genomeGenerate \
--runThreadN 40 \
--genomeDir ~/database/genome/Ensembl/Mus_musculus/GRCm39_release105/index/STAR \
--genomeFastaFiles ~/database/genome/Ensembl/Mus_musculus/GRCm39_release105/Mus_musculus.GRCm39.dna.primary_assembly.fa \
--sjdbGTFfile ~/database/genome/Ensembl/Mus_musculus/GRCm39_release105/Mus_musculus.GRCm39.105.gtf \
--sjdbOverhang 149
# --sjdbOverhang 这个值为你测序read的长度减1，是在注释可变剪切序列的时候使用的最大长度值
如果从fq文件开始运行rMATS，只需要索引文件，不需要得到比对后的bam文件。
4.2.2 STAR比对
cd ~/project/AS/Mapping/STAR
ls ../../data/cleandata/trim_galore/*.gz | awk -F'/' '{print$6}' | awk -F '_' '{print$1}'| uniq >sample.ID.txt
# 多个样本 vim STAR.sh，以下为sh的内容
# 目录改成自己的目录，否则会报错：permission deny.
index=~/database/genome/Ensembl/Mus_musculus/GRCm39_release105/index/STAR
cleandata=~/project/AS/data/cleandata/trim_galore
cat ~/project/AS/Mapping/STAR/sample.ID.txt | while read id
do
STAR --runThreadN 8 --genomeDir ${index} \
--readFilesIn ${cleandata}/${id}_1_val_1.fq.gz ${cleandata}/${id}_2_val_2.fq.gz \
--readFilesCommand zcat \
--outSAMtype BAM SortedByCoordinate \
--outFileNamePrefix ./${id} \
--outBAMsortingThreadN 8
done
## STAR程序线程数不能设置太大，不然会报错。
# 提交任务到后台
nohup sh STAR.sh >STAR.log &
5. 选择性剪切分析
5.1 从fastq文件开始
5.1.1 文件准备
首先将参考基因组、索引、过滤后fq文件软链接过来。
cd ~/miniconda3/envs/rmats/rMATS/path/to
ln -s ~/database/genome/Ensembl/Mus_musculus/GRCm39_release105/Mus_musculus.GRCm39.105.gtf ./
ln -s ~/database/genome/Ensembl/Mus_musculus/GRCm39_release105/index/STAR ./
ln -s ~/project/AS/Mapping/STAR/*.fq.gz ./
然后，制作b1.txt，存放FAD组文件目录。':'连接R1和R2文件，','连接生物重复。
./path/to/FAD_1.R1.fq.gz:./path/to/FAD_1.R2.fq.gz,./path/to/FAD_2.R1.fq.gz:./path/to/FAD_2.R1.fq.gz,./path/to/FAD_5.R1.fq.gz:./path/to/FAD_5.R1.fq.gz
其次，制作b2.txt，存放WT组文件目录
./path/to/WT_1.R1.fq.gz:./path/to/WT_1.R2.fq.gz,./path/to/WT_2.R1.fq.gz:./path/to/WT_2.R1.fq.gz,./path/to/WT_5.R1.fq.gz:./path/to/WT_5.R1.fq.gz
5.1.2 rMATS运行
cd ~/miniconda3/envs/rmats/rMATS/
mkdir -p path/to/output #创建输出目录
cd ~/miniconda3/envs/rmats/rMATS/path/to/
mkdir tmp_output #创建临时文件输出目录
cd ../../ #退回到rMATS目录
python rmats.py --s1 ./path/to/s1.txt --s2 ./path/to/s2.txt --gtf ./path/to/Mus_musculus.GRCm39.105.gtf --bi ./path/to/STAR -t paired --readLength 50 --nthread 8 --od ./path/to/output --tmp ./path/to/tmp_output
运行结果都存放在~/miniconda3/envs/rmats/rMATS/path/to/output目录里面。
5.2 从bam文件开始
5.2.1 文件准备
首先将参考基因组、STAR比对后的sorted BAM文件复制过来。
cd ~/miniconda3/envs/rmats/rMATS/path/to
##把一些文件软连接过来
cp ~/database/genome/Ensembl/Mus_musculus/GRCm39_release105/Mus_musculus.GRCm39.105.gtf ./
cp ~/project/AS/Mapping/STAR/*.out.bam ./
#从BAM开始不再需要STAR索引文件
然后，制作b1.txt，存放FAD组文件目录。
./path/to/FAD1Aligned.sortedByCoord.out.bam,./path/to/FAD2Aligned.sortedByCoord.out.bam,./path/to/FAD5Aligned.sortedByCoord.out.bam
其次，制作b2.txt，存放WT组文件目录。
./path/to/WT1Aligned.sortedByCoord.out.bam,./path/to/WT2Aligned.sortedByCoord.out.bam,./path/to/WT5Aligned.sortedByCoord.out.bam
5.2.2 rMATS运行
cd ~/miniconda3/envs/rmats/rMATS/
python rmats.py --b1 ./path/to/b1.txt --b2 ./path/to/b2.txt --gtf ./path/to/Mus_musculus.GRCm39.105.gtf -t paired --readLength 150 --nthread 8 --od ./path/to/output --tmp /path/to/tmp_output
image.png
5.2.3查看结果
cd ~/miniconda3/envs/rmats/rMATS/path/to/output
ll -h
总共有如下文件：
image.png
6. 数据可视化
6.1 rmats2sashimiplot安装
#还是采用conda安装
conda activate rmats
conda install -c bioconda rmats2sashimiplot
6.2 可视化
6.2.1 所有样本可视化
cd ~/miniconda3/envs/rmats/rMATS/path/to
mkdir test_coordinate_output
##这里用的是5.2步骤中的BAM文件。5.1步骤中的BAM文件在~/miniconda3/envs/rmats/rMATS/path/to/tmp_output里面子文件夹中。
rmats2sashimiplot \
--b1 ./FAD1Aligned.sortedByCoord.out.bam,./FAD2Aligned.sortedByCoord.out.bam,./FAD5Aligned.sortedByCoord.out.bam \ #第一组所有BAM文件
--b2 ./WT1Aligned.sortedByCoord.out.bam,./WT2Aligned.sortedByCoord.out.bam,./WT5Aligned.sortedByCoord.out.bam \ #第二组所有BAM文件
-t SE -e ./output/SE.MATS.JC.txt \ #这里选择SE event来可视化
--l1 FAD --l2 WT \ #标记可视化结果中每个组样本名字前缀
-o test_coordinate_output #输出目录
#复制代码时，需要将注释信息删除
如果一开始是以fastq文件运行rMATS的话，使用以下代码
cd ~/miniconda3/envs/rmats/rMATS/path/to
mkdir test_coordinate_output
rmats2sashimiplot \
--b1 ./tmp_output/2022-03-27-23_01_45_176979_bam1_1/Aligned.sortedByCoord.out.bam,./tmp_output/2022-03-27-23_01_45_176979_bam1_2/Aligned.sortedByCoord.out.bam,./tmp_output/2022-03-27-23_01_45_176979_bam1_3/Aligned.sortedByCoord.out.bam \
--b2 ./tmp_output/2022-03-27-23_01_45_176979_bam2_1/Aligned.sortedByCoord.out.bam,./tmp_output/2022-03-27-23_01_45_176979_bam2_2/Aligned.sortedByCoord.out.bam,./tmp_output/2022-03-27-23_01_45_176979_bam2_3/Aligned.sortedByCoord.out.bam \
-t SE -e ./output/SE.MATS.JC.txt \
--l1 FAD --l2 WT \
-o test_coordinate_output
查看结果
cd ~/miniconda3/envs/rmats/rMATS/path/to/test_coordinate_output/Sashimi_plot
可视化结果如下：
image.png
6.2.2 分组显示可视化
cd ~/miniconda3/envs/rmats/rMATS/path/to
mkdir test_coordinate_output_group
#先制作个分组文件grouping.gf，内容如下
#FAD: 1-3
#WT: 4-6
rmats2sashimiplot \
--b1 ./FAD1Aligned.sortedByCoord.out.bam,./FAD2Aligned.sortedByCoord.out.bam,./FAD5Aligned.sortedByCoord.out.bam \
--b2 ./WT1Aligned.sortedByCoord.out.bam,./WT2Aligned.sortedByCoord.out.bam,./WT5Aligned.sortedByCoord.out.bam \
-t SE -e ./output/SE.MATS.JC.txt \
--l1 FAD --l2 WT \
--group-info grouping.gf \
-o test_coordinate_output_group
查看结果
cd ~/miniconda3/envs/rmats/rMATS/path/to/test_coordinate_output_group/Sashimi_plot
image.png
可以发现分组显示和单个显示还是有些区别。
以上就是完整的基于rMATS进行转录组选择性剪切分析和可视化的流程了。
