1概述

SGSeq为分析RNA-seq数据中的剪接事件提供了一个框架。作为该框架的一部分，该软件实现了用于存储与拼接变量分析相关的信息的数据结构。概述SGSeq类以及它们之间的关系如图1所示，并总结如下:

• 的TxFeaturesclass存储在RNA转录本中观察到的离散转录本特征(外显子和剪接连接)。这些可以从注释中提取或从对齐的RNA-seq读取中预测。
• 的SGFeatures类存储定义拼接图的特性(Heber et al. 2002)．剪接图是一种有向无环图，其边对应于外显子区和剪接结点，节点对应于转录本的起始、结束和剪接位点。它是从5点开始的\ (^ \ ' \)3结束\ (^ \ ' \)基因的末端。
• 的SGVariants类存储拼接变量。如果拼接图中的两个节点由两条或多条路径连接，并且没有所有路径相交的中间节点，则认为备选路径为拼接变体。共享相同开始节点和结束节点的拼接变量一起形成一个拼接事件。

2预赛

库(SGSeq)

当开始一个新项目时，SGSeq需要被分析样品的信息。这些信息最初只获得一次，然后可以用于所有后续的分析。示例信息以数据帧的形式提供，包含以下列:

• sample_name每个样本具有唯一名称的字符向量
• file_bam字符向量或BamFileList指定使用支持拼接的对齐程序生成的BAM文件
• paired_end表示数据是成对端还是单端的逻辑向量
• read_length具有读取长度的数字向量
• frag_length具有平均片段长度的数字向量(用于成对端数据)
• lib_size数值向量，其中包含单端数据的对齐读总数，或成对端数据的一致对齐读对总数

示例信息可以存储在data.frame或DataFrame对象(如果BAM文件被指定为BamFileList，则必须存储在DataFrame)．具有自动获取样品信息的功能getBamInfo ()，它将一个带列的数据帧作为输入sample_name而且file_bam并从指定的BAM文件提取所需的信息。

为SGSeq为了正确地工作，读取必须使用一个感知拼接的对齐程序(如GSNAP)进行映射(t.d. Wu and Nacu 2010), HISAT(Kim, Langmead, and Salzberg 2015)或明星(Dobin et al. 2013)，生成带有自定义标记“XS”的SAM/BAM文件，用于拼接读取，指示转录方向。请注意,lib_size应该是对齐片段的总数，即使BAM文件是感兴趣的基因组区域的子集。为了精确计算FPKM值(每千碱基片段数和百万序列片段数)，需要总片段数。在这里，术语“片段”表示已测序的cDNA片段，它由单端数据中的一次读取表示，或成对端数据中的一对读取表示。

本插图说明了对来自四个肿瘤和四个正常结直肠样本的配对端RNA-seq数据的分析，这些数据是发表在(Seshagiri et al. 2012)．RNA-seq读数使用GSNAP映射到人类参考基因组(t.d. Wu and Nacu 2010)．分析基于BAM文件子集，以读取映射到感兴趣的单个基因(FBXO31)．一个data.frame如果从带函数的原始BAM文件中生成样本信息getBamInfo ()．注意这一栏lib_size反映原始BAM文件中对齐片段的总数。

如果

## sample_name file_bam paired_end read_length frag_length lib_size ## 1 N1 N1bam TRUE 75 293 12405197 ## 2 N2。bam TRUE 75 197 13090179 ## 3 N3 N3。bam TRUE 75 206 14983084 ## 4 N4 N4。bam TRUE 75 207 15794088 ## 5 T1 T1。bam TRUE 75 284 14345976 ## 6 T2 T2。bam TRUE 75 235 15464168 ## 7 T3 T3。bam TRUE 75 259 15485954 ## 8 T4 T4。bam TRUE 75 247 15808356

下面的代码块为当前文件设置正确的BAM文件路径SGSeq安装。

路径<- system. Path。file("extdata"， package = "SGSeq") si$file_bam <- file. file("extdata"， package = "SGSeq")路径(Path， "bams"， si$file_bam)

3.RNA转录本TxFeatures类

文本注释可以通过TxDb对象或使用函数从GFF格式导入importTranscripts ()．或者，可以将转录本指定为GRangesList外显子按转录本分组。在下面的代码块中，UCSC knownGene表作为TxDb对象。16号染色体上的转录本FBXO31基因定位)被保留，并且染色体名称被更改为与BAM文件中使用的命名约定相匹配。

