使用SUPPA进行可变剪切分析

前面我们的明码标价之普通转录组上游分析，受到了各大热心粉丝的吐槽，觉得太简单了我们居然还好意思收费。

额，其实呢，这些粉丝应该是“饱汉不知饿汉饥”，其实数据分析这个技能啊，难者不会，会者不难。而且普通转录组上游分析仅仅是一个引子，又不是说我们只能提供这些粗浅的分析，转录组相关的绝大部分高级分析我们都是有现成的代码，而且都分享过一些教程：

用LeafCutter探索转录组数据的可变剪切，需要star软件的
使用SGSeq探索可变剪切，一个读取bam文件的R包

最近看到了一个比较新，而且发表的杂志也还不错的工具：SUPPA2

首先认识 SUPPA 软件，它有5个子命令，分别如下:

generateEvents: Generates events from an annotation.
psiPerEvent : Quantifies event inclusion levels (PSIs) from multiple samples.
psiPerIsoform : Quantifies isoform inclusion levels (PSIs) from multiple samples.
diffSplice : Calculate differential splicing across multiple conditions with replicates.
clusterEvents : Cluster events according to PSI values across conditions.

可以根据已知的 gtf文件信息针对多种可变剪切形式进行差异分析，它还提供了一些独立的py脚本文件，如下所示:

准备工作

包括conda安装，以及SUPPA2环境安装，参考转录组fasta序列获取和构建索引文件，参考gtf文件处理。

统一使用gencode最新版： https://www.gencodegenes.org/human/

mkdir -p SUPPA2/ref
cd SUPPA2/ref

nohup wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_37/gencode.v37.annotation.gtf.gz &
nohup wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_37/gencode.v37.transcripts.fa.gz &
gunzip gencode.v37.annotation.gtf.gz
gunzip gencode.v37.transcripts.fa.gz
# 1.4G 2月 17 16:50 gencode.v37.annotation.gtf

使用conda创建rna小环境：

代码如下：

conda create -n rna
conda activate rna 
conda install -y trim-galore
# python-3.7.9 | 57.3 MB, openjdk-10.0.2 | 189.2 MB
conda install -y star hisat2 bowtie2 # 没有其它依赖
conda install -y subread bedtools deeptools
conda install -y salmon=1.4.0
conda install -y -c bioconda suppa

下载参考转录组fasta序列文件然后使用salmon构建index文件

conda activate rna 
nohup salmon index -t gencode.v37.transcripts.fa -i gencode.v37.transcripts.salmon.index &

得到的salmon的index是一个文件夹，内容如下:

$ ls -lh gencode.v37.transcripts.salmon.index/|cut -d" " -f 5-

913K 2月 17 17:19 complete_ref_lens.bin
 95M 2月 17 17:23 ctable.bin
4.3M 2月 17 17:23 ctg_offsets.bin
175K 2月 17 17:19 duplicate_clusters.tsv
1.1K 2月 17 17:24 info.json
 90M 2月 17 17:24 mphf.bin
479M 2月 17 17:24 pos.bin
 496 2月 17 17:24 pre_indexing.log
 23M 2月 17 17:23 rank.bin
1.8M 2月 17 17:23 refAccumLengths.bin
 16K 2月 17 17:24 ref_indexing.log
913K 2月 17 17:23 reflengths.bin
 92M 2月 17 17:23 refseq.bin
 45M 2月 17 17:23 seq.bin
 126 2月 17 17:24 versionInfo.json

下载基因注释gtf文件然后使用suppa构建可变剪切事件综合文件

统一使用gencode最新版： https://www.gencodegenes.org/human/

对gtf文件的全部转录本走 suppa.py generateEvents 命令：

# 使用 generateEvents 命令根据基因组的gtf注释文件生成所有的可变剪切事件，格式保存为ioe格式
# 针对 SE SS MX RI FL 可变剪切批量分析的
suppa.py generateEvents \
-i gencode.v37.annotation.gtf \
-o gencode.v37.events -e SE SS MX RI FL -f ioe

得到的文件如下:

$ wc -l *ioe
 20838 gencode.v37.events_A3_strict.ioe
 18921 gencode.v37.events_A5_strict.ioe
 112101 gencode.v37.events_AF_strict.ioe
 37317 gencode.v37.events_AL_strict.ioe
 8629 gencode.v37.events_MX_strict.ioe
 8474 gencode.v37.events_RI_strict.ioe
 50562 gencode.v37.events_SE_strict.ioe

可以看到是，基于 gencode.v37.annotation.gtf 文件的7种全部的可能的可变剪切事件，都被列出来了。

# 合并所有的ioe文件
awk 'FNR==1 && NR!=1 { while (/^<header>/) getline; } 1 {print}' *.ioe > gencode.v37.all.events.ioe
wc -l gencode.v37.all.events.ioe
# 256842 gencode.v37.all.events.ioe

上游分析拿到tpm矩阵

测序数据获取

得到文件如下:

ls -lh rawdata/*gz|cut -d" " -f 5-
1.7G 2月 8 14:33 rawdata/SRR1513329_1.fastq.gz
1.7G 2月 8 14:33 rawdata/SRR1513329_2.fastq.gz
1.7G 2月 8 14:33 rawdata/SRR1513330_1.fastq.gz
1.7G 2月 8 14:33 rawdata/SRR1513330_2.fastq.gz
1.8G 2月 8 14:33 rawdata/SRR1513331_1.fastq.gz
1.8G 2月 8 14:33 rawdata/SRR1513331_2.fastq.gz
1.8G 2月 8 14:33 rawdata/SRR1513332_1.fastq.gz
1.8G 2月 8 14:33 rawdata/SRR1513332_2.fastq.gz
1.9G 2月 8 14:33 rawdata/SRR1513333_1.fastq.gz
2.0G 2月 8 14:33 rawdata/SRR1513333_2.fastq.gz
1.7G 2月 8 14:33 rawdata/SRR1513334_1.fastq.gz
1.8G 2月 8 14:33 rawdata/SRR1513334_2.fastq.gz

然后质控，走fastqc软件，查看报告

fastqc -t 2 *.fastq.gz -o ../qc

通常情况下大家需要针对clean后的fq文件进行后续的salmon流程定量获得各个样品的TPM值。

走trim流程

因为是双端测序，所以:

cd rawData
 ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o ../cleanData ${id}*.gz &
done

得到的clean的fq文件如下:

 $ ls -lh *gz|cut -d" " -f 5-
1.6G 2月 18 10:00 SRR1513329_1_val_1.fq.gz
1.6G 2月 18 10:00 SRR1513329_2_val_2.fq.gz
1.6G 2月 18 10:00 SRR1513330_1_val_1.fq.gz
1.6G 2月 18 10:00 SRR1513330_2_val_2.fq.gz
1.7G 2月 18 10:02 SRR1513331_1_val_1.fq.gz
1.8G 2月 18 10:02 SRR1513331_2_val_2.fq.gz
1.7G 2月 18 10:02 SRR1513332_1_val_1.fq.gz
1.7G 2月 18 10:02 SRR1513332_2_val_2.fq.gz
1.9G 2月 18 11:15 SRR1513333_1_val_1.fq.gz
1.9G 2月 18 11:15 SRR1513333_2_val_2.fq.gz
1.7G 2月 18 18:04 SRR1513334_1_val_1.fq.gz
1.7G 2月 18 18:04 SRR1513334_2_val_2.fq.gz

项目介绍分组情况如下:

negative control siRNA ：SRR1513329,SRR1513330,SRR1513331 
TRA2A/B siRNA ：SRR1513332,SRR1513333,SRR1513334

salmon流程定量获得各个样品的TPM值

## 定量获得TPM值 
mkdir salmon_output
cd cleanData 
index=../ref/gencode.v37.transcripts.salmon.index/
 ls *gz|cut -d"_" -f 1|sort -u |while read id;do
 nohup salmon quant -i $index -l ISF --gcBias \
 -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -p 2 \
 -o ../salmon_output/${id}_output 1>${id}_salmon.log 2>&1 &
 done
 ## quant.sf文件很重要，要用于后续的分析
##ENST和ENSG的前三个字母（ENS），意思是“ENSENMBLE”。T是指转录本。G是指基因。当然，也会看到ENSP，P自然是指蛋白质了。

几分钟就可以走完全部的salmon流程，每一个样品虽然说会输出一个文件夹，但是最重要的都是各自的 quant.sf文件，行数都是统一的:

$ wc -l ../salmon_output/*/quant.sf
 233653 ../salmon_output/SRR1513329_output/quant.sf
 233653 ../salmon_output/SRR1513330_output/quant.sf
 233653 ../salmon_output/SRR1513331_output/quant.sf
 233653 ../salmon_output/SRR1513332_output/quant.sf
 233653 ../salmon_output/SRR1513333_output/quant.sf
 233653 ../salmon_output/SRR1513334_output/quant.sf

