Recent breakthroughs in bioinformatics now allow us to accurately reconstruct and quantify full-length gene isoforms from RNA-sequencing data via tools such as StringTie, Kallisto and Salmon. Alternatively long-read RNA-seq (third generation sequencing) now routinely provide us with full length transcripts. These full length transcripts makes it possible to analyze changes in 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 switches from RNA-seq derived quantification of novel and/or known full-length isoforms. IsoformSwitchAnalyzeR facilitates integration of many sources of annotation such as protein domains, signal peptides, Intrinsically Disordered Regions (IDR), subcellular localization, 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 of isoform switches of individual genes, the genome-wide consequences of isoform switches, and the associated changes in 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.
What To Cite (please remember)
Analyzing Alternative Splicing
Frequently Asked Questions, Problems and Errors
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. For an introduction to isoform switches please refere to Vitting-Seerup et al 2017.
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 StringTie, 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:
To solve all these problems, we developed IsoformSwitchAnalyzeR.
IsoformSwitchAnalyzeR is an easy-to-use R package that enables the user to import the result of the (reconstructed) full-length isoforms quantification an RNA-seq experiment into R. If annotated transcripts are analyzed, IsoformSwitchAnalyzeR offers integration with the multi-layer information stored in a GTF/GFF 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. Additionally, 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:
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.
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)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 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.
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 (also citations are what allow people to continue both maintaining and developing bioinformatic tools).
If you are using the:
Refrences:
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? and note that IsoformSwitchAnalyzeR both supports analysis of short- and long-read sequencing data.
From the isoform quantification IsoformSwitchAnalyzeR performs five high-level tasks:
Statistical identification of isoform switches.
Integration of a wide range of (predicted) annotations for the isoforms involved in the identified switches (e.g. protein domains).
Identification of which isoforms have a predicited functional consequnce (e.g. loss/gain of protein domain).
Visualization of predicted consequences of the isoform switches for individual genes
Analysis of genome wide patterns in both switch consequences and alternative splicing.
Please note that just like any other statistical tool IsoformSwitchAnalyzeR requires count data from 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. All of these analysis 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).
The resulting fasta files enables the usage of external sequence analysis tools. Currently the external sequence tools directly supported by IsoformSwitchAnalyzeR are:
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 (or is only interested in splicing), it is possible to run the full workflow using:
isoformSwitchAnalysisCombined()This corresponds to running isoformSwitchAnalysisPart1() and isoformSwitchAnalysisPart2() without adding the external results.
Back to Table of Content.
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
Also note this workflow have been written to be used with:
packageVersion('IsoformSwitchAnalyzeR')
#> [1] '2.6.0'Since IsoformSwitchAnalyzeR is (significantly) updated quite frequently please make sure you use the up to date version.
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/IsoQuant (for long-read data), IsoformSwitchAnalyzeR can be used with quantifications from any tool which produce isoform-level quantification.
For the purpose of illustrating importing data into R let us use some cufflinks data re-quantified quantified via Salmon which is included in the distribution of IsoformSwitchAnalyzeR.
Please note:
The approach illustrated below for Salmon data is identical if the data was quantified with Kallisto, RSEM or StringTie.
If you have used Salmon to quantify newer versions of the Human/Mouse/Dropsophila Gencode/Ensembl transcriptomes there is an even easier way to import the data (using tximeta) as described in Importing Data from Salmon via Tximeta.
For other sources of quantification (including StringTie and IsoQuant) 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.
# parentDir = "/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 abundance values (not counts) via edgeR...
#> DoneNote that:
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. Type “?importIsoformExpression” for more information.
When using the “parentDir” argument importIsoformExpression() can can import only a subset of the subfolders via the “pattern” argument (and its friends). Type “?importIsoformExpression” for details.
The result of using importIsoformExpression() is a list containing both count and abundance estimates for each isoform in each sample (count and abundance matrix):
head(salmonQuant$abundance, 2)
#> isoform_id Fibroblasts_0 Fibroblasts_1 hESC_0 hESC_1 iPS_0 iPS_1
#> 1 TCONS_00000001 7.481356 7.969927 0 0.000000 0.000000 5.320791
#> 2 TCONS_00000002 0.000000 0.000000 0 2.015807 6.627746 3.218813
head(salmonQuant$counts, 2)
#> isoform_id Fibroblasts_0 Fibroblasts_1 hESC_0 hESC_1 iPS_0 iPS_1
#> 1 TCONS_00000001 12.30707 14.02487 0 0.0000000 0.00000 18.13313
#> 2 TCONS_00000002 0.00000 0.00000 0 0.1116201 21.10248 10.96964Please note that (estimated) isoform level counts are nessesary for the statistical analysis.
Apart from the isoform quantification we need 2-3 additional sets of annotation in the IsoformSwitchAnalyzeR workflow:
A design matrix indicating which of the independent biological replicates belong to which condition (and if there are other covariates (incl batch effects) that should be taking into account - more info about that below).
The transcript structure of the isoforms (in genomic coordinates) as well as information about which isoforms originate from the same gene. This information 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. If you are using a refrence-only workflow (tools such as Kallisto/Salmon/RSEM etc) this argument should point to the refrence database GTF corresponding to the fasta file that you used to build the refrence index. If you use a guided/de-novo transcriptome assembly approach (tools like StringTie and Cufflinks) this argument should point to the GTF file created at the “merge” stage of the workflow. Please refere to [What Quantification Tool(s) Should I Use] for a more detailed description of the two different workflows. Please note that IsoformSwitchAnalyzeR also support RefSeq GFF files - for all other database you have to use a GTF file - see What Transcript Database Should I Use? for recomendations on databases and more info on how to obtain these.
If you have 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 via BSgenomes (more info in the Detailed Workflow) importing the transcript fasta file at this stage will provide a much faster (and potentially more accurate) analysis. Providing the sequences here will be (a lot!) easier for people analyzing non-model organisms. You can also find suggestions about how to obtain fasta files which are paired with the GTF/GFF files from some of the major databases in What Transcript Database Should I Use?.
Since the information in point 2-3 from above are just files stored on the computer which IsoformSwitchAnalyzeR itself imports into R, 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. For the example data this can be done from the file names as follows:
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 iPSPlease 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:
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.
# isoformExonAnnoation="/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"),
fixStringTieAnnotationProblem = TRUE,
showProgress = FALSE
)
#> Step 1 of 10: Checking data...
#> Step 2 of 10: Obtaining annotation...
#> importing GTF (this may take a while)...
#> Step 3 of 10: Fixing StringTie gene annoation problems...
#> There were no need to rescue any annotation
#> 302 genes_id were assigned their original gene_id instead of the StringTie gene_id.
#> This was only done when it could be done unambiguous.
#> Step 4 of 10: Calculating expression estimates from count data...
#> Skipped as user supplied expression via the "isoformRepExpression" argument...
#> Step 5 of 10: Testing for unwanted effects...
#> Added 1 batch/covariates to the design matrix
#> Step 6 of 10: Batch correcting expression estimates...
#> Step 7 of 10: Extracting data from each condition...
#> Step 8 of 10: Making comparisons...
#> Step 9 of 10: Making switchAnalyzeRlist object...
#> Step 10 of 10: Guestimating differential usage...
#> The GUESSTIMATED number of genes with differential isoform usage are:
#> comparison estimated_genes_with_dtu
#> 1 Fibroblasts vs hESC 34 - 57
#> 2 Fibroblasts vs iPS 47 - 78
#> 3 hESC vs iPS 23 - 39
#> Donesummary(aSwitchList)
#> This switchAnalyzeRlist list contains:
#> 1092 isoforms from 368 genes
#> 3 comparison from 3 conditions (in total 6 samples)
#>
#> Feature analyzed:
#> [1] "ORFs, ntSequence"Please 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.
By setting fixStringTieAnnotationProblem = TRUE (default) importRdata will automatically try and correct some of the annoation problems created when doing transcript assembly (unassigned transcripts and merged genes). For more information see documentation of importRdata (by typing ?importRdata in R).
The guesstimated provided is just that - a guess informed by an estimate. This is provided for getting a quick idea of what to expect. The final number will naturally not be known untill further down the workflow.
For more information about the switchAnalyzeRlist see IsoformSwitchAnalyzeR Background Information.
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:
### load example data
data("exampleSwitchList")
### Subset for quick runtime
exampleSwitchList <- subsetSwitchAnalyzeRlist(
switchAnalyzeRlist = exampleSwitchList,
subset = abs(exampleSwitchList$isoformFeatures$dIF) > 0.4
)
### Print summary of switchAnalyzeRlist
summary(exampleSwitchList)
#> This switchAnalyzeRlist list contains:
#> 20 isoforms from 12 genes
#> 1 comparison from 2 conditions (in total 6 samples)
#>
#> Feature analyzed:
#> [1] "ORFs, ntSequence"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,
# pathToOutput = 'path/to/where/output/should/be/'
outputSequences = FALSE, # change to TRUE whan analyzing your own data
prepareForWebServers = FALSE # change to TRUE if you will use webservers for external sequence analysis
)Note that:
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.
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().
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.
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 )
#> Comparison nrIsoforms nrSwitches nrGenes
#> 1 hESC vs Fibroblasts 10 5 5In a typical workflow, the user would here have to use produced fasta files to perform the external analysis tools (Pfam (for prediction of protein domains), SignalP (for prediction of signal peptides), CPAT or CPC2 (for prediction of coding potential, since they perform similar analysis so it is only necessary to run one of them)), Prediction of sub-cellular location (DeepLoc2) and isoform toplogy (DeepTmHmm), IUPred2A or NetSurfP-2 (for prediction of Intrinsically Disordered Regions (IDR), since they perform similar analysis so it is only necessary to run one of them. The main difference is that NetSurfP-2 is newer but IUPred2A also allows for prediction of IDBS (which are also integrated into the IsoformSwitchAnalyzeR workflow)). 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 some of these tools on the example data - these results are also included in the IsoformSwitchAnalyzeR package.
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,
n = 10, # if plotting was enabled, it would only output the top 10 switches
removeNoncodinORFs = TRUE,
pathToCPC2resultFile = system.file("extdata/cpc2_result.txt" , package = "IsoformSwitchAnalyzeR"),
pathToPFAMresultFile = system.file("extdata/pfam_results.txt" , package = "IsoformSwitchAnalyzeR"),
pathToIUPred2AresultFile = system.file("extdata/iupred2a_result.txt.gz" , package = "IsoformSwitchAnalyzeR"),
pathToSignalPresultFile = system.file("extdata/signalP_results.txt" , package = "IsoformSwitchAnalyzeR"),
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 : Predicting functional consequences...
#>
#> The number of isoform switches with functional consequences identified were:
#> Comparison nrIsoforms nrSwitches nrGenes
#> 1 hESC vs Fibroblasts 10 5 5Please 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_00000098 XLOC_000088:KIF1B KIF1B hESC Fibroblasts
#> 2 geneComp_00000205 XLOC_000177 LDLRAD2 hESC Fibroblasts
#> 3 geneComp_00000389 XLOC_001345 ARHGEF19 hESC Fibroblasts
#> gene_switch_q_value switchConsequencesGene Rank
#> 1 8.467463e-85 TRUE 1
#> 2 1.582351e-41 TRUE 2
#> 3 4.584799e-16 TRUE 3Now the data gathering and integration is complete (and some analysis have already been saved to your computer if you used outputPlots=TRUE. This means we are now ready to look at the results via the buld in visualizations as described in the next sections.
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. The fully annotated switchAnalyzeRlist (corresponding to what was created with isoformSwitchAnalysisPart1() and isoformSwitchAnalysisPart2() above) with 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 nrSwitches nrGenes
#> 1 COAD_ctrl vs COAD_cancer 653 372 349
#> 2 LUAD_ctrl vs LUAD_cancer 384 205 198
#> 3 Combined 946 531 497One 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_1','condition_2','gene_switch_q_value','Rank')],
gene_name == 'ZAK'
)
#> gene_name condition_1 condition_2 gene_switch_q_value Rank
#> 7 ZAK COAD_ctrl COAD_cancer 5.344466e-06 7Let us take a look at the isoform switch in the ZAK gene (7th most significant 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 (this is not necessary if there is only one condition in your data).
switchPlot(
exampleSwitchListAnalyzed,
gene='ZAK',
condition1 = 'COAD_ctrl',
condition2 = 'COAD_cancer',
localTheme = theme_bw(base_size = 13) # making text sightly larger for vignette
)