IsoformSwitchAnalyzeR

Enabling Identification and Analysis of Isoform Switches with Functional Consequences from RNA-sequencing data

Kristoffer Vitting-Seerup

2018-10-30

Abstract

Recent breakthroughs in bioinformatics now allow us to accurately reconstruct and quantify full-length gene isoforms from RNA-sequencing data (via tools such as Cufflinks, StringTie, Kallisto and Salmon). These tools make it possible to analyzing alternative isoform usage, but unfortunatly this is rarely done meaning RNA-seq data is typically not used to its full potential.

To solve this problem we developed IsoformSwitchAnalyzeR. IsoformSwitchAnalyzeR is an easy-to use-R package that enables statistical identification of isoform switching from RNA-seq derived quantification of novel and/or annotated full-length isoforms. IsoformSwitchAnalyzeR facilitates integration of many sources of (predicted) annotation such as Open Reading Frame (ORF), protein domains (via Pfam), signal peptides (via SignalP), coding potential (via CPAT) and sensitivity to Non-sense Mediated Decay (NMD) etc. The combination of identified isoform switches and their annotation also enables IsoformSwitchAnalyzeR to predict potential functional consequences of the identified isoform switches — such as loss of protein domains or coding potential — thereby identifying isoform switches of particular interest. Lastly, IsoformSwitchAnalyzeR provides article-ready visualization methods for isoform switches, and summary statistics describing the genome-wide occurences of isoform switches, their consequences as well as the associated alternative splicing.

In summary, IsoformSwitchAnalyzeR enables analysis of RNA-seq data with isoform resolution with a focus on isoform switching (with predicted consequences), thereby expanding the usability of RNA-seq data.


Table of Content

Preliminaries

What To Cite (please remember)

Quick Start

Detailed Workflow

Analyzing Alternative Splicing (new)

Other workflows

Frequently Asked Questions, Problems and Errors

Final Remarks

Sessioninfo


Preliminaries

Background and Package Description

The usage of Alternative Transcription Start sites (aTSS), Alternative Splicing (AS) and alternative Transcription Termination Sites (aTTS) are collectively collectively results in the production of different isoforms. Alternative isoforms are widely used as recently demonstrated by The ENCODE Consortium, which found that on average, 6.3 different transcripts are generated per gene; a number which may vary considerably per gene.

The importance of analyzing isoforms instead of genes has been highlighted by many examples showing functionally important changes. One of these examples is the pyruvate kinase. In normal adult homeostasis, cells use the adult isoform (M1), which supports oxidative phosphorylation. However, almost all cancer cells use the embryonic isoform (M2), which promotes aerobic glycolysis, one of the hallmarks of cancer. Such shifts in isoform usage are termed ‘isoform switching’ and cannot be detected at when only analyzing data on gene level.

On a more systematic level several recent studies suggest that isoform switches are quite common since they often identify hundres of switche events.

Tools such as Cufflinks, Salmon and Kallisto allows for reconstruction and quantification of full-length transcripts from RNA-seq data. Such data has the potential to facilitate genome-wide analysis of alternative isoform usage and identification of isoform switching — but unfortunately these types of analyses are still only rarely done; most analyses are on gene level only.

We hypothesize that there are multiple reasons why RNA-seq data is not used to its full potential:

  1. There is still a lack of tools that can identify isoform switches with isoform resolution
  2. Although there are many excellent tools to perform sequence analysis, there is no common framework which allows for integration of the analysis provided by these tools.
  3. There is a lack of tools facilitating easy and article-ready visualization of isoform switches.

To solve all these problems, we developed IsoformSwitchAnalyzeR.

IsoformSwitchAnalyzeR is an easy-to-use R package that enables the user to import (novel) full-length derived isoforms from an RNA-seq experiment into R. If annotated transcripts are analyzed, IsoformSwitchAnalyzeR offers integration with the multi-layer information stored in a GTF file including the annotated coding sequences (CDS). If transcript structures are predicted (either de-novo or guided) IsoformSwitchAnalyzeR offers an accurate tool for identifying the dominant ORF of the isoforms. The knowledge of isoform positions for the CDS/ORF allows for prediction of sensitivity to Nonsense Mediated Decay (NMD) — the mRNA quality control machinery that degrades isoforms with pre-mature termination codons (PTC).

