IsoformSwitchAnalyzeR

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.

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.

### 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

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

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 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"

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"

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
)

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

# 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

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

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

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

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

extractSwitchOverlap(
    exampleSwitchListAnalyzed,
    filterForConsequences=TRUE
)

extractConsequenceSummary(
    exampleSwitchListAnalyzed,
    consequencesToAnalyze='all',
    plotGenes = FALSE,           # enables analysis of genes (instead of isoforms)
    asFractionTotal = FALSE      # enables analysis of fraction of significant features
)

extractConsequenceEnrichment(
    exampleSwitchListAnalyzed,
    consequencesToAnalyze='all',
    analysisOppositeConsequence = TRUE,
    returnResult = FALSE # if TRUE returns a data.frame with the results
)

extractConsequenceEnrichmentComparison(
    exampleSwitchListAnalyzed,
    consequencesToAnalyze=c('domains_identified','intron_retention','coding_potential'),
    analysisOppositeConsequence = TRUE,
    returnResult = FALSE # if TRUE returns a data.frame with the results
)

extractSplicingEnrichment(
    exampleSwitchListAnalyzed,
    returnResult = FALSE # if TRUE returns a data.frame with the results
)

extractSplicingEnrichmentComparison(
    exampleSwitchListAnalyzed,
    splicingToAnalyze = c('A3','MES','ATSS','ATTS'), # the splicing highlighted above
    returnResult = FALSE # if TRUE returns a data.frame with the results
)

### Volcano like plot:
ggplot(data=exampleSwitchListAnalyzed$isoformFeatures, aes(x=dIF, y=-log10(isoform_switch_q_value))) +
     geom_point(
        aes( color=abs(dIF) > 0.1 & isoform_switch_q_value < 0.05 ), # default cutoff
        size=1
    ) +
    geom_hline(yintercept = -log10(0.05), linetype='dashed') + # default cutoff
    geom_vline(xintercept = c(-0.1, 0.1), linetype='dashed') + # default cutoff
    facet_wrap( ~ condition_2) +
    #facet_grid(condition_1 ~ condition_2) + # alternative to facet_wrap if you have overlapping conditions
    scale_color_manual('Signficant\nIsoform Switch', values = c('black','red')) +
    labs(x='dIF', y='-Log10 ( Isoform Switch Q Value )') +
    theme_bw()

### Switch vs Gene changes:
ggplot(data=exampleSwitchListAnalyzed$isoformFeatures, aes(x=gene_log2_fold_change, y=dIF)) +
    geom_point(
        aes( color=abs(dIF) > 0.1 & isoform_switch_q_value < 0.05 ), # default cutoff
        size=1
    ) + 
    facet_wrap(~ condition_2) +
    #facet_grid(condition_1 ~ condition_2) + # alternative to facet_wrap if you have overlapping conditions
    geom_hline(yintercept = 0, linetype='dashed') +
    geom_vline(xintercept = 0, linetype='dashed') +
    scale_color_manual('Signficant\nIsoform Switch', values = c('black','red')) +
    labs(x='Gene log2 fold change', y='dIF') +
    theme_bw()

Overview of Detailed Workflow

Before you start on this section, we recommend you to read though the Quick Start section, as it gives the large overview, introduces some basic concepts and gives a couple of tips not repeated here.

Compared to the workflow presented in Quick Start, a full workflow for analyzing isoform switches has many sub-steps which each can be customized/optimized. In this section, we will go more into depth with each of the steps as well as provide tips and shortcuts for working with IsoformSwitchAnalyzeR. Specifically, each of the main functions behind these steps will be described and illustrated. For a more comprehensive and detailed description of each individual function please refer to the individual function documentation build into the R package (easily accessed via ?functionName).

First we will start by getting an overview of the detailed workflow before diving into the individual parts. Then we will introduce the IsoformSwitchAnalyzeR Background Information - all you need to know before doing the detailed analysis. This is followed by a detailed step-by-step isoform switch analysis workflow and lastly an overview of useful Other Tools in IsoformSwitchAnalyzeR is provided.

The detailed workflow consists of the following steps (illustrated in Figure 2) which, just like before, can be divided into two parts:

Part 1) Extract Isoform Switches and Their Sequences.

• Importing Data Into R
  • Importing Data from Kallisto, Salmon, RSEM or StringTie
  • Importing Data from Cufflinks/Cuffdiff
  • Importing Data from other Full-length Transcript Assemblers
• Filtering
• Identifying Isoform Switches
• Analyzing Open Reading Frames
• Extracting Nucleotide and Amino Acid Sequences
• Advice for Running External Sequence Analysis Tools and Downloading Results

This corresponds to running the following functions in sequential order (which incidentally is just what the isoformSwitchAnalysisPart1() function does):

### Import data into R via one of:
myQantifications <- importIsoformExpression() # RSEM/Kallisto/Salmon/StringTie

### Create SwitchAnalyzeRlist
mySwitchList <- importRdata()         # OR
mySwitchList <- importCufflinksFiles()

### Filter
mySwitchList <- preFilter( mySwitchList )

### Test for isoform switches
mySwitchList <- isoformSwitchTestDEXSeq( mySwitchList )   # OR
mySwitchList <- isoformSwitchTestDRIMSeq( mySwitchList ) 

### If novel isoforms (else use CDS from ORF as explained in importRdata() )
mySwitchList <- analyzeORF( mySwitchList )

### Extract Sequences
mySwitchList <- extractSequence( mySwitchList )

### Summary
extractSwitchSummary(mySwitchList)

Part 2) Plot All Isoform Switches and Their annotation.

• Importing External Sequence Analysis
• Predicting Alternative Splicing
• Predicting Switch Consequences
• Post Analysis of Isoform Switches with Consequences
  • Analysis of Individual Isoform Switching
  • Genome-Wide Analysis of Isoform Switching

This corresponds to running the following functions in sequential order (just like the isoformSwitchAnalysisPart2() function does):

### Add annotation
mySwitchList <- analyzeCPAT( mySwitchList ) # OR
mySwitchList <- analyzeCPC2( mySwitchList )

mySwitchList <- analyzePFAM( mySwitchList )
mySwitchList <- analyzeSignalP( mySwitchList )
mySwitchList <- analyzeAlternativeSplicing( mySwitchList )

### Analyse consequences
mySwitchList <- analyzeSwitchConsequences( mySwitchList )

### Visual analysis
# Indiviudal switches
switchPlotTopSwitches( mySwitchList )

### Summary
extractSwitchSummary(mySwitchList, filterForConsequences = TRUE)

For a normal workflow, this would then typically be followed by the global analysis of alternative splicing and consequences:

# global consequence analysis
extractConsequenceSummary( mySwitchList )
extractConsequenceEnrichment( mySwitchList )
extractConsequenceGenomeWide( mySwitchList )

# global splicing analysis
extractSplicingSummary( mySwitchList )
extractSplicingEnrichment( mySwitchList )
extractSplicingGenomeWide( mySwitchList )

The combined workflow is visualized here:

Figure 2: Detailed workflow overview. The 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 along with a description and the main R functions for performing the steps.

IsoformSwitchAnalyzeR Background Information

data("exampleSwitchList")         # A newly created switchAnalyzeRlist + switch analysis
names(exampleSwitchList)
#> [1] "isoformFeatures"       "exons"                 "conditions"           
#> [4] "designMatrix"          "isoformCountMatrix"    "sourceId"             
#> [7] "isoformSwitchAnalysis" "orfAnalysis"           "ntSequence"

data("exampleSwitchListAnalyzed") # A fully analyzed switchAnalyzeRlist
names(exampleSwitchListAnalyzed)
#>  [1] "isoformFeatures"             "exons"                      
#>  [3] "conditions"                  "sourceId"                   
#>  [5] "orfAnalysis"                 "domainAnalysis"             
#>  [7] "signalPeptideAnalysis"       "AlternativeSplicingAnalysis"
#>  [9] "idrAnalysis"                 "switchConsequence"

### Preview
head(exampleSwitchList, 2)
#>             iso_ref          gene_ref     isoform_id     gene_id
#> 19 isoComp_00000007 geneComp_00000005 TCONS_00000007 XLOC_000005
#> 22 isoComp_00000008 geneComp_00000005 TCONS_00000008 XLOC_000005
#>    condition_1 condition_2 gene_name nearest_ref_id class_code
#> 19        hESC Fibroblasts      <NA>     uc009vjk.2          =
#> 22        hESC Fibroblasts      <NA>     uc001aau.2          =
#>    TSS_group_id length              locus gene_status gene_overall_mean
#> 19         TSS2   2750 chr1:322036-328580          OK          193.8339
#> 22         TSS3   4369 chr1:322036-328580          OK          171.6605
#>    gene_value_1 gene_value_2 gene_stderr_1 gene_stderr_2
#> 19      696.704      48.0566      3.592857      2.307488
#> 22      696.704      48.0566      3.592857      2.307488
#>    gene_log2_fold_change gene_p_value gene_q_value gene_significant
#> 19              -3.85774  2.66665e-09  3.20379e-08              yes
#> 22              -3.85774  2.66665e-09  3.20379e-08              yes
#>    iso_status iso_overall_mean iso_value_1 iso_value_2 iso_stderr_1
#> 19         OK         372.3803     358.383    29.28480     2.091049
#> 22         OK         372.3803     338.308     5.01291     1.322809
#>    iso_stderr_2 iso_log2_fold_change iso_p_value iso_q_value
#> 19    17.489950             -3.61328 4.83698e-05 0.000548629
#> 22     6.999295             -6.07655 2.64331e-03 0.015898000
#>    iso_significant IF_overall       IF1       IF2         dIF
#> 19             yes  0.5618896 0.5143978 0.6093814  0.09498364
#> 22             yes  0.2949481 0.4855835 0.1043126 -0.38127092
#>    isoform_switch_q_value gene_switch_q_value   PTC
#> 19                     NA                   1 FALSE
#> 22                     NA                   1 FALSE
# tail(exampleSwitchList, 2)

### Dimentions
dim(exampleSwitchList$isoformFeatures)
#> [1] 259  39

nrow(exampleSwitchList)
#> [1] 259
ncol(exampleSwitchList)
#> [1] 39
dim(exampleSwitchList)
#> [1] 259  39

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"

### Subset
subsetSwitchAnalyzeRlist(
    exampleSwitchList,
    exampleSwitchList$isoformFeatures$gene_name == 'ARHGEF19'
)
#> This switchAnalyzeRlist list contains:
#>  3 isoforms from 1 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"

