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

• 更詳細的說明可以參考https://github.com/comprna/SUPPA/wiki/SUPPA2-tutorial#differential-splicing-with-local-events
• 輸入輸出的格式參考https://github.com/comprna/SUPPA#output-files-1