IsoformSwitchAnalyzeR

Enabling Identification and Analysis of Isoform Switches with Functional Consequences and the Associated Alternative Splicing

Kristoffer Vitting-Seerup

2019-05-02

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 unfortunately 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/CDS), protein domains (via Pfam), signal peptides (via SignalP), Intrinsically Disordered Regions (IDR, via NetSurfP-2), coding potential (via CPAT or CPC2) and sensitivity to Non-sense Mediated Decay (NMD) and more. The combination of identified isoform switches and their annotation enables IsoformSwitchAnalyzeR to predict potential functional consequences of the identified isoform switches — such as loss of protein domains — thereby identifying isoform switches of particular interest. Lastly, IsoformSwitchAnalyzeR provides article-ready visualization methods for isoform switches for individual genes as well as both summary statistics and visualization of the genome-wide changes/consequences of isoform switches, their consequences and 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) and its associated alternative splicing, thereby expanding the usability of RNA-seq data.


Table of Content

Preliminaries

What To Cite (please remember)

Quick Start

Detailed Workflow

Analyzing Alternative Splicing

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 hundreds of switch events in both cancer and non-cancer studies.

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 (reconstructed) 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 tools for assessing coding potential of an isoform (both the CPAT and CPC2 tools are supported) which 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 amine acid sequence of the CDS/ORF. The amine 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 aTSS and aTTS sites - a functionality also implemented in IsoformSwitchAnalyzeR.

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 in both individual genes and on a genome wide level.

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 changes thereof.

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, although originally developed for model organisms now also directly support analysis of all organisms.

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 (or don’t know if you have used it), simply copy-paste the following into an R session to install the basic Bioconductor packages (will only done if you don’t already have them):

if (!requireNamespace("BiocManager", quietly = TRUE)){
    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 the developmental version of IsoformSwitchAnalyzeR this can be done from the GitHub. Please note that this is for advanced uses and should not be done unless you have good reason to. Installation from the GitHub can be done by copy/pasting the following into R:

if (!requireNamespace("devtools", quietly = TRUE)) {
    install.packages("devtools")
}
devtools::install_github("kvittingseerup/IsoformSwitchAnalyzeR", build_vignettes = TRUE)


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 of the individual function documentation. This vignette also contains a lot of information and will be continuously updated, so make sure to read both sources carefully as it contains the answers to the most questions (including 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 and ember to post i) the code that gave the error, ii) the error message as well as iii) the output of running sessionInfo().

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 code or error messages are included please post i) the code that gave the error, ii) the error message as well as iii) the output of running sessionInfo().

If you have suggestions for improvements, please put them on GitHub. This will allow other people to up vote 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) link.
  2. Vitting-Seerup et al. IsoformSwitchAnalyzeR: Analysis of changes in genome-wide patterns of alternative splicing and its functional consequences. Bioinformatics (2019) link.
  3. Nowicka et al. DRIMSeq: a Dirichlet-multinomial framework for multivariate count outcomes in genomics. F1000Research, 5(0), 1356. link.
  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. link.
  5. Weischenfeldt et al. Mammalian tissues defective in nonsense-mediated mRNA decay display highly aberrant splicing patterns. Genome Biol 2012, 13:R35 link.
  6. Huber et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat. Methods, 2015, 12:115-121. link.
  7. Wang et al. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 2013, 41:e74. link.
  8. Finn et al. The Pfam protein families database. Nucleic Acids Research (2012) link.
  9. _Almagro et al. SignalP 5.0 improves signal peptide predictions using deep neural networks.**. Nat. Biotechnol (2019)_ link
  10. Soneson et al. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research 4, 1521 (2015). link.
  11. Robinson et al. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology (2010) link.
  12. Ritchie et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research (2015) link.
  13. Anders et al. Detecting differential usage of exons from RNA-seq data. Genome Research (2012) link.
  14. Kang et al. CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res (2017) link.
  15. Klausen et al. NetSurfP-2.0: improved prediction of protein structural features by integrated deep learning. BioRxiv (2018) link


