Downloading and loading TCGA data

The quantification of each alternative splicing event is based on the proportion of junction reads that support the inclusion isoform, known as percent spliced-in or PSI (Wang et al., 2008).

To estimate this value for each splicing event, both alternative splicing annotation and junction quantification are required. While alternative splicing annotation is provided by the package, junction quantification may be retrieved from TCGA, GTEx, SRA or user-provided files.

Data is downloaded from Firebrowse, a service that hosts processed data from TCGA, as required to run the downstream analyses. Before downloading data, check the following options:

# Available tumour types
cohorts <- getFirebrowseCohorts()

# Available sample dates
date <- getFirebrowseDates()

# Available data types
dataTypes <- getFirebrowseDataTypes()

Note there is also the option for Gene expression (normalised by RSEM). However, we recommend to load the raw gene expression data instead, followed by filtering and normalisation as demonstrated afterwards.

After deciding on the options to use, download and load breast cancer data as follows:

# Set download folder
folder <- getDownloadsFolder()

# Download and load most recent junction quantification and clinical data from
# TCGA/Firebrowse for Breast Cancer
data <- loadFirebrowseData(folder=folder,
                           cohort="BRCA",
                           data=c("clinical", "junction_quantification",
                                  "RSEM_genes"),
                           date="2016-01-28")

# Select clinical and junction quantification dataset
clinical      <- data[[1]]$`Clinical data`
sampleInfo    <- data[[1]]$`Sample metadata`
junctionQuant <- data[[1]]$`Junction quantification (Illumina HiSeq)`
geneExpr      <- data[[1]]$`Gene expression`

Data is only downloaded if the files are not present in the given folder. In other words, if the files were already downloaded, the function will just load the files, so it is possible to reuse the code above just to load the requested files.

Windows limitations: If you are using Windows, note that the downloaded files have huge names that may be over Windows Maximum Path Length. A workaround would be to manually rename the downloaded files to have shorter names, move all downloaded files to a single folder and load such folder.

Filtering and normalising gene expression

As this package does not focuses on gene expression analysis, we suggest to read the RNA-seq section of limma’s user guide. Nevertheless, we present the following commands to quickly filter and normalise gene expression:

# Check genes where min. counts are available in at least N samples and filter
# out genes with mean expression and variance of 0
filter <- filterGeneExpr(geneExpr)
geneExprFiltered <- geneExpr[filter, ]

# Normalise gene expression and perform log2-transformation
geneExprNorm <- normaliseGeneExpression(geneExprFiltered)

Quantifying alternative splicing

After loading the clinical and alternative splicing junction quantification data from TCGA, quantify alternative splicing by clicking the green panel Alternative splicing quantification.

As previously mentioned, alternative splicing is quantified from the previously loaded junction quantification and an alternative splicing annotation file. To check current annotation files available:

# Available alternative splicing annotation
annotList <- listSplicingAnnotations()
annotList

##                        Human hg19/GRCh37 (2017-10-20) 
## "annotationHub_alternativeSplicingEvents.hg19_V2.rda" 
##                        Human hg19/GRCh37 (2016-10-11) 
##    "annotationHub_alternativeSplicingEvents.hg19.rda" 
##                               Human hg38 (2018-04-30) 
## "annotationHub_alternativeSplicingEvents.hg38_V2.rda"

Custom splicing annotation: Additional alternative splicing annotations can be prepared for psichomics by parsing the annotation from programs like VAST-TOOLS, MISO, SUPPA and rMATS. Note that SUPPA and rMATS are able to create their splicing annotation based on transcript annotation. Please read Preparing alternative splicing annotations.

To quantify alternative splicing, first select the junction quantification, alternative splicing annotation and alternative splicing event type(s) of interest:

# Load Human (hg19/GRCh37 assembly) annotation
hg19 <- listSplicingAnnotations(assembly="hg19")[[1]]
annotation <- loadAnnotation(hg19)

# Available alternative splicing event types (skipped exon, alternative 
# first/last exon, mutually exclusive exons, etc.)
getSplicingEventTypes()

