---
title: "Using proteasy to Retrieve and Analyze Protease Data"
author: "Martin Rydén"
package: proteasy
output:
  BiocStyle::html_document:
    toc_float: true
vignette: >
  %\VignetteIndexEntry{Using proteasy to Retrieve and Analyze Protease Data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteKeywords{Proteomics, Bioinformatics,
  Cleavage, Polypeptides, MEROPS, Proteases}
  %\VignetteEncoding{UTF-8}
  %\VignettePackage{proteasy}
---

```{r biocstyle, echo = FALSE, results = "asis", message = FALSE}
library(BiocStyle)
BiocStyle::markdown()

```

# Introduction

## Motivation

Protease cleavages affect many vital physiological processes, and dysregulation
of proteolytic activity is associated with a variety of pathological
conditions. The field of degradomics is committed to provide insights about
proteolytic events by incorporating techniques for identification and
functional characterization of proteases and their substrates and inhibitors.

`r BiocStyle::Biocpkg("proteasy")` allows for batch identification of possible
proteases for a set of substrates (protein IDs and peptide sequences), and may
be an important tool in peptide-centric analyses of endogenously cleaved
peptides.

This package utilizes data derived from the
[MEROPS database](https://www.ebi.ac.uk/merops/).
The database is a manually curated knowledgebase with information about
proteolytic enzymes, their inhibitors and substrates. 

This document illustrates the functionality of proteasy through some use cases.

## Package scope and limitations

Similarly to existing tools such as
[TopFind](https://topfind.clip.msl.ubc.ca/), and
[Proteasix](http://www.proteasix.org), `r BiocStyle::Biocpkg("proteasy")`
exists for the purpose of retrieving data about proteases by mapping peptide
termini positions to known sites where a protease cleaves. Unlike the
`r BiocStyle::Biocpkg("cleaver")` package, which allows for in-silico cleavage
of amino acid sequences by specifying an enzyme, the functions in
`r BiocStyle::Biocpkg("proteasy")` relies only on experimentally derived data
to find proteases annotated to cleavage sites.

The `r BiocStyle::Biocpkg("proteasy")` package is limited to entries annotated
in MEROPS. Moreover, `r BiocStyle::Biocpkg("proteasy")` currently does not
allow for retrieval of proteolytic details for multiple organisms at once, or
inter-organism events. The package does however provide annotation data for all
organisms available in MEROPS.

The function findProtease will automatically map a peptide's start- and end-
positions in its host-protein's sequence. When a protein sequence has repeated
instances of a peptide sequence, `r BiocStyle::Biocpkg("proteasy")` will 
only match to the first instance in the protein sequence.


# Installation

The package should be installed using as described below:

```{r installPackage, eval = FALSE}

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("proteasy")

```

The package can then be loaded with:

```{r loadPackage, eval = TRUE}

library("proteasy")

```

# Usage

## Find proteases by substrate search

A fast way to find which possible proteases, if any, are annotated as cleaving
actors for a substrate is by using the function *searchSubstrate*. Setting the 
parameter summarize = TRUE will return only a vector of reviewed proteases,
and summarize = FALSE will return a table with details about each cleaving
event. We will explore the two options for kininogen-1 (P01042).

```{r SearchSubstrate, eval = TRUE}

# Returns vector of reviewed proteases that cleave P01042
searchSubstrate(protein = "P01042", summarize = TRUE)

# Returns data.table with details on cleaving events
searchSubstrate(protein = "P01042", summarize = FALSE) |>
    head()

# The function also accepts multiple proteins as input. Let's inspect
# the last rows in the returned data.table.
searchSubstrate(protein = c("P01042", "P02461"), summarize = FALSE) |>
    tail()

# With summarize = FALSE we get both reviewed and unreviewed proteases.

```


## Find substrates by protease search

A corresponding function, *searchProtease*, exists to find which (if any)
substrates a protease cleaves. The function is demonstrated with
MMP-12 (P39900) as an example.

```{r SearchProtease, eval = TRUE}

# Returns vector of substrates that MMP-12 cleave
searchProtease(protein = "P39900", summarize = TRUE)

# Returns data.table with details on cleaving events that involve MMP-12
searchProtease(protein = "P39900", summarize = FALSE) |> head()

```

## Find possible proteases for cleaved peptides

The function *findProtease* automatically maps the peptide sequences to the
full-length protein sequence and obtains the start- and end-positions for the
peptide. Then, the positions are searched against the MEROPs database and
matches are returned.


```{r FindProtease, eval = TRUE}

# Create a vector of Uniprot IDs and corresponding peptide sequences

protein <- c("P02671", "P02671", "P68871",
             "P01011", "P68133", "P02461",
             "P0DJI8", "P0DJI8", "P0DJI8")
peptide <- c("FEEVSGNVSPGTR", "FVSETESR", "LLVVYPW",
             "ITLLSAL", "DSYVGDEAQS", "AGGFAPYYG",
             "FFSFLGEAFDGAR", "EANYIGSDKY", "GGVWAAEAISDAR")

# If we do not specify start position (start_pos) and end position
# (end_pos), the function will automatically assign these values by
# matching the provided peptide sequence against the full-length
# protein sequence.

res <- findProtease(protein = protein,
                    peptide = peptide,
                    organism = "Homo sapiens")

# The resulting S4 object can be inspected in three ways;
# (to save space we show only the first six rows)

# 1. Display sequence data for the provided input:

substrates(res) |> head()

# 2. Show which known proteases may have cleaved this protein:

proteases(res) |> head()

# 3. Display details of associated proteolytic events:

cleavages(res) |> head()

# We can find out what proportion of matching cleaving events by reviewed
# proteases occur at N- versus C-terminus

cl <- cleavages(res)[`Protease status` == "reviewed"]

cl$`Cleaved terminus` |> table() |> prop.table() |> round(digits = 2)

# And inspect the distribution of cleaved amino acids

cl$`Cleaved residue` |> table() |> barplot(cex.names = 2)


# Find which protease is involved in the greatest number of cleaving events

cl[!duplicated(Peptide), .(count = .N), by = `Protease (Uniprot)`]

# If start- and end-positions are specified, the function will not attempt
# to automatically look up sequence data for the specified protein/peptide.

cl_by_pos <- findProtease(
    protein = "P02671",
    peptide = "FEEVSGNVSPGTR",
    start_pos = 413,
    end_pos = 425
)

# However, this means sequence details for substrates is not available.

substrates(cl_by_pos)



```

## Look up a protease on MEROPS

We may want to read up on the details of an entry directly in MEROPS.
The function *browseProtease* takes a UniProt or MEROPS ID and opens the 
MEROPS summary page which corresponds to that ID in a web browser.

```{r Browse protease, eval = FALSE}

browseProtease("P07339", keytype = "UniprotID") # (opens web browser)

```

# Additional examples

## Plot cleavages as a protein-protein interaction network.

Here we visualize the cleaving events of two substrates as a protein-protein
interaction network between substrates (red) and associated proteases (blue).

```{r network, eval = TRUE, message = FALSE, fig.height = 12, fig.width = 12, crop = FALSE}
library(igraph)
library(data.table)
# Make a vector of substrates we want to investigate
proteins <- c('P01011','P02671')
# Look up known proteases which cleave these substrates, and associated data
# annotated to the cleaving events
res <- searchSubstrate(protein = proteins, summarize = FALSE)
# Filter to keep proteases with Uniprot status "reviewed"
res <- res[`Protease status` == "reviewed"]
# To create a network, we only need two columns of interactors
# i.e. the proteases with cleaving action on the substrates
res <- res[, c("Protease (Uniprot)", "Substrate (Uniprot)", "Cleavage type")]
# Construct the DAG
g <- igraph::graph_from_data_frame(res,
                                   directed = TRUE,
                                   vertices = unique(
                                       c(res$`Protease (Uniprot)`,
                                         res$`Substrate (Uniprot)`)))
# This will allow us to return to a layout we like
# (and where the legend fits nicely)
set.seed(104)
plot(g,
     vertex.label.family = "Helvetica",
     vertex.size = 14,
     vertex.color = ifelse(V(g)$name %in% res$`Protease (Uniprot)`,
                           "#377EB8", "#E41A1C"),
     vertex.label.cex = 1,
     vertex.label.color = "white",
     edge.arrow.size = .6,
     edge.color =  ifelse(E(g)$`Cleavage type` == "physiological",
                          "#01665E", "#8E0152"),
     layout = igraph::layout.davidson.harel)
legend(title = "Nodes", cex = 2, horiz = FALSE,
       title.adj = 0.0, inset = c(0.0, 0.2),
       "bottomleft", bty = "n",
       legend = c("Protease", "Substrate"),
       fill = c("#377EB8", "#E41A1C"), border = NA)
legend(title = "Edges", cex = 2, horiz = FALSE,
       title.adj = 0.0, inset = c(0.0, 0.0),
       "bottomleft", bty = "n",
       legend = c("Physiological", "Non-physiological"),
       fill = c("#01665E", "#8E0152"), border = NA)
```

## Annotated sequence similarity heatmaps

`r BiocStyle::Biocpkg("proteasy")` automatically finds protein sequences 
from protein IDs (thanks to `r BiocStyle::Biocpkg("ensembldb")` and
`r BiocStyle::Biocpkg("Rcpi")`). We can use *substrates()* to access them.
Here, we look study similarity matrices of some protein- and peptide-level
sequences and plot them as heatmaps, which we annotate with cleavage data.

```{r, message = FALSE, warning = FALSE, crop = FALSE}

# Load additional libraries

library(Rcpi)
library(viridis)
suppressPackageStartupMessages(library(ComplexHeatmap))

# Prepare input: protein and associated peptide sequences

protein <- c('P01011','P01011','P01034','P01034',
             'P01042','P02671','P02671','P68871',
             'P68871','P01042')

peptide <- c('LVETRTIVRFNRPFLMIIVPTDTQNIFFMSKVTNPK','ITLLSAL',
             'KAFCSFQIY','AFCSFQIY','DIPTNSPELEETLTHTITKL','FEEVSGNVSPGTR',
             'FVSETESR','LLVVYPW','VDEVGGEALGR','KIYPTVNCQPLGMISLM')

# Find cleaving data associated with above substrates

res <- findProtease(protein = protein,
                    peptide = peptide,
                    organism = "Homo sapiens")

# Get substrate info

ss <- substrates(res)

# Show only unique sequences

ss_uniq <- ss[!duplicated(`Substrate sequence`)]

# Calculate protein (substrate) sequence similarity

psimmat = Rcpi::calcParProtSeqSim(ss_uniq$`Substrate sequence`,
                                  type = 'global',
                                  submat = 'BLOSUM62')

rownames(psimmat) <- colnames(psimmat) <- ss_uniq$`Substrate (Uniprot)`

# Plot as a heatmap

ComplexHeatmap::Heatmap(psimmat, col = viridis::mako(100))

# We can do the same thing for peptide sequences,
# and annotate each row (cleaved peptide) with
# information about cleaved residue and termini

# Get cleavage details

cl <- cleavages(res)

# Calculate peptide sequence similarity

pep_psimmat = Rcpi::calcParProtSeqSim(cl$Peptide, type = 'global',
                                      submat = 'BLOSUM62')

# Right side annotation: cleaved residue

rsd <- cl$`Cleaved residue`
cols <- c("#8DD3C7", "#FFFFB3", "#BEBADA", "#FB8072")
names(cols) <- unique(rsd)
ha1 <- ComplexHeatmap::columnAnnotation(`cleaved residue` = rsd,
                         col = list(`cleaved residue` = cols))

# Right side annotation: Terminus

tn <- cl$`Cleaved terminus`
cols <- c("#B3E2CD", "#FDCDAC")
names(cols) <- unique(tn)
ha2 <- ComplexHeatmap::columnAnnotation(terminus = tn,
                     col = list(terminus = cols))

rownames(pep_psimmat) <- cl$`Substrate (Uniprot)`

# Plot as a heatmap

ComplexHeatmap::Heatmap(
    pep_psimmat,
    name = "sequence\nsimilarity",
    col = viridis::mako(100),
    show_column_names = FALSE,
    row_names_gp = grid::gpar(fontsize = 6.5),
    top_annotation = c(ha1, ha2)
)

```


# Session Information

```{r sessioninfo, echo = FALSE}
sessionInfo()
```
