【原】featureCounts和DEXseq做基于外显子定量的可变剪切

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

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

准备工作

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

统一使用gencode最新版：https://www./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)

文末友情推荐

与十万人一起学生信，你值得拥有下面的学习班：

• 数据挖掘（GEO,TCGA,单细胞）2021第5期
• 生信爆款入门-2021第5期