当前位置： 首页>前端>正文

转录组测序选择性剪切分析与可视化

前端2024-04-06 13:25:28

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
转录组测序选择性剪切分析与可视化,第1张

查看版本


转录组测序选择性剪切分析与可视化,第2张

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}' | awk -F '_' '{print}'| 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}' | awk -F '_' '{print}'| 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
转录组测序选择性剪切分析与可视化,第3张

5.2.3查看结果

cd ~/miniconda3/envs/rmats/rMATS/path/to/output
ll -h

总共有如下文件：


转录组测序选择性剪切分析与可视化,第4张

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

可视化结果如下：


转录组测序选择性剪切分析与可视化,第5张

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
转录组测序选择性剪切分析与可视化,第6张

可以发现分组显示和单个显示还是有些区别。
以上就是完整的基于rMATS进行转录组选择性剪切分析和可视化的流程了。


https://www.xamrdz.com/web/2ae1893302.html

相关文章：