##                                          Skipped exon 
##                                                  "SE" 
##                               Mutually exclusive exon 
##                                                 "MXE" 
##                            Alternative 5' splice site 
##                                                "A5SS" 
##                            Alternative 3' splice site 
##                                                "A3SS" 
##                                Alternative first exon 
##                                                 "AFE" 
##                                 Alternative last exon 
##                                                 "ALE" 
## Alternative first exon (exon-centred - less reliable) 
##                                            "AFE_exon" 
##  Alternative last exon (exon-centred - less reliable) 
##                                            "ALE_exon"

Afterwards, quantify alternative splicing using the previously defined parameters:

# Discard alternative splicing quantified using few reads
minReads <- 10 # default

psi <- quantifySplicing(annotation, junctionQuant, minReads=minReads)

# Check the identifier of the splicing events in the resulting table
events <- rownames(psi)
head(events)

## [1] "SE_3_+_13661331_13663275_13663415_13667945_FBLN2"   
## [2] "SE_3_+_57908750_57911572_57911661_57913023_SLMAP"   
## [3] "ALE_3_+_57908750_57911572_57913023_SLMAP"           
## [4] "SE_3_-_37136283_37133029_37132958_37125297_LRRFIP2" 
## [5] "SE_12_-_56558432_56558152_56558087_56557549_SMARCC2"
## [6] "AFE_4_+_56755098_56750094_56756389_EXOC1"

Note that the event identifier (for instance, SE_1_-_2125078_2124414_2124284_2121220_C1orf86) is composed of:

• Event type (SE stands for skipped exon)
• Chromosome (1)
• Strand (-)
• Relevant coordinates depending on event type (in this case, the first constitutive exon’s end, the alternative exon’ start and end and the second constitutive exon’s start)
• Associated gene (C1orf86)

Warning: all examples shown in this case study are performed using a small, yet representative subset of the available data. Therefore, values shown here may correspond to those when performing the whole analysis.

Data grouping

Let us create groups based on available samples types (i.e. Metastatic, Primary solid Tumor and Solid Tissue Normal) and tumour stages. As tumour stages are divided by sub-stages, we will merge sub-stages so as to have only tumour samples from stages I, II, III and IV (stage X samples are discarded as they are uncharacterised tumour samples).

# Group by normal and tumour samples
types  <- createGroupByAttribute("Sample types", sampleInfo)
normal <- types$`Solid Tissue Normal`
tumour <- types$`Primary solid Tumor`

# Group by tumour stage (I, II, III or IV) or normal samples
stages <- createGroupByAttribute(
    "patient.stage_event.pathologic_stage_tumor_stage", clinical)
groups <- list()
for (i in c("i", "ii", "iii", "iv")) {
    stage <- Reduce(union,
           stages[grep(sprintf("stage %s[a|b|c]{0,1}$", i), names(stages))])
    # Include only tumour samples
    stageTumour <- names(getSubjectFromSample(tumour, stage))
    elem <- list(stageTumour)
    names(elem) <- paste("Tumour Stage", toupper(i))
    groups <- c(groups, elem)
}
groups <- c(groups, Normal=list(normal))

# Prepare group colours (for consistency across downstream analyses)
colours <- c("#6D1F95", "#FF152C", "#00C7BA", "#FF964F", "#00C65A")
names(colours) <- names(groups)
attr(groups, "Colour") <- colours

# Prepare normal versus tumour stage I samples
normalVSstage1Tumour <- groups[c("Tumour Stage I", "Normal")]
attr(normalVSstage1Tumour, "Colour") <- attr(groups, "Colour")

# Prepare normal versus tumour samples
normalVStumour <- list(Normal=normal, Tumour=tumour)
attr(normalVStumour, "Colour") <- c(Normal="#00C65A", Tumour="#EFE35C")

Principal component analysis (PCA)

PCA is a technique to reduce data dimensionality by identifying variable combinations (called principal components) that explain the variance in the data (Ringnér, 2008). Use the following commands to perform PCA:

# PCA of PSI between normal and tumour stage I samples
psi_stage1Norm    <- psi[ , unlist(normalVSstage1Tumour)]
pcaPSI_stage1Norm <- performPCA(t(psi_stage1Norm))