head(exampleSwitchList$exons,2)
#> GRanges object with 2 ranges and 2 metadata columns:
#>       seqnames        ranges strand |     isoform_id     gene_id
#>          <Rle>     <IRanges>  <Rle> |    <character> <character>
#>   [1]     chr1 322037-322228      + | TCONS_00000007 XLOC_000005
#>   [2]     chr1 324288-324345      + | TCONS_00000007 XLOC_000005
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

createSwitchAnalyzeRlist()
preFilter()

Importing Data Into R

IsoformSwitchAnalyzeR can analyze the output from any tool that quantifies (annotated and/or de-novo/guided deconvoluted) full-length isoforms. All it requires are the following five sets of data:

• Isoform quantification (preferably both counts and abundance estimates (RPKM/FPKM/TPM)) in multiple samples from each condition (independent replicates are essential if differential isoform usage should be identified, see What constitute an independent biological replicate?)
• A design matrix indicating which samples belong to which conditions.
• The genomic coordinates of the isoform exon structure.
• Annotation of which isoform belongs to which gene.
• The nucleotide sequence of the isoforms quantified (can in some instances be circumvented - see below).

It is highly recommended to obtain the estimated replicate abundances and use them in the construction of the switchAnalyzeRlist as they better allow for effect size estimates.

From these data, a switchAnalyzeRlist (see The switchAnalyzeRlist for description) can be constructed. An important addition is that the design matrix (a data.frame with two mandatory columns: sampleID and condition) also allows for incorporation of additional cofactors (sources of variation not of interest for the analysis you are doing) in the experimental setup (see How to handle confounding effects (including batches) for more information).

To facilitate the usage of IsoformSwitchAnalyzeR, several dedicated wrappers for constructing a switchAnalyzeRlist from different data sources are included in IsoformSwitchAnalyzeR. See the following sections for details about the specific import methods:

• Importing Data from Kallisto, Salmon, RSEM or StringTie
• Importing Data from Cufflinks/Cuffdiff
• Importing Data from other Full-length Transcript Assemblers

For a (non-exhaustive) recommendations of bioinformatic tools which can perform the quantifications please refer to What Quantification Tool(s) Should I Use?

### 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.
# "/mySalmonQuantifications/" pointing to the parent directory (where 
# each sample is a seperate sub-directory) to the function.

### Import Salmon example data in R package
salmonQuant <- importIsoformExpression(
    parentDir = system.file("extdata/", package="IsoformSwitchAnalyzeR"),
    addIsofomIdAsColumn = TRUE
)
#> 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

### Make design matrix
myDesign <- data.frame(
    sampleID = colnames(salmonQuant$abundance)[-1],
    condition = gsub('_.*', '', colnames(salmonQuant$abundance)[-1])
)

### 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
)
#> Step 1 of 6: Checking data...
#> Step 2 of 6: Obtaining annotation...
#>     importing GTF (this may take a while)
#> Warning in importRdata(isoformCountMatrix = salmonQuant$counts, isoformRepExpression = salmonQuant$abundance, :    No CDS annotation was found in the GTF files meaning ORFs could not be annotated.
#>      (But ORFs can still be predicted with the analyzeORF() function)
#> Step 3 of 6: Calculating gene expression and isoform fraction...
#>      111 ( 9.23%) isoforms were removed since they were not expressed in any samples.
#> Step 4 of 6: Merging gene and isoform expression...
#> Step 5 of 6: Making comparisons...
#> Step 6 of 6: Making switchAnalyzeRlist object...
#> Done

head(aSwitchList$isoformFeatures,2)
#>            iso_ref          gene_ref     isoform_id     gene_id
#> 1 isoComp_00000001 geneComp_00000001 TCONS_00000001 XLOC_000001
#> 2 isoComp_00000002 geneComp_00000001 TCONS_00000002 XLOC_000001
#>   condition_1 condition_2 gene_name gene_biotype iso_biotype
#> 1 Fibroblasts        hESC      <NA>           NA          NA
#> 2 Fibroblasts        hESC      <NA>           NA          NA
#>   gene_overall_mean gene_value_1 gene_value_2 gene_stderr_1 gene_stderr_2
#> 1            6.4585     10.78392     1.007903      3.302564      1.007903
#> 2            6.4585     10.78392     1.007903      3.302564      1.007903
#>   gene_log2_fold_change gene_q_value iso_overall_mean iso_value_1
#> 1             -3.406546           NA         3.462012    7.725642
#> 2             -3.406546           NA         1.977061    0.000000
#>   iso_value_2 iso_stderr_1 iso_stderr_2 iso_log2_fold_change iso_q_value
#> 1    0.000000    0.2442856     0.000000            -9.595377          NA
#> 2    1.007903    0.0000000     1.007903             6.669457          NA
#>   IF_overall    IF1 IF2     dIF isoform_switch_q_value gene_switch_q_value
#> 1    0.43778 0.7829   0 -0.7829                     NA                  NA
#> 2    0.47538 0.0000   1  1.0000                     NA                  NA

head(aSwitchList$exons,2)
#> GRanges object with 2 ranges and 2 metadata columns:
#>       seqnames      ranges strand |     isoform_id     gene_id
#>          <Rle>   <IRanges>  <Rle> |    <character> <character>
#>   [1]     chr1 11874-12227      + | TCONS_00000001 XLOC_000001
#>   [2]     chr1 11874-12227      + | TCONS_00000002 XLOC_000001
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

head(aSwitchList$ntSequence,2)
#>   A DNAStringSet instance of length 2
#>     width seq                                          names               
#> [1]  1652 CTTGCCGTCAGCCTTTTCTTT...AAGCACACTGTTGGTTTCTG TCONS_00000001
#> [2]  1488 CTTGCCGTCAGCCTTTTCTTT...AAGCACACTGTTGGTTTCTG TCONS_00000002

# Etc...

importCufflinksFiles()

library(cummeRbund) # where the example data is stored

### Please note
# The way of importing file in this example with
# "system.file('pathToFile', package="cummeRbund") 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.
# "/myCuffdiffFiles/myFile.filetype" to the argument

testSwitchList <- importCufflinksFiles(
    pathToGTF                      = system.file('extdata/chr1_snippet.gtf',             package = "cummeRbund"),
    pathToGeneDEanalysis           = system.file('extdata/gene_exp.diff',                package = "cummeRbund"),
    pathToIsoformDEanalysis        = system.file('extdata/isoform_exp.diff',             package = "cummeRbund"),
    pathToGeneFPKMtracking         = system.file('extdata/genes.fpkm_tracking',          package = "cummeRbund"),
    pathToIsoformFPKMtracking      = system.file('extdata/isoforms.fpkm_tracking',       package = "cummeRbund"),
    pathToIsoformReadGroupTracking = system.file('extdata/isoforms.read_group_tracking', package = "cummeRbund"),
    pathToSplicingAnalysis         = system.file('extdata/splicing.diff',                package = "cummeRbund"),
    pathToReadGroups               = system.file('extdata/read_groups.info',             package = "cummeRbund"),
    pathToRunInfo                  = system.file('extdata/run.info',                     package = "cummeRbund"),
    fixCufflinksAnnotationProblem=TRUE,
    quiet=TRUE
)
testSwitchList

aSwitchList <- importGTF(pathToGTF = system.file("extdata/example.gtf.gz", package="IsoformSwitchAnalyzeR"))

aSwitchList
#> This switchAnalyzeRlist list contains:
#>  1203 isoforms from 400 genes
#>  1 comparison from 2 conditions (in total 2 samples)
#> 
#> Feature analyzed:
#> [1] "ORFs"

head(aSwitchList,2)
#>            iso_ref          gene_ref     isoform_id     gene_id
#> 3 isoComp_00000001 geneComp_00000001 TCONS_00000001 XLOC_000001
#> 2 isoComp_00000002 geneComp_00000001 TCONS_00000002 XLOC_000001
#>    condition_1  condition_2 gene_name gene_biotype iso_biotype class_code
#> 3 plaseholder1 plaseholder2      <NA>           NA          NA          =
#> 2 plaseholder1 plaseholder2      <NA>           NA          NA          =
#>   gene_overall_mean gene_value_1 gene_value_2 gene_stderr_1 gene_stderr_2
#> 3                 0            0            0            NA            NA
#> 2                 0            0            0            NA            NA
#>   gene_log2_fold_change gene_p_value gene_q_value iso_overall_mean
#> 3                     0            1            1                0
#> 2                     0            1            1                0
#>   iso_value_1 iso_value_2 iso_stderr_1 iso_stderr_2 iso_log2_fold_change
#> 3           0           0           NA           NA                    0
#> 2           0           0           NA           NA                    0
#>   iso_p_value iso_q_value IF_overall IF1 IF2 dIF isoform_switch_q_value
#> 3           1           1         NA  NA  NA  NA                     NA
#> 2           1           1         NA  NA  NA  NA                     NA
#>   gene_switch_q_value PTC
#> 3                  NA  NA
#> 2                  NA  NA
head(aSwitchList$conditions,2)
#>      condition nrReplicates
#> 1 plaseholder1            1
#> 2 plaseholder2            1

Filtering

Once you have a switchAnalyzeRlist, there is a good chance that it contains a lot of genes/isoforms that are irrelevant for an analysis of isoform switches. Examples of such could be single isoform genes or non-expressed isoforms. These extra genes/isoforms will make the downstream analysis take (much) longer than necessary. Therefore we have implemented a pre-filtering step to remove these features before continuing with the analysis. Importantly, filtering can enhance the reliability of the downstream analysis as described in detail below.

By using preFilter() it is possible to remove genes and isoforms from all aspects of the switchAnalyzeRlist by filtering on:

• Multi-isoform genes
• Gene expression
• Isoform expression
• Isoform Fraction (isoform usage)
• Unwanted isoform classes
• Unwanted gene biotypes
• Genes without differential isoform usage

Removal of single isoform genes is the default setting in preFilter() since these genes, per definition, cannot have changes in isoform usage.

Filtering on isoform expression allows removal of non-used isoforms that only appear in the switchAnalyzeRlist because they were in the isoform/gene annotation used. Furthermore, the expression filtering allows removal of lowly expressed isoforms where the expression levels might be untrustworthy. The filtering on gene expression allows for removal of lowly expressed genes which causes the calculations of the Isoform Fractions (IF) to become untrustworthy.

The filter on Isoform Fraction allows removal of isoforms that only contribute minimally to the gene expression thereby speeding up and simplifying the rest of the downstream analysis without losing important isoforms.

Filtering on unwanted isoform classes is only an option available if using data form Cufflinks/Cuffdiff since the Tuxedo workflow includes classification of transcript types. This allows for removal of, for example, transcripts classified as “Possible polymerase run-on fragment” or “Repeat”. See full description of Cufflinks/Cuffdiff transcript classification here

The preFilter() function is used as follows:

data("exampleSwitchList")