IsoformSwitchAnalyzeR facilitates identification of isoform switches via newly developed statistical methods that tests each individual isoform for differential usage and thereby identifies the exact isoforms involved in an isoform switch.

Since we know the exon structure of the full-length isoform, IsoformSwitchAnalyzeR can extract the underlying nucleotide sequence from a reference genome. This enables integration with the Coding Potential Assessment Tool (CPAT) which predicts the coding potential of an isoform and can be used to increase accuracy of ORF predictions. By combining the CDS/ORF isoform positions with the nucleotide sequence, we can extract the most likely amino acid sequence of the CDS/ORF. The amino acid sequence enables integration of analysis of protein domains (via Pfam) and signal peptides (via SignalP) — both supported by IsoformSwitchAnalyzeR. Lastly, since the structures of all expressed isoforms from a given gene are known, one can also annotate alternative splicing - including retentions - a functionality also implemented in IsoformSwitchAnalyzeR.

Thus, in summary, IsoformSwitchAnalyzeR enables annotation of isoforms with intron retention, ORF, NMD sensitivity, coding potential, protein domains and signal peptides (and many more), resulting in the ability to predict important functional consequences of isoform switches.

IsoformSwitchAnalyzeR contains tools that allow the user to create article-ready visualizations of:

  1. Individual isoform switches
  2. Genome-wide analysis of isoform switches and their predicted consequences
  3. Genome-wide analysis of alternative splicing and isoform switches and their predicted consequences.

These visualizations are easy to understand and integrate all information gathered throughout the workflow. Example of visualizations can be found in the Examples Visualizations section.

Lastly IsoformSwitchAnalyzeR is based on standard Bioconductor classes such as GRanges and BSgenome. Thus, it supports all species- and annotation versions facilitated by the Bioconductor annotation packages.

Back to Table of Content.

Installation

IsoformSwitchAnalyzeR is part of the Bioconductor repository and community which means it is distributed with, and dependent on, Bioconductor. Installation of IsoformSwitchAnalyzeR is easy and can be done from within the R terminal. If it is the first time you use Bioconductor, simply copy-paste the following into an R session to install the basic Bioconductor packages:

install.packages("BiocManager")
BiocManager::install()

If you already have installed Bioconductor, running these two commands will check whether updates for installed packages are available.


After you have installed the basic Bioconductor packages you can install IsoformSwitchAnalyzeR by copy-pasting the following into an R session:

BiocManager::install("IsoformSwitchAnalyzeR")

This will install the IsoformSwitchAnalyzeR package as well as other R packages that are needed for IsoformSwitchAnalyzeR to work.


If you need to install from the developmental branch of Bioconductor it is nesseary to explicitely specify that. Please note that this is for advanced uses and should not be done unless you have good reason to. Installation from the developmental branch can be done by copy/pasting:

BiocManager::install("IsoformSwitchAnalyzeR", version = "devel")


How To Get Help

This R package comes with plenty of documentation. Much information can be found in the R help files (which can easily be accessed by running the following command in R “?functionName”, for example “?importRdata”) - here a lot of information can be found in the individual argument description as well as in the details section. This vignette contains a lot of information and will be continously updated, so make sure to read both sources carefully as it contains the answers to the most Frequently Asked Questions, Problems and Errors.

If you want to report a bug/error (found in the newest version of the R package!) please make an issue with a reproducible example at github.

If you have unanswered questions or comments regarding IsoformSwitchAnalyzeR or how to run it please post them on the associated google group (after make sure the question was not already answered).

If you have suggestions for improvements, please put them on github. This will allow other people to upvote your idea, thereby showing us there is wide support of implementing your idea.

Back to Table of Content.


What To Cite

The analysis performed by IsoformSwitchAnalyzeR is only possible due to a string of other tools and scientific discoveries — please read this section thoroughly and cite the appropriate articles to give credit where credit is due (and allow people to continue both maintaining and developing bioinformatic tools).