As PCA cannot be performed on data with missing values, missing values need to be either removed (thus discarding data from whole splicing events or genes) or impute them (i.e. attributing to missing values the median of the non-missing ones). Use the argument missingValues within function performPCA to select the number of missing values that are tolerable per event (i.e. if a splicing event or gene has less than N missing values, those missing values will be imputed; otherwise, the event is discarded from PCA).

# Explained variance across principal components
plotVariance(pcaPSI_stage1Norm)

# Score plot (clinical individuals)
plotPCA(pcaPSI_stage1Norm, groups=normalVSstage1Tumour)

# Loading plot (variable contributions)
plotPCA(pcaPSI_stage1Norm, loadings=TRUE, individuals=FALSE)

# Table of variable contributions (as used to plot PCA, also)
table <- calculateLoadingsContribution(pcaPSI_stage1Norm)
knitr::kable(head(table, 5))

For performance reasons, the loading plot is able to exclusively render the top variables that most contribute to the select principal components by using the argument nLoadings within function plotPCA.

Hint: As most plots in psichomics, PCA plots can be zoomed-in by clicking-and-dragging within the plot (click Reset zoom to zoom-out). To toggle the visibility of the data series represented in the plot, click its respective name in the plot legend.

To perform PCA using alternative splicing quantification and gene expression data (both using all samples and only Tumour Stage I and Normal samples):

# PCA of PSI between all samples (coloured by tumour stage and normal samples)
pcaPSI_all <- performPCA(t(psi))
plotPCA(pcaPSI_all, groups=groups)
plotPCA(pcaPSI_all, loadings=TRUE, individuals=FALSE)

# PCA of gene expression between all samples (coloured by tumour stage and 
# normal samples)
pcaGE_all <- performPCA(t(geneExprNorm))
plotPCA(pcaGE_all, groups=groups)
plotPCA(pcaGE_all, loadings=TRUE, individuals=FALSE)

# PCA of gene expression between normal and tumour stage I samples
ge_stage1Norm    <- geneExprNorm[ , unlist(normalVSstage1Tumour)]
pcaGE_stage1Norm <- performPCA(t(ge_stage1Norm))
plotPCA(pcaGE_stage1Norm, groups=normalVSstage1Tumour)
plotPCA(pcaGE_stage1Norm, loadings=TRUE, individuals=FALSE)

NUMB exon 12 inclusion and correlation with QKI gene expression

One of the splicing events that most contribute the separation between tumour stage I and normal samples is NUMB exon 12 inclusion, whose protein is crucial for cell differentiation as a key regulator of the Notch pathway. The RNA-binding protein QKI has been shown to repress NUMB exon 12 inclusion in lung cancer cells by competing with core splicing factor SF1 for binding to the branch-point sequence, thereby repressing the Notch signalling pathway, which results in decreased cancer cell proliferation (Zong et al., 2014).

# Find the right event
ASevents <- rownames(psi)
(tmp     <- grep("NUMB", ASevents, value=TRUE))

## [1] "SE_14_-_73749067_73746132_73745989_73744001_NUMB"

NUMBskippedExon12 <- tmp[1]

# Plot the representation of NUMB exon 12 inclusion
plotSplicingEvent(NUMBskippedExon12)

