轉(zhuǎn)錄組測序選擇性剪切分析與可視化

1. 創(chuàng)建小環(huán)境

rMATS是常用的選擇性剪切（alternative splicing）分析軟件。目前最新版是rMATS 4.1.2，已經(jīng)可以通過conda安裝，步驟如下：

1.1 創(chuàng)建rmats小環(huán)境
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. 數(shù)據(jù)質(zhì)控

2.1 質(zhì)控軟件安裝

這里主要安裝質(zhì)控軟件

conda install -y fastqc trim-galore multiqc
2.2 fastqc質(zhì)控

使用fastqc對原始數(shù)據(jù)進行質(zhì)量評估

# 激活conda環(huán)境
conda activate rmats

# 上傳數(shù)據(jù)到自己的文件夾~/project/AS/data/rawdata
cd ~/project/AS/data

# 使用FastQC軟件對單個fastq文件進行質(zhì)量評估，結(jié)果輸出到qc/文件夾下
# 目錄改成自己的目錄，否則會報錯：permission deny.
qcdir=~/project/AS/data/qc
fqdir=~/project/AS/data/rawdata

# 多個數(shù)據(jù)質(zhì)控
fastqc -t 20 -o $qcdir $fqdir/*.fq.gz

# 使用MultiQc整合FastQC結(jié)果
cd ~/project/AS/data/qc
multiqc *.zip

3. 數(shù)據(jù)過濾

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的內(nèi)容
# 目錄改成自己的目錄，否則會報錯：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

# 提交任務(wù)到后臺
nohup sh trim_galore.sh >trim_galore.log &

4. 數(shù)據(jù)比對

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 數(shù)據(jù)比對
4.2.1 建立索引

先進入?yún)⒖蓟蚰夸?/p> 

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的內(nèi)容
# 目錄改成自己的目錄，否則會報錯：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程序線程數(shù)不能設(shè)置太大，不然會報錯。

# 提交任務(wù)到后臺
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文件，','連接生物重復(fù)。

./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 #創(chuàng)建輸出目錄
cd ~/miniconda3/envs/rmats/rMATS/path/to/
mkdir tmp_output #創(chuàng)建臨時文件輸出目錄
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

運行結(jié)果都存放在~/miniconda3/envs/rmats/rMATS/path/to/output目錄里面。

5.2 從bam文件開始

5.2.1 文件準備

首先將參考基因組、STAR比對后的sorted BAM文件復(fù)制過來。

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查看結(jié)果

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

總共有如下文件：


image.png

6. 數(shù)據(jù)可視化

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 \ #標記可視化結(jié)果中每個組樣本名字前綴
-o test_coordinate_output #輸出目錄

#復(fù)制代碼時，需要將注釋信息刪除

如果一開始是以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

查看結(jié)果

cd ~/miniconda3/envs/rmats/rMATS/path/to/test_coordinate_output/Sashimi_plot

可視化結(jié)果如下：


image.png

6.2.2 分組顯示可視化

cd ~/miniconda3/envs/rmats/rMATS/path/to
mkdir test_coordinate_output_group
#先制作個分組文件grouping.gf，內(nèi)容如下
#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

查看結(jié)果

cd ~/miniconda3/envs/rmats/rMATS/path/to/test_coordinate_output_group/Sashimi_plot
image.png

可以發(fā)現(xiàn)分組顯示和單個顯示還是有些區(qū)別。
以上就是完整的基于rMATS進行轉(zhuǎn)錄組選擇性剪切分析和可視化的流程了。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容