Functionality

The main functions of OutSplice achieve the following for either user provided data or data provided from the TCGA.

1. Junction normalization
2. Outlier analysis
3. Determination of a junctional outlier as a skipping, insertion, or deletion
4. Calculation of splicing burden
5. Plotting expression levels

Minimum required packages

OutSplice will import the below packages and their dependencies so please ensure they are installed before using the software:

1. AnnotationDbi
2. GenomicRanges
3. GenomicFeatures
4. IRanges
5. org.Hs.eg.db
6. Repitools
7. TxDb.Hsapiens.UCSC.hg19.knownGene
8. TxDb.Hsapiens.UCSC.hg38.knownGene
9. S4Vectors

Installation from Bioconductor

The OutSplice package is available at https://bioconductor.org and can be installed via BiocManager::install:

A package only needs to be installed once. Load the package into an R session with:

library(OutSplice)

outspliceAnalysis/outspliceTCGA

Full Sample data for OutSplice can be downloaded from The Cancer Genome Atlas (TCGA).

The outspliceAnalysis function has 4 required inputs and 13 optional Inputs. The outspliceTCGA function has 3 required inputs and 9 optional Inputs and should only be used if the data is in TCGA format.

In the below examples we run outspliceAnalysis and outspliceTCGA on a subset of Head and Neck squamous cell carcinoma data obtained from The Broad Institute TCGA GDAC Firehose database.

junction <- system.file("extdata", "HNSC_junctions.txt.gz", package = "OutSplice")
gene_expr <- system.file("extdata", "HNSC_genes_normalized.txt.gz", package = "OutSplice")
rawcounts <- system.file("extdata", "Total_Rawcounts.txt", package = "OutSplice")
sample_labels <- system.file("extdata", "HNSC_pheno_table.txt", package = "OutSplice")
output_file_prefix <- "OutSplice_Example"
TxDb_hg19 <- "TxDb.Hsapiens.UCSC.hg19.knownGene"
dir <- paste0(tempdir(), "/")
message("Output is located at: ", dir)

## Output is located at: /tmp/RtmpVH4mPw/

results <- outspliceAnalysis(junction, gene_expr, rawcounts, sample_labels, saveOutput = TRUE, output_file_prefix, dir, filterSex = TRUE, annotation = "org.Hs.eg.db", TxDb = TxDb_hg19, offsets_value = 0.00001, correction_setting = "fdr", p_value = 0.05)

## Loading data

## pheno
## Normal  Tumor 
##     44     50

## convert to RPM

## filter the putative junctions

## run the ogsa function for pre filtering

## get the genomic information for all the junctions

##   403 genes were dropped because they have exons located on both strands
##   of the same reference sequence or on more than one reference sequence,
##   so cannot be represented by a single genomic range.
##   Use 'single.strand.genes.only=FALSE' to get all the genes in a
##   GRangesList object, or use suppressMessages() to suppress this message.

## 'select()' returned 1:1 mapping between keys and columns

## filtering for putative junctions, skipping, insertion and deletion events based on known exons.

## skipping

## insertions

## deletions

## remove all that map to 'NA' no gene name, and assign gene expression from gene_expr

## align with gene_expr data

## subset removing any genes without normalization

## Perform normalization using gene_expr values

## run the ogsa function

junction <- system.file("extdata", "TCGA_HNSC_junctions.txt.gz", package = "OutSplice")
gene_expr <- system.file("extdata", "TCGA_HNSC_genes_normalized.txt.gz", package = "OutSplice")
rawcounts <- system.file("extdata", "Total_Rawcounts.txt", package = "OutSplice")
output_file_prefix <- "TCGA_OutSplice_Example"
dir <- paste0(tempdir(), "/")
message("Output is located at: ", dir)

## Output is located at: /tmp/RtmpVH4mPw/

results_TCGA <- outspliceTCGA(junction, gene_expr, rawcounts, saveOutput = TRUE, output_file_prefix, dir, filterSex = TRUE, annotation = "org.Hs.eg.db", TxDb = "TxDb.Hsapiens.UCSC.hg19.knownGene", offsets_value = 0.00001, correction_setting = "fdr", p_value = 0.05)

## Loading data

## pheno
## Normal  Tumor 
##     44     50

## convert to RPM

## filter the putative junctions

## run the ogsa function for pre filtering

## get the genomic information for all the junctions

##   403 genes were dropped because they have exons located on both strands
##   of the same reference sequence or on more than one reference sequence,
##   so cannot be represented by a single genomic range.
##   Use 'single.strand.genes.only=FALSE' to get all the genes in a
##   GRangesList object, or use suppressMessages() to suppress this message.

## 'select()' returned 1:1 mapping between keys and columns

## filtering for putative junctions, skipping, insertion and deletion events based on known exons.

## skipping

## insertions

## deletions

## remove all that map to 'NA' no gene name, and assign gene expression from gene_expr

## align with gene_expr data

