SUPPA 可变剪切分析

  • SUPPA是一款通过转录本定量来获取可变剪切定量结果的软件。转录本的定量方式有很多，例如count，FPKM, TPM等，作者建议使用TPM，因为先均一化了基因的长度，然后均一化了测序的深度。同时建议使用salmon软件进行定量
  • 软件的下载与安装

首先下载salmon二进制版本

wget https://github.com/COMBINE-lab/salmon/releases/download/v0.14.0/salmon-0.14.0_linux_x86_64.tar.gz
tar zxvf salmon-0.14.0_linux_x86_64.tar.gz -C ../

下面使用小鼠的转录本来建立索引

wget ftp://ftp.ensembl.org/pub/release-96/fasta/mus_musculus/cds/Mus_musculus.GRCm38.cds.all.fa.gz
gunzip Mus_musculus.GRCm38.cds.all.fa.gz
perl -lane 'if(/^>/){$id=(split/\./,$_)[0];print $id}else{print}' Mus_musculus.GRCm38.cds.all.fa >Mus_musculus.GRCm38.cds.all.format.fa

使用salmon建立索引

mkdir Mus_musculus.GRCm38.cds.index
export LD_LIBRARY_PATH=/media/sdb/user/yueyao/software/salmon-latest_linux_x86_64/lib:$LD_LIBRARY_PATH
/media/sdb/user/yueyao/software/salmon-latest_linux_x86_64/bin/salmon index -t Mus_musculus.GRCm38.cds.all.format.fa -i Mus_musculus.GRCm38.cds.index

SUPPA2的安装，需要注意的是SUPPA2是使用python3.4编写的，所以不要使用python2来进行安装了

#使用pip安装
pip install SUPPA==2.2.1

#使用conda进行安装
conda install -c bioconda suppa

#从github下载
wget https://github.com/comprna/SUPPA/archive/master.zip -O suppa2.zip

salmon进行定量

/media/sdb/user/yueyao/software/salmon-latest_linux_x86_64/bin/salmon quant -i Mus_musculus.GRCm38.cds.index -l ISF --gcBias -1 /media/sdb/user/yueyao/CircosRNA/00.data/107190A/107190A_1.fq.gz -2 /media/sdb/user/yueyao/CircosRNA/00.data/107190A/107190A_2.fq.gz -p 12 -o 107190A

-l 参数表示

提取salmon结果中的TPM值

python /media/sdb/user/yueyao/software/SUPPA-master/multipleFieldSelection.py -i 107190A/quant.sf -k 1 -f 4 -o 107190A_iso_tpm.txt

-k 表示转录本id在第一列
-f 表示TPM值在第四列
使用suppa进行生成一个ioe文件，需要注意的是gtf注释文件是进行一些简单的处理，这里将小鼠的ensemble下载的gtf处理成如下格式

chr14 Ensembl exon  73741918  73744001  0.0 - . gene_id "ENSG00000000001"; transcript_id "ENST00000000001.1"; 
chr14 Ensembl exon  73749067  73749213  0.0 - . gene_id "ENSG00000000001"; transcript_id "ENST00000000001.1";  
chr14 Ensembl exon  73750789  73751082  0.0 - . gene_id "ENSG00000000001"; transcript_id "ENST00000000001.1"; 
chr14 Ensembl exon  73753818  73754022  0.0 - . gene_id "ENSG00000000001"; transcript_id "ENST00000000001.1";

使用如下命令

perl -F"\t" -lane 'if(/^[1-9]+/){$F[8]=~/(gene_id\s"ENSMUSG\d+";)\s/;$gene=$1;$F[8]=~/(transcript_id\s"ENSMUST\d+";)\s/;$trans=$1;if($F[2] eq "exon"){print join("\t",@F[0..7])."\t".qq{$gene $trans}}}' Mus_musculus.GRCm38.96.gtf  >Mus_musculus.GRCm38.96.format.gtf

