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

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
)