# Preparing an Alternative Splicing Annotation for psichomics

#### Nuno Saraiva Agostinho

#### 2018-10-02

psichomics currently quantifies alternative splicing based on an alternative splicing annotation for Human (hg19). This annotation was created based on the annotations from MISO, SUPPA, VAST-TOOLS and rMATS.

The preparation of a new alternative splicing annotation for usage in psichomics starts by parsing the alternative splicing events. Currently, events from SUPPA, rMATS, MISO and VAST-TOOLS are parsable. Support can be extended for alternative splicing events from other programs.

Also, start by loading the following packages:

    library(psichomics)

    ## Loading required package: shiny

    ## Loading required package: shinyBS

    ## Setting options('download.file.method.GEOquery'='auto')

    ## Setting options('GEOquery.inmemory.gpl'=FALSE)

    ## No methods found in package 'IRanges' for request: 'subset' when loading 'derfinder'

    ## Start the visual interface by running the function psichomics()

    ##
    ## Attaching package: 'psichomics'

    ## The following object is masked from 'package:base':
    ##
    ##     rowMeans

    library(plyr)

## SUPPA annotation

SUPPA generates alternative splicing events based on a transcript annotation. Start by running SUPPA’s `generateEvents` script with a transcript file (GTF format) for all event types, if desired. See SUPPA’s page for more information.

The resulting output will include a directory containing tab-delimited files with alternative splicing events (one file for each event type). Hand the path of this directory to the function `parseSuppaAnnotation()`, prepare the annotation using `prepareAnnotationFromEvents()` and save the output to a RDS file:

    # Replace `genome` for the string with the identifier before the first
    # underscore in the filenames of that directory (for instance, if one of your
    # filenames of interest is "hg19_A3.ioe", the string would be "hg19")
    suppa <- parseSuppaAnnotation(suppaOutput, genome="hg19")

    ##
    ## Retrieving SUPPA annotation...
    ##
    ## Parsing SUPPA annotation...

    annot <- prepareAnnotationFromEvents(suppa)

    ##
    ## Sorting coordinates...
    ##
    ## A3SS SUPPA
    ##
    ## A5SS SUPPA
    ##
    ## AFE SUPPA
    ##
    ## ALE SUPPA
    ##
    ## MXE SUPPA
    ##
    ## Joining events per event type...
    ##
    ## A3SS
    ##
    ## A5SS
    ##
    ## AFE
    ##
    ## ALE
    ##
    ## MXE
    ##
    ## RI
    ##
    ## SE
    ##
    ## Cleaning the annotation...

    # suppaFile <- "suppa_hg19_annotation.RDS"
    saveRDS(annot, file=suppaFile)

## rMATS annotation

Just like SUPPA, rMATS also allows to generate alternative splicing events based on a transcript annotation, although two BAM or FASTQ files are required to generate alternative splicing events. Read rMATS’s page for more information.

The resulting output of rMATS is then handed out to the function `parseMatsAnnotation()`:

    mats <- parseMatsAnnotation(
        matsOutput,         # Output directory from rMATS
        genome = "fromGTF", # Identifier of the filenames
        novelEvents=TRUE)   # Parse novel events?

    ##
    ## Retrieving rMATS annotation...
    ##
    ## Parsing rMATS annotation...

    annot <- prepareAnnotationFromEvents(mats)

    ##
    ## Sorting coordinates...
    ##
    ## A3SS MATS
    ##
    ## A5SS MATS
    ##
    ## AFE MATS
    ##
    ## ALE MATS
    ##
    ## MXE MATS
    ##
    ## Joining events per event type...
    ##
    ## A3SS
    ##
    ## A5SS
    ##
    ## AFE
    ##
    ## ALE
    ##
    ## MXE
    ##
    ## RI
    ##
    ## SE
    ##
    ## Cleaning the annotation...

    # matsFile <- "mats_hg19_annotation.RDS"
    saveRDS(annot, file=matsFile)

## MISO annotation

Simply retrieve MISO’s alternative splicing annotation and give the path to the downlaoded folder as input.

    miso <- parseMisoAnnotation(misoAnnotation)

    ##
    ## Retrieving MISO annotation...
    ##
    ## Parsing MISO annotation...

    annot <- prepareAnnotationFromEvents(miso)

    ##
    ## Sorting coordinates...
    ##
    ## A3SS MISO
    ##
    ## A5SS MISO
    ##
    ## AFE MISO
    ##
    ## ALE MISO
    ##
    ## MXE MISO
    ##
    ## TandemUTR MISO
    ##
    ## Joining events per event type...
    ##
    ## A3SS
    ##
    ## A5SS
    ##
    ## AFE
    ##
    ## ALE
    ##
    ## MXE
    ##
    ## RI
    ##
    ## SE
    ##
    ## TandemUTR
    ##
    ## Cleaning the annotation...

    # misoFile <- "miso_AS_annotation_hg19.RDS"
    saveRDS(annot, file=misoFile)

## VAST-TOOLS annotation