## $`SE_14_-_73749067_73746132_73745989_73744001_NUMB`
## [1] "<svg height=\"50px\" width=\"262\">\n  <rect class=\"diagram\" x=\"1\" y=\"10\" width=\"60\" height=\"20\" style=\"stroke-width: 1.5px;; fill: lightgray; stroke: darkgray\"></rect>\n  <text class=\"diagram outside\" x=\"61\" y=\"8\" style=\"font-size: 10px; font-family: Helvetica, sans-serif; text-anchor: end; dominant-baseline: auto\">73749067</text>\n  <rect class=\"diagram\" x=\"76\" y=\"10\" width=\"110\" height=\"20\" style=\"stroke-width: 1.5px;; fill: #ffb153; stroke: #faa000\"></rect>\n  <text class=\"diagram outside\" x=\"76\" y=\"8\" style=\"font-size: 10px; font-family: Helvetica, sans-serif; text-anchor: start; dominant-baseline: auto\">73746132</text>\n  <text class=\"diagram outside\" x=\"186\" y=\"8\" style=\"font-size: 10px; font-family: Helvetica, sans-serif; text-anchor: end; dominant-baseline: auto\">73745989</text>\n  <text class=\"diagram\" x=\"131\" y=\"20\" style=\"font-size: 10px; font-family: Helvetica, sans-serif; text-anchor: middle; dominant-baseline: middle\">143 nts</text>\n  <rect class=\"diagram\" x=\"201\" y=\"10\" width=\"60\" height=\"20\" style=\"stroke-width: 1.5px;; fill: lightgray; stroke: darkgray\"></rect>\n  <text class=\"diagram outside\" x=\"201\" y=\"8\" style=\"font-size: 10px; font-family: Helvetica, sans-serif; text-anchor: start; dominant-baseline: auto\">73744001</text>\n  <path class=\"diagram\" style=\"fill: none; stroke: darkgray; stroke-width: 1.5px\" d=\"M 61 30 C 108 45 154 45 201 30\"></path>\n  <path class=\"diagram\" style=\"fill: none; stroke: #faa000; stroke-width: 1.5px\" d=\"M 61 30 C 66 45 71 45 76 30\"></path>\n  <path class=\"diagram\" style=\"fill: none; stroke: #faa000; stroke-width: 1.5px\" d=\"M 186 30 C 191 45 196 45 201 30\"></path>\n</svg>"
## attr(,"class")
## [1] "splicingEventPlot" "character"        
## 
## attr(,"class")
## [1] "splicingEventPlotList" "list"

# Plot its PSI distribution
plotDistribution(psi[NUMBskippedExon12, ], normalVStumour)

# Find the right gene
genes <- rownames(geneExprNorm)
(tmp  <- grep("QKI", genes, value=TRUE))

## [1] "QKI|9444"

QKI   <- tmp[1] # "QKI|9444"

# Plot its gene expression distribution
plotDistribution(geneExprNorm[QKI, ], normalVStumour, psi=FALSE)

plotCorrelation(correlateGEandAS(
    geneExprNorm, psi, QKI, NUMBskippedExon12, method="spearman"))

## $`SE_14_-_73749067_73746132_73745989_73744001_NUMB`
## $`SE_14_-_73749067_73746132_73745989_73744001_NUMB`$`QKI|9444`

Differential splicing analysis

To analyse alternative splicing between normal and tumour stage I samples:

diffSplicing <- diffAnalyses(psi, normalVSstage1Tumour)

# Filter based on |∆ Median PSI| > 0.1 and q-value < 0.01
deltaPSIthreshold <- abs(diffSplicing$`∆ Median`) > 0.1
pvalueThreshold   <- diffSplicing$`Wilcoxon p-value (BH adjusted)` < 0.01

# Plot results
library(ggplot2)
ggplot(diffSplicing, aes(`∆ Median`, 
                         -log10(`Wilcoxon p-value (BH adjusted)`))) +
    geom_point(data=diffSplicing[deltaPSIthreshold & pvalueThreshold, ],
               colour="orange", alpha=0.5, size=3) + 
    geom_point(data=diffSplicing[!deltaPSIthreshold | !pvalueThreshold, ],
               colour="gray", alpha=0.5, size=3) + 
    theme_light(16) +
    ylab("-log10(q-value)")

# Events already tested which have prognostic value
events <- c(
    "SE_9_+_6486925_6492303_6492401_6493826_UHRF2",
    "SE_4_-_87028376_87024397_87024339_87023185_MAPK10",
    "SE_2_+_152324660_152324988_152325065_152325155_RIF1",
    "SE_2_+_228205096_228217230_228217289_228220393_MFF",
    "MXE_15_+_63353138_63353397_63353472_63353912_63353987_63354414_TPM1",
    "SE_2_+_173362828_173366500_173366629_173368819_ITGA6",
    "SE_1_+_204957934_204971724_204971876_204978685_NFASC")