## subset removing any genes without normalization

## Perform normalization using gene_expr values

## run the ogsa function

-junction: filepath to a matrix of junction raw read counts

-gene_expr: filepath to a matrix of normalized gene expression data. RSEM quartile normalized data is recommended.

-rawcounts: filepath to a matrix of raw read counts for each gene. Can either be per gene, or a summed total for each sample.

-sample_labels: filepath to a matrix of phenotypes (in Section 2.2)

-saveOutput: saves results output to a user specified directory. Default is FALSE [OPTIONAL]

-output_file_prefix: user defined string for what the prefix of the output file should be named. Default is NULL [OPTIONAL]

-dir: string containing the desired output path. Ex) “~/my_outsplice_output/”. Default is NULL [OPTIONAL]

-filtersex: ignores sex chromosomes when generating results. Default is TRUE [OPTIONAL]

-annotation: the bioconductor package containing the annotations to use. Uses the human genome by default (in Section 3.4) [OPTIONAL]

-TxDb: the bioconductor package containing the transcript annotations to use. The outspliceAnalysis function uses the hg38 genome by default, and outspliceTCGA uses hg19 by default (in Section 3.4) [OPTIONAL]

-offsets_value: the normalized junction expression threshold. Uses 0.00001 by Default (in Section 3.6) [OPTIONAL]

-correction_setting: the correction value to be used for p-value adjustment during Fisher analyses. The available options are: “holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, and “none”. Uses ‘fdr’ by Default. [OPTIONAL]

-p_value: the significance threshold to use during Fisher analyses. Uses 0.05 by default. [OPTIONAL]

-use_junc_col: an integer indicating which column in the junction matrix contains the junction regions in your matricies. Uses the first column (1) by default. outSpliceAnalysis only. [OPTIONAL]

-use_gene_col: which column in the gene_expr matrix contains the entrez ids of your genes. Uses the first column (1) by default. outSpliceAnalysis only. [OPTIONAL]

-use_rc_col: which column in the rawcounts matrix contains the row names. Uses the first column (1) by default. outSpliceAnalysis only. [OPTIONAL]

-use_labels_col: which column in the sample_labels matrix contains the sample names. Uses the first column (1) by default. outSpliceAnalysis only. [OPTIONAL]

This algorithm is compatible with any organism provided the genome object, genome wide annotation, and transcript annotations exist as packages on Bioconductor.

Phenotype matrix

If using the outspliceAnalysis function, users must provide a phenotype matrix, designating which samples in the junction file belong to the tumor group (labeled as “T”) and the normal group (labeled as “F”). Please ensure the matrix file contains a header row with the first column designating the sample names, and the second column designating the phenotype. If using TCGA data, the two phenotypes are “tumor” and “normal.” outspliceTCGA can automatically infer the phenotype of TCGA data using the sample names. Only if the phenotype matrix has more than 10 control samples, the program proceeds.

plotJunctionData

After running outspliceAnalysis or outspliceTCGA, the user can create plots of specified junctions using the plotJunctionData function.

The function plotJunctionData has 1 required input and 10 optional inputs

In the below example we can plot the expression levels for a “skipping” junction event on the ECM1 gene based on our example data from the TCGA.

data_file <- system.file("extdata", "OutSplice_Example_2023-01-06.RDa", package = "OutSplice")
ecm1_junc <- "chr1:150483674-150483933"
pdf <- "ecm1_expression.pdf"
pdf_output <- paste0(tempdir(), "/", pdf)
message("Output is located at: ", pdf_output)

## Output is located at: /tmp/RtmpVH4mPw/ecm1_expression.pdf

plotJunctionData(data_file, NUMBER = 1, junctions = ecm1_junc, tail = NULL, p_value = 0.05, GENE = FALSE, SYMBOL = NULL, makepdf = TRUE, pdffile = pdf_output, tumcol = "red", normcol = "blue")

## chr1:150483674-150483933

-data_file: filepath to an R Data file containing output from the outspliceAnalysis or outspliceTCGA functions.

-NUMBER: number of junctions to plot. This can be top number of junctions (over or under expressed), or can be specific junctions in a list. Default is 1 [OPTIONAL]

-junctions: you can input the specific junction you want to graph (or vector of junctions). Default is NULL [OPTIONAL]

-tail: you can specify if you want top over or under expressed with tail=‘RIGHT’ for junctions overexpressed in tumors, or tail=‘LEFT’ for junctions underexpressed in tumors. Default is NULL [OPTIONAL]

-p_value: p-value threshold to use when plotting the top over or under-expressed junctions with “tail”. Default is 0.05 [OPTIONAL]

-GENE: Pick junctions based on a specific gene. TRUE means you will pick all the junctions mapping to a certain gene. FALSE means you do not pick based on the gene. Default is False. [OPTIONAL]

-SYMBOL: HGNC SYMBOL of gene you want to graph [OPTIONAL]

