理论篇
一、什么是可变剪切
二、可变剪切的类型
分析篇
一、搜官网
二、逐个拆解
当需要分析一个新的测序数据或者对一种数据做新的分析的时候(按我一脸懵逼的心情,我给它取了一个名字,“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
联系客服