If you are using the


Refrences:

  1. Vitting-Seerup et al. The Landscape of Isoform Switches in Human Cancers. Cancer Res. (2017)
  2. Vitting-Seerup et al. IsoformSwitchAnalyzeR: Analysis of changes in genome-wide patterns of alternative splicing and its functional consequences. bioRxiv (2018)
  3. Nowicka et al. DRIMSeq: a Dirichlet-multinomial framework for multivariate count outcomes in genomics. F1000Research, 5(0), 1356.
  4. Vitting-Seerup et al. spliceR: an R package for classification of alternative splicing and prediction of coding potential from RNA-seq data. BMC Bioinformatics 2014, 15:81.
  5. Weischenfeldt et al. Mammalian tissues defective in nonsense-mediated mRNA decay display highly aberrant splicing patterns. Genome Biol 2012, 13:R35
  6. Huber et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat. Methods, 2015, 12:115-121.
  7. Wang et al. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 2013, 41:e74.
  8. Finn et al. The Pfam protein families database. Nucleic Acids Research (2014) Database Issue 42:D222-D230
  9. Petersen et al. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nature Methods, 8:785-786, 2011
  10. Soneson et al. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research 4, 1521 (2015).
  11. Robinson et al. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology (2010)
  12. Ritchie et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research (2015)
  13. Anders et al. Detecting differential usage of exons from RNA-seq data. Genome Research (2012)


Quick Start

Overview of Isoform Switch Workflow

The idea behind IsoformSwitchAnalyzeR is to make it easy to do advanced post-analysis of full-length RNA-seq derived transcripts quantification with a focus on finding, annotating and visualizing isoform switches with functional consequences. Furthermore IsoformSwitchAnalyzeR also allows for analysis of alternative splicing. Here we will go though the isoform switch workflow. For the workflow of alternative splicing see Analyzing Alternative Splicing. If you want to know more about recomendations for bioinformatic tools which can perform the initial isoform quantifications please refere to What Quantification Tool(s) Should I Use?.

From the isoform quantifications IsoformSwitchAnalyzeR performs three high-level tasks:

  • Statistical identification of isoform switches.
  • Annotation of the transcripts involved in the isoform switches.
  • Visualization of predicted consequences of the isoform switches, for individual switches and globally.

Please note that just like any other statistical tool IsoformSwitchAnalyzeR requires independent biological replicates (see What constitute an independent biological replicate?) and that recent benchmarks highlights that at least three independent replicates are needed for good statistical performance - most tools have a hard time controling the False Discovery Rate (FDR) with only two replicates so extra caution is needed for interpreting results based on few replicates.

A normal workflow for identification and analysis of isoform switches with functional consequences can be divided into two parts (also illustrated below in Figure 1).


Part 1) Extract Isoform Switches and Their Sequences. This part includes importing the data into R, identifying isoform switches, annotating those switches with open reading frames (ORF) and extracting the nucleotide and amino acid (peptide) sequences. The latter step enables the usage of external sequence analysis tools such as

  • CPAT : The Coding-Potential Assessment Tool, which can be run either locally or via their webserver.
  • Pfam : Prediction of protein domains, which can be run either locally or via their webserver.
  • SignalP : Prediction of Signal Peptides, which can be run either locally or via their webserver.

All of the above steps are performed by the high-level function:

isoformSwitchAnalysisPart1()

See below for example of usage, and Detailed Workflow for details on the individual steps.


Part 2) Plot All Isoform Switches and Their annotation. This part involves importing and incorporating the results of the external sequence analysis, identifying intron retention, predicting functional consequences and plotting i) all genes with isoform switches and ii) summaries of general consequences of switching.

All of this can be done using the function:

isoformSwitchAnalysisPart2()

See below for usage example, and Detailed Workflow for details on the individual steps.


Alternatively if one does not plan to incorporate external sequence analysis, it is possible to run the full workflow using:

isoformSwitchAnalysisCombined()