-makepdf: Save graphs to a pdf? Default is no/FALSE [OPTIONAL]

-pdffile: if you want to save to a pdf, you will need to write the file path [OPTIONAL when makepdf = F]

-tumcol: color for tumors on graph. Default is red [OPTIONAL]

-normcol: color for normals on graph. Default is blue [OPTIONAL]

Junction RPM normalization

The program automatically normalizes the junction counts by dividing the junction counts by the total raw counts and then dividing each count by 10^6 to generate RPM junction data.

OGSA initial filtering

The dotheogsa function from the Bioconductor package OGSA is sourced to remove junctions that may not be biologically relevant due to low expression or that have any difference between tumor and normal. In this package, we set a 0.1 RPM expression threshold for pre-filtering.

OGSA outlier analysis

The dotheogsa function is again employed to determine splicing events as outliers, which are defined as any normalized junctions that are two standard deviations above or below the median of the normal distribution. A Fisher exact test is used to determine which junctions are significantly over or under expressed in tumors compared to the normal samples.

Genomic references

The Bioconductor GenomicRanges packages are used to assign each junction to a known gene. The user has the option in the main function to input which genome and its associated Bioconductor packages to use as the reference. These Bioconductor packages should include the genome object, annotations, and transcript annotations.

Ex) For mouse genomes aligned to mm39, install and load: Mus.musculus (genome object), org.Mm.eg.db (annotations), and TxDb.Mmusculus.UCSC.mm39.refGene (transcript annotations) using the library() function. Then, when using the OutSplice functions, specify org.Mm.eg.db for the annotation argument, and TxDb.Mmusculus.UCSC.mm39.refGene for the TxDb argument.

This process will also work for any genome provided the correct annotation package, its corresponding transcript annotation package, and their dependencies from Bioconductor are provided and loaded.

Using this genomic assignment, the dotheogsa function determines insertion, skipping, or deletion events based on the following criteria:

insertion: junction that starts outside a known exon

skipping: junction that skips over a known exon

deletion: junction that is inside a known exon but not as its start or end

Junction expression normalization

Junction expression is normalized based on its corresponding gene expression from the gene_expr input. This is achieved by dividing the junction RPM data by the normalized gene expression counts from a junction’s corresponding gene. If a junction is aligned to more than one gene, then the gene with the lower entrez id will be the one selected for the normalization.

Filter by expression via offsets

Offsets, which the user can specify, sets a minimum value relative to the normal samples in order to call a junction an outlier. The goal is to remove data with low expression that may not be biologically relevant. In this example example, an outlier junction must have a normalized expression greater than 0.00001 in order to be called an outlier. Any outliers with expressions below this value are too low to be relevant for the analysis in this example.

Splice Burden Calculation

Sums the number of splicing events in each sample that were marked as a TRUE outlier for both over-expressed and under-expressed events. The total number of outliers is then calculated as the sum of the over and under-expressed outliers.

Junction Plotting

Creates bar and waterfall plots of junction expression in both the tumor and normal samples. The data for these plots comes from the raw junction input, the gene expression values to reflect overall gene expression, and the junction expression normalized by gene expression.

outSpliceAnalysis/outSpliceTCGA

Returns a list and, if specified, an R data file and tab deliminated text files with the following data:

-ASE.type: significant junction events labeled by type (skipping, insertion, or deletion)

-FisherAnalyses: Data Frame of junction events containing the number of under/over-expressed outliers in the tumor group (Num_UE_Outliers/Num_OE_Outliers), the Fisher p-value for under/over-expressed events (FisherP1/FisherP2), and a ranking of the under/over expressed events (UE_Rank/OE_Rank)

-geneAnnotations: object containing gene names corresponding to each junction region

junc.Outliers: list containing the logical matrices TumorOverExpression and TumorUnderExpression. “True” indicates an over-expressed event in TumorOverExpression, or an under-expressed event in TumorUnderExpression.

-junc.RPM: junction counts in reads per million following a division of the junction counts input by the total rawcounts for each sample

-junc.RPM.norm: junction counts normalized by each event’s total gene expression value

-gene_expr: gene expression values for each junction event

-splice_burden: matrix containing the number of Fisher-P significant over-expressed, under-expressed, and total number of outliers per sample

-NORM.gene_expr.norm: Median of junction data normalized by gene expression for normal samples only (Used for Junction Plotting Only)

-pheno: Phenotypes of Samples (Tumor or Normal)

-pvalues: Junction Fisher P-values

If file output was specified, the files are:

-Data File: _.RDa

-ASE.type: event_types.txt

-FisherAnalyses: FisherAnalyses.txt

-geneAnnotations: gene_annotations.txt

-junc.Outliers: TumorOverExpression.txt, TumorUnderExpression.txt

-splice_burden: splice_burden.txt

plotJunctionData

Outputs junction expression plots of user specified junctions as defined in Section 3.8. Plots can be saved to a user defined pdf file.