exampleSwitchListFiltered <- preFilter(
    switchAnalyzeRlist = exampleSwitchList,
    geneExpressionCutoff = 1,
    isoformExpressionCutoff = 0,
    removeSingleIsoformGenes = TRUE
)
#> The filtering removed 52 ( 20.08% of ) transcripts. There is now 207 isoforms left

exampleSwitchListFilteredStrict <- preFilter(
    switchAnalyzeRlist = exampleSwitchList,
    geneExpressionCutoff = 10,
    isoformExpressionCutoff = 3,
    removeSingleIsoformGenes = TRUE
)
#> The filtering removed 87 ( 33.59% of ) transcripts. There is now 172 isoforms left

Back to Table of Content

Identifying Isoform Switches

As identification of isoform switches are essential for IsoformSwitchAnalyzeR, multiple ways of identifying isoform switches are supported. Currently IsoformSwitchAnalyzeR directly supports three different approaches:

• DEXSeq : Using DEXSeq (state-of-the art) to test for differential isoform usage. DEXSeq was originally designed for testing for differential exon usage but have recently been shown to perform exceptionally well for differential isoform usage. The default and recommend in IsoformSwitchAnalyzeR. See Testing for Isoform Switches via DEXSeq.
• DRIMSeq : Using the DRIMSeq package. See Testing for Isoform Switches via DRIMSeq.
• Cufflinks/Cuffdiff : the splicing analysis performed by Cuffdiff is automatically imported when importing Cufflinks/Cuffdiff data. See Importing Data from Cufflinks/Cuffdiff.

However, results of other tools which can test for differential isoform usage at gene and/or isoform level can also be used. See below in Testing for Isoform Switches with other Tools.

No matter which method is used, all the downstream functionality in IsoformSwitchAnalyzeR (e.g. not only the test) uses two parameters to define a significant isoform switch:

1. The alpha argument indicating the FDR corrected P-value (Q-value) cutoff.
2. The dIFcutoff argument indicating the minimum (absolute) change in isoform usage (dIF).

This combination is used since a Q-value is only a measure of the statistical certainty of the difference between two groups and thereby does not reflect the effect size which is measured by the dIF values.

IsoformSwitchAnalyzeR supports tests for differential isoform usage performed both at isoform and gene resolution. Note, however, that for gene-level analysis you rely on cutoffs on dIF values for identifying the isoforms involved in a switch — something that, when compared to the isoform-level analysis, could give false positives. We therefore recommend the use of isoform-level analysis whenever possible.

# Load example data and pre-filter
data("exampleSwitchList")
exampleSwitchList <- preFilter(exampleSwitchList) # preFilter to remove lowly expressed features
#> The filtering removed 52 ( 20.08% of ) transcripts. There is now 207 isoforms left

# Perform test
exampleSwitchListAnalyzed <- isoformSwitchTestDEXSeq(
    switchAnalyzeRlist = exampleSwitchList,
    reduceToSwitchingGenes=TRUE
)
#> Step 1 of 3: Calculating isoform abundances (from counts)...
#> Step 2 of 3: Calculating Isoform Fraction matrix...
#> Step 3 of 3: Testing each pairwise comparisons with DEXSeq (this might be a bit slow)...
#>     Estimated time (for dataset with ~30.000 isoforms): 1.4 min
#> Step 4 of 3: Integrating result into switchAnalyzeRlist...
#>     Isoform switch analysis was performed for 62 gene comparisons (100%).
#> Done

# Summarize switching features
extractSwitchSummary(exampleSwitchListAnalyzed)
#>            Comparison nrIsoforms nrGenes
#> 1 hESC vs Fibroblasts         87      40

# Load example data and pre-filter
data("exampleSwitchList")
exampleSwitchList <- preFilter(exampleSwitchList, geneExpressionCutoff = 4) # preFilter for fast runtime
#> The filtering removed 69 ( 26.64% of ) transcripts. There is now 190 isoforms left

# Perform test (takes ~10 sec)
exampleSwitchListAnalyzed <- isoformSwitchTestDRIMSeq(
    switchAnalyzeRlist = exampleSwitchList,
    testIntegration='isoform_only',
    reduceToSwitchingGenes=TRUE
)
#> Step 1 of 6: Creating DM data object...
#> Step 2 of 6: Filtering DM data object...
#> Step 3 of 6: Estimating precision paramter (this may take a while)...
#> Step 4 of 6: Fitting linear models (this may take a while)...
#> Step 5 of 6: Testing pairwise comparison(s)...
#> Step 6 of 6: Preparing output...
#> Result added switchAnalyzeRlist
#> An isoform switch analysis was performed for 56 gene comparisons (98.2%).
#> Done

# Summarize switching features
extractSwitchSummary(exampleSwitchListAnalyzed)
#>            Comparison nrIsoforms nrGenes
#> 1 hESC vs Fibroblasts         76      35

Analyzing Open Reading Frames

Once the isoform switches have been found, the next step is to annotate the isoforms involved in the isoform switches. The first step in such annotation is to obtain Open Reading Frames (ORFs).

If known annotated isoforms were quantified (i.e. you did not perform a (guided) de novo isoform reconstruction) we recommend that you use the annotated CDS (Coding Sequence) from the annotation — which can be imported and integrated though one of the implemented import methods, see Importing Data Into R.

If you performed (guided) de-novo isoform reconstruction (isoform deconvolution), or just want to re-do them, prediction of ORFs can be done with high accuracy from the transcript sequence alone (see supplementary data in Vitting-Seerup et al 2017 for benchmark). Such prediction can be done with:

analyzeORF()

This function relies on analyzing the transcript nucleotide sequences which can be obtained in two ways:

1. If you quantified known annotated isoforms were quantified (i.e. you did not perform a (guided) de novo isoform reconstruction) we recommend that you supply the transcript fasta file while constructing the switchAnalyzeRlist (see above in Importing Data Into R). This will result in more accurate and much faster analysis in the entire workflow.

2. Alternatively the analyzeORF() function can utilizes the genomic coordinates of each transcript to extract the transcript nucleotide sequence from a reference genome (supplied via the genomeObject argument). Such genomic sequences are available in Bioconductor for most model organisms as BSgenome object and can be found here. If necessary check out the BSGenomes for non-model organisms section. Importantly this option should only be used if the isoforms quantified include guide/de-novo reconstructed transcripts (such as produced by Cufflinks/StringTie), else we recommend the first option.

If sequences are extracted from a BSgenome object they will also be added to the switchAnalyzeRlist to facilitate internal sequence analysis and used in all downstream analysis. An example of the internal sequence analysis is a pairwise comparison of the predicted amino acid sequence of the CDS/ORFs in two switching isoforms (see Predicting Switch Consequences below).

In analyzeORF four different methods for predicting the ORF are implemented, suitable for different purposes and circumstances. The four methods are:

1. The ‘longest’ method. This method identifies the longest ORF based on finding the canonical start and stop codons in the transcript nucleotide sequence. This approach is what the CPC2 and CPAT tool uses in its analysis of coding potential. This is the default as it is the most common use case and has the highest accuracy in benchmarks against known annotated ORFs (>90% accuracy against GENCODE data, Vitting-Seerup et al 2017).
2. The ‘mostUpstream’ method. This method identifies the most upstream ORF based on finding the canonical start and stop codons in the transcript nucleotide sequence.
3. The ‘longestAnnotated’ method. This method identifies the longest ORF downstream of an annotated translation start site. It requires known translation start sites are supplied to the cds argument.
4. The ‘mostUpstreamAnnoated’ method. This method identifies ORF downstream of the most upstream overlapping annotated translation start site. It requires known translation start sites are supplied to the cds argument.

One important argument is the minORFlength, which will ensure that only ORFs longer than minORFlength nucleotides are annotated. Besides predicting the ORF, information about the most likely stop codon also allows for prediction of Pre-mature Termination Codon (PTC) and thereby sensitivity to degradation via the Nonsense Mediated Decay (NMD) pathway. This analysis is also implemented in the analyzeORF() function and is controlled by the PTCDistance argument.

Note that the ORF prediction can be integrated with both the CPC2/CPAT results (via the removeNoncodinORFs parameter, see ?analyzeCPAT or ?analyzeCPC2) as well as Pfam results (see Augmenting ORF Predictions with Pfam Results)

The analyzeORF() function can be used as follows:

### This example relies on the example data from the 'Usage of The Statistical Test' section above 

exampleSwitchListAnalyzed <- analyzeORF(
    exampleSwitchListAnalyzed,
    orfMethod = "longest",
    # genomeObject = Hsapiens, # not necessary since sequences are already stored in the switchAnalyzeRlist
    showProgress=FALSE
)
#> Step 1 of 3 : Extracting transcript sequences...
#> Step 2 of 3 : Locating ORFs of interest...
#> Step 3 of 3 : Scanning for PTCs...
#> 120 putative ORFs were identified and analyzed
#> Done

head(exampleSwitchListAnalyzed$orfAnalysis, 3)
#>       isoform_id orfTransciptStart orfTransciptEnd orfTransciptLength
#> 1 TCONS_00000007               248            1414               1167
#> 2 TCONS_00000008               225            1391               1167
#> 3 TCONS_00000009               225            1391               1167
#>   orfStarExon orfEndExon orfStartGenomic orfEndGenomic
#> 1           2          3          324343        325602
#> 2           2          3          324343        325602
#> 3           2          3          324343        325602
#>   stopDistanceToLastJunction stopIndex   PTC
#> 1                      -1164         0 FALSE
#> 2                      -1164         0 FALSE
#> 3                       1336        -1  TRUE

As seen above, genomic and transcript coordinates for ORF beginning and end, as well as PTC status, are added to the ‘orfAnalysis’ entry of the switchAnalyzeRlist.

If the user wants to use the longestAnnotated or mostUpstreamAnnoated methods, the analyzeORF() function requires known CDS to be supplied as described above. The CDS must be stored as a CDSSeq (see ?CDSSet) and a wrapper for downloading the CDS of the most frequently used datasets from UCSC genome browser is implemented in:

getCDS()

Extracting Nucleotide and Amino Acid Sequences

Since we know the CDS/ORF of a transcript, we can obtain the corresponding protein amino acid sequence of the ORF simply by translating the nucleotide sequence of the ORF into amino acids. This opens the possibility for performing both internal and external sequence analysis which enables us to annotate the isoforms involved in isoform switches even further. To facilitate this we have implemented

extractSequence()

which can write fasta files (one file per sequence type) to the hard disk from where they can be used with the external analysis tools.

Please note:

• If you are planning on using the webserver version of the external sequence analysis tools you should set filterAALength=TRUE and alsoSplitFastaFile=TRUE. These arguments exist 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 those arguments are true IsoformSwitchAnalyzeR will generate multiple filtered amino acid fasta files which can then be supplied as multiple individual analysis at the websites.
• If you have quantified known annotated isoforms were quantified (i.e. you did not perform a (guided) de novo isoform reconstruction) we recommend that you supply the transcript fasta file while constructing the switchAnalyzeRlist (see above in Importing Data Into R). This will result in more accurate and much faster analysis in the entire workflow.
• Alternatively the analyzeORF() function can utilizes the genomic coordinates of each transcript to extract the transcript nucleotide sequence from a reference genome (supplied via the genomeObject argument). Such genomic sequences are available in Bioconductor for most model organisms as BSgenome object and can be found here. In necessary check out BSGenomes for non-model organisms. Importantly this option should only be used if the isoforms quantified include guide/de-novo reconstructed transcripts (such as produced by Cufflinks/StringTie), else we recommend the first option.

If sequences are extracted from a BSgenome object they will also be added to the switchAnalyzeRlist to facilitate internal sequence analysis and used in all downstream analysis. An example of the internal sequence analysis is a pairwise comparison of the CDS/ORFs in two switching isoforms (see Predicting Switch Consequences below).

### This example relies on the example data from the 'Analyzing Open Reading Frames' section above 

exampleSwitchListAnalyzed <- extractSequence(
    exampleSwitchListAnalyzed, 
    pathToOutput = '<insert_path>',
    writeToFile=FALSE # to avoid output when running this example data
)
#> Step 1 of 3 : Extracting transcript nucleotide sequences...
#> Step 1 of 3 : Extracting ORF AA sequences...
#> Done

Now you are ready for the next step - running the external sequence analysis tools.

Back to Table of Content

Advice for Running External Sequence Analysis Tools and Downloading Results

The two fasta files that are the output by extractSequence() (if writeToFile=TRUE) can be used as input, to among, others:

• CPAT : The Coding-Potential Assessment Tool which is a tool for predicting whether an isoform is coding or not. CPAT can be run either locally or via their webserver. This is an alternative to CPC2.
• CPC2 : The Coding Potential Calculator 2 which is a tool for predicting whether an isoform is coding or not. CPC2 can be run either locally or via their webserver. This is an alternative to CPAT.
• Pfam : Prediction of protein domains, which can be run either locally (using the pfam_scan.pl script as described in the readme found here or via their webserver.
• SignalP : Prediction of Signal Peptides — a short N-terminal sequence of a peptide indicating where a protein should be membrane bound or secreted. SignalP can be run either locally or via their webserver (V5 is supported).
• NetSurfP-2 : Prediction of Intrinsically Disordered Regions (IDR) — regions of a protein which does not have a fixed three-dimensional structure (opposite protein domains). NetSurfP-2 can be run either locally or via their webserver.

Please note that:

• CPAT and CPC2 performs the same type of analysis so it is only necessary to run one of them.
• Some webservers (Pfam and SignalP) have (quite) strict limits on the number (and length) of amino acid sequences one can submit at the time. As indicated above extractSequence() can generate multiple filtered amino acid fasta files that can be uploaded one at the time (aka as multiple runs) to circumvent this.
• These five tools are the ones currently supported, but if you have additional ideas please do not hesitate to contact us as described in How To Get Help.

We suggest that the external sequence analysis tools are run as follows:

• CPAT : Use default parameters and the nucleotide fasta file (_nt.fasta). If the webserver was used, download the tab-delimited result file (available at the bottom of the result page). If a stand-alone version was used, just supply the path to the result file.
• CPC2 : Use default parameters and if required select the most similar species. If the webserver (batch submission) was used, download the tab-delimited result file (via the “Download the result” button). If a stand-alone version was just just supply the path to the result file.
• Pfam : Use default parameters and the amino acid fasta file (_AA.fasta). If the webserver is used you need to copy/paste the result part of the mail you receive into an empty plain text document (notepad, sublimetext, TextEdit or similar (not Word)) and save that to a plain text (txt) file. The path to that file should be supplied. If a stand-alone version was used, just supply the path to the result file. A more detailed walkthrough is found under details in the documentation of the analyzePFAM() function (?analyzePFAM).
• SignalP : Use the amino acid fasta file (_AA.fasta). If using the webserver SignalP should be run with the parameter “Short output (no figures)” under “Output format” and one should select the appropriate “Organism group”. When using a stand-alone version SignalP should be run with the ‘-f summary’ option. If using the webserver the results can be downloaded using the “Downloads” button in the top-right corner where the user should select “Prediction summary” and supply the path to the resulting file to the “pathToSignalPresultFile” argument. If a stand-alone version was just supply the path to the summary result file.
• NetSurfP-2 : 1) Go to webserver. 2) Upload the amino acid file (_AA) created with extractSequence() function. 3) Submit your job. 4) Wait till job is finished (if you submit your email you will receive a notification). 5) In the top-right corner of the result site use the “Export All” button to download the results as a CNV file. 6) Supply a string indicating the path to the downloaded csv file directly to the “pathToNetSurfP2resultFile” argument. If you run NetSurfP-2 locally just use the “–csv” argument and provide the resulting csv file to the pathToNetSurfP2resultFile argument.

Back to Table of Content

Importing External Sequence Analysis

After the external sequence analysis with CPC2/CPAT, Pfam, and SignalP have been performed (please remember to cite as describe in What To Cite), the results can be extracted and incorporated in the switchAnalyzeRlist via (respectively)

analyzeCPAT() # OR
analyzeCPC2()
analyzePFAM()
analyzeSignalP()
analyzeNetSurfP2()

Please note that CPC2 and CPAT performs the same type of analysis so only one of them is necessary.

The functions are used as:

### Load test data (matching the external sequence analysis results)
data("exampleSwitchListIntermediary")
exampleSwitchListIntermediary
#> This switchAnalyzeRlist list contains:
#>  162 isoforms from 40 genes
#>  1 comparison from 2 conditions (in total 6 samples)
#> 
#> Switching features:
#>            Comparison switchingIsoforms switchingGenes
#> 1 hESC vs Fibroblasts                83             40
#> 
#> Feature analyzed:
#> [1] "Isoform Switch Identification, ORFs, ntSequence, aaSequence"

### Add CPAT analysis
exampleSwitchListAnalyzed <- analyzeCPAT(
    switchAnalyzeRlist   = exampleSwitchListIntermediary,
    pathToCPATresultFile = system.file("extdata/cpat_results.txt", package = "IsoformSwitchAnalyzeR"),
    codingCutoff         = 0.725, # the coding potential cutoff we suggested for human
    removeNoncodinORFs   = TRUE   # because ORF was predicted de novo
)
#> Added coding potential to 162 (100%) transcripts

### Add CPC2 analysis
exampleSwitchListAnalyzed <- analyzeCPC2(
    switchAnalyzeRlist   = exampleSwitchListAnalyzed,
    pathToCPC2resultFile = system.file("extdata/cpc2_result.txt", package = "IsoformSwitchAnalyzeR"),
    codingCutoff         = 0.725, # the coding potential cutoff we suggested for human
    removeNoncodinORFs   = TRUE   # because ORF was predicted de novo
)
#> Added coding potential to 162 (100%) transcripts

### Add PFAM analysis
exampleSwitchListAnalyzed <- analyzePFAM(
    switchAnalyzeRlist   = exampleSwitchListAnalyzed,
    pathToPFAMresultFile = system.file("extdata/pfam_results.txt", package = "IsoformSwitchAnalyzeR"),
    showProgress=FALSE
)
#> Converting AA coordinats to transcript and genomic coordinats...
#> Added domain information to 124 (76.54%) transcripts

### Add SignalP analysis
exampleSwitchListAnalyzed <- analyzeSignalP(
    switchAnalyzeRlist       = exampleSwitchListAnalyzed,
    pathToSignalPresultFile  = system.file("extdata/signalP_results.txt", package = "IsoformSwitchAnalyzeR")
)
#> Added signal peptide information to 17 (10.49%) transcripts

### Add NetSurfP2 analysis
exampleSwitchListAnalyzed <- analyzeNetSurfP2(
    switchAnalyzeRlist        = exampleSwitchListAnalyzed,
    pathToNetSurfP2resultFile = system.file("extdata/netsurfp2_results.csv.gz", package = "IsoformSwitchAnalyzeR"),
    showProgress = FALSE
)
#> Step 1 of 3: Reading results into R...
#> Step 2 of 3: Analyzing data to extract IDRs...
#> Step 3 of 3: Converting AA coordinats to transcript and genomic coordinats...
#> Added IDR information to 94 (58.02%) transcripts

exampleSwitchListAnalyzed
#> This switchAnalyzeRlist list contains:
#>  162 isoforms from 40 genes
#>  1 comparison from 2 conditions (in total 6 samples)
#> 
#> Switching features:
#>            Comparison switchingIsoforms switchingGenes
#> 1 hESC vs Fibroblasts                83             40
#> 
#> Feature analyzed:
#> [1] "Isoform Switch Identification, ORFs, ntSequence, aaSequence, Protein Domains, Signal Peptides, IDR, Coding Potential"

Here, the pathTo<tool_name>resultFile points to the result files constructed as suggested in Advice for Running External Sequence Analysis Tools and Downloading Results and in the detailed documentation of each of the individual functions. The part of the examples using system.file() is not necessary when using your own data - just supply the path as a string (e.g. “/myFiles/myCPATresults.txt”)

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 to the functions above and IsoformSwitchAnalyzeR will read them all and integrate them.

Of particular interest is the removeNoncodinORFs argument in analyzeCPAT() and analyzeCPC2() since it allows integration of the ORF analysis with the CPAT/CPC2 analyses by removing ORFs from isoforms not predicted to be coding by CPAT/CPC2. This can be particularly useful if isoforms and ORFs have been predicted de novo (both guided or non-guided). Note that if enabled (by setting to TRUE), it will affect all downstream analyses and plots as both analysis of protein domains and signal peptides require that ORFs are annotated (e.g. analyzeSwitchConsequences will not consider the domains (potentially) found by Pfam if the ORFs have been removed).

Back to Table of Content

Predicting Alternative Splicing

Another type of annotation we easily can obtain, since we know the exon structure of all isoforms in a given gene (with isoform switching), is alternative splicing. Here intron retention events are of particular interest as a consequence in isoform switches since they represent the largest changes in isoforms. Identification of alternative splicing, alternative transcription start and termination sites can be obtained via the analyzeAlternativeSplicing().

### This example relies on the example data from the 'Importing External Sequence Analysis' section above 

exampleSwitchListAnalyzed <- analyzeAlternativeSplicing(exampleSwitchListAnalyzed, quiet=TRUE)

