前面我们在《生信技能树》已经是多次分享了变剪切相关教程:

• rMATS这款差异可变剪切分析软件的使用体验
• 用LeafCutter探索转录组数据的可变剪切
• 用Expedition来分析单细胞转录组数据的可变剪切
• 使用SGSeq探索可变剪切
• 用DEXSeq分析可变剪切，外显子差异表达

但是部分教程因为时间关系，有点过时了，我们已经是大量招募小伙伴重新系统性整理梳理我们过去七年做的工作。也就是这个求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，走大运结识了几位优秀小伙伴！其中featureCounts和DEXseq做基于外显子定量的可变剪切就需要更新了：

准备工作

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

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

mkdir -p $HOME/rna/SUPPA2
mkdir -p leanData gtf rawData

mkdir -p $HOME/rna/SUPPA2/gtf
cd $HOME/rna/SUPPA2/gtf

nohup wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_37/gencode.v37.annotation.gtf.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序列文件然后使用hisat2构建index文件

这里我们直接使用服务器现成的hisat2的index文件 （其实就是hisat2官网下载的，无需自己构建）

$ ls -lh /home/data/server/reference/index/hisat/hg38/ |cut -d" " -f 5-

974M 12月 31 17:42 genome.1.ht2
728M 12月 31 17:41 genome.2.ht2
 15K 12月 31 17:42 genome.3.ht2
728M 12月 31 17:42 genome.4.ht2
1.3G 12月 31 17:42 genome.5.ht2
741M 12月 31 17:42 genome.6.ht2
 8 12月 31 17:42 genome.7.ht2
 8 12月 31 17:42 genome.8.ht2
1.3K 12月 31 17:42 make_hg38.sh

测序数据获取

这里就不赘述数据下载细节了，就是SUPPA软件文章里面举例的数据，走：使用ebi数据库直接下载fastq测序数据 , 这个教程来直接下载fastq文件啦。

得到文件如下:

$ ls -lh rawData/*gz|cut -d" " -f 5-
1.7G 2月 17 19:56 rawData/SRR1513329_1.fastq.gz
1.7G 2月 17 19:56 rawData/SRR1513329_2.fastq.gz
1.7G 2月 17 19:57 rawData/SRR1513330_1.fastq.gz
1.7G 2月 17 19:57 rawData/SRR1513330_2.fastq.gz
1.8G 2月 17 19:57 rawData/SRR1513331_1.fastq.gz
1.8G 2月 17 19:57 rawData/SRR1513331_2.fastq.gz
1.8G 2月 17 19:56 rawData/SRR1513332_1.fastq.gz
1.8G 2月 17 19:56 rawData/SRR1513332_2.fastq.gz
1.9G 2月 17 19:57 rawData/SRR1513333_1.fastq.gz
2.0G 2月 17 19:57 rawData/SRR1513333_2.fastq.gz
1.7G 2月 17 19:57 rawData/SRR1513334_1.fastq.gz
1.8G 2月 17 19:57 rawData/SRR1513334_2.fastq.gz

项目介绍 分组情况如下:

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

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

nohup fastqc -t 4 *.fq.gz 1> fastq.log 2>&1 &

走trim流程

因为是双端测序，所以:

cd $HOME/rna/SUPPA2/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 
# 如果是单端数据，需要去除 --paired 参数

得到的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

走hisat比对流程

cd $HOME/rna/SUPPA2/cleanData

indexPrefix=/home/data/server/reference/index/hisat/hg38/genome 
# 单端测序
ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup hisat2 -p 2 -x $indexPrefix -U ${id}*.gz -S ${id}.hisat.sam & 
done

# 双端测序
ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup hisat2 -p 2 -x $indexPrefix -1 ${id}*_1_val_1.fq.gz -2 ${id}*_2_val_2.fq.gz -S ${id}.hisat.sam & 
done

# sam文件转bam
ls *.sam|while read id ;do (nohup samtools sort -O bam -@ 2 -o $(basename ${id} ".sam").bam ${id} & );done
# 这个过程会输出大量中间文件
rm *.sam 
# 为bam文件建立索引
ls *.bam |xargs -i samtools index {}

上面的区分不同步骤，每次都是并行处理全部的样品，另外一个并行处理的策略 这里先不讲解，因为这个项目 的样品数量太少了。

得到全部的bam文件如下:

$ ls -lh *bam |cut -d" " -f 5-
2.6G 2月 18 20:53 SRR1513329.hisat.bam
2.6G 2月 18 20:52 SRR1513330.hisat.bam
2.8G 2月 18 20:53 SRR1513331.hisat.bam
2.8G 2月 18 20:53 SRR1513332.hisat.bam
3.0G 2月 18 20:53 SRR1513333.hisat.bam
2.6G 2月 18 20:53 SRR1513334.hisat.bam

走featureCounts定量流程

featureCounts命令的参数多种多样，可以设置基于基因的ID或者名字来进行计算表达量，或者基于外显子来。

gtf=$HOME/rna/SUPPA2/gtf/gencode.v37.annotation.gtf 
ls -lh $gtf # 1.4G

# 针对bam文件进行qc 
ls *bam |while read id;do
nohup qualimap rnaseq --java-mem-size=20G -gtf $gtf -bam $id -pe -oc ${id%%.*} & 
done

# 走featureCounts定量流程
nohup featureCounts -T 5 -p -t exon -g gene_id -a $gtf \
-o all.id.txt *.bam 1>counts.id.log 2>&1 &
nohup featureCounts -T 5 -p -t exon -g gene_name -a $gtf \
-o all.name.txt *.bam 1>counts.name.log 2>&1 &

# https://github.com/vivekbhr/Subread_to_DEXSeq
nohup featureCounts -f -O -s 2 -p -T 5 \
-F GTF -a $gtf \
-o all.exon.txt *.bam 1>counts.exon.log 2>&1 &

仔细看了看以前的教程，但是被辣鸡的《简书》平台都封杀了：曾经我给你带来了十万用户，但现在祝你倒闭，

• https://www.jianshu.com/p/6e6c4ff9f595
• https://www.jianshu.com/p/b8f480a68cae

原来是把人家写好的轮子看错了，https://github.com/vivekbhr/Subread_to_DEXSeq

修改后的gtf走featureCounts

首先需要格式化我们的 gencode.v37.annotation.gtf 文件，代码如下:

cd $HOME/rna/SUPPA2/ref 
git clone https://github.com/vivekbhr/Subread_to_DEXSeq
pip install HTSeq 
python Subread_to_DEXSeq/dexseq_prepare_annotation2.py \
-f gencode.v37.for_dexseq.gtf gencode.v37.annotation.gtf gencode.v37.for_dexseq.gff

# 在 umac 服务器，指定全路径 
gtf="/home/yb77613/reference/gtf/gencode/gencode.v32.annotation.gtf"
gtf=/home/yb77613/reference/gtf/gencode/gencode.v37.annotation.gtf 
python Subread_to_DEXSeq/dexseq_prepare_annotation2.py \
-f gencode.v37.for_dexseq.gtf $gtf gencode.v37.for_dexseq.gff

# 这个脚本把1.4G的gencode.v37.annotation.gtf格式化成为了如下的两个文件：
# 152M 2月 19 10:45 gencode.v37.for_dexseq.gff
# 143M 2月 19 10:45 gencode.v37.for_dexseq.gtf

cd ../align/ 
gtf="$HOME/rna/SUPPA2/ref/gencode.v37.for_dexseq.gtf" 
ls -lh $gtf 
nohup featureCounts -f -O -s 2 -p -T 5 \
-F GTF -a $gtf \
-o all.for_DEXSeq.txt ../align_star/*bam 1>counts.for_DEXSeq.log 2>&1 &

下游的DEXSeq分析

项目介绍 分组情况如下:

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

还是要借助 https://github.com/vivekbhr/Subread_to_DEXSeq 里面的代码

source("Subread_to_DEXSeq/load_SubreadOutput.R")
library(tidyverse)
samp <- data.frame(row.names = c('SRR1513329','SRR1513330','SRR1513331',
 'SRR1513332','SRR1513333','SRR1513334'), 
 condition = rep(c("control","trt"),each=3))
samp
sampleData = samp
countfile="featurecounts/all.for_DEXSeq.txt"
flattenedfile = "gtf/gencode.v37.for_dexseq.gff"
dxd <- DEXSeqDataSetFromFeatureCounts(countfile,
 flattenedfile = flattenedfile,
 sampleData = samp)

# 重难点就是构建 dxd 这个对象。 
colData(dxd) 
head( counts(dxd), 5 ) 
split( seq_len(ncol(dxd)), colData(dxd)$exon ) 
head( featureCounts(dxd), 5 ) 
head( rowRanges(dxd), 3 )
sampleAnnotation( dxd )
dxd = estimateSizeFactors( dxd )
# 差异分析方法
dxd = estimateDispersions( dxd )
plotDispEsts( dxd )
dxd = testForDEU( dxd )
dxd = estimateExonFoldChanges( dxd, fitExpToVar="condition")
dxr1 = DEXSeqResults( dxd )
dxr1
mcols(dxr1)$description
table ( dxr1$padj < 0.01 )
table ( tapply( dxr1$padj < 0.01, dxr1$groupID, any ) )

## ----MvsA, fig.small=TRUE, fig.cap="MA plot. Mean expression versus $log{2}$ fold change plot. 
# Significant hits at an $FDR=0.1$ are coloured in red."----
plotMA( dxr1, cex=0.8 )
sampleAnnotation(dxd)

save(dxr1,file = 'output_of_DEXSeq.Rdata')

统计结果存储在一个非常复杂的对象里面哦！

可视化差异可变剪切结果

首先进行简单的统计：

> mcols(dxr1)$description
 [1] "group/gene identifier" 
 [2] "feature/exon identifier" 
 [3] "mean of the counts across samples in each feature/exon" 
 [4] "exon dispersion estimate" 
 [5] "LRT statistic: full vs reduced" 
 [6] "LRT p-value: full vs reduced" 
 [7] "BH adjusted p-values" 
 [8] "exon usage coefficient" 
 [9] "exon usage coefficient" 
[10] "relative exon usage fold change" 
[11] "GRanges object of the coordinates of the exon/feature" 
[12] "matrix of integer counts, of each column containing a sample"
[13] "list of transcripts overlapping with the exon" 
> table ( dxr1$padj < 0.01 )

FALSE TRUE 
227203 2333 
> table ( tapply( dxr1$padj < 0.01, dxr1$groupID, any ) )

FALSE TRUE 
 796 1206

任何可以选取特定的基因去看可变剪切情况，前提是拿到了全部的统计学显著的基因：

dat=dxr1@listData
kp=dat$padj < 0.01
table(kp)
kp[is.na(kp)]='FALSE'
table(kp)
allG=dat$groupID
sg=allG[as.logical(kp)];head(sg)
sg=unique(sg);length(sg)
head(sg)

> head(sg)
[1] "ENSG00000001497.18" "ENSG00000002586.20_PAR_Y" 
[3] "ENSG00000002834.18" "ENSG00000003393.16" 
[5] "ENSG00000004455.17" "ENSG00000004534.15+ENSG00000003756.17"

我们就可以随机选取一个进行可视化啦！

plotDEXSeq( dxr1, "ENSG00000003393.16", 
 legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 )

这应该是展现了一个很明显的外显子跳跃吧：

而其它软件，比如 SUPPA通常是会对可变剪切进行分类，如下所示：

• SE (exon skipping)
• RI (retained intron)
• MXE (mutually exclusive exons)
• A5SS (alternative 5’ splice site)
• A3SS (alternative 3’ splice site)
• ALE (alternative last exon)
• AFE (alternative first exon)

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

## 定量获得TPM值 
mkdir salmon_output
cd cleanData 
index=$HOME/rna/SUPPA2/ref/gencode.v37.transcripts.salmon.index/
ls -lh $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文件，行数都是统一的!

如果是几百个样品就需要批量操作啦：

 ls *gz|cut -d"_" -f 1|sort -u |while read id;do
 echo salmon quant -i $index -l ISF --gcBias \
 -1 ${id}_rmrRNA.fq.1.gz -2 ${id}_rmrRNA.fq.2.gz -p 2 \
 -o ../salmon_output/${id}_output
 done > run_salmon.sh 
 # for i in {0..19};do ( nohup bash ../submit.sh run_salmon.sh 20 $i 1>log.salmon.$i.txt 2>&1 & ) ;done

后续走 SUPPA2流程需要参考： https://mp.weixin.qq.com/s/kG2uI3me8Xo69mClThRutA

multipleFieldSelection.py \
-i quants/*/quant.sf -k 1 -f 4 \
-o iso_tpm.txt

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

ioe_merge_file=$HOME/rna/SUPPA2/ref/gencode.v37.all.events.ioe
ioe_merge_file=SUPPA2/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 &

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

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 # 前缀

ls -lh project_diffSplice* |cut -d" " -f 5-
cat project_diffSplice.dpsi|perl -alne '{print if $F[2] <0.01}' |wc