In pkimes/spliceclust: Visualization and Clustering of Spliced Components

spliceclust

Introduction
SpliceGraHM Examples
SplicePCA Examples
SplicePCP Examples
SpliceGraHM2 Examples

Introduction

This package may be used to plot exon and splice junction coverage across a large cohort of RNA-seq samples. The plotting approach is based on the idea of transposing expression heatmaps on splicing diagrams. The plots are generated using the ggplot2 and ggbio packages.

## load annotations packages
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
library("BSgenome.Hsapiens.UCSC.hg19")

## load current package
devtools::load_all()

SpliceGraHM Examples

First we must construct a concomp object from a GRangesList object with exon and junction boundaries and corresponding coverage values. For illustration purposes, we use a previously constructed GRanges object, klk12_lusc_gr, containing annotation information for the KLK12 gene locus along with coverage for 177 lung squamous cell carcinoma samples.

data("klk12_lusc_gr")

The exon and junction boundaries, and the corresponding kind label are given by:

ranges(klk12_lusc_gr)
mcols(klk12_lusc_gr)$kind

In addition to kind, the klk12_lusc_gr metadata columns also include the gene name (gIdx), gene boundaries (gStart, gStop), and coverage values for the 177 samples.

mcols(klk12_lusc_gr)[1:5, 1:6]

To construct the concomp (connected component) object for analysis, we first convert the GRanges object into a GRangesList object of length 2, corresponding to exon and junction information.

klk12_gl <- split(klk12_lusc_gr, mcols(klk12_lusc_gr)$kind)
klk12_cc <- concomp(klk12_gl)

We first demonstrate the default and basic SpliceGraHM (Splice Graph Heat Map) plotting procedure.

splicegrahm(klk12_cc)

In the plot above, each box arranged horizontally corresponds to a contiguous exonic region along the genome. Each box is colored by 177 horizontal lines, showing the expression level for the 177 samples being analyzed. Note that the exons are plotted along genomic coordinates, and log expression is shown in color.

In addition to the boxes, the plot contains arrows which correspond to a splicing events with sufficient support in the data (e.g. at least 8 samples, each with at least 5 reads spanning the splice junction). The arrows are colored such that darker arrows were present in a higher proportion of the 177 samples (scale shown on right).

The SpliceGraHM name is in reference to the fact that after removing the spacing between each exon the plot simply reduces to a standard heatmap of expression along a single gene. The corresponding figure after removing the intronic gaps is shown below. The rectangle of color is a heatmap with rows and columns corresponding to samples and exons, respectively.

splicegrahm(klk12_cc, genomic = FALSE, ex_use = 1, log_base = 2)

It is possible to show the coverage of each splice junction using a similar convention with rows corresponding to samples, and with color being used for coverage.

splicegrahm(klk12_cc, genomic = FALSE, j_incl = TRUE)

If annotation information is available, e.g. from UCSC KnownGenes, these can be passed to the function to add an additional track to the plot. The appropriate GRangesList object to be passed is illustrated in the following example.

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
isActiveSeq(txdb)[seqlevels(txdb)] <- FALSE
isActiveSeq(txdb)[paste0("chr", 1:22)] <- TRUE
exbytx <- exonsBy(txdb, "tx")

splicegrahm(klk12_cc, genomic = TRUE, j_incl = TRUE, txlist = exbytx)

If gene names are desired, the following can be used to match the transcript ID in txdb against gene symbols (e.g. in org.Hs.eg.db).

suppressPackageStartupMessages(library("org.Hs.eg.db"))

splicegrahm(klk12_cc, genomic = TRUE, j_incl = TRUE, txlist = exbytx,
            txdb = txdb, orgdb = org.Hs.eg.db)

The splicegrahm function can now also plot gene models in non-genomic space with an additional parameter eps. The eps parameter determines how far up/down from the connected component to lookfor overlapping gene models. If eps = NULL, all overlapping gene models are included. If eps = 1000, only overlapping gene models which are fully contained within 1000bp of the connected component are included.