### overview of number of intron retentions (IR)
table(exampleSwitchListAnalyzed$AlternativeSplicingAnalysis$IR)
#> 
#>   0   1   2   3   5 
#> 135  20   3   3   1

Meaning 20 isoforms contain one intron retention (IR) and 7 isoforms each contain two or more intron retentions.

IsoformSwitchAnalyzeR also supports many types of genome wide analysis of alternative splicing - for a more thorough walkthrough of the analysis of alternative splicing see the Analyzing Alternative Splicing workflow.

Back to Table of Content

Predicting Switch Consequences

If an isoform has a significant change in its contribution to gene expression, there must per definition be reciprocal changes in one (or more) isoforms in the opposite direction, compensating for the change in the first isoform. We utilize this by extracting the isoforms that are significantly differentially used and compare them to the isoforms that are compensating. Using all the information gathered through the workflow described above, the annotation of the isoform(s) used more (positive dIF) can be compared to the isoform(s) used less (negative dIF) and by systematically identify differences annotation we can identify potential function consequences of the isoform switch.

Specifically, IsoformSwitchAnalyzeR contains a function analyzeSwitchConsequences() which extracts the isoforms with significant changes in their isoform usage (defined by the alpha and dIFcutoff parameters, see Identifying Isoform Switches for details) and the isoform, with a large opposite change in isoform usage (also controlled via the dIFcutoff parameters) that compensate for the changes. Note that if an isoform-level test was not used, the gene is require to be significant (defined by the alpha parameter); but, isoforms are then selected purely based on their changes in dIF values.

These isoforms are then divided into the isoforms that increase their contribution to gene expression (positive dIF values larger than dIFcutoff) and the isoforms that decrease their contribution (negative dIF values smaller than -dIFcutoff). The isoforms with increased contribution are then (in a pairwise manner) compared to the isoform with decreasing contribution. In each of these comparisons the isoforms compared are analyzed for differences in their annotation (controlled by the consequencesToAnalyze parameter). Currently 22 different features of the isoforms can be compared, which include features such as intron retention, coding potential, NMD status, protein domains and the sequence similarity of the amino acid sequence of the annotated ORFs. For a full list, see “Details” on the documentation page for analyzeSwitchConsequences() (by running ?analyzeSwitchConsequences).

A more strict analysis can be performed by enabling the onlySigIsoforms argument, which causes analyzeSwitchConsequences() to only consider significant isoforms (defined by the alpha and dIFcutoff parameters) meaning the compensatory changes in isoform usage are ignored unless they themselves are significant.

The analysis of consequences can be performed as follows:

### This example relies on the example data from the 'Predicting Alternative Splicing' section above 

# the consequences highlighted in the text above
consequencesOfInterest <- c('intron_retention','coding_potential','NMD_status','domains_identified','ORF_seq_similarity')

exampleSwitchListAnalyzed <- analyzeSwitchConsequences(
    exampleSwitchListAnalyzed,
    consequencesToAnalyze = consequencesOfInterest, 
    dIFcutoff = 0.4, # very high cutoff for fast runtimes
    showProgress=FALSE
)
#> Step 1 of 4: Extracting genes with isoform switches...
#> Step 2 of 4: Analyzing 5 pairwise isoforms comparisons...
#> Step 3 of 4: Massaging isoforms comparisons results...
#> Step 4 of 4: Preparing output...
#> Identified  genes with containing isoforms switching with functional consequences...

extractSwitchSummary(exampleSwitchListAnalyzed, dIFcutoff = 0.4, filterForConsequences = FALSE)
#>            Comparison nrIsoforms nrGenes
#> 1 hESC vs Fibroblasts         12       7
extractSwitchSummary(exampleSwitchListAnalyzed, dIFcutoff = 0.4, filterForConsequences = TRUE)
#>            Comparison nrIsoforms nrGenes
#> 1 hESC vs Fibroblasts         10       5

Meaning that ~5/7 genes with switches (with dIF > 0.4) have isoform switches with functional consequences. For larger dataset with default parameters this is typical 1/3.

Please note that the analyzeSwitchConsequences() function contains many parameters and cutoffs for deciding when a difference in the annotation is sufficiently different (e.g. supplying ‘ntCutoff=50’ ensures that differences in nucleotide lengths of two ORF are larger than 50 nucleotides — not just 3 nucleotides (one codon)).

Back to Table of Content

Post Analysis of Isoform Switches with Consequences

Now that we have a genome-wide characterization of isoform switches with potential consequences there are a lot of possibilities for analyzing these. IsoformSwitchAnalyzeR currently directly supports two different types of post-analysis (of switches - see Genome-Wide Analysis of Alternative Splicing below for analysis of alternative splicing):

1. Analysis of Individual Isoform Switching which aims to identify the top genes and provide a visual overview summarizing all the relevant information about an isoform switch single gene.
2. Genome-Wide Analysis of Isoform Switching which summarizes and analyzes larger patterns of isoform switches with predicted functional consequences

### Load already analyzed data
data("exampleSwitchListAnalyzed")

### Let's reduce the switchAnalyzeRlist to only one condition
exampleSwitchListAnalyzedSubset <- subsetSwitchAnalyzeRlist(
    exampleSwitchListAnalyzed, 
    exampleSwitchListAnalyzed$isoformFeatures$condition_1 == 'COAD_ctrl'
)
exampleSwitchListAnalyzedSubset
#> This switchAnalyzeRlist list contains:
#>  1000 isoforms from 369 genes
#>  1 comparison from 2 conditions (in total 322 samples)
#> 
#> Switching features:
#>                 Comparison switchingIsoforms switchingGenes
#> 1 COAD_ctrl vs COAD_cancer               692            369
#> 
#> Feature analyzed:
#> [1] "Isoform Switch Identification, ORFs, Protein Domains, Signal Peptides, Alternative splicing, IDR, Switch Consequences, Coding Potential"

### Extract top switching genes (by q-value)
extractTopSwitches(
    exampleSwitchListAnalyzedSubset, 
    filterForConsequences = TRUE, 
    n = 2, 
    sortByQvals = TRUE
)
#>                               gene_ref                     gene_id
#> 1 chr1:244998639..245008357,+COAD_ctrl chr1:244998639..245008357,+
#> 2      chr12:569543..671058,+COAD_ctrl      chr12:569543..671058,+
#>   gene_name condition_1 condition_2 gene_switch_q_value
#> 1    FAM36A   COAD_ctrl COAD_cancer        5.344466e-06
#> 2  B4GALNT3   COAD_ctrl COAD_cancer        5.344466e-06
#>   switchConsequencesGene Rank
#> 1                   TRUE    1
#> 2                   TRUE    2

### Extract top 2 switching genes (by dIF values)
extractTopSwitches(
    exampleSwitchListAnalyzedSubset, 
    filterForConsequences = TRUE, 
    n = 2, 
    sortByQvals = FALSE
)
#>                              gene_ref                    gene_id gene_name
#> 1 chr10:70980059..71027314,+COAD_ctrl chr10:70980059..71027314,+     HKDC1
#> 2  chr3:13573824..13679921,+COAD_ctrl  chr3:13573824..13679921,+     FBLN2
#>   condition_1 condition_2 gene_switch_q_value combinedDIF
#> 1   COAD_ctrl COAD_cancer        7.463277e-06       1.216
#> 2   COAD_ctrl COAD_cancer        5.344466e-06       1.167
#>   switchConsequencesGene Rank
#> 1                   TRUE    1
#> 2                   TRUE    2

### Extract data.frame with all switching isoforms
switchingIso <- extractTopSwitches( 
    exampleSwitchListAnalyzedSubset, 
    filterForConsequences = TRUE, 
    n = NA,                  # n=NA: all features are returned
    extractGenes = FALSE,    # when FALSE isoforms are returned
    sortByQvals = TRUE
)

subset(switchingIso, gene_name == 'ZAK')
#>                 iso_ref                             gene_ref isoform_id
#> 9   uc002uhz.2COAD_ctrl chr2:173940565..174132736,+COAD_ctrl uc002uhz.2
#> 237 uc002uhy.2COAD_ctrl chr2:173940565..174132736,+COAD_ctrl uc002uhy.2
#>                         gene_id gene_name condition_1 condition_2   IF1
#> 9   chr2:173940565..174132736,+       ZAK   COAD_ctrl COAD_cancer 0.176
#> 237 chr2:173940565..174132736,+       ZAK   COAD_ctrl COAD_cancer 0.181
#>       IF2    dIF isoform_switch_q_value switchConsequencesGene
#> 9   0.408  0.232           5.344466e-06                   TRUE
#> 237 0.024 -0.157           5.387221e-04                   TRUE
#>                   comparison Rank
#> 9   COAD_ctrl_vs_COAD_cancer    9
#> 237 COAD_ctrl_vs_COAD_cancer  237

switchPlot(exampleSwitchListAnalyzedSubset, gene = 'ZAK')

switchPlotTopSwitches(
    switchAnalyzeRlist = exampleSwitchListAnalyzed, 
    n = 10,
    filterForConsequences = FALSE, 
    splitFunctionalConsequences = TRUE
)

Other Tools in IsoformSwitchAnalyzeR

Apart from the analysis presented above, IsoformSwitchAnalyzeR also contains a few other useful tools which will be presented in this section.

The central visualization is the switch plot, created with switchPlot() function as shown above. This plot is made from four subplots, which each can be created individually using:

switchPlotTranscript() # Visualizes the transcripts and their annotation
switchPlotGeneExp()    # Visualizes the gene expression
switchPlotIsoExp()     # Visualizes the isoform expression
switchPlotIsoUsage()   # Visualizes the isoform usage

As illustrated here:

### exampleSwitchListAnalyzedSubset is created above

switchPlotTranscript(exampleSwitchListAnalyzedSubset, gene = 'ZAK')

switchPlotGeneExp (exampleSwitchListAnalyzedSubset, gene = 'ZAK')

switchPlotIsoExp(exampleSwitchListAnalyzedSubset, gene = 'TNFRSF1B')

switchPlotIsoUsage(exampleSwitchListAnalyzedSubset, gene = 'ZAK')

Back to Table of Content

Overview of Alternative Splicing Workflow

Analysis of alternative splicing, here made possible by a liftover of functionality from the now deprecated R package spliceR, is an important area of research. IsoformSwitchAnalyzeR is well suited for performing alternative splicing analysis since alternative splicing literally will result in isoform switches. Please note that alternative splicing here is used as a term which covers both Alternative Splicing (AS), Alternative Transcription Start Sites (ATSS, sometimes called alternative first exon) as well as Alternative Transcription Start Sites (ATTS, sometimes called alternative last exon): although this is technically incorrect it is here used for convenience.