Quick Start

Overview of Isoform Switch Workflow

In this section, we will outline the IsoformSwitchAnalyzeR workflow - in the following sections we will show you how to do it including example data and code which can be copy/pasted to do your own analysis.

The idea behind IsoformSwitchAnalyzeR is to make it easy to do advanced post-analysis of full-length RNA-seq derived isoform/transcript level quantification of RNASeq data with a focus on finding, annotating and visualizing isoform switches with functional consequences as well as its associated alternative splicing. If you want to know more about considerations and recommendations for bioinformatic tools which can perform the initial isoform quantification, please refer to What Quantification Tool(s) Should I Use?.

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

  • Statistical identification of isoform switches.
  • Integration of a wide range of annotations for the isoforms involved in the identified switches.
  • Visualization of predicted consequences of the isoform switches, for individual genes as well as for analyzing the genome wide patterns.

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 (if not more) independent replicates are needed for good statistical performance - most tools have a hard time controlling the False Discovery Rate (FDR) with only two replicates so extra caution is needed for interpreting results based on three or fewer replicates.


The first step of a IsoformSwitchAnalyzeR workflow is to import and integrate the isoform quantification with its basic annotation. Once you have all the relevant data imported into R (IsoformSwitchAnalyzeR will also help you with that), the 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 filtering out lowly expressed genes/isoforms, statistical analysis to identify isoform switches, annotating those switches with open reading frames (ORF) and writing the nucleotide and amino acid (peptide) sequences to fasta files. These fasta files enables the usage of external sequence analysis tools such as

  • Pfam: Prediction of protein domains, which can be run either locally or via their webserver.
  • Tools for Coding Potential Analysis:
    • CPC2: The Coding Calculator 2, which can be run either locally or via their webserver.
    • CPAT: The Coding-Potential Assessment Tool, 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 (V5 is supported).
  • NetSurfP-2: Prediction of Intrinsically Disordered Regions (IDR), 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 analysis step (and what function actually perform them).


Part 2) Plot All Isoform Switches and Their annotation. This part involves importing and incorporating the results of all the external sequence analysis, analyzing alternative splicing, predicting functional consequences of isoform switches and plotting i) individual genes with isoform switches and ii) genome wide patterns in the consequences of isoform switching.

All of this can be done using the function:

isoformSwitchAnalysisPart2()

See below for example of usage, and Detailed Workflow for details on the individual analysis step (and what function actually perform them).


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()). Please note the boxes marked with “external analysis” means the user have to run tools outside of R either locally or on a web-server - more info can be found below.

Back to Table of Content.


Example Workflow

As indicted above, a comprehensive, 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. For the detailed, fully customizable workflow please refer to the Detailed Workflow (though we still recommend reading this section as it will provide a nice overview of the workflow).

Let us start by loading IsoformSwitchAnalyzeR into the R environment.

library(IsoformSwitchAnalyzeR)

Note that newer versions of RStudio (a wrapper for R which makes it much easier to work with R) support auto-completion of package names so you just have to type “Iso” and press then RStudio will suggest IsoformSwitchAnalyzeR.


Importing the Data

The first step is to import all data needed for the analysis into R and then concatenate them as a switchAnalyzeRlist object. Although IsoformSwitchAnalyzeR has different functions for importing data from tools such as Salmon/Kallisto/RSEM/StringTie/Cufflinks, IsoformSwitchAnalyzeR can be used with quantifications from all isoform-level quantification tools.

For the purpose of illustrating importing data into R let us use some data quantified via Salmon which is included in the distribution of IsoformSwitchAnalyzeR. Note the approach for Kallisto/RSEM/StringTie would be identical - for other sources of quantification (including Cufflinks/Cuffdiff) see Importing Data Into R. The import of the example quantification can be done as follows:

### Please note
# The way of importing files in the following example with
# "system.file('pathToFile', package="IsoformSwitchAnalyzeR") is
# specialized way of accessing the example data in the IsoformSwitchAnalyzeR package
# and not something  you need to do - just supply the string e.g.
# "/path/to/mySalmonQuantifications/" pointing to the parent directory (where 
# each sample is a separate 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
#>     Found 6 quantification file(s) of interest
#> Step 2 of 3: Reading data...
#> reading in files with read_tsv
#> 1 2 3 4 5 6 
#> Step 3 of 3: Normalizing FPKM/TxPM values via edgeR...
#> Done

Note that:

  1. When using the “parentDir” argument importIsoformExpression() assumes each quantification files are stored in seperate sub-directory. To accomodate cases where that is not the case the “sampleVector” argument should be used instead.

  2. When using the “parentDir” argument importIsoformExpression() can can import only a subset of the subfolders via the “pattern” argument (and its friends). See documentation for importIsoformExpression() for more details.

The result of using importIsoformExpression is 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      7.481356      7.969927      0 0.000000 0.000000
#> 2 TCONS_00000002      0.000000      0.000000      0 2.015807 6.627746
#>      iPS_1
#> 1 5.320791
#> 2 3.218813

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 2-3 additional sets of annotation.

  1. 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).

  2. The transcript structure of the isoforms (in genomic coordinates) as well as information about which isoforms originate from the same gene. All this is typically stored in a GTF file - a file format which IsoformSwitchAnalyzeR can directly import and integrate (as described below) you just have to tell IsoformSwitchAnalyzeR where the file is located on your computer/server.

  3. if you have “only” quantified the transcriptome via tools such as Salmon/Kallisto/RSEM (aka you did NOT do any guided/de-novo transcriptome assembly (using tools such as Cufflinks/StringTie)) we highly recommended that the transcript nucleotide sequences (the fasta file used to make the index for the quantification tool) is also supplied to IsoformSwitchAnalyzeR. If a fasta file is provided IsoformSwitchAnalyzeR which will import and integrate these sequences into the workflow. Although options are available to extracting transcript sequences later (more info in the Detailed Workflow) importing the transcript fasta file at this stage will provide a more accurate (and much faster) analysis. Providing the sequences here will be (a lot!) easier for people analyzing non-model organisms.

Since the information in point 2-3 from above are just files stored on the computer which IsoformSwitchAnalyzeR itself imports, the only thing you need to make yourself is the design matrix. The design matrix will have to be generated manually to ensure correct grouping of samples:

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


Please also note that if you have additional covariates in the experimental setup 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 conditions 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 might influence the result (e.g. sex or age). See How to handle confounding effects (including batches) for more info.

Now that we have the isoform quantifications, the design matrix, the isoform annotation and the isoform nucleotide sequence (respectively in a GTF file and a fasta file stored on the computer, both of which IsoformSwitchAnalyzeR will import itself) we can now integrate it all into a switchAnalyzeRlist.

Please note that:

  • It is highly recommended to both supply count and abundance expression matrices to the importRdata() as count are better for the statistical analysis but the abundance estimates are better to measure effect sizes. If necessary IsoformSwitchAnalyzeR can calculate FPKM values from the counts but the values estimated by the quantification tools are much better!
  • It is essential that the quantification, the GTF file and the fasta file all contain information on the same isoforms. Although this is not always the case with Ensemble data (which is why we recommend using files from GENCODE genes - where The “Comprehensive gene annotation” GTF and the “Transcript sequences” fasta file is a perfect pair) IsoformSwitchAnalyzeR can handle it using its build in arguments - but it might require a couple of tries to get those arguments correct. For more information (including instructions for how to use Ensemble annotation) refer to this FAQ: The error “The annotation and quantification (count/abundance matrix and isoform annotation) seems to be different (Jaccard similarity < 0.925)”.
  • It is highly recommended to also supply the nucleotide fasta file as IsoformSwitchAnalyzeR will import and propagate the sequences in the entire workflow making it more accurate and much faster.


Since we now have all the data we can use the importRdata() function to import and integrate it all into a switchAnalyzeRlist:

### 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"),
    isoformNtFasta       = system.file("extdata/example_isoform_nt.fasta.gz", package="IsoformSwitchAnalyzeR"),
    showProgress = FALSE
)
aSwitchList
#> This switchAnalyzeRlist list contains:
#>  1092 isoforms from 362 genes
#>  3 comparison from 3 conditions (in total 6 samples)
#> 
#> Feature analyzed:
#> [1] "ntSequence"

Note that by supplying the GTF file to the “isoformExonAnnoation” argument IsoformSwitchAnalyzeR will automatically import and integrate CoDing Sequence (CDS, aka the part that is predicted to give rise to a protein) regions from in the GTF file as the ORF regions used by IsoformSwitchAnalyzeR.

For more information about the switchAnalyzeRlist see IsoformSwitchAnalyzeR Background Information.



Part 1

Before we can run the analysis, it is necessary to know that IsoformSwitchAnalyzeR measures isoform usage vian isoform fraction (IF) values which quantifies the fraction of the parent gene expression originating from a specific isoform (calculated as isoform_exp / gene_exp). Consequently, the difference in isoform usage is quantified 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).

To illustrate the IsoformSwitchAnalyzeR workflow we will use a smaller example dataset equal to what you would get from following the import steps described above (used for convenience). The example 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 (in total 6 samples)
#> 
#> Switching features:
#>            Comparison switchingIsoforms switchingGenes
#> 1 hESC vs Fibroblasts                 0              0
#> 
#> Feature analyzed:
#> [1] "Isoform Switch Identification, ORFs, ntSequence"

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

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 separate fasta files (here turned off to make the example run - so a user should use outputSequences=TRUE).

