SpliceR：一个用RNA-Seq数据进行可变剪接分类和预测潜在编码区域的R包

SpliceR：一个用RNA-Seq数据进行可变剪接分类和预测潜在编码区域的R包

Kristoffer  Knudsen,  Johannes  Waage 5 Dec 2013

翻译：斑斑 <23920620>

2016年7月14日

欢迎加入生物信息QQ群78750864讨论相关问题

1 简介

SpliceR是一个可以对转录本完整isoform（剪接模式）进行分类的R包。这些转录本可以是由Cufflinks等工具通过组装RNA-seq数据产生的。SpliceR输出的是关于可变剪接的类别及转录本对无义介导的降解机制的敏感性。SpliceR还提供导出GTF格式的文件用以在基因组浏览器中浏览，以及少量的绘图功能。spliceR位于典型RNA-Seq分析流程的末端，使用由组装工具生成的full devoncoluted isoforms，并支持输出剪接元件的基因组位置，促进了序列元件的下游分析。

SpliceR在SpliceRList对象中存储转录本信息，由两个GRanges对象构成，‘转录本特征’和‘外显子特征’。作为GRanges对象的‘转录本特征’ 包含关于每个isoform的额外信息，涉及产生它的基因和表达值。作为GRanges对象的‘GRanges’包含所有的外显子和一个能回溯到‘转录本特征’的isoform ID列。SpliceRList对象可以使用内建的prepareCuff 函数通过Cufflinks的工作流程来生成，然后直接从一个Cufflins的辅助包——cummeRbund创建的cuffDB对象中读取。可选的，SpliceRList还可以从用户提供的由RNA-Seq组装工具生成的转录本和外显子信息中手动产生。

此外，spliceR通过annotatePTC()提供了关于无义介导的降解敏感性的转录本注释，通过spliceRPlot()进行可视化，通过generateGTF()生成用于基因组浏览器的GTF格式文件，以及一些辅助功能。

 

1.1 可变剪接的类型

 

spliceR使用以下剪接类型注释‘transcript features’GRanges对象。

 

?          单外显子跳跃/包含（ESI, Single exon skipping/inclusion）

?          多外显子跳跃/包含（ESI, Multiple exon skipping/inclusion）

?          内含子滞留/包含（ISI, Intron retention/inclusion）

?          可变5'/3'剪接位点（A5/A3, Alternative 5’/3’ splice sites）

?          可变转录起始位点/可变转录终止位点（ATSS/ATTS, Alternative transcription start site/ transcription termination site）

?          多外显子互斥（MEE, Mutually exclusive exons）

 

2 Cufflinks工作流程

Cufflinks是Tuxido套装的一部分，有着完善的文档和技术支持，用来对RNA-Seq数据进行映射、组装和定量。Tuxido工作流程的末端分析通常依赖cummeRbund包，它能利用R的分析和可视化能力将cufflinks的输出结果进行转换。SpliceR提供了prepareCuff()函数来将cummeRbund中的cuffSet类轻松的转换为spliceR中的SpliceRList类。SpliceR需要关于isoform及构成它的外显子的基因组坐标的全部信息。当生成cuffSet对象时，必须提供Cufflinks输出的转录本GTF文件。这里，我们使用由cummeRbund提供的内置数据来生成一个SpliceRList：

 

>   library("spliceR")

> dir <- tempdir()

>  extdata  <-  system.file("extdata",  package="cummeRbund")

>   file.copy(file.path(extdata,   dir(extdata)),   dir)

 

[1] TRUE TRUE TRUE TRUE  TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE  TRUE