The workflow implemented in IsoformSwitchAnalyzeR is mostly focused on identifying global changes in alternative splicing but also allows for analysis of individual isoforms/switches. A typical workflow consists of five steps. Please note that step 1-4 are included in the standard switch analysis workflow, meaning if you have already analyzed isoform switches with consequences you can skip 1-4. The five steps correspond to using the following functions in succession:

### Step (1) Import data and create SwitchAnalyzeRlist
mySwitchList <- importCufflinksData()      # OR
mySwitchList <- importRdata()

### Step (2) Prefilter
mySwitchList <- preFilter( mySwitchList )

### Step (3) Identify differentially used isoforms
mySwitchList <- isoformSwitchTestDEXSeq(  mySwitchList, reduceToSwitchingGenes=FALSE ) # OR
mySwitchList <- isoformSwitchTestDRIMSeq( mySwitchList, reduceToSwitchingGenes=FALSE ) 

### Step (4) Analyze alternative splicing
mySwitchList <- analyzeAlternativeSplicing( mySwitchList )

### Step (5) Global splicing analysis
extractSplicingSummary( mySwitchList )
extractSplicingEnrichment( mySwitchList )
extractSplicingEnrichmentComparison( mySwitchList )
extractSplicingGenomeWide( mySwitchList )

Step 1-4 are shared with the standard switch analysis workflow which is describe above in the Detailed Workflow (specifically each step is described in Importing Data Into R, Filtering, Identifying Isoform Switches and Predicting Alternative Splicing) and will therefore not be repeated here. The only difference is that in some cases one might not want reduce the switchAnalyzeRList to only contain the genes which contain differentially used isoforms — as discussed below. This is done by setting reduceToSwitchingGenes=FALSE in the switch test.

Genome-Wide Analysis of Alternative Splicing

The genome-wide analysis of splicing patterns is facilitated by the following four functions:

• extractSplicingSummary()
• extractSplicingEnrichment()
• extractSplicingEnrichmentComparison()
• extractSplicingGenomeWide()

The extractSplicingSummary function does exactly what it says on the tin - it summarizes the number of events found in each comparison. To showcase it, let’s look at some already analyzed example data (meaning step 1-4 have already been done).

data("exampleSwitchListAnalyzed")
exampleSwitchListAnalyzed
#> This switchAnalyzeRlist list contains:
#>  1469 isoforms from 530 genes
#>  2 comparison from 4 conditions (in total 891 samples)
#> 
#> Switching features:
#>                 Comparison switchingIsoforms switchingGenes
#> 1 COAD_ctrl vs COAD_cancer               692            369
#> 2 LUAD_ctrl vs LUAD_cancer               422            218
#> 3                 combined              1010            530
#> 
#> Feature analyzed:
#> [1] "Isoform Switch Identification, ORFs, Protein Domains, Signal Peptides, Alternative splicing, IDR, Switch Consequences, Coding Potential"

Note that via the in “Feature analyzed” part of the summary we can see that the analysis of alternative splicing has been done and added to the switchAnalyzeRlist.

extractSplicingSummary(
    exampleSwitchListAnalyzed,
    asFractionTotal = FALSE,
    plotGenes=FALSE
)

Note that it is possible to both summarize per gene or per isoform via the plotGenes argument and that the asFractionTotal argument enables analysis the fraction (if TRUE) or Number (if FALSE) of significant isoforms/genes which have a particular splice type.

From the summary above, it is quite clear that some of the alternative splicing events are not equally used - most prominently illustrated by the use of ATSS in COAD, where there is (a lot) more gain than losses. To formally analyze this type of uneven alternative splicing usage we have implemented extractSplicingEnrichment(). This function summarizes the uneven usage within each comparison by for each alternative splicing type calculate the fraction of events being gains (as opposed to loss) and perform a statistical analysis of this fraction.

splicingEnrichment <- extractSplicingEnrichment(
    exampleSwitchListAnalyzed,
    splicingToAnalyze='all',
    returnResult=TRUE,
    returnSummary=TRUE
)

From this plot it is clear to see that some alternative splicing types are preferentially used in the two cancer types, and from the statistical summary also created by extractSplicingEnrichment() (due to the returnSummary=TRUE argument):

subset(splicingEnrichment, splicingEnrichment$AStype == 'ATSS')
#>    condition_1 condition_2 AStype nUp nDown    propUp propUpCiLo
#> 3    COAD_ctrl COAD_cancer   ATSS 229    35 0.8674242  0.8190962
#> 11   LUAD_ctrl LUAD_cancer   ATSS  83    53 0.6102941  0.5226713
#>    propUpCiHi   propUpPval   propUpQval Significant
#> 3   0.9047325 1.534035e-32 2.147649e-31        TRUE
#> 11  0.6915895 1.289234e-02 2.087173e-02        TRUE
#>                    Comparison                             AStype2
#> 3  COAD_ctrl\nvs\nCOAD_cancer ATSS gain \n(paried with ATSS loss)
#> 11 LUAD_ctrl\nvs\nLUAD_cancer ATSS gain \n(paried with ATSS loss)

We can indeed see that ATSS has a very significant skew in COAD. Note that the “splicingToAnalyze” argument implemented in all these summary functions enables analysis of only a subset of the splice types.

Another observation that can be made from the enrichment analysis above is that although the overall trend in usage of alternative splicing is the same the two cancer types there seem to be some differences. These differences can be analyzed using extractSplicingEnrichmentComparison() function which for each splice type, in a pairwise manner contrast the individual condition comparison, to assess whether the ratio of gains are indeed different across comparisons.

extractSplicingEnrichmentComparison(
    exampleSwitchListAnalyzed,
    splicingToAnalyze=c('A3','MES','ATSS'), # Those significant in COAD in the previous analysis
    returnResult=FALSE # Preventing the summary statistics to be returned as a data.frame
)

And indeed only ATSS usage is significantly different between the two cancer types.

All of the above splicing analysis relies on the identification of pairs of isoforms involved in a switch, but the analysis of alternative splicing can also be done from the perspective of the individual splice types. This can be done by asking the question: How does the isoform usage of all isoforms utilizing a particular splicing type change - in other words is all isoforms or only only a subset of isoforms that are affected. This can be analyzed with the extractSplicingGenomeWide() function.

extractSplicingGenomeWide(
    exampleSwitchListAnalyzed,
    featureToExtract = 'all',                 # all isoforms stored in the switchAnalyzeRlist
    splicingToAnalyze = c('A3','MES','ATSS'), # Splice types significantly enriched in COAD
    plot=TRUE,
    returnResult=FALSE  # Preventing the summary statistics to be returned as a data.frame
)

From which it can be seen that all isoforms using ATSS is generally used more, meaning it is a global phenomenon, whereas none of the other splice types are general (i.e. they seem to only happen in a specific subset of the data).

Please note that:

• The the dots in the violin plots above indicate 25th, 50th (median) and 75th percentiles (just like the box in a boxplot would).
• The analysis provided by extractSplicingGenomeWide() should only be expected to give a result when splicing is (severely) affected.
• When using the featureToExtract=‘all’ one should not use the reduceToSwitchingGenes=TRUE during the testing of isoform switches (as then it is not all isoforms but only those in genes containing at least one differentially used isoform) - then it will not be ‘all’ isoforms that can be analyzed.

Back to Table of Content

Augmenting ORF Predictions with Pfam Results

The workflow described here is an extension of the workflow enabled in the analyzeCPAT() function via the “removeNoncodinORFs” argument which if TRUE removes ORF information from the isoforms where the CPAT analysis classifies them as non-coding. Please note this can be done exactly the same way if coding potential were analyzed with analyzeCPC2(). Here we will “rescue” the ORF of isoforms predicted to be noncoding by CPAT/CPC2 but which have a predicted protein domain. This workflow is an alternative approach to setting “removeNoncodinORFs=TRUE” and should be used instead of that option.

Note that this is only recommended if ORFs were predicted (i.e. not imported from a GTF file). After this procedure, isoforms with ORFs will be isoforms with an ORF longer than minORFlength (if specified in analyzeORF) which are predicted to be coding by CPAT/CPC2 OR have a predicted protein domain (by Pfam).

Since the ORF information is stored in the ‘orfAnalysis’ analysis entry of the switchAnalyzeRlist we can remove it (by replacing it with NAs) as follows:

### Test data
data("exampleSwitchListAnalyzed")
exampleSwitchListAnalyzed

### Extract coding isoforms
nonCodingIsoforms <- unique(exampleSwitchListAnalyzed$isoformFeatures$isoform_id[
    which( ! exampleSwitchListAnalyzed$isoformFeatures$codingPotential )
    ])

### Rescue those with protein domains
nonCodingIsoformsRescued <- setdiff(nonCodingIsoforms, exampleSwitchListAnalyzed$domainAnalysis$isoform_id)

# nr rescued
length(nonCodingIsoforms) - length(nonCodingIsoformsRescued)


### Remove noncoding isoforms ORF annotation
sum(is.na(exampleSwitchListAnalyzed$orfAnalysis$orfTransciptStart))

exampleSwitchListAnalyzed$orfAnalysis[
    which( exampleSwitchListAnalyzed$orfAnalysis$isoform_id %in% nonCodingIsoformsRescued), 2:ncol(exampleSwitchListAnalyzed$orfAnalysis)
    ] <- NA
    
exampleSwitchListAnalyzed$isoformFeatures$PTC[which(exampleSwitchListAnalyzed$isoformFeatures$isoform_id %in% nonCodingIsoforms)] <- NA

sum(is.na(exampleSwitchListAnalyzed$orfAnalysis$orfTransciptStart))

Back to Table of Content

Analyze Small Upstream ORFs

Recent research suggests that small upstream ORFs are far more frequent that previously assumed. It is therefore of particular interest to start analyzing these and here we have indirectly presented a tool which can do just so: analyzeORF().

Here we show how one could start such an analysis of small upstream ORFs.

# run ORF analysis on longest ORF
exampleSwitchListAnalyzed <- analyzeORF(exampleSwitchListAnalyzed, genomeObject = Hsapiens, method='longest')
mean(exampleSwitchListAnalyzed$orfAnalysis$orfTransciptLength)

# run ORF analysis on most upstream ORF
exampleSwitchListAnalyzed2 <- analyzeORF(exampleSwitchListAnalyzed, orfMethod = 'mostUpstream', minORFlength = 50)
mean(exampleSwitchListAnalyzed2$orfAnalysis$orfTransciptLength)

# calculate pairwise difference
summary(
    exampleSwitchListAnalyzed2$orfAnalysis$orfTransciptLength -
        exampleSwitchListAnalyzed$orfAnalysis$orfTransciptLength[
            match(
                exampleSwitchListAnalyzed2$orfAnalysis$isoform_id,
                exampleSwitchListAnalyzed$orfAnalysis$isoform_id
            )
        ]
)

