作者：钱胜

原载于本实验室微信公众号EvoDevo好好玩（http://mp.weixin.qq.com/s/XY5u88NikEO7WQ42F6zLZg ）

   理论篇

       一、什么是可变剪切

       二、可变剪切的类型

   分析篇

       一、搜官网

       二、逐个拆解

当需要分析一个新的测序数据或者对一种数据做新的分析的时候（按我一脸懵逼的心情，我给它取了一个名字，“de novo"分析），网上教程很多，也有很多软件，我选择了SGSeq。

理论篇

一、什么是可变剪切

可变剪切又叫选择性剪切（Alternative splicing, AS），生物的基因序列包含了外显子（exon）和内含子（intron），两者相互间隔。在mRNA前体的剪接过程中，参加剪接的外显子可以不按其线性次序剪接，内含子也可以不被切除而保留，即一个外显子或内含子是否出现在成熟mRNA中是可以选择的，这种剪接方式称为选择性剪接。AS也是转录本复杂性的一个主要来源。

二、可变剪切的类型

关于可变剪切的类型，有几种不同的说法，但都主要包含下面五种：

http://www.cnblogs.com/daimakun/p/5079689.html

1、外显子跳跃（Exon Skipping）

2、内含子保留（Intron Retention）

3、5'端可变剪接（Alternative 5' splice Site）

4、3'端可变剪接（Alternative 3' splice Site）

5、最后一个外显子可变剪接（Alternative Last Exon）

6、第一个外显子可变剪接（Alternative First Exon）

最后两个（5、6）可以概括为可变外显子（Alternative Exon），这张图的好处就是有表达量，可以很清晰地看出到底哪些地方转录了，哪些地方没有转录。

还不过瘾？那看看这个详细的介绍吧：http://mp.weixin.qq.com/s/2EC2NRPNxj_ZeNWv7qCb_w

分析篇

一、搜官网

首先到bioconductor里面找到这个包，传送门https://www.bioconductor.org/packages/release/bioc/html/SGSeq.html

这里告诉了我们SGSeq是用来分析可变剪切的R package，输入是RNA-seq比对后的BAM file，下面是安装方法，以及PDF/HTML两种格式的manual和一个版本更新信息。点开HTML/PDF其中的一种，即可看到详细的用法，共14个版块。

二、逐个拆解

①转录本的特征和可变剪切事件，SGSeq分析步骤

②下载调用SGSeq

source("https://bioconductor.org/biocLite.R")
biocLite("SGSeq")
library(SGSeq)
si   #SGSeq中的测试数据
path <- system.file("extdata", package ="SGSeq")
si$file_bam <- file.path(path,"bams", si$file_bam)   #添加路径
si   #和前面的si不一样

③注释信息：TxDb.Hsapiens.UCSC.hg19.knownGene

biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <-TxDb.Hsapiens.UCSC.hg19.knownGene
txdb <- keepSeqlevels(txdb,"chr16")
seqlevelsStyle(txdb)<-"NCBI"

用convertToTxFeatures函数转换类型对象

txf_ucsc <- convertToTxFeatures(txdb)
txf_ucsc <- txf_ucsc[txf_ucsc %over% gr]
head(txf_ucsc)

这里每个字母代表的剪切类型

• J (splice junction)

• I (internal exon)

• F (first/5′′-terminal exon)

• L (last/5′′-terminal exon)

• U (unspliced transcript)

④再次转换

sgf_ucsc <- convertToSGFeatures(txf_ucsc)
head(sgf_ucsc)

字母代表的含义

• J (splice junction)

• E (disjoint exon bin)

• D (splice donor site)

• A (splice acceptor site

⑤根据注释的转录本信息转换成剪切信息

sgfc_ucsc <- analyzeFeatures(si, features = txf_ucsc)
sgfc_ucsc
colData(sgfc_ucsc)
rowRanges(sgfc_ucsc)
head(counts(sgfc_ucsc))
head(FPKM(sgfc_ucsc))
df <- plotFeatures(sgfc_ucsc, geneID =1)

历经千辛万苦，终于出来一张图！每一行表示一个bam file，每一列表示外显子的表达，表达程度用染色表示。

⑥用从RNAseq数据中预测的信息作注释，instead of relying on existing anntation

sgfc_pred <- analyzeFeatures(si, which = gr)
head(rowRanges(sgfc_pred))
sgfc_pred <- annotate(sgfc_pred, txf_ucsc)
head(rowRanges(sgfc_pred))
df <- plotFeatures(sgfc_pred, geneID =1, color_novel ="red")

将可变剪切的外显子高亮

最后一部分，可视化

plotFeatures(sgfc_pred, geneID =1)
plotFeatures(sgfc_pred, geneName ="79791")
plotFeatures(sgfc_pred, which = gr)
plotFeatures(sgfc_pred, geneID =1, include ="junctions")
plotFeatures(sgfc_pred, geneID =1, include ="exons")
plotFeatures(sgfc_pred, geneID =1, include ="both")
plotFeatures(sgfc_pred, geneID =1, toscale ="gene")
plotFeatures(sgfc_pred, geneID =1, toscale ="exon")
plotFeatures(sgfc_pred, geneID =1, toscale ="none")
par(mfrow = c(5,1), mar = c(1,3,1,1))
plotSpliceGraph(rowRanges(sgfc_pred), geneID =1, toscale ="none", color_novel ="red")
for(j in 1:4){
  plotCoverage(sgfc_pred[, j], geneID =1, toscale ="none")
}

这里使用的都是SGSeq的测试数据，那我们自产自销的数据呢，只需要修改一下si即可，例：

si =read.table('~/data/project/scRNA/bam/tumor/SC1/config',stringsAsFactors = F)
head(si)
si <- data.frame("sample_name"= si$V1,
                 "file_bam"=  si$V2,
                 "paired_end"=TRUE ,
                 "read_length"=50,
                 "frag_length"=300,
                 "lib_size"= si$V3,stringsAsFactors = F
)
rownames(si)=si$sample_name
str(si)

快去试试吧

参考资料：

https://www.bioconductor.org/packages/devel/bioc/vignettes/SGSeq/inst/doc/SGSeq.html

http://mp.weixin.qq.com/s/XZX5pGapOMQ7EXMiJVHE4w