# Survival curves based on optimal PSI cutoff
library(survival)

# Assign alternative splicing quantification to patients based on their samples
samples <- colnames(psi)
match <- getSubjectFromSample(samples, clinical, sampleInfo=sampleInfo)

survPlots <- list()
for (event in events) {
    # Find optimal cutoff for the event
    eventPSI <- assignValuePerSubject(psi[event, ], match, clinical,
                                      samples=unlist(tumour))
    opt <- optimalSurvivalCutoff(clinical, eventPSI, censoring="right", 
                                 event="days_to_death", 
                                 timeStart="days_to_death")
    (optimalCutoff <- opt$par)    # Optimal exon inclusion level
    (optimalPvalue <- opt$value)  # Respective p-value
    
    label     <- labelBasedOnCutoff(eventPSI, round(optimalCutoff, 2), 
                                    label="PSI values")
    survTerms <- processSurvTerms(clinical, censoring="right",
                                  event="days_to_death", 
                                  timeStart="days_to_death",
                                  group=label, scale="years")
    surv <- survfit(survTerms)
    pvalue <- testSurvival(survTerms)
    plotSurvivalCurves(surv, pvalue=pvalue, mark=FALSE)
}

Differential gene expression

Detected alterations in alternative splicing may simply be a reflection of changes in gene expression levels. Therefore, to disentangle these two effects, differential expression analysis between tumour stage I and normal samples should also be performed. In order to do so:

# Prepare groups of samples to analyse and further filter unavailable samples in
# selected groups for gene expression
ge           <- geneExprNorm[ , unlist(normalVSstage1Tumour), drop=FALSE]
isFromGroup1 <- colnames(ge) %in% normalVSstage1Tumour[[1]]
design       <- cbind(1, ifelse(isFromGroup1, 0, 1))

# Fit a gene-wise linear model based on selected groups
library(limma)
fit <- lmFit(ge, design)

# Calculate moderated t-statistics and DE log-odds using limma::eBayes
ebayesFit <- eBayes(fit, trend=TRUE)

# Prepare data summary
pvalueAdjust <- "BH" # Benjamini-Hochberg p-value adjustment (FDR)
summary <- topTable(ebayesFit, number=nrow(fit), coef=2, sort.by="none",
                    adjust.method=pvalueAdjust, confint=TRUE)
names(summary) <- c("log2 Fold-Change", "CI (low)", "CI (high)", 
                    "Average expression", "moderated t-statistics", "p-value", 
                    paste0("p-value (", pvalueAdjust, " adjusted)"),
                    "B-statistics")
attr(summary, "groups") <- normalVSstage1Tumour

# Calculate basic statistics
stats <- diffAnalyses(ge, normalVSstage1Tumour, "basicStats", 
                      pvalueAdjust=NULL)

final <- cbind(stats, summary)

# Differential gene expression between breast tumour stage I and normal samples
library(ggplot2)
library(ggrepel)
cognateGenes <- unlist(parseSplicingEvent(events)$gene)
logFCthreshold  <- abs(final$`log2 Fold-Change`) > 1
pvalueThreshold <- final$`p-value (BH adjusted)` < 0.01

final$genes <- gsub("\\|.*$", "\\1", rownames(final))
ggplot(final, aes(`log2 Fold-Change`, 
                  -log10(`p-value (BH adjusted)`))) +
    geom_point(data=final[logFCthreshold & pvalueThreshold, ],
               colour="orange", alpha=0.5, size=3) + 
    geom_point(data=final[!logFCthreshold | !pvalueThreshold, ],
               colour="gray", alpha=0.5, size=3) + 
    geom_text_repel(data=final[cognateGenes, ], aes(label=genes),
                    box.padding=0.4, size=5) +
    theme_light(16) +
    ylab("-log10(q-value)")

# UHRF2 skipped exon 10's PSI values per tumour stage I and normal samples
UHRF2skippedExon10 <- events[1]
plotDistribution(psi[UHRF2skippedExon10, ], normalVSstage1Tumour)

# Find optimal cutoff for the event
UHRF2skippedExon10 <- events[1]
eventPSI <- assignValuePerSubject(psi[UHRF2skippedExon10, ], match, clinical,
                                  samples=unlist(tumour))
