2.1 Creating the maser object

We illustrate the workflow using the hypoxia dataset from the previous vignette.

library(maser)
library(rtracklayer)

# Creating maser object using hypoxia dataset
path <- system.file("extdata", file.path("MATS_output"),
                    package = "maser")
hypoxia <- maser(path, c("Hypoxia 0h", "Hypoxia 24h"))

# Remove low coverage events
hypoxia_filt <- filterByCoverage(hypoxia, avg_reads = 5)

2.2 Query available protein features at Uniprot

Available UniprotKB protein annotation can be queried using availableFeaturesUniprotKB(). Currently there are 30 distinct features grouped into broader categories, including Domain and Sites, PTM (Post-translational modifications), Molecule Processing, Topology, Mutagenesis and Structural Features.

Name	Description	Category
act-site	UniProtKB active site sequence annotations	Domain_and_Sites
binding	UniProtKB binding site sequence annotations	Domain_and_Sites
Ca-binding	UniProtKB Calcium binding site sequence annotations	Domain_and_Sites
coiled	UniProtKB coiled coil sequence annotations	Domain_and_Sites
DNA-bind	UniProtKB DNA binding site sequence annotations	Domain_and_Sites
domain	UniProtKB domain sequence annotations	Domain_and_Sites
metal	UniProtKB metal ion binding site sequence annotations	Domain_and_Sites
motif	UniProtKB motif of interest sequence annotations	Domain_and_Sites
NP bind	UniProtKB Nucleotide Phosphate binding sequence annotations	Domain_and_Sites
region	UniProtKB region of interest sequence annotations	Domain_and_Sites

2.3 Steps for annotation

Protein feature annotation of splicing events is performed in two steps.

1. Use mapTranscriptsToEvents() to add both transcript and protein IDs to all events in the maser object.
2. Use mapProteinFeaturesToEvents() for annotation specifying UniprotKB features or categories.

mapTranscriptsToEvents() identifies transcripts compatible with the splicing event by overlapping exons involved in splicing to the gene models provided in the Ensembl GTF. Each type of splice event applies a specific overlapping rule (described in the Introduction vignette). The function also maps transcripts to their corresponding protein identifiers in Uniprot when available.

mapTranscriptsToEvents() requires an Ensembl or Gencode GTF using the hg38 build of the human genome. Ensembl GTFs can be retrieved using AnnotationHub or imported using import.gff() from the rtracklayer package. Several GTF releases are available, and maser is compatible with any version using the hg38 build.

We are using reduced GTF extracted from Ensembl Release 85 for running examples.

## Ensembl GTF annotation
gtf_path <- system.file("extdata", file.path("GTF","Ensembl85_examples.gtf.gz"),
                        package = "maser")
ens_gtf <- rtracklayer::import.gff(gtf_path)

In the second step, mapProteinFeaturesToEvents() retrieves data from UniprotKB and overlaps splicing events to genomic coordinates of protein features.

2.4 SRSF6 example

The splicing factor SRSF6 undergoes splicing during hypoxia by expressing an alternative exon. We will annotate the exon skipping event with domain, sites and topology information. The first step is to obtain a maser object containing SRSF6 splicing information, and then map transcripts to splicing events.

# Retrieve gene specific splicing events
srsf6_events <- geneEvents(hypoxia_filt, "SRSF6")
srsf6_events
#> A Maser object with 1 splicing events.
#> 
#> Samples description: 
#> Label=Hypoxia 0h     n=3 replicates
#> Label=Hypoxia 24h     n=3 replicates
#> 
#> Splicing events: 
#> A3SS.......... 0 events
#> A5SS.......... 0 events
#> SE.......... 1 events
#> RI.......... 0 events
#> MXE.......... 0 events

# Map transcripts to splicing events
srsf6_mapped <- mapTranscriptsToEvents(srsf6_events, ens_gtf)

If transcript mapping worked correctly, Ensembl and Uniprot identifiers will be added to splicing events. Possible NA values indicates non-protein coding transcripts. In this case, the splicing involves two Ensembl transcripts coding for the Q13247 isoform of SRSF6.

head(annotation(srsf6_mapped, "SE"))

ID	GeneID	geneSymbol	txn_3exons	txn_2exons	list_ptn_a	list_ptn_b
33209	ENSG00000124193.14	SRSF6	ENST00000483871	ENST00000244020	Q13247	Q13247

Now we are ready to call mapProteinFeaturesToEvents() for annotation. Feature annotation can be interactively displayed in a web browser using display() or retrieved as a data.frame using annotation().

mapProteinFeaturesToEvents() will add extra columns describing the feature name, feature description and protein identifiers for which the annotation has been assigned. Possible NA values indicate the particular feature is not annotated for the splice event.

# Annotate splicing events with protein features
srsf6_annot <- mapProteinFeaturesToEvents(srsf6_mapped, c("Domain_and_Sites", "Topology"), by="category")

head(annotation(srsf6_annot, "SE"))

ID	GeneID	geneSymbol	txn_3exons	txn_2exons	list_ptn_a	list_ptn_b	act-site	binding	Ca-binding	coiled	DNA-bind	domain	metal	motif	NP bind	region	repeat	site	Zn-fing	intramem	topo-dom	transmem
33209	ENSG00000124193.14	SRSF6	ENST00000483871	ENST00000244020	Q13247	Q13247	NA	NA	NA	NA	NA	Q13247:M1-G72;RRM1,Q13247:Y110-P183;RRM2	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA

By inspecting the results, we see that the SRSF6 exon skipping event is annotated with the Uniprot features domain, chain and mod-res (modidifed residue). Visualization of the splice event, transcripts and protein features is performed with plotUniprotKBFeatures(). In this example, exons in the splice event overlap the Serine/arginine-rich splicing factor 6 region of the protein, while the upstream exon and downstream exons are overlapping the RRM1 and RRM2 domains of SRSF6, respectively.

# Plot splice event, transcripts and protein features
plotUniprotKBFeatures(srsf6_mapped, "SE", event_id = 33209, gtf = ens_gtf, 
   features = c("domain", "chain"), show_transcripts = TRUE)

2.5 RIPK2 example

RIPK2 has an exon skipping event in the hypoxia dataset. Following the example above, we map transcripts to splicing events and annotate protein features overlapping the splice event. We find out that the alternative exon overlaps the kinase domain of the protein, thus possibly changing the configuration of this domain during hypoxia. The ATP and proton acceptor binding sites are overlapping exons flanking the alternative exon.

ripk2_events <- geneEvents(hypoxia_filt, "RIPK2")
ripk2_mapped <- mapTranscriptsToEvents(ripk2_events, ens_gtf)
ripk2_annot <- mapProteinFeaturesToEvents(ripk2_mapped, 
                                          tracks = c("Domain_and_Sites"), 
                                          by = "category")

plotUniprotKBFeatures(ripk2_annot, type = "SE", event_id = 14319, 
                      features = c("domain", "binding", "act-site"), gtf = ens_gtf, 
                      zoom = FALSE, show_transcripts = TRUE)