Simply retrieve VAST-TOOLS’s alternative splicing annotation and give the path to the downlaoded folder as input. Note that, however, complex events (i.e. alternative coordinates for the exon ends) are not yet parseable.

    vast <- parseVastToolsAnnotation(vastAnnotation)

    ##
    ## Retrieving VAST-TOOLS annotation...
    ##
    ## Parsing VAST-TOOLS annotation...
    ##
    ## ALT3
    ##
    ## ALT5
    ##
    ## COMBI
    ##
    ## IR
    ##
    ## MERGE3m
    ##
    ## MIC
    ##
    ## EXSK
    ##
    ## MULTI

    annot <- prepareAnnotationFromEvents(vast)

    ##
    ## Sorting coordinates...
    ##
    ## A3SS VAST-TOOLS
    ##
    ## A5SS VAST-TOOLS
    ##
    ## Joining events per event type...
    ##
    ## A3SS
    ##
    ## A5SS
    ##
    ## RI
    ##
    ## SE
    ##
    ## Cleaning the annotation...

    # vastFile <- "vast_AS_annotation_hg19.RDS"
    saveRDS(annot, file=vastFile)

## Combining annotation from different sources

To combine the annotation from different sources, provide the parsed annotations of interest simultasneously to the function `prepareAnnotationFromEvents`:

    # Let's combine the annotation from SUPPA, MISO, rMATS and VAST-TOOLS
    suppa <- parseSuppaAnnotation(suppaOutput, genome="hg19")

    ##
    ## Retrieving SUPPA annotation...
    ##
    ## Parsing SUPPA annotation...

    miso <- parseMisoAnnotation(misoAnnotation)

    ##
    ## Retrieving MISO annotation...
    ##
    ## Parsing MISO annotation...

    mats <- parseMatsAnnotation(matsOutput, genome = "fromGTF", novelEvents=TRUE)

    ##
    ## Retrieving rMATS annotation...
    ##
    ## Parsing rMATS annotation...

    vast <- parseVastToolsAnnotation(vastAnnotation)

    ##
    ## Retrieving VAST-TOOLS annotation...
    ##
    ## Parsing VAST-TOOLS annotation...
    ##
    ## ALT3
    ##
    ## ALT5
    ##
    ## COMBI
    ##
    ## IR
    ##
    ## MERGE3m
    ##
    ## MIC
    ##
    ## EXSK
    ##
    ## MULTI

    annot <- prepareAnnotationFromEvents(suppa, vast, mats, miso)

    ##
    ## Sorting coordinates...
    ##
    ## A3SS MATS
    ##
    ## A3SS MISO
    ##
    ## A3SS SUPPA
    ##
    ## A3SS VAST-TOOLS
    ##
    ## A5SS MATS
    ##
    ## A5SS MISO
    ##
    ## A5SS SUPPA
    ##
    ## A5SS VAST-TOOLS
    ##
    ## AFE MATS
    ##
    ## AFE MISO
    ##
    ## AFE SUPPA
    ##
    ## ALE MATS
    ##
    ## ALE MISO
    ##
    ## ALE SUPPA
    ##
    ## MXE MATS
    ##
    ## MXE MISO
    ##
    ## MXE SUPPA
    ##
    ## TandemUTR MISO
    ##
    ## Joining events per event type...
    ##
    ## A3SS
    ##
    ## A5SS
    ##
    ## AFE
    ##
    ## ALE
    ##
    ## MXE
    ##
    ## RI
    ##
    ## SE
    ##
    ## TandemUTR
    ##
    ## Cleaning the annotation...

    # annotFile <- "AS_annotation_hg19.RDS"
    saveRDS(annot, file=annotFile)

## Quantifying alternative splicing using the created annotation

The created alternative splicing annotation can then be used when quantifying alternative splicing events. To do so, when using the GUI version of psichomics, be sure to select the Load annotation from file… option, click the button that appears below and select the recently created RDS file. Otherwise, if you are using the CLI version, perform the following steps:

    annot <- readRDS(annotFile) # "file" is the path to the annotation file
    junctionQuant <- readFile("ex_junctionQuant.RDS") # example set

    psi <- quantifySplicing(annot, junctionQuant,
                            c("SE", "MXE", "A3SS", "A5SS", "AFE", "ALE"))
    psi # may have 0 rows because of the small junction quantification set

    ##                                              Normal 1  Normal 2 Normal 3
    ## SE_1_+_23385660_23385840_23385851_23395032_ 0.7777778 0.4444444      0.6
    ##                                             Cancer 1  Cancer 2  Cancer 3
    ## SE_1_+_23385660_23385840_23385851_23395032_      0.4 0.4285714 0.5555556

# Feedback

All feedback on the program, documentation and associated material (including this tutorial) is welcome. Please send any suggestions and comments to:

> Nuno Saraiva-Agostinho (nunoagostinho@medicina.ulisboa.pt)

> [Disease Transcriptomics Lab, Instituto de Medicina Molecular (Portugal)][iMM]