Package version: SGSeq 1.4.3
Contents
Preliminaries
library(SGSeq)
This vignette illustrates an analysis of paired-end RNA-seq data from four tumor and four normal colorectal samples, which are part of a data set published in (Seshagiri et al. 2012). For the purpose of this vignette, we created BAM files that only include reads mapping to a single gene of interest (FBXO31).
When starting a new project, SGSeq requires information about the samples to be analyzed. This information can be provided as a data.frame, which must include columns sample_name (specificying a unique name for each sample) and file_bam (specifying the location of the BAM file). Function getBamInfo can be used to extract additional required information from BAM files, including paired-end status, median read length, median insert size and the total number of aligned reads. This information must be obtained once initially and can then be used for all subsequent analyses. It is essential that BAM files are generated using a splice-aware alignment program that generates the custom tag ‘XS’ indicating the direction of transcription for spliced reads. In the following, we work with a data.frame si previously generated from the original (complete) BAM files with function getBamInfo.
si
## sample_name file_bam paired_end read_length frag_length lib_size
## 1 N1 N1.bam TRUE 75 293 12405197
## 2 N2 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
The following code block sets the correct BAM file paths in the sample information for this vignette.
path <- system.file("extdata", package = "SGSeq")
si$file_bam <- file.path(path, "bams", si$file_bam)
Transcript features and the TxFeatures class
We use the UCSC knownGene table as reference annotation, which is available as a Bioconductor annotation package TxDb.Hsapiens.UCSC.hg19.knownGene. We retain transcripts on chromosome 16, where the FBXO31 gene is located, and change chromosome names in the annotation to match chromosome names in the BAM files.
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
txdb <- keepSeqlevels(txdb, "chr16")
seqlevelsStyle(txdb) <- "NCBI"
To work with transcript annotation in the SGSeq framework, we first extract exons and splice junctions from the TxDb object using function convertToTxFeatures. We only retain features overlapping the FBXO31 gene (genomic coordinates of the FBXO31 gene are stored in the GRanges object gr).
txf_ucsc <- convertToTxFeatures(txdb)
## Warning in convertToTxFeatures(txdb): Merged adjacent exons
txf_ucsc <- txf_ucsc[txf_ucsc %over% gr]
txf_ucsc
## TxFeatures object with 23 ranges and 0 metadata columns:
## seqnames ranges strand type
## <Rle> <IRanges> <Rle> <factor>
## [1] 16 [87362942, 87365116] - L
## [2] 16 [87365116, 87367492] - J
## [3] 16 [87367492, 87367892] - I
## [4] 16 [87367892, 87368910] - J
## [5] 16 [87368910, 87369063] - I
## ... ... ... ... ...
## [19] 16 [87417011, 87417394] - F
## [20] 16 [87417628, 87417700] - U
## [21] 16 [87423343, 87423454] - I
## [22] 16 [87423454, 87425689] - J
## [23] 16 [87425689, 87425708] - F
## txName geneName
## <CharacterList> <CharacterList>
## [1] uc002fjv.3,uc002fjw.3,uc010vot.2 79791
## [2] uc002fjv.3,uc002fjw.3,uc010vot.2 79791
## [3] uc002fjv.3,uc002fjw.3,uc010vot.2 79791
## [4] uc002fjv.3,uc002fjw.3,uc010vot.2 79791
## [5] uc002fjv.3,uc002fjw.3,uc010vot.2 79791
## ... ... ...
## [19] uc002fjw.3 79791
## [20] uc021tmi.1
## [21] uc010vot.2 79791
## [22] uc010vot.2 79791
## [23] uc010vot.2 79791
## -------
## seqinfo: 1 sequence from hg19 genome
SGSeq makes extensive use of the Bioconductor infrastructure for genomic ranges (Lawrence et al. 2013). The TxFeatures class shown above extends the GRanges class with additional columns. Column type can take values
- J (splice junction)
- I (internal exon)
- F (first/5\(^\prime\)-terminal exon)
- L (last/3\(^\prime\)-terminal exon)
- U (unspliced).
Columns txName and geneName indicate the transcript and the gene that each feature derives from. Note that a feature can belong to more than one transcript. Accordingly these columns can store multiple values for each feature.
columns can be accessed using accessor functions named after the columns they access (e.g. use function type to obtain feature type).
If transcript annotation is not available as a TxDb object, function importTranscripts can be used to import annotation in GFF format. Function convertToTxFeatures can construct TxFeatures from a GRangesList of exons grouped by transcript.
Splice graph features and the SGFeatures class
Exons stored as TxFeatures can be overlapping (e.g. due to alternative splice sites). Overlapping exons can result in ambiguities when attempting to assign reads to individual exons. We therefore partition exonic regions into disjoint exon bins. Splice junctions and disjoint exon bins uniquely determine a genome-wide splice graph (Heber et al. 2002). To store splice graph features, SGSeq implements the SGFeatures class.
sgf_ucsc <- convertToSGFeatures(txf_ucsc)
sgf_ucsc
## SGFeatures object with 42 ranges and 0 metadata columns:
## seqnames ranges strand type splice5p splice3p
## <Rle> <IRanges> <Rle> <factor> <logical> <logical>
## [1] 16 [87362942, 87365116] - E TRUE FALSE
## [2] 16 [87365116, 87365116] - A <NA> <NA>
## [3] 16 [87365116, 87367492] - J <NA> <NA>
## [4] 16 [87367492, 87367492] - D <NA> <NA>
## [5] 16 [87367492, 87367892] - E TRUE TRUE
## ... ... ... ... ... ... ...
## [38] 16 [87423343, 87423454] - E TRUE TRUE
## [39] 16 [87423454, 87423454] - A <NA> <NA>
## [40] 16 [87423454, 87425689] - J <NA> <NA>
## [41] 16 [87425689, 87425689] - D <NA> <NA>
## [42] 16 [87425689, 87425708] - E FALSE TRUE
## featureID geneID txName
## <integer> <integer> <CharacterList>
## [1] 1 1 uc002fjv.3,uc002fjw.3,uc010vot.2
## [2] 2 1 uc002fjv.3,uc002fjw.3,uc010vot.2
## [3] 3 1 uc002fjv.3,uc002fjw.3,uc010vot.2
## [4] 4 1 uc002fjv.3,uc002fjw.3,uc010vot.2
## [5] 5 1 uc002fjv.3,uc002fjw.3,uc010vot.2
## ... ... ... ...
## [38] 38 1 uc010vot.2
## [39] 39 1 uc010vot.2
## [40] 40 1 uc010vot.2
## [41] 41 1 uc010vot.2
## [42] 42 1 uc010vot.2
## geneName
## <CharacterList>
## [1] 79791
## [2] 79791
## [3] 79791
## [4] 79791
## [5] 79791
## ... ...
## [38] 79791
## [39] 79791
## [40] 79791
## [41] 79791
## [42] 79791
## -------
## seqinfo: 1 sequence from hg19 genome
Similar to TxFeatures, SGFeatures extends the GRanges class with additional columns. Column type for an SGFeatures object takes values
- J (splice junction)
- E (disjoint exon bin)
- D (splice donor site)
- A (splice acceptor site).
By convention, splice donor and acceptor sites correspond to exonic positions immediately upstream and downstream of the intron, respectively. Note that splice sites are redundant in the sense that they are determined by the splice junctions included in the SGFeatures object. When assigning read counts to each feature (see below), counts for exons and splice junctions are based on structurally compatible reads. In the case of splice donor and acceptor sites, counts indicate the number of reads that extend across the spliced boundary (i.e. overlapping the splice site, as well as the flanking intronic position). Splice sites are included in the SGFeatures object as their counts are subsequently used for splice variant quantification.
SGFeatures includes additional columns not included in TxFeatures. spliced5p and spliced3p indicate whether exon bins have a mandatory splice at the 5\(^\prime\) and 3\(^\prime\) boundary, respectively. This information is used to determine whether a read is structurally compatible with an exon bin, as well as to determine whether an exon bin is consistent with an annotated transcript.
Column featureID provides a unique identifier for each feature, while columnn geneID indicates the unique connected component of the splice graph a feature belongs to.
Both TxFeatures and SGFeatures objects can be exported to BED files using function exportFeatures.
Analysis based on annotated transcripts
We can now start analyzing the RNA-seq data at the FBXO31 gene locus. We first perform an analysis based on annotated transcripts. The following example converts the transcript features into splice graph features and obtains counts of compatible RNA-seq reads for each feature and each sample.
sgfc_ucsc <- analyzeFeatures(si, features = txf_ucsc)
## Process features...
## Obtain counts...
## N1 complete.
## N2 complete.
## N3 complete.
## N4 complete.
## T1 complete.
## T2 complete.
## T3 complete.
## T4 complete.
sgfc_ucsc
## class: SGFeatureCounts
## dim: 42 8
## metadata(0):
## assays(2): counts FPKM
## rownames: NULL
## rowRanges metadata column names(0):
## colnames(8): N1 N2 ... T3 T4
## colData names(6): sample_name file_bam ... frag_length lib_size
analyzeFeatures returns an object of class SGFeatureCounts, which extends the RangedSummarizedExperiment class from the SummarizedExperiment package. SGFeatureCounts contains the sample information as colData, splice graph features as rowRanges and assays counts and FPKM, which store structurally compatible counts and FPKMs, respectively. Accessor functions colData, rowRanges, counts and FPKM can be used to access the data. For example, counts and FPKMs can be extracted from an SGFeatureCounts object as shown below.
head(counts(sgfc_ucsc))
## N1 N2 N3 N4 T1 T2 T3 T4
## [1,] 248 424 385 477 398 459 325 394
## [2,] 0 0 2 1 1 2 1 1
## [3,] 4 10 8 13 7 19 6 13
## [4,] 0 1 1 0 0 0 0 0
## [5,] 28 41 40 48 33 81 33 63
## [6,] 0 0 0 0 0 0 0 0
head(FPKM(sgfc_ucsc))
## N1 N2 N3 N4 T1 T2
## [1,] 8.103616 13.6611958 10.796489 12.6842394 11.2868069 12.3210956
## [2,] 0.000000 0.0000000 0.901918 0.4278029 0.4709862 0.8738597
## [3,] 2.178686 5.1616993 3.607672 5.5614378 3.2969034 8.3016673
## [4,] 0.000000 0.5161699 0.450959 0.0000000 0.0000000 0.0000000
## [5,] 3.257025 5.2464307 4.405408 5.0067741 3.3630067 8.2486853
## [6,] 0.000000 0.0000000 0.000000 0.0000000 0.0000000 0.0000000
## T3 T4
## [1,] 8.6258775 10.2947244
## [2,] 0.4363152 0.4274168
## [3,] 2.6178911 5.5564183
## [4,] 0.0000000 0.0000000
## [5,] 3.2336318 6.1595582
## [6,] 0.0000000 0.0000000
Compatible FPKMs for splice graph features can be visualized with function plotFeatures. plotFeatures generates a two-panel figure with a splice graph shown in the top panel and a heatmap illustrating expression levels of individual features in the bottom panel. For customization of plotFeatures output, see section Visualization. The plotting function invisibly returns a data.frame with information on splice graph features, including genomic coordinates.
df <- plotFeatures(sgfc_ucsc, geneID = 1)
df
## id name type featureID color
## 1 E1 E:16:87425689-87425708:- E 42 black
## 2 E2 E:16:87423343-87423454:- E 38 black
## 3 E3 E:16:87417011-87417394:- E 35 black
## 4 E4 E:16:87393973-87394561:- E 33 black
## 5 E5 E:16:87393901-87393972:- E 29 black
## 6 E6 E:16:87380780-87380856:- E 25 black
## 7 E7 E:16:87377204-87377371:- E 21 black
## 8 E8 E:16:87376483-87376557:- E 17 black
## 9 E9 E:16:87369761-87369870:- E 13 black
## 10 E10 E:16:87368910-87369063:- E 9 black
## 11 E11 E:16:87367492-87367892:- E 5 black
## 12 E12 E:16:87362942-87365116:- E 1 black
## 13 J1 J:16:87423454-87425689:- J 40 black
## 14 J2 J:16:87393972-87423343:- J 32 black
## 15 J3 J:16:87393972-87417011:- J 31 black
## 16 J4 J:16:87380856-87393901:- J 27 black
## 17 J5 J:16:87377371-87380780:- J 23 black
## 18 J6 J:16:87376557-87377204:- J 19 black
## 19 J7 J:16:87369870-87376483:- J 15 black
## 20 J8 J:16:87369063-87369761:- J 11 black
## 21 J9 J:16:87367892-87368910:- J 7 black
## 22 J10 J:16:87365116-87367492:- J 3 black
Note that the splice graph derived from annotated transcripts includes three alternative transcript start sites (TSSs). However, the heatmap indicates that the first TSS is not used in the samples in our data set.
Analysis based on de novo prediction
Instead of relying on existing annotation, SGSeq can predict features from BAM files directly. The following code block predicts splice graph features with read evidence in our data set.
sgfc_pred <- analyzeFeatures(si, which = gr)
## Predict features...
## N1 complete.
## N2 complete.
## N3 complete.
## N4 complete.
## T1 complete.
## T2 complete.
## T3 complete.
## T4 complete.
## Process features...
## Obtain counts...
## N1 complete.
## N2 complete.
## N3 complete.
## N4 complete.
## T1 complete.
## T2 complete.
## T3 complete.
## T4 complete.
For interpretability, we annotate predicted features with respect to transcripts included in the UCSC knownGene table. The annotate function assigns compatible transcripts to each feature and stores them in column txName. column geneName behaves transitively, meaning all features belonging to the same connected component of the splice graph (with identical geneID) have the same value for geneName. This behavior makes it easy to identify unannotated features (with empty txName) that belong to an annotated gene (non-empty geneName).
sgfc_pred <- annotate(sgfc_pred, txf_ucsc)
Predicted splice graph features and compatible FPKMs can be visualized as previously. Splice graph features with missing annotation can be highlighted using argument color_novel.
df <- plotFeatures(sgfc_pred, geneID = 1, color_novel = "red")
df
## id name type featureID color
## 1 E1 E:16:87417011-87417348:- E 38 black
## 2 E2 E:16:87393901-87393972:- E 34 black
## 3 E3 E:16:87392017-87392103:- E 30 black
## 4 E4 E:16:87380780-87380856:- E 25 black
## 5 E5 E:16:87377204-87377371:- E 21 black
## 6 E6 E:16:87376483-87376557:- E 17 black
## 7 E7 E:16:87369761-87369870:- E 13 black
## 8 E8 E:16:87368910-87369063:- E 9 black
## 9 E9 E:16:87367492-87367892:- E 5 black
## 10 E10 E:16:87362930-87365116:- E 1 black
## 11 J1 J:16:87393972-87417011:- J 36 black
## 12 J2 J:16:87392103-87393901:- J 32 black
## 13 J3 J:16:87380856-87393901:- J 28 black
## 14 J4 J:16:87380856-87392017:- J 27 black
## 15 J5 J:16:87377371-87380780:- J 23 black
## 16 J6 J:16:87376557-87377204:- J 19 black
## 17 J7 J:16:87369870-87376483:- J 15 black
## 18 J8 J:16:87369063-87369761:- J 11 black
## 19 J9 J:16:87367892-87368910:- J 7 black
## 20 J10 J:16:87365116-87367492:- J 3 black
Note that most exons and splice junctions predicted from the RNA-seq data are consistent with transcripts in the UCSC knownGene table (shown in gray). However, in contrast to the previous figure, the predicted gene model does not include parts of the splice graph that are not expressed in our data set. Also, an unannotated exon (E3, shown in red) was discovered from the RNA-seq data, which is expressed in three of the four normal colorectal samples (N2, N3, N4).
Analysis of predicted splice variants
Instead of considering the complete splice graph of a gene, we can focus the analysis on individual splice events. In the SGSeq framework, the splice graph is a directed acyclic graph with nodes corresponding to transcript starts, ends and splice sites, and edges corresponding to disjoint exon bins and splice junctions, directed from 5\(^\prime\) to the 3\(^\prime\) end. A splice event is defined by a start node and an end node connected by two or more paths and no intervening nodes with all paths intersecting. SGSeq identifies splice events recursively from the graph, and estimates relative usage of splice variants based on compatible reads spanning the event boundaries. The following example identifies splice events from the splice graph and obtains representative counts for each splice variant.
sgvc_pred <- analyzeVariants(sgfc_pred)
## Find segments...
## Find variants...
## Annotate variants...
sgvc_pred
## class: SGVariantCounts
## dim: 2 8
## metadata(0):
## assays(5): countsVariant5p countsTotal5p countsVariant3p
## countsTotal3p variantFreq
## rownames(2): 1 2
## rowRanges metadata column names(16): from to ... variantType
## variantName
## colnames(8): N1 N2 ... T3 T4
## colData names(6): sample_name file_bam ... frag_length lib_size
analyzeVariants returns an SGVariantCounts object. Similar to SGFeatureCounts, SGVariantCounts extends the RangedSummarizedExperiment class. SGVariantCounts contains sample information as colData and SGVariants as rowRanges. Assay variantFreq stores estimates of relative usage for each splice variant and sample. Accessor functions colData, rowRanges and variantFreq can be used to access the data. For example, estimates of relative usage can be extracted from an SGVariantCounts object as shown below.
variantFreq(sgvc_pred)
## N1 N2 N3 N4 T1 T2 T3
## 1 0.88 0.5769231 0.6153846 0.5925926 0.93333333 0.90196078 0.8484848
## 2 0.12 0.4230769 0.3846154 0.4074074 0.06666667 0.09803922 0.1515152
## T4
## 1 0.95833333
## 2 0.04166667
Information on splice variants is stored in metadata columns and can be accessed as follows.
mcols(sgvc_pred)
## DataFrame with 2 rows and 16 columns
## from to type featureID segmentID
## <character> <character> <character> <character> <character>
## 1 D:16:87393901:- A:16:87380856:- J 28 4
## 2 D:16:87393901:- A:16:87380856:- JEJ 32,30,27 2
## closed3p closed5p geneID eventID variantID featureID5p
## <logical> <logical> <integer> <integer> <integer> <IntegerList>
## 1 TRUE TRUE 1 1 1 28
## 2 TRUE TRUE 1 1 2 32
## featureID3p txName geneName
## <IntegerList> <CharacterList> <CharacterList>
## 1 28 uc002fjv.3,uc002fjw.3,uc010vot.2 79791
## 2 27 79791
## variantType variantName
## <CharacterList> <character>
## 1 SE:S 79791_1_1/2_SE
## 2 SE:I 79791_1_2/2_SE
- from and to indicate the variant start and end node, respectively, from nodes are splice donor sites (D) or transcript starts (S), while to nodes are splice acceptor sites (A) or transcript ends (E)
- columns type and featureID describe the variant in terms of the splice graph features that make up the variant
- segmentID are unique identifiers labelling unbranched segments of the splice graph (internal use only)
- closed5p indicates whether the nodes belonging to a splice variant can be reached from nodes outside of the splice event exclusively through the from node
- closed3p indicates whether the nodes belonging to a splice variant can reach nodes outside of the splice event exclusively through the to node
- eventID and variantID are unique identifiers for each splice event and splice variant, respectively
- featureID5p and featureID3p indicate representative features that are used for obtaining relative usage estimates at the variant start and end, respectively
- variantType indicates whether a splice variant is consistent with a canonical splice event (for a list of possible values, see the manual page for annotateSGVariants)
- variantName provides a unique identifier for each splice variant that is intended to be more human-readable than the numeric identifier stored in variantID (for details, see the manual page for makeVariantNames)
Splice variants and estimates of relative usage can be visualized with function plotVariants.
plotVariants(sgvc_pred, eventID = 1, color_novel = "red")
plotVariants generates a two-panel figure similar to plotFeatures. The splice graph in the top panel illustrates the selected splice event. In this example, the splice event consists of two splice variants that correspond to a skip or inclusion of the unannotated exon. The heatmap illustrates estimates of relative usage for each splice variant. We observe that samples N2, N3 and N4 show evidence for both transcripts that include the exon as well as transcripts that skip the exon. The remaining samples show little evidence for exon inclusion.
Visualization
Functions plotFeatures and plotVariants support many options for customizing figures. Note that the splice graph in the top figure panel is plotted by function plotSpliceGraph, which can be called directly.
plotFeatures includes multiple alternative arguments for selecting features to be displayed. The following code illustrates three different ways of selecting and plotting the splice graph and expression levels for FBXO31 (Entrez ID 79791).
plotFeatures(sgfc_pred, geneID = 1)
plotFeatures(sgfc_pred, geneName = "79791")
plotFeatures(sgfc_pred, which = gr)
By default, the heatmap generated by plotFeatures displays splice junctions. Alternatively, exon bins, or both exon bins and splice junctions can be displayed.
plotFeatures(sgfc_pred, geneID = 1, include = "junctions")
plotFeatures(sgfc_pred, geneID = 1, include = "exons")
plotFeatures(sgfc_pred, geneID = 1, include = "both")
Argument toscale controls which parts of the gene model are drawn to scale.
plotFeatures(sgfc_pred, geneID = 1, toscale = "gene")
plotFeatures(sgfc_pred, geneID = 1, toscale = "exon")
plotFeatures(sgfc_pred, geneID = 1, toscale = "none")
Heatmaps allow the visualization of expression values summarized for splice junctions and exon bins. Alternatively, per-base read coverages and splice junction counts can be visualized with function 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")
}
Advanced use
Functions analyzeFeatures and analyzeVariants wrap multiple analysis steps for convenience. Alternatively, the functions performing individual steps can be called directly. For example, the previous analysis using de novo prediction can be performed as follows.
txf <- predictTxFeatures(si, gr)
## N1 complete.
## N2 complete.
## N3 complete.
## N4 complete.
## T1 complete.
## T2 complete.
## T3 complete.
## T4 complete.
sgf <- convertToSGFeatures(txf)
sgf <- annotate(sgf, txf_ucsc)
sgfc <- getSGFeatureCounts(si, sgf)
## N1 complete.
## N2 complete.
## N3 complete.
## N4 complete.
## T1 complete.
## T2 complete.
## T3 complete.
## T4 complete.
sgv <- findSGVariants(sgf)
## Find segments...
## Find variants...
## Annotate variants...
sgvc <- getSGVariantCounts(sgv, sgfc)
predictTxFeatures and getSGFeatureCounts can be run on individual samples (e.g. for distribution across a high-performance computing cluster). predictTxFeatures predicts features for each sample, merges features across samples and finally performs filtering and processing of predicted terminal exons. When using predictTxFeatures for individual samples, with predictions intended to be merged at a later point in time, run predictTxFeatures with argument min_overhang = NULL to suppress processing of terminal exons. Then predictions can subsequently be merged and processed with functions mergeTxFeatures and processTerminalExons, respectively.