使用generateEvents命令根据基因组的gtf注释文件生成所有的可变剪切事件，格式保存为ioe格式

python /media/sdb/user/yueyao/software/miniconda3/envs/suppa2/bin/suppa.py generateEvents -i Mus_musculus.GRCm38.96.format.gtf -o 107190A.events -e SE SS MX RI FL -f ioe

-i 输入的gtf文件
-o 输出的文件前缀
-e 输出可变剪切的类型
-f 设置输出格式
将不同的可变剪切事件合并成一个结果

awk '
    FNR==1 && NR!=1 { while (/^<header>/) getline; }
    1 {print}
' *.ioe > ensembl_mm10.events.ioe

需要注意的时如果使用的是RefSeq annotation注释文件，要使用--pooled-genes参数，因为RefSeq基因是根据它们所包含的mRNA序列来识别，而不是基因位置信息，这可能导致同一个基因会比对到基因组两个不同的地方去，或者是两个同源异构体具有相同的exon或者剪切位点被标记为不同的基因。pooled-genes重新定义了这种情况

得到一个PSI值，需要注意的是使用的转录本ID和gtf的转录本ID应该是一致，数目不一样可能会有错误提示:

python /media/sdb/user/yueyao/software/miniconda3/envs/suppa2/bin/suppa.py psiPerEvent -i ensembl_mm10.events.ioe -e 107190A_iso_tpm.txt -o TRA2_events

提取某一个基因的可变剪切事件进行不同分组的作图

python /media/sdb/user/yueyao/software/SUPPA-master/scripts/generate_boxplot_event.py -e "ENSMUSG00000000244;MX:7:143007005-143011034:143011108-143015581:143007005-143014959:143015060-143015581:+" -i TRA2_events.psi.psi -g 1-2,3-4,5-6 -c A1,A2,A3 -o /media/sdb/user/yueyao/Test/suppa2/result/

-i 输入PSI矩阵
-e 某一个基因的某种类型的可变剪切事件
-g 设置分组的样品
-c 设置组名，与前面的分组对应

根据PSI文件，将可变剪切的结果和表达量的结果按照分组分成两个文件
需要注意的是样品的名称命名不能以数字开头，因为R脚本默认会对读入的title的数字前面加一个X

#注意样品的命名，可能会影响这一步
# Split the PSI and TPM files between the 2 conditions:
Rscript /media/sdb/user/yueyao/software/SUPPA-master/scripts/split_file.R iso_tpm.txt S107190A,S107192A,S107194A S107195A,S107196A,S107197A /media/sdb/user/yueyao/Test/suppa2/Group1_iso_tmp.txt /media/sdb/user/yueyao/Test/suppa2/Group2_iso_tmp.txt
Rscript /media/sdb/user/yueyao/software/SUPPA-master/scripts/split_file.R TRA2_events.psi S107190A,S107192A,S107194A S107195A,S107196A,S107197A /media/sdb/user/yueyao/Test/suppa2/Group1.psi /media/sdb/user/yueyao/Test/suppa2/Group2.psi

计算两个不同分组的可变剪切差异

python /media/sdb/user/yueyao/software/SUPPA-master/suppa.py diffSplice -m empirical -gc -i ensembl_mm10.events.ioe -p Group1.psi Group2.psi -e Group1_iso_tmp.txt Group2_iso_tmp.txt -o result/Group1_vs_Group2

根据差异可变剪切的结果，可以选区某一个类型的可变剪切事件进行聚类分析,我这里用了所有的事件，好像结果不正常

python /media/sdb/user/yueyao/software/SUPPA-master/suppa.py clusterEvents --dpsi Group1_vs_Group2.dpsi.temp.0 --psivec Group1_vs_Group2.psivec --sig-threshold 0.05 --eps 0.2 --separation 0.11 -dt 0.2 --min-pts 10 --groups 1-3,4-6 -c OPTICS -o cluster

 

posted @ 2019-07-01 14:09  raisok  阅读(4362)  评论(0)  编辑  收藏  举报