splicegrahm(klk12_cc, genomic = FALSE, j_incl = TRUE, txlist = exbytx,
            txdb = txdb, orgdb = org.Hs.eg.db, eps = NULL)

splicegrahm(klk12_cc, genomic = FALSE, j_incl = TRUE, txlist = exbytx,
            txdb = txdb, orgdb = org.Hs.eg.db, eps = 1)

splicegrahm(klk12_cc, genomic = FALSE, j_incl = TRUE, txlist = exbytx,
            txdb = txdb, orgdb = org.Hs.eg.db, eps = 1000)

SplicePCA Examples

Principal Component Analysis (PCA) is a popular exploratory analysis tool for studying low rank structure in high-dimensional datasets by projecting the data along the directions of greatest variation. These directions of greatest variation are often referred to as the PC "loading" directions. The splicepca function is written to visualize the loading vectors of splicing data.

As an example, the following code can be used to visualize the first 3 PC loadings in the KLK12 dataset from above.

splicepca(klk12_cc, npc = 3)

In the above PC loadings plot, red and blue are used to denote the magnitude of the PC loadings for each exon and splice junction. From the plot above, we see that the greatest source of variation in expression at the KLK12 gene locus is the overall expression of the gene. However, the second PC loading shows interesting behavior where the primary source of variation is attributable to the expression of the central exon. Scrolling up a little to the UCSC KnownGenes shown included above, we see that the presence or absence of the central exon actually corresponds to differing isoforms annotated to this gene. Note that above the PCs are computed separately for exons and junctions information. As such, we obtain separate PC "scores" for the exon and junction PC loadings.

splicepca(klk12_cc, npc = 3, scores = TRUE)

It is also possible to perform the PCA analysis using the concatenated exon and junction information by setting the pc_sep parameter to FALSE, and specifying the relative "weight" of each with ej_w. The exon and junction data are rescaled to each have sum-of-squares equal to the values specified by ej_w. In the following example, we use equal weights for the two data sources.

splicepca(klk12_cc, npc = 3, pc_sep = FALSE, ej_w = c(1, 1))
splicepca(klk12_cc, npc = 3, pc_sep = FALSE, ej_w = c(1, 1), scores = TRUE)

SplicePCP Examples

All above plots have used color to represent expression. However, low frequency or outlier events may become lost in this particular view. To handle this problem, we also provide the option to plot the splicing objects with vertical height (the y-axis) corresponding to expression. The name of the function is a reference to parallel coordinates plots. The same KLK12 data is shown below using the splicepcp function.

splicepcp(klk12_cc, genomic = TRUE, txlist = exbytx, txdb = txdb, orgdb = org.Hs.eg.db)

The plot includes 3 tracks: 1. log-expression values for each exon or region of an exon 2. the complete gene model and splice junctions 3. annotated transcripts from the UCSC KnownGene database.

Currently, the function is being rewritten to also include junction coverage using a parallel coordinates plot in a separate track.

SpliceGraHM2 Examples

To make direct comparison of two populations possible, we have created splicegrahm2 which draws two splicegrahm plots in a single figure. Consider, for example, the task of comparing the behavior of splicing between the LUSC samples and head and neck squamous cell carcinoma (HNSC) samples. We have loaded a connected component for 279 HNSC samples as klk12_hnsc_cc.

data("klk12_hnsc_cc")

splicegrahm2(klk12_cc, klk12_hnsc_cc, genomic=FALSE)

When gene models are also desired, they are placed in a central track for easy comparison to the two splicegrahms.

splicegrahm2(klk12_cc, klk12_hnsc_cc, genomic=FALSE,
             txlist = exbytx, txdb = txdb, orgdb = org.Hs.eg.db)

Alternatively, the bottom plot can be drawn with the same orientation as the top plot by setting mirror=FALSE.

splicegrahm2(klk12_cc, klk12_hnsc_cc, genomic=FALSE, mirror=FALSE)

While not modified above, the option same_scale can be used to specify whether the two plots should be drawn using the same or separate scale along the y-axis. This can be helpful when the two populations are of substantially different sizes.