This corresponds to running isoformSwitchAnalysisPart1() and isoformSwitchAnalysisPart2() without adding the external results.

Figure 1: Workflow overview. Grey transparent boxes indicate the two parts of a normal workflow for analyzing isoform switches. The individual steps in the two sub-workflows are indicated by arrows. Speech bubbles summarize how this full analysis can be done in a two-step process using the high-level functions (isoformSwitchAnalysisPart1() and isoformSwitchAnalysisPart2()).

Back to Table of Content.


Example Workflow

As indicted above, a full, but less customizable, analysis of isoform switches can be done using the two high-level functions isoformSwitchAnalysisPart1() and isoformSwitchAnalysisPart2(). This section aims to show how these functions are used as well as illustrate what IsoformSwitchAnalyzeR can be used for.

Lets start by loading the R package

library(IsoformSwitchAnalyzeR)


Note that newer versions of RStudio support auto-completion of package names.

Importing the Data

The first step is to import all data needed for the analysis and store them in a switchAnalyzeRlist object. IsoformSwitchAnalyzeR has different functions for importing data from tools such as Salmon/Kallisto/RSEM/StringTie/Cufflinks, but can be used with all isoform-level expression data using implemented general-purpose functions. For the purpose of illustrating data import lets use some data quantified via Salmon. Note the approach for Kallisto/RSEM/StringTie would be identical - for other sources of quantification (including Cufflinks/Cuffdiff) see Importing Data Into R.

To illustrate the data import lets use some Salmon data included in the IsoformSwitchAnalyzeR package.

### Please note
# The way of importing files in the following example with
# "system.file('pathToFile', package="IsoformSwitchAnalyzeR") is
# specialiced way of accessing the example data in the IsoformSwitchAnalyzeR package
# and not somhting you need to do - just supply the string e.g.
# "/mySalmonQuantifications/" pointing to the parent directory (where 
# each sample is a seperate sub-directory) to the function.

### Import quantifications
salmonQuant <- importIsoformExpression(
    parentDir = system.file("extdata/",package="IsoformSwitchAnalyzeR")
)
#> Step 1 of 3: Identifying which algorithm was used...
#>     The quantification algorithm used was: Salmon
#> Step 2 of 3: Reading data...
#>     Found 6 quantification file(s) of interest
#> reading in files with read_tsv
#> 1 2 3 4 5 6 
#> Step 3 of 3: Normalizing FPKM/TxPM values via edgeR...
#> Done

Which results in a list containing both count and abundance estimates for each isoform:

head(salmonQuant$abundance, 2)
#>       isoform_id Fibroblasts_0 Fibroblasts_1 hESC_0    hESC_1    iPS_0
#> 1 TCONS_00000001      11.36932       14.0215      0 0.0000000  0.00000
#> 2 TCONS_00000002       0.00000        0.0000      0 0.6904076 13.07866
#>       iPS_1
#> 1 11.303557
#> 2  6.838089

head(salmonQuant$counts, 2)
#>       isoform_id Fibroblasts_0 Fibroblasts_1 hESC_0    hESC_1    iPS_0
#> 1 TCONS_00000001      12.30707      14.02487      0 0.0000000  0.00000
#> 2 TCONS_00000002       0.00000       0.00000      0 0.1116201 21.10248
#>      iPS_1
#> 1 18.13313
#> 2 10.96964


Apart from the isoform quantification we need three additional pices of information.

  • The transcript structure of the isoforms (in genomic coordinats). This is typically stored in a GTF file - a file format directly supported by IsoformSwitchAnalyzeR (as described below).
  • The information about which isoforms belongs to the same gene (also in the GTF file)
  • A design matrix indicating which of the independent biological replicates belong to which condition (and if there are other covariates that should be taking into account).

The design matrix will have to be generated manually to ensure correct grouping:

myDesign <- data.frame(
    sampleID = colnames(salmonQuant$abundance)[-1],
    condition = gsub('_.*', '', colnames(salmonQuant$abundance)[-1])
)
myDesign
#>        sampleID   condition
#> 1 Fibroblasts_0 Fibroblasts
#> 2 Fibroblasts_1 Fibroblasts
#> 3        hESC_0        hESC
#> 4        hESC_1        hESC
#> 5         iPS_0         iPS
#> 6         iPS_1         iPS