opt <- optimalSurvivalCutoff(clinical, eventPSI, censoring="right", 
                             event="days_to_death", timeStart="days_to_death")
(optimalCutoff <- opt$par)    # Optimal exon inclusion level

## [1] 0.1436954

(optimalPvalue <- opt$value)  # Respective p-value

## [1] 0.0358

label     <- labelBasedOnCutoff(eventPSI, round(optimalCutoff, 2), 
                                label="PSI values")
survTerms <- processSurvTerms(clinical, censoring="right",
                              event="days_to_death", timeStart="days_to_death",
                              group=label, scale="years")
surv <- survfit(survTerms)
pvalue <- testSurvival(survTerms)
plotSurvivalCurves(surv, pvalue=pvalue, mark=FALSE)

plotDistribution(geneExprNorm["UHRF2", ], normalVSstage1Tumour)

UHRF2ge <- assignValuePerSubject(geneExprNorm["UHRF2", ], match, clinical, 
                                 samples=unlist(tumour))

# Survival curves based on optimal gene expression cutoff
opt <- optimalSurvivalCutoff(clinical, UHRF2ge, censoring="right",
                             event="days_to_death", timeStart="days_to_death")
(optimalCutoff <- opt$par)    # Optimal exon inclusion level

## [1] 10.42619

(optimalPvalue <- opt$value)  # Respective p-value

## [1] 0.176

# Process again after rounding the cutoff
roundedCutoff <- round(optimalCutoff, 2)
label     <- labelBasedOnCutoff(UHRF2ge, roundedCutoff, label="Gene expression")
survTerms <- processSurvTerms(clinical, censoring="right",
                              event="days_to_death", timeStart="days_to_death",
                              group=label, scale="years")
surv   <- survfit(survTerms)
pvalue <- testSurvival(survTerms)
plotSurvivalCurves(surv, pvalue=pvalue, mark=FALSE)

Load GTEx data

GTEx data (subject phenotype, sample attributes, gene expression and junction quantification) can be automatically retrieved and loaded by following these commands:

# Check GTEx tissues available based on the sample attributes
getGtexTissues(sampleAttr)

tissues <- c("blood", "brain")
gtex <- loadGtexData("~/Downloads", tissues=tissues)

Junction quantification can be loaded for all GTEx tissues by typing:

gtex <- loadGtexData("~/Downloads", tissues=NULL)

Load SRA project data using recount

recount2 is a resource of pre-processed data for thousands of SRA projects (including gene read counts, splice junction quantification and sample metadata). psichomics supports automatic downloading and loading of SRA data from recount2, as exemplified below:

library(recount)
View(recount_abstract)
sra <- loadSRAproject("SRP053101")

Please refer to our methods article for more information (the code for performing the analysis can be found at GitHub):

Nuno Saraiva-Agostinho and Nuno L. Barbosa-Morais (2020). Interactive Alternative Splicing Analysis of Human Stem Cells Using psichomics. In: Kidder B. (eds) Stem Cell Transcriptional Networks. Methods in Molecular Biology, vol 2117. Humana, New York, NY

Load other SRA and user-provided local files

Although only select SRA projects are available to be automatically downloaded (based on pre-processed data from the recount2 project), other SRA projects and user-provided data can be manually aligned using a splice-aware aligner and loaded by following the instructions in Loading SRA and user-provided RNA-seq data.

Local files can be loaded by indicating their containing folder Any files located in this folder and sub-folders will be loaded.

For instance, to load GTEx files via local files, create a directory called GTEx, put all GTEx files within that folder and type these commands:

folder <- "~/Downloads/GTEx/"
ignore <- c(".aux.", ".mage-tab.")
data <- loadLocalFiles(folder, ignore=ignore)

# Select clinical and junction quantification dataset
clinical      <- data[[1]][["Clinical data"]]
sampleInfo    <- data[[1]][["Sample metadata"]]
geneExpr      <- data[[1]][["Gene expression"]]
junctionQuant <- data[[1]][["Junction quantification"]]