各自的 quant.sf文件内容如下:

==> ../salmon_output/SRR1513334_output/quant.sf <==
Name Length EffectiveLength TPM NumReads
ENST00000456328.2|ENSG00000223972.5|OTTHUMG00000000961.2|OTTHUMT00000362751.1|DDX11L1-202|DDX11L1|1657|processed_transcript| 1657 1587.157 0.000000 0.000
ENST00000450305.2|ENSG00000223972.5|OTTHUMG00000000961.2|OTTHUMT00000002844.2|DDX11L1-201|DDX11L1|632|transcribed_unprocessed_pseudogene| 632 452.000 0.000000 0.000
ENST00000488147.1|ENSG00000227232.5|OTTHUMG00000000958.1|OTTHUMT00000002839.1|WASH7P-201|WASH7P|1351|unprocessed_pseudogene| 1351 1273.856 3.098744 27.492

tpm表达量矩阵的获取

这个时候，可以直接使用作者的 multipleFieldSelection.py 脚本，对我们的全部样品的 quant.sf文件进行信息提取即可。其实自己写一个脚本也是很简单的：

-k indicates the row used as index
-f TPM值在第四列，在这里要提取每个样本的TPM值

## 提取所有样本的TPM值并合并为一个文件
 cd ../salmon_output/
multipleFieldSelection.py \
-i SRR*/quant.sf -k 1 -f 4 \
-o iso_tpm.txt

得到 iso_tpm.txt的内容，还需要格式化一下，因为它跟gtf文件里面的基因信息不一致。

# gencode.v37.transcripts.fa 内容如下
>ENST00000450305.2|ENSG00000223972.5|OTTHUMG00000000961.2|OTTHUMT00000002844.2|DDX11L1-201|DDX11L1|632|transcribed_unprocessed_pseudogene|
# gencode.v37.annotation.gtf 内容如下
chr1 HAVANA gene 11869 14409 . + . gene_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; level 2; hgnc_id "HGNC:37102"; havana_gene "OTTHUMG00000000961.2";

各自的 quant.sf文件里面的ID是gencode.v37.transcripts.fa里面的转录本ID，有点过于复杂了，所以需要进行简化，需要与 gencode.v37.annotation.gtf 里面的转录本ID保持一致。

 perl -alne '{/(\|.*\|)\t/; ;s/$1//g;s/\|//g;print}' iso_tpm.txt > iso_tpm_formatted.txt

输出文件是 iso_tpm_formatted.txt ，是 233653 行的普通文本，可以看到简洁很多：

SRR1513329_output SRR1513330_output SRR1513331_output SRR1513332_output SRR1513333_output SRR1513334_output
ENST00000456328.2 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
ENST00000450305.2 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000

拿到了这样的以转录本ID为行的tpm矩阵，再加上 gencode.v37.annotation.gtf 的可变剪切形式文件，就可以进行下游分析获取统计学显著的可变剪切事件。

下游分析获取统计学显著的可变剪切事件

项目介绍分组情况如下:

negative control siRNA ：SRR1513329,SRR1513330,SRR1513331 
TRA2A/B siRNA ：SRR1513332,SRR1513333,SRR1513334

对tpm表达量矩阵，根据可变剪切事件计算psi值矩阵：

ioe_merge_file=../ref/gencode.v37.all.events.ioe
ls $ioe_merge_file
suppa.py psiPerEvent \
-i $ioe_merge_file -e iso_tpm_formatted.txt -o project_events 1>psiPerEvent_log.txt 2>&1 &
# 会输出 project_events.psi 文件