Note that if you have additional covariates in the data of interest these can also be added to the design matrix to ensure they are taking into account during the statistical analysis. ( Covariates are (unwanted) sources of variation not due to experimental groups you are interested in. This can be anything from a batch effects (e.g. data produced at different points in time) to an effect you are not interested in but which migth influence the result (e.g. sex or age). See How to handle cofounding effects (including batches) for more info ).

Now have the isoform quantifications, the design matrix and the isoform annotation (in a GTF file which IsoformSwitchAnalyzeR will import itsef) we can now combine and store all the relevant data in a switchAnalyzeRlist.

Please note that:

  • It is highly recommended to both supply count and abundance expression matrixes to the importRdata() but a switchAnalyzeRlist can also be created with only one of them.
  • It is essential that the isoforms quantified (with Salmon, Kallisto etc) and the annotation stored in the GTF file is identical. We highly recomend using files from GENCODE genes - where The “Comprehensive gene annotation” GTF and the “Transcript sequences” fasta file is a perfect pair. For more information (including instructions for how to use Ensemble annotation) refere to the FAQ The error “The annotation does not fit the expression data”.


Once we have all the nessesary data the easiest way to create a switchAnalyzeRlist is with the importRdata() function - which also directly can handle import and integration of annotation from a GTF file:

### Please note
# The way of importing files in the following example with
# "system.file("extdata/example.gtf.gz", package="IsoformSwitchAnalyzeR")"" is
# specialiced way of accessing the example data in the IsoformSwitchAnalyzeR package
# and not somhting you need to do - just supply the string e.g.
# "/myAnnotation/myQuantified.gtf" to the isoformExonAnnoation argument

### Create switchAnalyzeRlist
aSwitchList <- importRdata(
    isoformCountMatrix   = salmonQuant$counts,
    isoformRepExpression = salmonQuant$abundance,
    designMatrix         = myDesign,
    isoformExonAnnoation = system.file("extdata/example.gtf.gz", package="IsoformSwitchAnalyzeR"),
    showProgress = FALSE
)
aSwitchList
#> This switchAnalyzeRlist list contains:
#>  1092 isoforms from 362 genes
#>  3 comparison from 3 conditions

Note that by supplying the GTF file to the “isoformExonAnnoation” argument IsoformSwitchAnalyzeR will automatically import and integrate CDS regions from in the GTF file as the ORF regions used by IsoformSwitchAnalyzeR (if addAnnotatedORFs=TRUE (default)).

For more information about the switchAnalyzeRlist see IsoformSwitchAnalyzeR Background Information.



To illustrate the IsoformSwitchAnalyzeR workflow we will use a smaller example dataset originally quantified with Cufflinks/CuffDiff. The corresponding switchAnalyzeRlist is included in IsoformSwitchAnalyzeR and can loaded as follows:

data("exampleSwitchList")
exampleSwitchList
#> This switchAnalyzeRlist list contains:
#>  259 isoforms from 84 genes
#>  1 comparison from 2 conditions
#> 
#> Switching features:
#>            Comparison switchingIsoforms switchingGenes
#> 1 hESC vs Fibroblasts                 0              0
#> 
#> Feature analyzed:
#> [1] "Isoform Swich Identification"

Note that although a isoform switch identification analysis was performed (by CuffDiff), no genes with differential usage were found.


Part 1

Before we can run the analysis it is necessary to know that IsoformSwitchAnalyzeR measures isoform usage via isoform fraction (IF) values which quantifies the fraction of the parent gene expression originating from a specific isoform (calculated as / ). Consequently the difference in isoform usage is quantifed as the difference in isoform fraction (dIF) calculated as IF2 - IF1, and these dIF are used to measure the effect size (like fold changes are in gene/isoform expression analysis).