txdb <- txdb . hsapiens . ucsc .hg19. knowngeneknownGene txdb <- keepSeqlevels(txdb， "chr16") seqlevelsStyle(txdb) <- "NCBI"

中的注释SGSeq框架，从文本特征中提取TxDb使用函数的对象convertToTxFeatures ()．只有基因位点重叠的特征FBXO31基因被保留。基因组坐标FBXO31存储在农庄对象gr．

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

## seqnames范围链类型##     ## [1] 16 [87362942, 87365116] - l# # [2] 16 [87365116, 87367492] - j# # [4] 16 [87367892,87367892] - j# # [5] 16 [87368910, 87369063] - i# # [6] 16 [87369063, 87369761] - j# # txName geneName ##   ## [1] uc002fjv.3,uc002fjw.3,uc010vot。2 79791 ## [2] uc002fjv.3,uc002fj .3,uc010vot。2 79791 ## [3] uc002fjv.3,uc002fj .3,uc010vot。2 79791 ## [4] uc002fjv.3,uc002fj .3,uc010vot。2 79791 ## [5] uc002fjv.3,uc002fj .3,uc010vot。2 79791 ## [6] uc002fjv.3,uc002fj .3,uc010vot。279791## ------- ## seqinfo: 1 sequence from hg19 genome

convertToTxFeatures ()返回一个TxFeatures对象，它是a农庄类对象，具有额外的列。列类型特征类型，可以取值