exampleSwitchList <- isoformSwitchAnalysisPart1(
    switchAnalyzeRlist   = exampleSwitchList,
    dIFcutoff            = 0.3, # Cutoff for finding switches - set high for short runtime in example data
    # pathToOutput = 'path/to/where/output/should/be/'
    outputSequences      = FALSE, # prevents creation of fasta files used for external sequence analysis (you should change this)
    prepareForWebServers = FALSE  # change to TRUE if you will use webservers 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 should 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.

  4. If you are planning on using the webserver version of the external sequence analysis tools you should set prepareForWebServers=TRUE. This argument exists because some webservers (Pfam and SignalP) have (quite) strict limits on the number (and length) of amino acid sequences one can submit at the time. When prepareForWebServers=TRUE (default) IsoformSwitchAnalyzeR will generate multiple filtered amino acid fasta files which can then be supplied as multiple individual analysis at the websites as well as trim the sequences which are too long (will not cause problems downstream since IsoformSwitchAnalyzeR have been optimized to take this into account).

The number of switching features is easily summarized as follows:

extractSwitchSummary(
    exampleSwitchList,
    dIFcutoff = 0.3 # the same cutoff to the summary function as used above
)
#>            Comparison nrIsoforms nrGenes
#> 1 hESC vs Fibroblasts         26      13


In a typical workflow, the user would here have to use produced fasta files to perform the external analysis tools (Pfam (protein domains), SignalP (signal peptides), CPAT or CPC2 (coding potential, since they perform similar analysis so it is only necessary to run one of them)), NetSurfP-2 (Intrinsically Disordered Regions (IDR)). For more information on how to run those tools refer to Advice for Running External Sequence Analysis Tools and Downloading Results. To illustrate the workflow, we will here use the result of running those tools on the example data which 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()" command not nesseary when using your own
# data - just supply the path as a string
# (e.g. pathToCPC2resultFile = "/myFiles/myCPC2results.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
  pathToCPC2resultFile      = system.file("extdata/cpc2_result.txt"         , package = "IsoformSwitchAnalyzeR"),
  pathToPFAMresultFile      = system.file("extdata/pfam_results.txt"        , package = "IsoformSwitchAnalyzeR"),
  pathToNetSurfP2resultFile = system.file("extdata/netsurfp2_results.csv.gz", 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
)
#> Step 1 of 3 : Importing external sequence analysis...
#> Step 2 of 3 : Analyzing alternative splicing...
#> Step 3 of 3 : Prediciting functional consequences...
#> 
#> The number of isoform swithces with functional consequences identified were:
#>            Comparison nrIsoforms nrGenes
#> 1 hESC vs Fibroblasts         22      11

Please note that if you had to submit your data as multiple runs at the Pfam or SignalP website you can just supply a vector of strings indicating the path to each of the resulting files and IsoformSwitchAnalyzeR will read them all and integrate them.

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

The top switching genes can be extracted as follows:

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

To illustrate visualizations, both for individual genes and for the genome-wide analysis, we will use a larger dataset which is a subset of two of the TCGA Cancer types analyzed in Vitting-Seerup et al 2017. Please note that all analysis describe below can be done with any data analyzed as describe above. The TCGA example data can be loaded as follows:

data("exampleSwitchListAnalyzed")

The number of isoform switches with predicted functional consequences are extracted by setting “filterForConsequences = TRUE” in the extractSwitchSummary() function:

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 switchPlot

One isoform switch of particular interest is found in the ZAK gene (as highlighted in Vitting-Seerup et al 2017). ZAK is a signaling gene involved in the response to radiation (one of the most common treatments for cancers). It is ranked as the 7th most significant isoform switch in the COAD as can easily be seen by subsetting on the top isoform switches (extracted with extractTopSwitches()):

subset(
    extractTopSwitches(
        exampleSwitchListAnalyzed,
        filterForConsequences = TRUE,
        n=10,
        inEachComparison = TRUE
    )[,c('gene_name','condition_2','gene_switch_q_value','Rank')],
    gene_name == 'ZAK'
)
#>   gene_name condition_2 gene_switch_q_value Rank
#> 7       ZAK COAD_cancer        5.344466e-06    7

Let us take a look at the isoform switch in the ZAK gene using the switchPlot functionality. Note that since the switchAnalyzeRlist contains two different comparison (LUAD and COAD) we need to specify the conditions to plot (not necessary if there is only one).

switchPlot(
    exampleSwitchListAnalyzed,
    gene='ZAK',
    condition1 = 'COAD_ctrl',
    condition2 = 'COAD_cancer'
)

From this plot, we first note that no difference in the gene expression of ZAK between normal and cancer samples (bottom left) meaning that a standard gene level analysis would not pick up any changes to ZAK. Next, we see the large highly 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 inferred that in the healthy samples it is primary the sort isoforms that are used. In the COAD patients, there is however a large increase in the use of the long isoform. Interestingly the long isoform contains a SAM_1 protein domain, a domain which facilitates interaction with both other protein and RNA molecules. This leads to the hypothesis that this isoform, use a lot more in the tumors, results in a protein which can interact with a new/different set of interaction partners in the cancer. Furthermore, it is noteworthy all 3 isoforms share an Intrinsically Disordered Regions (IDR) in the N-terminal even though they are encoded by different regions of the chromosome.

Note that if you want to save this plot as a pdf file via the pdf command, you need to specify “onefile = FALSE”. The following 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:

  • The differential gene/isoform expression analysis is not a part of the standard IsoformSwitchAnalyzeR workflow but can easily be added as described in Adding differential gene expression.
  • the switchPlot() function has an argument called IFcutoff which requires the Isoform Fraction of an Isoform to be larger than IFcutoff (default 0.05) to be included in the plot. Increasing this cutoff can result in “cleaner” plots as less used isoforms will be removed.
  • Some regions of transcripts might be marked as “Not Analyzed” if you have been using webservers - for more info on that see The switchPlot contains regions “Not Annotated”


Genome-wide Summaries

To analyze the global usage of isoform switches the next step would typically be to look at the number, and overlap, between conditions - which can be visualized with the extractSwitchOverlap() function:

extractSwitchOverlap(
    exampleSwitchListAnalyzed,
    filterForConsequences=TRUE
)