同时可以检查看看 psiPerEvent_log.txt 文件。

project_events.psi 文件，内容如下所示：

SRR1513329_output SRR1513330_output SRR1513331_output SRR1513332_output SRR1513333_output SRR1513334_output
ENSG00000000003.15;A5:chrX:100635746-100636191:100635746-100636608:- 0.20449107275352194 0.20506877462332357 0.18431225084229827 0.32174184236220055 0.27519806289263166 0.318119876931685
ENSG00000000003.15;A5:chrX:100635746-100636608:100635746-100636793:- 1.0 1.0 0.9009096685944438 1.0 0.9145937251391483 1.0

根据分组做差异可变剪切事件寻找

把前面的全部转录本可变剪切事件的psi文件以及 iso_tpm_formatted.txt ，根据分组信息，拆分为两个psi矩阵和两个tpm矩阵。

# negative control siRNA ：SRR1513329,SRR1513330,SRR1513331 
cut -f 1-4 project_events.psi > control.psi
cut -f 1-4 iso_tpm_formatted.txt > control.tpm
# TRA2A/B siRNA ：SRR1513332,SRR1513333,SRR1513334 
 cut -f 1,5-7 project_events.psi > treat.psi
 cut -f 1,5-7 iso_tpm_formatted.txt > treat.tpm
 # 很有意思的事情是，这个时候的表头似乎是并不是很重要

然后走 suppa.py diffSplice 命令，对两个分组的psi矩阵和tpm矩阵进行差异可变剪切事件寻找

ioe_merge_file=../ref/gencode.v37.all.events.ioe
ls $ioe_merge_file
suppa.py diffSplice \
-m empirical -gc -i $ioe_merge_file \
--save_tpm_events \
-p treat.psi control.psi \
-e treat.tpm control.tpm \
-o project_diffSplice # 前缀

差异分析生成下列3个文件：

ls -lh project_diffSplice* |cut -d" " -f 5-

7.9M 2月 18 09:16 project_diffSplice_avglogtpm.tab
 23M 2月 18 09:38 project_diffSplice.dpsi
 29M 2月 18 09:38 project_diffSplice.psivec

到此为止，统计学显著的可变剪切事件就获取完毕，后面就是各式各样的可视化啦！

当然了，一般来说，统计学显著的才有意义:

$ cat project_diffSplice.dpsi|perl -alne '{print if $F[2] <0.01}' |wc
 118 353 11552
$ cat project_diffSplice.dpsi|perl -alne '{print if $F[2] <0.05}' |wc
 750 2249 74204

查看差异比较大的可变剪切事件：

treat-control_dPSI treat-control_p-val
ENSG00000008517.18;AF:chr16:3065297:3065494-3065784:3065677:3065710-3065784:+ -1.0 0.0
ENSG00000009307.16;SE:chr1:114738072-114749821:114750207-114757925:- 0.6712142721 0.008991009
ENSG00000009307.16;SE:chr1:114739890-114749821:114750207-114757925:- 0.1488639077 0.008991009

这些差异比较大的可变剪切事件可以进行一些可视化，可以是igv里面的，举例如下:

b0d439d54469f5ef4cb9fc3f6fbdb6b

网页工具绝对是完成不了这样的命令行数据分析哦

这个是基于Linux的ngs数据的上游处理，目前没有成熟的网页工具支持这样的分析。其实呢，如果你有时间请务必学习编程基础，自由自在的探索海量的公共数据，辅助你的科研，那么：

如果你没有时间从头开始学编程，也可以委托专业的团队付费拿到同样的数据分析，这样的20个样品以内的简单分组转录组测序数据的表达矩阵获取，首先需800元，然后如果要进行高级分析，比如上面讲解的可变剪切，再需要800元，仍然是可以拿到全部的数据和代码哦！

需要自己读文献筛选合适的数据集
提供半个小时左右的一对一讲解转录组数据处理背景知识。

如果需要委托，直接在我们《生信技能树》公众号留言即可，我们会安排合适的生信工程师对接具体的项目。

生信菜鸟团

欢迎去论坛biotrainee.com留言参与讨论，或者关注同名微信公众号biotrainee