We can now run the first part of the isoform switch analysis workflow which filters for non-expressed genes/isoforms, identifies isoform switches, annotates open reading frames (ORF), switches with and extracts both the nucleotide and peptide (amino acid) sequences and output them as two seperate fasta files.

### isoformSwitchAnalysisPart1 needs the genomic sequence to predict ORFs

# Genome sequences are available from Bioconductor as BSgenome objects: 
# http://bioconductor.org/packages/release/BiocViews.html#___BSgenome
# Here we use the hg19 reference genome - which can be downloaded by
# copy/pasting the following two lines into the R terminal:
# BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")

library(BSgenome.Hsapiens.UCSC.hg19)

exampleSwitchList <- isoformSwitchAnalysisPart1(
    switchAnalyzeRlist = exampleSwitchList,
    genomeObject = Hsapiens, # the reference to the human BS genome
    dIFcutoff = 0.3,         # Cutoff for finding switches - set high for short runtime in example data
    # pathToOutput = 'path/to/where/ouput/should/be/'
    outputSequences = FALSE # prevents outputting of the fasta files used for external sequence analysis
)

Note that:

  1. In the example above we set outputSequences=FALSE only to make the example data run nicely (e.g. to prevent the function from out the two fasta files of the example data to your computer). When you analyze your own data you want to set outputSequences=TRUE and use the pathToOutput to specify where the files sould be saved. The two files outputted are needed to do the external analysis as described below.

  2. The isoformSwitchAnalysisPart1() function has an argument, overwritePvalues, which overwrites the result of an already existing switch test (such as those imported by cufflinks) with the result of running isoformSwitchTestDEXSeq().

  3. The switchAnalyzeRlist returned by isoformSwitchAnalysisPart1() has been reduced to only contain genes where an isoform switch (as defined by the alpha and dIFcutoff arguments) was identified. This enables much faster runtimes for the rest of the pipeline, as only isoforms from a gene with a switch are analyzed.

The number of switching features is easily summarized as follows:

extractSwitchSummary(
    exampleSwitchList,
    dIFcutoff = 0.3 # supply the same cutoff to the summary function
)
#>            Comparison nrIsoforms nrGenes
#> 1 hESC vs Fibroblasts         31      18


In a typical workflow the user would here have to use produced fasta files to perfrom the external analysis tools (Pfam (protein domains), SignalP (signal peptides), CPAT (coding potential)). For more information on how to run those tools refere to Advice for Running External Sequence Analysis Tools. To illustrate the workflow we will here use the result of running those tools on the example data wich we have also included in our R package.


Part 2

The second part of the isoform switch analysis workflow, which includes importing and incorporating external sequence annotation, analysis of alternative splicing, predicting functional consequences and visualizing both the general effects of isoform switches and the individual isoform switches. The combined analysis can be done by:

# Please note that in the following the part of the examples using
# the "system.file()" commandis not nesseary when using your own
# data - just supply the path as a string
# (e.g. pathToCPATresultFile = "/myFiles/myCPATresults.txt" )

exampleSwitchList <- isoformSwitchAnalysisPart2(
    switchAnalyzeRlist      = exampleSwitchList, 
    dIFcutoff               = 0.3,   # Cutoff for finding switches - set high for short runtime in example data
    n                       = 10,    # if plotting was enabled, it would only output the top 10 switches
    removeNoncodinORFs      = TRUE,  # Because ORF was predicted de novo
    pathToCPATresultFile    = system.file("extdata/cpat_results.txt"   , package = "IsoformSwitchAnalyzeR"),
    pathToPFAMresultFile    = system.file("extdata/pfam_results.txt"   , package = "IsoformSwitchAnalyzeR"),
    pathToSignalPresultFile = system.file("extdata/signalP_results.txt", package = "IsoformSwitchAnalyzeR"),
    codingCutoff            = 0.725, # the coding potential cutoff we suggested for human 
    outputPlots             = FALSE  # keeps the function from outputting the plots from this example
)

The exampleSwitchList now contains all the information needed to analyze isoform switches and alterantive splicing both for individual genes as well as genome-wide summary statistics or analysis.

The number of isoform switches with predicted functional consequences are extracted by setting “filterForConsequences = TRUE”:

