|
理论篇 一、什么是可变剪切 二、可变剪切的类型 分析篇 一、搜官网 二、逐个拆解 当需要分析一个新的测序数据或者对一种数据做新的分析的时候（按我一脸懵逼的心情，我给它取了一个名字，“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./packages/release/bioc/html/SGSeq.html 
这里告诉了我们SGSeq是用来分析可变剪切的R package，输入是RNA-seq比对后的BAM file，下面是安装方法，以及PDF/HTML两种格式的manual和一个版本更新信息。点开HTML/PDF其中的一种，即可看到详细的用法，共14个版块。 
二、逐个拆解①转录本的特征和可变剪切事件，SGSeq分析步骤 
②下载调用SGSeq source('https:///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)
这里每个字母代表的剪切类型 ④再次转换 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./packages/devel/bioc/vignettes/SGSeq/inst/doc/SGSeq.html http://mp.weixin.qq.com/s/XZX5pGapOMQ7EXMiJVHE4w
|