[16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

 

>  cuffDB  <-  readCufflinks(dir=dir,+  gtf=system.file("extdata/chr1_snippet.gtf", package="cummeRbund"),genome="hg19" , rebuild = T)

>  cuffDB_spliceR  <-  prepareCuff(cuffDB)

 

SpliceRList是一个列表，包含两个GRanges对象，转录本特征和外显子特征，包含的信息分别为基因组坐标，转录本和转录本的基因ID，以及构成它们的外显子。使用prepareCuff()在转录本特征中能够创建附加的元列，这可以在spliceR稍后的工作流程中做进一步过滤。此外，SpliceRList对象还包含关于样品名称、数据来源、应用过的过滤条件以及被spliceR使用的其它数据。

我们还提供了一些有用的函数来提取转录本和外线子的GRanges数据，以及样品的名称：

 

>  myTranscripts  <-  transcripts(cuffDB_spliceR)

>  myExons  <-  exons(cuffDB_spliceR)

>    conditions(cuffDB_spliceR)

 

[1] "iPS"              "hESC"            "Fibroblasts"

 

一旦SpliceRList对象准备完毕，它可以用spliceR()函数，或者基于每个转录本的潜在编码区域信息的annotatePTC()函数来进行进一步的可变剪接类型的注释。

 

2.1 Cufflinks基因注释的问题

Cufflinks使用一个内部的（自定义的）‘gene id’。根据官方手册，这个基因ID 是一个用于描述基因的唯一标识。但是根据调查，我们发现许多Cufflinks的基因ID实际上包含多个注释的基因。据我们所知，这个问题在用Ensembl/Gencode作为参考转录组，对unstranded RNA-seq数据执行Cufflinks/Cuffdiff时尤为明显。在这些条件下，我们发现200-400个（总体的2-4%）cuffgenes包含多个注释了的基因。我们的假设为，这一问题的出现是由于这些基因间存在跨越基因间区的长的Ensembl转录本，因此欺骗了Cufflinks使其认为所有源自这一区域的转录本都属于同一个基因。

Cufflinks注释问题的全局效应是使所有这类基因的表达水平错误，因为这里的表达水平是几个基因的平均值。这反过来又影响基因的差异表达分析。Cufflinks注释问题只是影响基因的表达水平，这意味着它们转录本的表达水平和差异表达检验并不受这个问题的影响。

我们对这一问题粗略的解决方案为：1）识别所有受影响的基因。2）为cuffgene中的每个真实注释的基因创建一个新的cuffgene。3）根据新的转录本排布重新计算基因表达。4）设定显著性检验为负。

这会校正基因的表达并移除假阳性的差异表达基因。这一设定默认为开（prepareCuff()中的fixCufflinksAnnotationProblem选项）

 

3 GRanges 工作流程

对从其他非Cufflinks的转录本组装工具得到的RNA-seq数据，SpliceRList可以简单的通过转录本和外显子GRanges类的创建函数来手动创建。在下面的例子中，我们下载了一个公用的转录本和外显子模型数据，虚拟出样本间的表达量和p值来模拟常见的RNA-seq组装数据。

首先，从UCSC储存仓库中下载UCSC “knownGene”的转录本列表：

 

> session <- browserSession("UCSC")

>  genome(session)  <-  "hg19"

> query <- ucscTableQuery(session, "knownGene")

> tableName(query) <- "knownGene"

>  cdsTable  <-  getTable(query)

 

接下来，我们从UCSC下载kgXref表，它包含了基因ID和转录本ID间的关系。

 

>  tableName(query)  <-  "kgXref"

> kgXref <- getTable(query)

 

创建包含转录本特征的GRanges（设定假的基因表达，倍数变化和p值），从kgXref表中获取基因名字：

 

> knownGeneTranscripts <- GRanges(

+                   seqnames=cdsTable$"chrom",

+                 ranges=IRanges(

+                               start=cdsTable$"txStart",

+                              end=cdsTable$"txEnd"),

+                  strand=cdsTable$"strand",

+                    spliceR.isoform_id  =  cdsTable$"name",

+                  spliceR.sample_1="placeholder1",

+                  spliceR.sample_2="placeholder2",

+                  spliceR.gene_id=kgXref[match(cdsTable$"name",  kgXref$"kgID"),

+                              "geneSymbol"],

+                   spliceR.gene_value_1=1,

+                   spliceR.gene_value_2=1,

+                  spliceR.gene_log2_fold_change=log2(1/1),

+                   spliceR.gene_p_value=1,

+                   spliceR.gene_q_value=1,

+                 spliceR.iso_value_1=1,

+                 spliceR.iso_value_2=1,

+                    spliceR.iso_log2_fold_change=log2(1/1),

+                 spliceR.iso_p_value=1,

+                 spliceR.iso_q_value=1

+               )

 

创建包含外显子特征的GRanges：

 

> startSplit              <- strsplit(as.character(cdsTable$"exonStarts"), split=",")

> endSplit                  <- strsplit(as.character(cdsTable$"exonEnds"), split=",")

> startSplit   <- lapply(startSplit, FUN=as.numeric)

> endSplit             <- lapply(endSplit, FUN=as.numeric)

> knownGeneExons <- GRanges(

+             seqnames=rep(cdsTable$"chrom", lapply(startSplit, length)),

+                 ranges=IRanges(

+                                start=unlist(startSplit)+1,

+                               end=unlist(endSplit)),

+                 strand=rep(cdsTable$"strand", lapply(startSplit, length)),

+                     spliceR.isoform_id=rep(knownGeneTranscripts$"spliceR.isoform_id",

+                lapply(startSplit, length)),

+                  spliceR.gene_id=rep(knownGeneTranscripts$"spliceR.gene_id",

+                lapply(startSplit, length))

+               )

 

最后，从两个GRanges来创建SpliceRList，设置和assembly，coditions source和coditions条件：

 

>  knownGeneSpliceRList  <-  SpliceRList(

+                  transcript_features=knownGeneTranscripts,

+                   exon_features=knownGeneExons,

+                  assembly_id="hg19",

+                source="granges",

+                  conditions=c("placeholder1", "placeholder2")

+               )

 

4 preFiltering

数据集由e.g产生。在多个样本中，考虑到所有样本的组合，Cufflinks（占用的内存）会很快变的非常大。Cufflinks报告的许多基因和转录本是基于运行时提供的GTF参考文件得到的，（这些基因和转录本）可能在任意的或全部样品中都没有表达。

类似的，cufflinks可能并不信任所有报告得到的转录本，这些转录本被标记为‘FAIL’或 ‘LOW DATA’。在运行spliceR()和下游分析前将这些基因和转录本从SpliceRList中过滤出去将减少相当多的运行时间。

这里给出一个使用标准过滤器进行预处理的例子。更多其它的过滤条件，请参考spliceR() 的帮助页面。注意：preSpliceRFilter()报告的是总的成对比较的数量，而spliceR()报告的是唯一的isoform的总数量。

 

>  cuffDB_spliceR_filtered  <-  preSpliceRFilter(

+                cuffDB_spliceR,

+               filters=c("expressedIso", "isoOK", "expressedGenes", "geneOK")

+               )

 

 

5 SpliceR

spliceR()获取一个SpliceRList对象并对isoform注释可变剪接类型，以及一些其它的元列，数据可以是由prepareCuff()通过GRanges对象生成，也可以手动生成。设定给preTranscript 或一个已知样本的compareTo参数会告诉spliceR()当存在分类事件时，如何挑选出isoform的比较结果。使用preTranscript参数，spliceR()聚集全部的外显子信息，为每一个唯一的gene ID创建一个基于全部isoform的pre-RNA。然后当分类时，每个isoform同这个结构进行比较。通过给出一个样本名来使spliceR()把这个样本的主要isoform（根据isoform相对于整个基因的表达）作为主要的参考RNA。当分类时，将每个isoform同这个参考进行比较。

spliceR()可以进行一定量的过滤来减少总数据的大小和运行时间。大部分过滤器只在运行cufflinks产生的数据时才有意义，如同这个软件对每个转录本和基因产生的附加信息，这些信息表明模型中整体的稳定性和可信性。

下面展示了一个典型的对spliceR()的调用，设置常用的过滤器和pre-RNA的参考。全部过滤器的详细内容请参考spliceR()的帮助页面。

 

> #Commented out due to problems with vignettes and progress bars.

>  mySpliceRList  <-  spliceR(

+                cuffDB_spliceR,

+                 compareTo="preTranscript",

+                 filters=c("expressedGenes","geneOK", "isoOK", "expressedIso", "isoClass"),

+                  useProgressBar=F

+ )

 

在主流的工作站上，运行4个样本，每个样本不超过100，000转录本的数据大概需要1.5-5个小时的时间。spliceR()结束后，transcriptFeatures GRanges被一些额外的列填充。它们包含了isoform fraction（IF）值和delta IF。每个剪接类型占一列，表示每个转录本在对应的类型中有多少剪接事件被探测到。最后，“analyzed”列表示这个转录本是否通过了筛选。

 

6 PTC-注释

annotatePTC()获取一个SpliceRList对象，一个BSGenome序列对象和一个CDSList对象。BSGenome序列对象可以从Bioconductor下载并安装。此处的组装名和染色体名应当同用于cufflinks或其它组装工具进行映射的组装名相一致（如：BSgenome.Hsapiens.UCSC.hg19或BSgenome.Mmusculus.UCSC.mm9）。CDSList对象由getCDS()生成，getCDS()获取一个组装名、基因追踪名以及从UCSC基因组浏览器表中获取的CDS信息（需要联网）。可选的，我们可以通过一个选定的CDS表手动创建一个CDSList对象（参考？CDSList获得更多信息）。通过CDSList和BSGenome对象，annotatePTC()为 SpliceRList exon_features GRanges中包含的所有外显子获取其基因组序列，历遍transcript_features GRanges中的每个转录本，连接它们的外显子并翻译含有一个CDS的转录本。如果存在一个以上的CDS，则采用最上游的。transcript_features GRanges通过终止密码子的位置（在Mrna坐标系中），终止密码子所在外显子（0是最后一个外显子，-1是倒数第二个外显子，以此类推），最后的外显子-外显子连接距离和使用的CDS的基因组位置等信息来对每个转录本进行注释。如果没有发现CDS，这些列被设为NA。最终，转录本被注释为TURE或FALSE及关于无义介导的降解敏感性（NMD）。这个参数默认值为50（表示转录本终止位置距最后的外显子-外显子剪接位点超过50nt，且不在最后的外显子上），但可能会在程序运行时被改变。下面是一个典型的annotatePTC()工作流程：

 

> ucscCDS <- getCDS(selectedGenome="hg19", repoName="UCSC")

>    require("BSgenome.Hsapiens.UCSC.hg19",character.only    =    TRUE)

> #Commented out due to problems with vignettes and progress bars.

> #PTCSpliceRList <- annotatePTC(cuffDB_spliceR, cds=ucscCDS, Hsapiens,

>  #                PTCDistance=50)

 

7 GTF 输出

在cufflinks或其它组装工具装配过后，又或者spliceR() 及 annotatePTC()处理过后，将转录本在基因组浏览器中可视化是非常有帮助的。为了更加便利，spliceR提供了generateGTF函数。generateGTF()为每个样本生成GTF文件，并允许采用同spliceR()和annotatePTC()中一样的过滤条件来进行转录本的筛选。此外，generateGTF()颜色根据表达量编码转录本。通过参数scoreMethod设定精确打分方法，“local”打分是转录本相对于唯一基因ID中表达量最高的转录本，而“global”打分是转录本相对于整个样本。对generateGTF()函数的典型调用看起来是这样的：

 

> generateGTF(mySpliceRList, filters=

+  c("geneOK",  "isoOK",  "expressedGenes",  "expressedIso"),

+                 scoreMethod="local",

+                 useProgressBar=F)

 

GTF文件输出在当前的工作目录下，它的名字同样本的名字相一致，并在每个文件中添加了filePrefix标记。

 

8 可视化

运行spliceR()后，使用可变剪接事件的相关信息注释spliceRlist对象。spliceRPlot()函数可以用于对原始数据的探索。spliceRPlot()生成一个用圆代表各种情况的韦恩图来分析数据的各个方面。

spliceRPlot()包含两个控制层次。第一个层次，用户可以通过expressionCutoff和 filters两个参数来控制用于绘图的数据。由于韦恩图的基本目的是区分不同情况中表达了的转录本，因此expressionCutoff参数用来控制界定“表达”的阈值。接下来类似spliceR的其它核心函数，spliceRPlot()含有一个过滤参数。通过这个过滤参数，可以选择数据的子集。

过滤数据的准备和转换过程可能需要几分钟的时间，之后spliceRPlot()返回一个spliceRlist对象，中间的数据结构会被存储，这可以使用户在生成了其它图的情况下跳过这步。如果用户改变了expressionCutoff或filters参数，reset参数需要被设置为TRUE。

第二个层次的控制决定实际的绘图。有一些选项可用于估计参数，在图形中安排每种情况中表达基因的数目，及表达的转录本中发现的可变剪接的平均数量。此外，“asType”参数允许用户来指定可变剪接的子类型来进行评估。

 

> #Plot the average number of transcripts pr gene

>   mySpliceRList   <-   spliceRPlot(mySpliceRList,

+                   evaluate="nr_transcript_pr_gene")

 

 

> #Plot the average number of Exon skipping/inclucion evens pr gene

> mySpliceRList <- spliceRPlot(mySpliceRList, evaluate="mean_AS_transcript",

+     asType="ESI")

 

 

> #Plot the average number of Exon skipping/inclucion evens pr gene,

> #but only using transcripts that are significantly differntially expressed

> mySpliceRList <- spliceRPlot(mySpliceRList, evaluate="mean_AS_transcript",

+     asType="ESI",      filters="sigIso",reset=TRUE)

 

 

>  sessionInfo()

 

R version 3.0.2 Patched (2013-12-18 r64488) Platform:  x86_64-unknown-linux-gnu  (64-bit)

 

locale:

[1] LC_CTYPE=en_US.UTF-8    LC_NUMERIC=C

[3] LC_TIME=en_US.UTF-8 LC_COLLATE=C

[5] LC_MONETARY=en_US.UTF-8        LC_MESSAGES=en_US.UTF-8

[7] LC_PAPER=en_US.UTF-8    LC_NAME=C

[9] LC_ADDRESS=C LC_TELEPHONE=C

[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

 

attached base	packages:
[1] grid	parallel       stats	graphics     grDevices utils	datasets
[8] methods	base

 

other attached packages:

[1] BSgenome.Hsapiens.UCSC.hg19_1.3.19 BSgenome_1.30.0

[3] Biostrings_2.30.1        spliceR_1.3.1

 

 

[5] plyr_1.8     RColorBrewer_1.0-5

[7] VennDiagram_1.6.5        cummeRbund_2.4.1

[9] Gviz_1.6.0      rtracklayer_1.22.2

[11] GenomicRanges_1.14.4    XVector_0.2.0

[13] IRanges_1.20.6  fastcluster_1.1.13

[15] reshape2_1.2.2   ggplot2_0.9.3.1

[17] RSQLite_0.11.4    DBI_0.2-7

[19] BiocGenerics_0.8.0

 

loaded via a namespace (and not attached):

[1] AnnotationDbi_1.24.0   Biobase_2.22.0 Formula_1.1-1 [4] GenomicFeatures_1.14.2 Hmisc_3.14-0   MASS_7.3-29 [7] RCurl_1.95-4.1      Rsamtools_1.14.2    XML_3.98-1.1 [10] biomaRt_2.18.0      biovizBase_1.10.7    bitops_1.0-6

[13] cluster_1.14.4     colorspace_1.2-4   dichromat_2.0-0 [16] digest_0.6.4   gtable_0.1.2     labeling_0.2 [19] lattice_0.20-24    latticeExtra_0.6-26 munsell_0.4.2 [22] proto_0.3-10     scales_0.2.3     splines_3.0.2 [25] stats4_3.0.2      stringr_0.6.2     survival_2.37-7 [28] tools_3.0.2      zlibbioc_1.8.0