• J(接头连接)
• 我(内部外显子)
• F(第一个/ 5\ (^ \ ' \)终端外显子)
• l(最后一个/ 5\ (^ \ ' \)终端外显子)
• U(unspliced成绩单)。

列txName而且geneName指出每个特征所属于的转录本和基因。请注意，一个特性可以属于多个转录本，因此列可以为每个特性存储多个值。为TxFeatures中定义的其他数据结构SGSeq，列可以使用以相关列命名的函数访问，如下面的代码块所示。

类型(txf_ucsc)

## [1] L J I J I J I J I J I J I J I J I F J J F U I J F F ##级别:J I F L U

头(txName (txf_ucsc))

##长度为6的字符列表## [[1]]uc002fjv。3 uc002fjw。3 uc010vot。2## [[2]] uc002fjv.3 uc002fjw.3 uc010vot.2 ## [[3]] uc002fjv.3 uc002fjw.3 uc010vot.2 ## [[4]] uc002fjv.3 uc002fjw.3 uc010vot.2 ## [[5]] uc002fjv.3 uc002fjw.3 uc010vot.2 ## [[6]] uc002fjv.3 uc002fjw.3 uc010vot.2

头(geneName (txf_ucsc))

长度6 # # # # CharacterList [[1]] 79791 # # 79791 # # [[2]] [[3]] 79791 # # 79791 # # [[4]] [[5]] 79791 # # [[6]] 79791

4拼接图和SGFeatures类

外显子在TxFeatures对应于RNA分子中拼接在一起的外显子。在拼接图的上下文中，外显子由唯一的非重叠外显子区域表示。函数convertToSGFeatures ()转换TxFeatures来SGFeatures．在这个过程中，重叠的外显子被分离成非重叠的外显子箱。

sgf_ucsc <- convertToSGFeatures(txf_ucsc) head(sgf_ucsc)

## SGFeatures对象，包含6个范围和0个元数据列:## seqnames range strand type splice5p splice3p ##       ## [1] 16 [87362942,87365116] - E TRUE FALSE ## [2] 16 [87365116, 87365116] - A   ## [3] 16 [87365116, 87367492] - J   ## [4] 16 [87367492, 87367492] - D   ## [5] 16 [87367492, 87367892] - E TRUE TRUE ## [6] 16 [87367892, 87367892，87367892] - A   ## featureID geneID txName geneName ##     ## [1] 1 1 uc002fjv.3,uc002fjw.3,uc010vot. 3,uc010vot. 32 79791 ## [2] 2 1 uc002fjv.3,uc002fjw.3,uc010vot。2 79791 ## [3] 3 1 uc002fjv.3,uc002fjw.3,uc010vot。2 79791 ## [4] 4 1 uc002fjv.3,uc002fjw.3,uc010vot。2 79791 ## [5] 5 1 uc002fjv.3,uc002fjw.3,uc010vot。2 79791 ## [6] 6 1 uc002fjv.3,uc002fjw.3,uc010vot。279791## ------- ## seqinfo: 1 sequence from hg19 genome

类似于TxFeatures,一个SGFeatures对象是农庄类对象，具有额外的列。列类型对于一个SGFeatures对象接受值

• J(接头连接)
• E(分离外显子仓)
• D(拼接供体部位)
• 一个(剪接受体位点)。

按照惯例，剪接的供体和受体位点对应于内含子侧翼的外显子位置。SGFeatures是否没有包含其他列TxFeatures：spliced5p而且spliced3p指示外显子箱是否在5处有强制剪接\ (^ \ ' \)和3\ (^ \ ' \)分别结束。该信息用于确定读取的内容在结构上是否与外显子箱兼容，以及外显子箱是否与注释的转录本一致。列featureID为每个特性提供唯一标识符，而列geneID表示一个特征所属的拼接图的唯一连接组件。

5基于注释文本的拼接图分析

本节说明基于注释文本的分析。函数analyzeFeatures ()将转录本特征转换为拼接图特征，并获得每个特征和每个样本的兼容片段计数。

sgfc_ucsc <- analyzeFeatures(si, features = txf_ucsc) sgfc_ucsc

##类:sgfeaturecrets ## dim: 42 8 ##元数据(0):## assays(2):计数FPKM ## rownames: NULL ## rowData names(0): ## colnames(8): N1 N2…T3 T4 ## colData names(6): sample_name file_bam…frag_length lib_size

analyzeFeatures ()返回一个SGFeatureCounts对象。SGFeatureCounts包含样例信息为colData，拼接图特征为rowRanges和化验计数而且FPKM，分别存储兼容的片段计数和fpkm。可以使用同名的访问器函数访问不同的数据类型。

colData(sgfc_ucsc) rowRanges(sgfc_ucsc) head(counts(sgfc_ucsc)) head(FPKM(sgfc_ucsc))

外显子和剪接连接的计数是基于结构兼容的片段。在剪接供体和受体的情况下，计数表示跨越剪接边界(重叠剪接位点和侧翼内含子位置)的读片段的数量。

FPKM值计算为\ \(压裂{x}{问}10 ^ 6 \),在那里\ \ (x)是兼容片段的数量，\ (N \)库大小(存储在lib_size),l有效特征长度(兼容片段可能位置的数量)。对于成对端数据，假设片段长度等于frag_length．

用于拼接图特征的fpkm可以用函数进行可视化plotFeatures．plotFeatures生成一个双面板图，其中显示在顶部面板中的拼接图和底部面板中单个特征的表达水平热图。用于定制plotFeatures输出，参见节可视化．绘图函数无形地返回data.frame在图中包含有关拼接图特征的详细信息。

df <- plotFeatures(sgfc_ucsc, geneID = 1)

请注意，拼接图包括三个可选的转录起始点(tss)。然而，热图表明在这个数据集中的样本中没有使用第一个TSS。

6基于拼接图分析新创预测

而不是依赖于现有的注释，注释可以通过RNA-seq数据的预测来增强，或者可以在不使用注释的情况下从RNA-seq数据构建拼接图。下面的代码块预测RNA-seq读取支持的转录本特征，将其转换为拼接图特征，然后获得兼容的片段计数。有关如何获得预测的详细信息，请参阅(Goldstein et al. 2016)．

sgfc_pred <- analyzeFeatures(si, which = gr) head(rowRanges(sgfc_pred))

## SGFeatures对象，包含6个范围和0个元数据列:## seqnames range链类型splice5p splice3p ##       ## [1] 16 [87362930,87365116] - E TRUE FALSE ## [2] 16 [87365116, 87365116] - A   ## [3] 16 [87365116, 87367492] - J   ## [4] 16 [87367492, 87367492] - D   ## [5] 16 [87367492, 87367892] - E TRUE TRUE ## [6] 16 [87367892, 87367892，87367892] - A   ## featureID geneID txName geneName ##     ## [1] 1 1 ## [2] 2 1 ## [3] 3 1 ## [4] 4 1 ## [5] 5 1 ## [6] 6 1 ## ------- ## seqinfo: 84 sequences from an unspecified genome

对于解释，预测的特征可以根据已知的转录本进行注释。的注释()函数为每个特征分配兼容的转录本，并将相应的转录本和基因名存储在列中txName而且geneName,分别。

sgfc_pred <-注释(sgfc_pred, txf_ucsc) head(rowRanges(sgfc_pred))

## SGFeatures对象，包含6个范围和0个元数据列:## seqnames range链类型splice5p splice3p ##       ## [1] 16 [87362930,87365116] - E TRUE FALSE ## [2] 16 [87365116, 87365116] - A   ## [3] 16 [87365116, 87367492] - J   ## [4] 16 [87367492, 87367492] - D   ## [5] 16 [87367492, 87367892] - E TRUE TRUE ## [6] 16 [87367892, 87367892，87367892] - A   ## featureID geneID txName geneName ##     ## [1] 1 1 uc002fjv.3,uc002fjw.3,uc010vot. 3,uc010vot. 32 79791 ## [2] 2 1 uc002fjv.3,uc002fjw.3,uc010vot。2 79791 ## [3] 3 1 uc002fjv.3,uc002fjw.3,uc010vot。2 79791 ## [4] 4 1 uc002fjv.3,uc002fjw.3,uc010vot。2 79791 ## [5] 5 1 uc002fjv.3,uc002fjw.3,uc010vot。2 79791 ## [6] 6 1 uc002fjv.3,uc002fjw.3,uc010vot。279791## ------- ## seqinfo: 84 sequences from an unspecified genome

预测的拼接图和fpkm可以像前面一样可视化。使用参数突出显示缺少注释的拼接图特征color_novel．

df <- plotFeatures(sgfc_pred, geneID = 1, color_novel = "red")

请注意，从RNA-seq数据预测的大多数外显子和剪接连接与UCSC knownGene表中的转录本一致(灰色显示)。然而，与前图相比，预测的基因模型不包括拼接图中未在数据中表达的部分。此外，从RNA-seq数据中发现了一个未注释的外显子(E3，红色部分)，该外显子在4个正常结肠直肠样本中的3个(N2, N3, N4)中表达。

7剪接变体识别

而不是考虑一个基因的完整剪接图，分析可以集中在个别剪接事件。函数analyzeVariants ()从图中递归地标识剪接事件，获得每个剪接变体的代表性计数，并估计它们的相对使用情况。

sgvc_pred <- analyzeVariants(sgfc_pred) sgvc_pred

##类:SGVariantCounts ## dim: 2 8 ##元数据(0):## assays(5): countsVariant5p countsVariant3p countsev5p ## countsev3p variantFreq ## rownames: NULL ## rowData names(20): from to…## colnames(8): N1 N2…T3 T4 ## colData names(6): sample_name file_bam…frag_length lib_size

analyzeVariants ()返回一个SGVariantCounts对象。示例信息存储为colData,SGVariants作为rowRanges．分析variantFreq存储每个拼接变体和样本的相对使用估计。与前面一样，可以使用同名的访问器函数访问不同的数据类型。关于拼接变量的信息存储在SGVariants元数据列，可以通过函数访问mcols ()如下图所示。有关列的详细描述，请参阅的手册页SGVariants．

mcols (sgvc_pred)

# # DataFrame与2 # # 20行和列类型featureID segmentID # # <人物> <人物> <人物> <人物> <人物> # # 1 D: 16:87393901:——:16:87380856:28 - J 4 # # 2 D: 16:87393901:——:16:87380856:——JEJ 32岁,30日,27日2 # # closed5p closed3p closed5pEvent closed3pEvent geneID eventID # # <逻辑> <逻辑> <逻辑> <逻辑> <整数> <整数> # # 1真的真的真的真的1 1 # # 2真的真的真的真的1 1 # # variantID featureID5p featureID3p featureID5pEvent featureID3pEvent # #     ## 1 1 28 28 28,32 28,27 ## 2 2 32 27 28,32 28,27 ## txName geneName variantType ##    ## 1 uc002fjv.3,uc002fjw.3,uc010vot。2 79791 SE:S ## 2 79791 SE:I ## variantName ## <字符> ## 1 79791_1_1/2_SE ## 2 79791_1_2/2_SE

8拼接变量量化

拼接变体是局部量化的，基于结构兼容的片段，这些片段重叠在每个变体的开始或结束。相对使用量的本地估计\ \ psi_i \ ()的变体我\ \ ()得到的片段数与我\ \ ()除以与属于同一事件的任何变体兼容的片段数。对于变型启动\ (\)变量末端\ \ (E)\(\psi_i^S = x_i^S / x_.^S\)而且\(\psi_i^E = x_i^E / x_.^E\),分别。对于具有有效估计的变量\ (\ psi_i ^ \)而且\ \ (\ psi_i ^ E)单个估计值是局部估计值的加权平均值\(\psi_i = x_.^S/(x_.)^S + x_.^E) \psi_i^S + x_.^E/(x_。^S + x_.^E) \psi_i \ E\．

对于读计数较低的事件，相对使用率的估计可能不可靠。如果参数min_denominator为函数指定。analyzeVariants ()或getSGVariantCounts ()，估计值设置为NA除非至少有一个\(间。^ \)或\(间E ^ \)。等于或大于指定的值。

请注意,SGVariantCounts对象还存储原始计数数据。这些数据可以用于统计建模，例如小节中所建议的测试不同拼接变量的使用情况．

variantFreq (sgvc_pred)

## n1 n2 n3 n4 t1 t2 t3 ## [1，] 0.88 0.56 0.6153846 0.5925926 0.93333333 0.90196078 0.8484848 ## [2，] 0.12 0.44 0.3846154 0.4074074 0.06666667 0.09803922 0.1515152 ## t4 ## [1，] 0.95833333 ## [2，] 0.04166667

剪接变量和相对使用的估计是可视化的功能plotVariants．

plotvariables (sgvc_pred, eventID = 1, color_novel = "red")

plotVariants生成类似于的双面板图形plotFeatures．顶部面板中的拼接图说明了所选的拼接事件。在本例中，事件由两个变体组成，对应于跳过或包含未注释的外显子。热图说明了每个拼接变体的相对使用估计。样本N2, N3和N4显示了包括外显子以及跳过外显子的转录本的证据。其余样本显示很少有外显子包含的证据。

9拼接变体解释

剪接变异的功能后果可以通过预测其对蛋白质编码潜力的影响来评估。函数predictVariantEffects ()作为输入SGVariants对象，其中包含感兴趣的拼接变体、一组带注释的转录本和存储为BSgenome对象。

library(BSgenome.Hsapiens.UCSC.hg19) seqlevelsStyle(Hsapiens) <- "NCBI" vep <- predictvarianteeffects (sgv_pred, txdb, Hsapiens) vep

## variantID txName geneName RNA_change ## 1 uc002fjv。2 2 uc002fjw. 1;3.79791r.412_413ins412+1798_412+1884 ## 3 2 uc010vot.2 79791 r.-105_-104ins-105+1798_-105+1884 ## RNA_variant_type protein_change ## 1 insertion p.K29_L30insRINPRVKSGRFVKILPDYEHMAYRDVYTC ## 2 insertion p.K137_L138insRINPRVKSGRFVKILPDYEHMAYRDVYTC ## 3 insertion p.= ## protein_variant_type ## 1 in-frame_insertion ## 2 in-frame_insertion ## 3 no_change

输出是一个数据帧，每一行都描述了一个特定的剪接变体对注释的蛋白质编码转录本的影响。变量的影响如下所述HGVS发布建议．在其目前的实现中，变异效应预测相对较慢，建议运行predictVariantEffects ()仅在选定的变量上。

10可视化

功能plotFeatures ()而且plotVariants ()支持许多自定义图形的选项。顶部面板中的拼接图是按函数绘制的plotSpliceGraph，可直接调用。

plotFeatures ()包含多个参数，用于选择要显示的特性。下面的代码块说明了用于绘制的拼接图和表达式级别的三个不同选项FBXO31(Entrez ID 79791)。

plotFeatures(sgfc_pred, geneID = 1) plotFeatures(sgfc_pred, geneName = "79791") plotFeatures(sgfc_pred, which = gr)

缺省情况下，由plotFeatures ()显示拼接连接。或者，可以显示外显子箱，或者外显子箱和剪接结都可以显示。

plotFeatures(sgfc_pred, geneID = 1, include = " connections ") plotFeatures(sgfc_pred, geneID = 1, include = "exons")

论点toscale控制按比例绘制基因模型的哪些部分。

plotFeatures(sgfc_pred, geneID = 1, toscale = "gene") plotFeatures(sgfc_pred, geneID = 1, toscale = "exon")

热图可以可视化剪接结和外显子箱的表达值。或者，每个基读覆盖和拼接结计数可以用函数可视化plotCoverage．

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")}

11测试不同拼接变量的使用情况

SGSeq不实现对差异拼接变量使用情况的统计测试。然而，现有的软件包等DEXSeq(安德斯，雷耶斯，Huber 2012)而且limma(Ritchie et al. 2015)可用于此目的。这些包允许在基因内(在样本组之间)推断差异外显子的使用。在SGSeq框架，这种方法可以用来测试在剪接事件中的差异剪接变体的使用，其中剪接变体和剪接事件分别类似于外显子和基因。

为了使这些方法适用，每个拼接变体都需要一个计数。SGVariantCounts如上所述的对象为每个拼接变量存储两个计数，一个用于5\ (^ \ ' \)一个是3\ (^ \ ' \)变体结束。这些数字可以很容易地从SGFeatureCounts对象，但对于基于计数的差分测试是不切实际的。每个变体的单个计数(基于变体两端兼容的片段)可以使用函数从BAM文件中获得getSGVariantCounts ()．输出是SGVariantCounts带有附加分析的对象countsVariant5pOr3p．

sgv <- rowRanges(sgvc_pred) sgvc <- getSGVariantCounts(sgv, sample_info = si) sgvc

##类:SGVariantCounts ## dim: 2 8 ##元数据(0):## assays(6): countsVariant5p countsVariant3p…countsVariant5pOr3p ## variantFreq ## rownames: NULL ## rowData names(20): from to…## colnames(8): N1 N2…T3 T4 ## colData names(6): sample_name file_bam…frag_length lib_size

执行差异测试需要每个变量计数，每个变量的唯一标识符，以及指示变量如何按事件分组的变量。这三种方法都可以从SGVariantCounts对象。

x <- counts(sgvc) vid <- variantID(sgvc) eid <- eventID(sgvc)

将这三个对象类比为每外显子计数、外显子标识符和基因标识符，它们可以用来构建一个DEXSeqDataSet对象，用于DEXSeq或者作为函数的输入diffSplice ()结合轰()用于limma．

12高级用法

功能analyzeFeatures ()而且analyzeVariants ()包装多个分析步骤方便。或者，可以直接调用执行单独步骤的函数。例如，基于的分析新创可以进行如下预测。

txf <- predictTxFeatures(si, gr) sgf <- convertToSGFeatures(txf) sgf <- annotate(sgf, txf_ucsc) sgfc <- getsgfeatucounts (si, sgf) sgv <- findsgvariables (sgf) sgvc <- getSGVariantCounts(sgv, sgfc)

predictTxFeatures ()预测每个样本的特征，合并样本之间的特征，并对预测的末端外显子进行过滤和处理。predictTxFeatures ()而且getSGFeatureCounts ()也可以在单个示例上运行(例如在高性能计算集群上分布)。当使用predictTxFeatures ()对于单独的样本，其预测将在稍后合并，运行predictTxFeatures ()与参数min_overhang = NULL抑制末端外显子的加工。然后，预测随后可以与函数合并和处理mergeTxFeatures ()而且processTerminalExons (),分别。

13练习

练习1在SGSeq，已知的转录本可以指定为GRangeList的每一项GRangesList对应一个副本。(a)构造一个GRangesList一个有三个外显子的转录本。(b)转换GRangesList到一个TxFeatures对象。(c)重复(a)和(b)以创建与第一转录本中的外显子共享或重叠的第二转录本。(d)将两份抄本合并为一份GRangesList把它转换成aTxFeatures对象。输出如何变化?(e)转换TxFeatures对象从(d)转换为SGFeatures对象。如何TxFeatures对象与TxFeatures对象?(f)使用plotSpliceGraph ()函数可视化与两个转录本对应的拼接图。

tx_1 <- GRangesList(tx_1 = GRanges("1"， IRanges(c(101, 301, 501)， c(200, 400, 600))， "+") txf_1 <- convertToTxFeatures(tx_1) tx_2 <- GRangesList(tx_2 = GRanges("1"， IRanges(c(101, 351, 701)， c(200, 400, 800))， "+") txf_2 <- converttoxfeatures (tx_2) txf <- converttoxfeatures (c(tx_1, tx_2)) sgf <- convertToSGFeatures(txf) par(mfrow = c(1,1)) plotSpliceGraph(sgf)

其余的练习是基于从16个正常人体组织的配对端RNA-seq数据中获得的全基因组预测，这些数据是Illumina Body Map 2.0的一部分。数据的处理如下面的代码块所示(该代码将无法运行，因为BAM文件不能作为研讨会的一部分)。

sgvc <- analyzeFeatures(si, alpha = 2, beta = 0.2, gamma = 0.2) sgvc <- analyzvariants (sgfc, min_分母= 10)

加载先前生成的对象sgfc而且sgvc，并恢复seqlevels为TxDb对象(因此所有染色体都可用)。请注意sgfc而且sgvc已经使用UCSC knownGene表中的转录本特征进行了注释。

data(sgfc, sgvc, package = "SGSeqBioC2016") txdb <- restoreSeqlevels(txdb) seqlevelsStyle(txdb) <- "NCBI" txdb <- keepSeqlevels(txdb, c(1:22， "X"， "Y"))

练习2(a)使用功能plotSpliceGraph ()绘制基因的基因模型KIFAP3(Entrez ID 22920)。(b)预测基因模型的哪些部分是没有注释的?(c)使用功能plotFeatures ()检查体图组织中剪接连接或外显子的表达水平。

sgf <- rowRanges(sgfc) plotSpliceGraph(sgf, geneName = "22920") plotSpliceGraph(sgf, geneName = "22920"， color_novel = "red") plotFeatures(sgfc, geneName = "22920"， color_novel = "red") plotFeatures(sgfc, geneName = "22920"， color_novel = "red"， include = "exons")

练习3(a)检查中拼接事件的细节KIFAP3基因。(b)绘制未注释的拼接事件。

mcols(sgvc)[any(geneName(sgvc) == "22920")，] plotvariant (sgvc, eventID = 2872, color_novel = "red")

练习4(a)从SGVariantCounts目的，提取出与新型盒式外显子相对应的变体KIFAP3基因。(b)评估这一变种对KIFAP3编码蛋白质的潜力。

variant <- rowRanges(sgvc)[variantID(sgvc) == 6796] variant <- keepSeqlevels(variant, c(1:22， "X"， "Y")) vep <- predictvarianteeffects (variant, txdb, Hsapiens) vep

练习5按照练习2-4中概述的方法，分析基因ABCD3(entrezid5825)，ESYT2(entrezid57488)，ROCK2(Entrez ID 9475)或其他感兴趣的基因。

14会话信息

sessionInfo ()

## R版本3.3.0(2016-05-03)##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 14.04.2 LTS下## ## 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] stats4 parallel stats graphics grDevices utils datasets ## [8] methods base ## ## other attached packages: ## [1] BSgenome.Hsapiens.UCSC.hg19_1.4.0 ## [2] BSgenome_1.41.2 ## [3] rtracklayer_1.33.7 ## [4] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 ## [5] GenomicFeatures_1.25.12 ## [6] AnnotationDbi_1.35.3 ## [7] SGSeq_1.7.6 ## [8] SummarizedExperiment_1.3.5 ## [9] Biobase_2.33.0 ## [10] Rsamtools_1.25.0 ## [11] Biostrings_2.41.4 ## [12] XVector_0.13.2 ## [13] GenomicRanges_1.25.4 ## [14] GenomeInfoDb_1.9.1 ## [15] IRanges_2.7.6 ## [16] S4Vectors_0.11.5 ## [17] BiocGenerics_0.19.1 ## [18] knitr_1.13 ## [19] BiocStyle_2.1.8 ## ## loaded via a namespace (and not attached): ## [1] igraph_1.0.1 Rcpp_0.12.5 ## [3] magrittr_1.5 zlibbioc_1.19.0 ## [5] GenomicAlignments_1.9.3 BiocParallel_1.7.4 ## [7] stringr_1.0.0 tools_3.3.0 ## [9] DBI_0.4-1 htmltools_0.3.5 ## [11] yaml_2.1.13 digest_0.6.9 ## [13] formatR_1.4 bitops_1.0-6 ## [15] RUnit_0.4.31 biomaRt_2.29.2 ## [17] RCurl_1.95-4.8 RSQLite_1.0.0 ## [19] evaluate_0.9 rmarkdown_0.9.6 ## [21] stringi_1.1.1 XML_3.98-1.4