Back to Table of Content

Remove Sequences Stored in SwitchAnalyzeRlist

The sequences stored in the SwitchAnalyzeRlist are not needed after consequences have been predicted and can thereby be removed, reducing the object size as well as loading/saving times. This has been done for the example dataset ‘exampleSwitchListAnalyzed’. This is simply done as follows

summary(exampleSwitchListIntermediary)

exampleSwitchListIntermediary$ntSequence <- NULL
exampleSwitchListIntermediary$aaSequence <- NULL

summary(exampleSwitchListIntermediary)

Back to Table of Content

Adding Uncertain Category to Coding Potential Predictions

There has been considerable debate about whether the default parameters for the coding potential calculated via CPAT are too lenient (too low). Although CPC2 is a more recent tool it suffers from the same problem. In fact this will always be a problem when having a single cutoff on a coding potential value. One possible solution is to introduce an “unknown” class with medium size coding potential which we can then disregard.

Let’s start by looking at the distribution of coding potential values.

data("exampleSwitchListAnalyzed")
hist(exampleSwitchListAnalyzed$isoformFeatures$codingPotentialValue)

These coding potential values are summarized by the cutoff supplied in the ‘codingPotential’ column — which is what IsoformSwitchAnalyzeR uses in the downstream analysis.

### These are summarized by the cutoff supplied in the 'codingPotential' column --- which is what is used by IsoformSwitchAnalyzeR in the downstream analysis.
table(exampleSwitchListAnalyzed$isoformFeatures$codingPotential, exclude = NULL)

By simply setting the mid-range values to NA it will cause IsoformSwitchAnalyzeR to ignore them, thereby removing the isoforms with “unknown” coding potential. This can be done as follows:

exampleSwitchListAnalyzed$isoformFeatures$codingPotential <- NA
exampleSwitchListAnalyzed$isoformFeatures$codingPotential[which(exampleSwitchListAnalyzed$isoformFeatures$codingPotentialValue > 0.75)] <- TRUE
exampleSwitchListAnalyzed$isoformFeatures$codingPotential[which(exampleSwitchListAnalyzed$isoformFeatures$codingPotentialValue < 0.25)] <- FALSE

table(exampleSwitchListAnalyzed$isoformFeatures$codingPotential, exclude = NULL)

Back to Table of Content

Quality control of ORF of known annotation

As we have shown there are many problems with known CDS annotations (see supplementary data in Vitting-Seerup et al 2017) we have implemented two ways to ensure a high quality of the CDS imported from a GTF annotation file:

1. When you import the GTF file (via importRdata() or importGTF()) you can enable the ‘onlyConsiderFullORF’ argument which makes sure you only add the ORF stored in the GTF file if it is annotated with both a start and stop codon.
2. When you extract the biological sequences (nucleotide and amino acid) with extractSequence() you can enable the argument ‘removeORFwithStop’ which will remove ORFs that contain un-annotated stop codons, defined as * when the ORF nucleotide sequences are translated to the amino acid sequence. If enabled, the ORFs will both be removed from the switchAnalyzeRList and from the sequences outputted by extractSequence().

Analyzing the Biological Mechanisms Behind Isoform Switching

The difference between the isoforms involved in an isoform switch can arise by changes in three distinct biological mechanisms:

1. Alternative Transcription Start Site (aTSS)
2. Alternative Splicing (AS)
3. Alternative Transcription Termination Site (aTTS)

Since we how the structure of the isoforms involved in an isoform switch we can also analyze which (combination) of these biological mechanisms gives rise to the difference between the two isoforms involved in an isoform switch.

This is simply done by, in addition to running analyzeSwitchConsequences with the consequences you find interesting, making a separate consequence analysis of consequences (also with analyzeSwitchConsequences) where the consequences you analyze (supplied to the consequencesToAnalyze argument) are:

1. ‘tss’ — which will analyze the isoforms for aTSS
2. ‘intron_structure’ — which will analyze the isoforms for AS
3. ‘tts’ — which will analyze the isoforms for aTSS

Then we can simply compare the result of this analysis to the isoform switches with consequences we already have identified to be of interest, and thereby identify which biological mechanisms give rise to the isoform switches with consequence you are interested in. One suggestion for such an analysis is illustrated here:

### Load example data
data("exampleSwitchListAnalyzed")

### Reduce datasize for fast runtime
selectedGenes <- unique(exampleSwitchListAnalyzed$isoformFeatures$gene_id)[50:100]
exampleSwitchListAnalyzedSubset <- subsetSwitchAnalyzeRlist(
    exampleSwitchListAnalyzed,
    exampleSwitchListAnalyzed$isoformFeatures$gene_id %in% selectedGenes
)

### analyze the biological mechanisms
bioMechanismeAnalysis <- analyzeSwitchConsequences(
    exampleSwitchListAnalyzedSubset, 
    consequencesToAnalyze = c('tss','tts','intron_structure'),
    showProgress = FALSE
)$switchConsequence # only the consequences are interesting here

### subset to those with differences
bioMechanismeAnalysis <- bioMechanismeAnalysis[which(bioMechanismeAnalysis$isoformsDifferent),]

### extract the consequences of interest already stored in the switchAnalyzeRlist
myConsequences <- exampleSwitchListAnalyzedSubset$switchConsequence
myConsequences <- myConsequences[which(myConsequences$isoformsDifferent),]
myConsequences$isoPair <- paste(myConsequences$isoformUpregulated, myConsequences$isoformDownregulated) # id for specific iso comparison

### Obtain the mechanisms of the isoform switches with consequences
bioMechanismeAnalysis$isoPair <- paste(bioMechanismeAnalysis$isoformUpregulated, bioMechanismeAnalysis$isoformDownregulated)
bioMechanismeAnalysis <- bioMechanismeAnalysis[which(bioMechanismeAnalysis$isoPair %in% myConsequences$isoPair),]  # id for specific iso comparison

This result is best summarized in a Venn diagram:

### Create list with the isoPair ids for each consequence
AS   <- bioMechanismeAnalysis$isoPair[ which( bioMechanismeAnalysis$featureCompared == 'intron_structure')]
aTSS <- bioMechanismeAnalysis$isoPair[ which( bioMechanismeAnalysis$featureCompared == 'tss'             )]
aTTS <- bioMechanismeAnalysis$isoPair[ which( bioMechanismeAnalysis$featureCompared == 'tts'             )]

mechList <- list(
    AS=AS,
    aTSS=aTSS,
    aTTS=aTTS
)

### Create Venn diagram
library(VennDiagram)
#> Loading required package: grid
#> Loading required package: futile.logger
myVenn <- venn.diagram(mechList, col='transparent', alpha=0.4, fill=RColorBrewer::brewer.pal(n=3,name='Dark2'), filename=NULL)

### Plot the venn diagram
grid.newpage() ; grid.draw(myVenn)

From this, we can see the relative importance of each of the three mechanisms, as well as the combination of these.

Back to Table of Content

Analyzing experiments without replicates

Generally we do not recommend analysis of experiments without replicates as it is impossible to assess technical and biological variation. On the other hand, single-replicate pilot studies can be useful for planning larger experiments: hence this small guide. For running IsoformSwitchAnalyzeR without replicates, we have to solely rely on the dIF values to call isoform switches. This can be done by following the steps outlined in Detailed Workflow with a few modifications:

• Omit the isoform switch test step
• Manually set the switch q-value to 1, which can be done as follows:

aSwitchList$isoformFeatures$gene_switch_q_value <- 1

• In all subsequent steps of the Detailed Workflow, make sure IsoformSwitchAnalyzeR ignores the q-value by setting the alpha argument to something larger than 1 (e.g. 1.1). You will have to ignore the resulting warnings IsoformSwitchAnalyzeR produces as it tests the validity of the arguments supplied by the user.

Back to Table of Content

FAQ Table of Content

• What Quantification Tool(s) Should I Use?
• How to handle confounding effects (including batches)
• What constitute an independent biological replicate?
• Adding differential gene expression
• The error “The annotation and quantification (count/abundance matrix and isoform annotation) seems to be different (Jaccard similarity < 0.925)”
• The warning “The annotation and quantification (count/abundance matrix and isoform annotation) seem to be slightly different”
• The error “The supplied design matrix will result in a model matrix that is not full rank”
• How to prevent visualization of NMD in the switchPlot
• BSGenomes for non-model organisms
• The switchPlot contains regions “Not Annotated”
• Can I use a GFF file instead of a GTF file?

What Quantification Tool(s) Should I Use?

This section is about considerations and recommendations for quantification of RNASeq data and are relevant both if you are:

1. Interested in isoform/transcript level quantification (as necessary for using with IsoformSwitchAnalyzeR (see Abstract))
2. Interest in any gene analysis (as described in the end of this paragraph)

But why not do both?

Isoform/transcript quantification

Generally speaking there are two different ways to obtain gene/transcript quantification and which to choose boils down to whether you A) trust the available transcript annotation or B) You think you RNA-seq data contains important additional isoforms not included in the available annotations.

Option A): If you generally trust the transcript annotation, state-of-the-art quantifications can be obtained directly from the sequence files (FASTQ) via Salmon or Kallisto. These methods give much more accurate quantifications (and are a lot faster) than for example first mapping to the genome and afterwards quantifying with Cufflinks. You can find more information in the Kallisto paper. For comparison here is the link to the Salmon paper but please not the comparison in figure 1 is a special case - under normal circumstances Kallisto and Salmon performs quite similarly and choice is up to personal preference. IsoformSwitchAnalyzeR directly supports the quantifications of both tools (and more) via the importIsoformExpression() function - see Importing Data Into R. Kallisto can be downloaded from here and the manual for running Kallisto can be found here. Salmon can be downloaded from here and a manual for running Salmon can be found here.

Option B): If you on the other hand do not have any (trustworthy) annotation, or want to improve it, you need to do a transcript deconvolution aka predicting novel transcripts from the RNAseq data. A transcript deconvolution can be described as a 4-step process:

1. First you map reads to the genome (with tools such as STAR, HISAT2 etc.).
2. Then you assemble transcripts (typically guided by existing annoation) with tools such as Cufflinks or StringTie. Please note that many other tools for solving this problem and new methods are still published regularly. Since the transcript assembly is the central step in the analysis, you probably want to read up on these tools yourself to make sure you use the tools that currently are the better one.
3. The third step is to merge the assembled transcripts from all individual Cufflinks/StringTie/other runs into one combined reference transcriptome with tools such as Cuffmerge, StringTie or TACO.
4. In step 4 this combined transcriptome can be used to quantify all the original FASTQ files with Cufflinks/CuffDiff, StringTie, Salmon or Kallisto to ensure the same transcripts were quantified in all samples (aka that samples are comparable).