extractSwitchSummary(
    exampleSwitchList,
    filterForConsequences = TRUE,
    dIFcutoff = 0.3               # supply the same cutoff to the summary function
) 
#>            Comparison nrIsoforms nrGenes
#> 1 hESC vs Fibroblasts         24      12

For each of these top genes, a switch plot will be generated if “outputPlots=TRUE” were used. Let’s take a closer look at the top candidates:

The top genes with isoform switches are:

extractTopSwitches(exampleSwitchList, filterForConsequences = TRUE, n=3)
#>            gene_ref           gene_id gene_name condition_1 condition_2
#> 1 geneComp_00000237 XLOC_000202:SRRM1     SRRM1        hESC Fibroblasts
#> 2 geneComp_00000098 XLOC_000088:KIF1B     KIF1B        hESC Fibroblasts
#> 3 geneComp_00000389       XLOC_001345  ARHGEF19        hESC Fibroblasts
#>   gene_switch_q_value switchConsequencesGene Rank
#> 1       1.553471e-185                   TRUE    1
#> 2        1.879647e-87                   TRUE    2
#> 3        1.312385e-83                   TRUE    3


Examples Visualizations

Let’s take a look at the isoform switch in the SRRM1 gene via the switch plot produced by IsoformSwitchAnalyzeR:

switchPlot(exampleSwitchList, gene='SRRM1')

From this plot, we first note that the gene expression is not significantly changed (bottom left). Next, we see the large significant switch in isoform usage across conditions (bottom right). By comparing which isoforms are changing (bottom right) to the isoform structure (top) it can be deduced that in hESC it is primarly the short isoforms that is used, but as the cells are differentiated to fibroblasts, there is a change towards mainly using the long isoform. Interestingly, this isoform switch seems to result in a truncation of the PWI domain which is important for DNA/RNA binding. SRRM1 is a splice factor and this switch could lead to the hypothesis that SRRM1 is nonfunctional in hESC but functional in Fibroblasts even though the total output from the gene is the same — a hypothesis which would naturallt have to be confirmed experimentally.

Note that if you want to save this plot as a pdf file via the pdf command, you need to specify “onefile = FALSE”. The folloing code chunk will produce a nicely-sized figure:

pdf(file = '<outoutDirAndFileName>.pdf', onefile = FALSE, height=5, width = 8)
switchPlot(exampleSwitchList, gene='SRRM1')
dev.off()

Furthermore, note that:

  • If the switchAnalyzeRList contains multiple comparison, you will also need to specify the ‘condition1’ and ‘condition2’ arguments in the switchPlot() function to indicate specifically which comparison you want to plot.
  • The differential isoform/gene expression analysis is not a part of the IsoformSwitchAnalyzeR workflow but can easily be added as described in Adding differential gene expression.
  • the switchPlot() function have a argument called IFcutoff which requires the Isoform Fraction of an isforom to be larger than IFcutoff (default 0.01) to be included in the plot. Increasing this cutoff can result in “cleaner” plots as minor isoforms will be removed.


To illustrate the genome-wide analysis of consequences of isoform switching in the different comparisons, we will use a larger dataset which is a subset of two of the TCGA Cancer types analyzed in Vitting-Seerup et al 2017.

data("exampleSwitchListAnalyzed")

The first step is to take a look at the number and overlap of isoform switches. The numbers can be extracted via the extractSwitchSummary() function and by setting “filterForConsequences=TRUE” we only extract the features involved in a switch which is predicted to have a functional consequence:

extractSwitchSummary(
    exampleSwitchListAnalyzed,
    filterForConsequences=TRUE
)
#>                 Comparison nrIsoforms nrGenes
#> 1 COAD_ctrl vs COAD_cancer        690     368
#> 2 LUAD_ctrl vs LUAD_cancer        422     218
#> 3                 combined       1008     529


The number, and overlap, between condtions can be viusalized with the extractSwitchOverlap() function:

extractSwitchOverlap(
    exampleSwitchListAnalyzed,
    filterForConsequences=TRUE
)