As an alternative, tools such as Bayesembler ( paper ) or BIISQ ( paper) performs step 2-4 described above on all samples simultaneously.

You will need to figure out which combination of tools you choose on your own - but IsoformSwitchAnalyzeR can analyze the data from all of them. For more information see Importing Data Into R.

Gene quantification

With regards to gene-leve quantification it is now widely recognized that performing transcript quantification and afterwards obtaining the gene expression by adding up the expression from the individual transcripts will result in improved gene-level analysis. Prime examples illustrating concept is this blogpost and this article. The “To Long - Did not Read” is that compared to the old-school gene-quantification (genome-mapping + featureCounts/HTSeq) the main advantages of the transcript-level gene quantification are:

1. The inclusion of multi-mappers in the quantification (ignored by featureCounts/HTSeq) which typically correspind to 20-50% of all reads in a sample that would otherwise be ignored!.
2. The more accurate quantification algorithm (which among others include correction of various types of sequence-biases).
3. The ability to correctly quantify genes even if they contain isoform switches.
4. The runtime. Lightweight algorithms are MUCH faster!

In summary: unless you have single end RNASeq data and don’t know the fragment length there is no good reason for not doing transcript level quantification and afterwards summarizing to gene level. The gene level summaries can easily be obtainedwith tools such as tximport or IsoformSwitchAnalyzeR’s importIsoformExpression() (which acutally is just a tximport wrapper with extra functionalities).

Back to Table of Content

How to handle confounding effects (including batches)

“Confounding effect” is a common description of all sources of variation not of interest for the analysis you are doing. This can be anything from a batch effects (data/sample produced/analyzed at different locations, time points, personnel etc.) to an effect you are not interested in but which might influence the result (e.g. sex or age). Taking all this into account to avoid false positive of false negative results can easily be done when testing for isoform switches (or anything else really) when using tools based on linear models.

Naturally taking confounding effects into account is also supported in IsoformSwitchAnalyzeR via the isoformSwitchTestDRIMSeq() and isoformSwitchTestDEXSeq() functions. Both these functions automatically take additional confounding effects into account if they have been annotated in the switchAnalyzeRlist. isoformSwitchTestDEXSeq() furthermore also supports extract of batch corrected effect sizes (IF and dIFs).

Confounding effects should be annotated in ‘designMatrix’ used when creating the switchAnalyzeRList (via createSwitchAnalyzeRlist() or importRdata()). This is simply done by adding extra columns to the design matrix before creating the switchAnalyzeRlist:

So if the design matrix looked like:

myDesign <- data.frame(
    sampleID=1:6,
    condition=c('wt','wt','wt','ko','ko','ko') # WT = wildtype, KO = Knock out
)

myDesign
#>   sampleID condition
#> 1        1        wt
#> 2        2        wt
#> 3        3        wt
#> 4        4        ko
#> 5        5        ko
#> 6        6        ko

Adding a co-founder is simply done as follows:

myDesign$batch <- factor(c('a','a','b','a','b','b'))
myDesign
#>   sampleID condition batch
#> 1        1        wt     a
#> 2        2        wt     a
#> 3        3        wt     b
#> 4        4        ko     a
#> 5        5        ko     b
#> 6        6        ko     b

Please note that confounding effects can both be categorical (e.g. sex) (indicated by type character or factor) and continuous variables (e.g. age) (indicated by numeric or integers) - both will be taken into account. Due to the different interpretation it is recommended to use strings/factors for categorical variables (aka not use 1 and 2 to indicate batch in the example above (instead of “a” and “b”)). The conditions analyzed (as indicated by the “condition” column in the design matrix) must however always be categorical (character or factor).

Back to Table of Content

What constitute an independent biological replicate?

A very nice discussion, which I highly recommend, can be found in Vaux et al.

As a TL;DR here are a couple of guidelines:

• Science can under NO circumstances be made from sample sizes of 1 (n=1). To call anything science proof that a result is reproducible is needed.
• If all replicates were made at the same time no extrapolations can be made! Running an experiment in parallel does not result in independent samples - it is still only n=1. The experiment will need to be done multiple times.
• If all experiments were all done on the same cell line, extrapolation to anything other than that cell line cannot be made. The same actually goes for mouse strains.
• If a knock-down is performed multiple different siRNA/shRNAs are needed (both to be independent and due to off target effects).
• If a knock-out is created is not enough to only investigate a single cell line originating from single clone. Multiple independent clones are needed.
• If you want to investigate a disease, analyzing a single mouse or human is not enough - you would need data from multiple patients/mice to extrapolate to other humans/mice.

All these replicates which does not allow for extrapolation would be termed independent replicates - and if you ask me the work replicate should never be allowed to be used without specifying whether it is independent or not.

Back to Table of Content

Adding differential gene expression

Gene differential expression is currently not build into IsoformSwitchAnalyzeR since there are so many other good tools to do this. Some frequently used tools are edgeR, DESeq2, limma and Sleuth (for Kallisto users). This list is not exhaustive but clearly shows the variety of well-maintained tools out there. For a recent review of these tool we refer you to Love et al 2018. No matter which tool you choose we support integration of the result by showing the result in the plot created by the switchPlot() function.

To add the result of a differential expression analysis all you need to do is overwrite the following columns:

exampleSwitchList$isoformFeatures$gene_q_value
exampleSwitchList$isoformFeatures$iso_q_value

With the corresponding FDR corrected P-values of a differential expression analysis. We recommend the base function match() to do this - and do remember to consider conditions (indicated by the “condition_1” and “condition_2” columns) if you have more than one comparison.

Back to Table of Content

The error “The annotation and quantification (count/abundance matrix and isoform annotation) seems to be different (Jaccard similarity < 0.925)”

This problem can arise from multiple sources but essentially means that we determined that not all isoforms quantified are found in the annotation or vice versa. Ensuring that it is the exact same isoforms that are prescient in the annotation and quantification is essential because it is the foundation for all downstream analysis with IsoformSwitchAnalyzeR - therefore we need to be quite strict with this.

If there is no overlap at all (Jaccard similarity is zero) there are two possibilities:

1. It most likely has something to do with how the isoform_ids are stored in the different data sources. This most likely can be solved using some combination of the “ignoreAfterBar”, “ignoreAfterSpace” and “ignoreAfterPeriod” in the import* functions.

2. When using RSEM for the quantification some older versions of tximport could cause this problem. To fix this problem you need to update the underlying package tximport. This can be done by copy/pasting this command into R:

   devtools::install_github(“mikelove/tximport”)

If there is a large overlap but still far from perfect, there are two possibilities:

1. If you are using Ensembl data you have supplied GTF with haplotypes ( .chr_patch_hapl_scaff.gtf, - NOT the .chr.gtf). This in combination with the default status of the “ignoreAfterBar”, “ignoreAfterSpace” and “ignoreAfterPeriod” in the import* functions should do the trick.
2. One data source could contain non-canonical chromosomes while the other do not - this might be solved using the removeNonConvensionalChr argument.

Back to Table of Content

The warning “The annotation and quantification (count/abundance matrix and isoform annotation) seem to be slightly different”

This warning most likely occurs because there are differences in which isoforms are found in the annotation and which isoforms have been quantified. This problem has (so far) only been observed in data quantified with Salmon. It typically occurs because during the index building Salmon collapses completely identical transcripts (based on the nucleotide sequence). A discussion with the Salmon authors of pros and cons of this behavior can be found here.

We ask the user to make sure why this problem occurs in the user-specific setting, and that removing the non-quantified isoforms is the desired solution, before continuing. This can be done by comparing the number of sequences to the number of unique sequences.

Back to Table of Content

The error “The supplied design matrix will result in a model matrix that is not full rank”

If a design matrix is not of full rank it simply means that one or more columns (or combinations of columns) in the design matrix describe the same thing. Imaging a case when you had two parameters describing the same person - say fingerprint and social security number. If you then try to estimate the effect of these factors you can have an infinite number of solutions by having them being opposites (-1 and 1, -2 and 2 etc.) which means you cannot estimate a linear model. If you have problems with this the easiest thing is to write a post about it on Bioconductor support site.

How to prevent visualization of NMD in the switchPlot

In some biological systems the role of NMD is not well established whereby one potentially would avoid visualizing it. This can be done simply by removing the PTC annotation from the switchAnalyzeRlist as follows:

aSwitchList$isoformFeatures$PTC <- NULL

BSGenomes for non-model organisms

If you are analyzing a non-model organism the easiest way of incorporating the underlying sequence analysis is to use the isoformNtFasta argument in importRdata. This is the recommended analysis workflow even for model organism and will make the entire analysis more accurate as well as speed up specific steps of the workflow.

Alternatively one can build a custom BSgenome as describe here.

The switchPlot contains regions “Not Annotated”

The problem is the increasingly strict limitations on some of the external sequences analysis tools (primary EBI’s Pfam) specifically their websites which among others have limitations on the max length of the sequences which can be analyzed. To extract the maximum information IsoformSwitchAnalyzeR supports trimming of the amino acid sequences so the fit the length allowed to be analyzed. This naturally mean that some parts of a sequence cannot be analyzed - and that is what is visualized in the switchPlot. Furthermore please note that un-analyzed regions are ignored during the consequence predictions meaning that the trimming result in no bias. Specifically this is controlled with the “filterAALength” argument of extractSequence() function (but the websites will not run them if they are too long).

Complaints about these limitations can be send to EBI via their website or on twitter via @hmm3r

Can I use a GFF file instead of a GTF file?

Unfortunately not. But there is a simple work around using the following pearl script:

#!/usr/bin/perl

use warnings;
use strict;
while (<>) {
chomp;
my @parts = split /\t/;
if ( !(($parts[2] eq 'three_prime_UTR')|($parts[2] eq 'five_prime_UTR')|($parts[2] eq 'exon')|($parts[2] eq 'CDS')) ) {
next;
}
elsif ( $parts[2] eq 'three_prime_UTR') {
$parts[2] = '3UTR';
}
elsif ($parts[2] eq 'five_prime_UTR') {
$parts[2] = '5UTR';
}
$parts[8] =~ s/^Parent=((\w+)\.\w*).*/gene_id \"$2\"; transcript_id \"$1\";/sg;
$_ = join "\t" , @parts;
print "$_\n";
}

Simply save the script as “converter_gff_to_gtf.pl”"

and run

perl converter_gff_to_gtf.pl input.gff > output.gtf

Generating additional links for sharing online

• The switchPlot
• Genome-wide Summaries
• Consequence Enrichment Analysis
• Splicing Enrichment Analysis
• Overview Plots