---
title: "Cytoscape StringApp"
author: "Kozo Nishida, Kristina Hanspers and Alex Pico"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 4
    df_print: paged
---
```{r, echo = FALSE}
knitr::opts_chunk$set(
  eval=FALSE
)
```

*The R markdown is available from the pulldown menu for* Code *at the upper-right, choose "Download Rmd", or [download the Rmd from GitHub](https://raw.githubusercontent.com/cytoscape/cytoscape-automation/master/for-scripters/R/notebooks/stringApp.Rmd).*

<hr />

In these exercises, we will use the [stringApp](http://apps.cytoscape.org/apps/stringApp) for [Cytoscape](http://cytoscape.org/) to retrieve molecular networks from the [STRING](https://string-db.org/) and [STITCH](http://stitch-db.org/) databases. The exercises will teach you how to:

- retrieve networks for proteins or small-molecule compounds of interest
- retrieve networks for a disease or arbitrary topics in PubMed
- layout and visually style the resulting networks
- import external data and map them onto a network
- perform enrichment analyses and visualize the results
- merge and compare networks
- select proteins by attributes
- identify functional modules through network clustering

The original version of this tutorial was developed by Lars Juhl Jensen of the Novo Nordisk Center for Protein Research at the University of Copenhagen. We thank professor Jensen for his gracious willingness to allow us to repackage the content for delivery as a Cytoscape tutorial.

<hr />

# Installation
```{r, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

if(!"RCy3" %in% installed.packages())
  BiocManager::install("RCy3")
library(RCy3)

if(!"RColorBrewer" %in% installed.packages()){
    install.packages("RColorBrewer")
}
library(RColorBrewer)
```

# Getting started
First, launch Cytoscape and keep it running whenever using RCy3. Confirm that you have everything installed and running:
```{r}
cytoscapePing()
cytoscapeVersionInfo()
```

The exercises require you to have certain Cytoscape apps and R packages installed.

```{r}
installApp('stringApp')
installApp('enhancedGraphics')
installApp('clusterMaker2')
```

If you are not already familiar with the STRING database, we highly recommend that you go through the short [STRING exercises](https://jensenlab.org/training/string/) provided by the [Jensen lab](https://jensenlab.org/) to learn about the underlying data before working with them in these exercises.

# Exercise 1

In this exercise, we will perform some simple queries to retrieve molecular networks based on a protein, a small molecule compound, a disease, and a topic in PubMed.

## Protein queries

You can query **STRING** as protein data source (in the following query the **protein name** is **SORCS2**).
You can select the appropriate organism by setting the **species** option (e.g. **Homo sapiens**).
(And the following chunk also imports the result into Cytoscape.)

```{r}
string.cmd = 'string protein query query="SORCS2" species="Homo sapiens"'
commandsRun(string.cmd)
```

The limit (Maximum number of interactors) option determines how many interaction partners of your protein(s) of interest will be added to the network.
By default, if you enter only one protein name, the resulting network will contain 10 additional interactors. If you enter more than one protein name, the network will contain only the interactions among these proteins, unless you explicitly ask for additional proteins.

```{r}
string.cmd = 'string protein query query="SORCS2" species="Homo sapiens" limit=100'
commandsRun(string.cmd)
```

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_SORCS2.png){width=60%}
</center>

How many nodes are in the resulting network?
How does this compare to the maximum number of interactors you specified?
What types of information does the **Node Table** provide?

## Compound queries

You can query **STITCH** as protein/compound data source (in the following query the **compound name** is **imatinib**).
You can select the organism and number of additional interactors just like for the protein query above.

```{r}
string.cmd = 'string compound query query="imatinib" species="Homo sapiens"'
commandsRun(string.cmd)
```

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_imatinib.png){width=60%}
</center>

How is this network different from the protein-only network with respect to node types and the information provided in the Node Table?

## Disease queries

You can query **STRING** as disease data source (in the following query the **disease term** is **Alzheimer**).
The stringApp will retrieve a STRING network for the top-N proteins (by default 100) associated with the disease.

```{r}
string.cmd = 'string disease query disease="Alzheimer"'
commandsRun(string.cmd)
```

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_alzheimer.png){width=60%}
</center>

Which additional attribute column do you get in the Node Table for a disease query compared to a protein query? (Hint: check the last column.)

## PubMed queries

You can query **STRING** as **PubMed** data source (in the following query the query representing a topic of interest is **jet-lag**).
You can use any query that would work on the PubMed website, but it should obviously a topic with related genes or proteins. The stringApp will query PubMed for the abstracts, find the top-N proteins (by default 100) associated with these abstracts, and retrieve a STRING network for them.

```{r}
string.cmd = 'string pubmed query pubmed="jet-lag"'
commandsRun(string.cmd)
```
<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_jetlag.png){width=60%}
</center>

Which attribute column do you get in the Node Table for a PubMed query compared to a disease query? (Hint: check the last columns.)

# Exercise 2
In this exercise, we are going to use the stringApp to query the DISEASES database for proteins associated with ovarian cancer, retrieve a STRING network for them, and explore the resulting network.

Close the current session in Cytoscape:
```{r}
closeSession(save.before.closing=FALSE)
```

## Retrieve disease network
First, we will run a disease query for *ovarian cancer*:

```{r}
string.cmd = 'string disease query disease="ovarian cancer" limit=250'
commandsRun(string.cmd)
```

The retrieved network contains a lot of additional information associated with the nodes and edges, such as the protein sequence, tissue expression data, subcellular localization, disease score (Node Table) as well as the confidence scores for the different interaction evidences (Edge Table). In the next few steps, we will explore these data.

Create a dataframe containing the *disease score* and sort it descending by values to see the highest disease scores:

```{r}
disease.score <- getTableColumns('node', "stringdb::disease score")
disease.score.sorted <- disease.score[order(-disease.score$`stringdb::disease score`),,drop=FALSE]
head(disease.score.sorted)
```

You can highlight the top nodes by selecting the corresponding rows in the table. Here we are selecting the top few nodes based on disease score:

```{r}
selectNodes(rownames(head(disease.score.sorted)))
```

## Continuous color mapping
The stringApp automatically retrieves information about which compartments the proteins are located from the COMPARTMENTS database. Cytoscape allows you to map attributes of the nodes and edges to visual properties such as node color and edge width. Here, we will map the subcellular localization data for nucleus to the node color.

First, let's remove the String style:
```{r}
deleteStyleMapping(style.name = 'STRING style v1.5 - ovarian cancer', visual.prop = 'NODE_CUSTOMGRAPHICS_1')
```

For each cellular compartment, a score is defined for each node. Next, let's define a dataframe corresponding to the score for nucleus:
```{r}
nucleus.nodes <- getTableColumns('node', 'compartment::nucleus')
```

Define the bounds of the values in the *compartment::nucleus* column:
```{r}
nucleus.min <- min(nucleus.nodes, na.rm = T)
nucleus.max <- max(nucleus.nodes, na.rm = T)
nucleus.center <- nucleus.min + (nucleus.max - nucleus.min)/2
```

```{r}
data.values = c(nucleus.min, nucleus.center, nucleus.max)
node.colors <- c(brewer.pal(length(data.values), "YlOrRd"))
setNodeColorMapping('compartment::nucleus', data.values, node.colors, style.name = "STRING style v1.5 - ovarian cancer")
```

Because many proteins are located in the nucleus, we will identify the proteins with highest confidence of 5 and create a subnetwork.

```{r}
top.nodes <- row.names(nucleus.nodes)[which(nucleus.nodes[,1]>=5)]
createSubnetwork(top.nodes,subnetwork.name ='nucleus score 5')
```

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_nucleus-sub.png){width=60%}
</center>

# Exercise 3
In this exercise, we will work with a list of 541 proteins associated with epithelial ovarian cancer (EOC) as identified by phosphoproteomics in the study by [Francavilla et al](https://doi.org/10.1016/j.celrep.2017.03.015). An adapted table with the data from this study is available [here](https://raw.githubusercontent.com/cytoscape/cytoscape-automation/master/for-scripters/R/notebooks/stringApp/Francavilla2017CellRep.tsv).

## Protein network retrieval

Here, we run a query with the first column (UniProt IDs) in the table:

```{r}
eoc.df <- read.table("https://raw.githubusercontent.com/cytoscape/cytoscape-automation/master/for-scripters/R/notebooks/stringApp/Francavilla2017CellRep.tsv", header = TRUE, sep = "\t", quote="\"", stringsAsFactors = FALSE, check.names = FALSE)
```

```{r}
string.cmd = paste('string protein query query="', paste(eoc.df$UniProt, collapse = '\n'), '"', sep = "")
commandsRun(string.cmd)
```

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_eoc.png){width=60%}
</center>

How many nodes and edges are there in the resulting network? Do the proteins all form a connected network? Why?

Since we will need to refer to this network later, let's give it a more specific name:
```{r}
renameNetwork("EOC - data")
```

Cytoscape provides several network layout options.
For example, you can try the **Degree Sorted Circle Layout**

```{r}
layoutNetwork('degree-circle')
```

and the **Prefuse Force Directed Layout** with **score** as edge weight

```{r}
layoutNetwork('force-directed edgeAttribute="score"')
```

*Can you find a layout that allows you to easily recognize patterns in the network? Try the Edge-weighted Spring Embedded (Kamada-Kawai) Layout with the attribute ‘score’, which is the combined STRING interaction score.*

```{r}
layoutNetwork('kamada-kawai edgeAttribute="score"')
```

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_eoc-kamada-kawai.png){width=60%}
</center>

**Note that [yFiles Layout Algorithms App](http://apps.cytoscape.org/apps/yfileslayoutalgorithms) does not support any automation.**

## Discrete color mapping

Cytoscape allows you to map attributes of the nodes and edges to visual properties such as node color and edge width.
Here, we will map the target family data to the node color.

Lets change the node **Fill Color** with **target family** column in the node table. We can also remove the String styling for nodes to better see our data:

```{r}
deleteStyleMapping(style.name = 'STRING style v1.5', visual.prop = 'NODE_CUSTOMGRAPHICS_1')
column <- "target::family"
values <- c('Kinase', 'GPCR')
colors <- c('#FF0000', '#0000FF')
setNodeColorDefault('#CCCCCC', style.name = "STRING style v1.5")
setNodeColorMapping(column, values, colors, mapping.type = "d", style.name = "STRING style v1.5")
```

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_style.png){width=60%}
</center>
*How many of the proteins in the network are kinases?*

Note that the retrieved network contains a lot of additional information associated with the nodes and edges, such as the protein sequence, tissue expression data (Node Table) as well as the confidence scores for the different interaction evidences (Edge Table). In the following steps, we will explore these data using Cytoscape.

## Data import

Network nodes and edges can have additional information associated with them that we can load into Cytoscape and use for visualization.
We already imported the data from an Excel spreadsheet derived from data provided in the paper mentioned above. Here we check it again with:

```{r}
head(eoc.df)
```

Now we need to map unique identifiers between the entries in the data and the nodes in the network. The key point of this is to identify which nodes in the network are equivalent to which entries in the table. This enables mapping of data values into visual properties like Fill Color and Shape. This kind of mapping is typically done by comparing the unique Identifier attribute value for each node (Key Column for Network) with the unique Identifier value for each data value (key symbol).

The **Key Column** for Network allows you to set the node attribute column that is to be used as key to map to.
In this case it is **query term** because this attribute contains the UniProt accession numbers you entered when retrieving the network.

To import the node attributes file into Cytoscape, run

```{r}
loadTableData(as.data.frame(eoc.df), data.key.column = "UniProt", table = "node", table.key.column = "query term")
```

If there is a match between the value of a Key in the dataset and the value the Key Column for Network field in the network, all attribute–value pairs associated with the element in the dataset are assigned to the matching node in the network. 
You will find the imported columns at the end of the Node Table.

## Continuous color mapping

Now, we want to color the nodes according to the quantitative phosphorylation data (log ratio) between disease and healthy tissues for the most significant site for each protein:

```{r}
logratio.score.table <- getTableColumns('node', "EOC vs FTE&OSE")
```

Since this is a numeric value, we will use the Continuous Mapping as the Mapping Type, and set a color gradient for how abundant each protein is.

```{r}
logratio.min <- min(logratio.score.table, na.rm = T)
logratio.max <- max(logratio.score.table, na.rm = T)
```

```{r}
logratio.data.values = c(logratio.min, 0, logratio.max)
logratio.node.colors <- c(brewer.pal(length(logratio.data.values), "RdBu"))
```

```{r}
setNodeColorMapping('EOC vs FTE&OSE', logratio.data.values, logratio.node.colors, style.name = "STRING style v1.5")
```
The color gradient blue-white-red gives a visualization of the negative-to-positive abundance ratio.

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_eoc-final.png){width=60%}
</center>

*Are the up-regulated nodes grouped together?*

## Functional enrichment

Next, we will retrieve functional enrichment for the proteins in our network.

The STRING app has built-in enrichment analysis functionality, which includes enrichment for GO Process, GO Component, GO Function, InterPro, KEGG Pathways, and PFAM.

First, we will run the enrichment on the whole network, against the genome: 

```{r}
string.cmd = 'string retrieve enrichment allNetSpecies="Homo sapiens", background=genome  selectedNodesOnly="false"'
commandsRun(string.cmd)
string.cmd = 'string show enrichment'
commandsRun(string.cmd)
```

A new STRING Enrichment tab will appear in the Table Panel on the bottom. 
It contains a table of enriched terms and corresponding information for each enrichment category.

*Which are the three most statistically significant terms?*

To explore only specific types of terms, e.g. GO terms, and to remove redundant terms from the table, we are going to filter the table to only show GO Process:

```{r}
string.cmd = 'string filter enrichment categories="GO Process", overlapCutoff = "0.5", removeOverlapping = "true"'
commandsRun(string.cmd)
```

In this way, you will see only the statistically significant GO terms that do not represent largely the same set of proteins within the network.

*Do the functional terms assigned to a protein correlate with whether it is up- or down-regulated?*

Next, we will visualize the top-3 enriched terms in the network by using split charts.

```{r}
string.cmd = 'string show charts'
commandsRun(string.cmd)
```

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_eoc-enrich.png){width=60%}
</center>

*Do these terms give good coverage of the proteins in network?*

Finally, you can save the list of enriched terms and associated p-values as a text file locally. Note that this will export a cvs file to your current working directory.

```{r}
commandsPOST(paste0('table export table="STRING Enrichment: All" options="CSV" outputFile="',paste(getwd(),"string-enrichment-all.csv",sep = "/"),'"'))
```

## Network overlap

Cytoscape provides functionality to merge two or more networks, building either their union, intersection or difference. We will now merge the "ovarian cancner" network we have from the DISEASES query with the one we have from the data, so that we can identify the overlap between them. 

```{r}
mergeNetworks(c("String Network", "String Network - ovarian cancer"), "Merged Network", operation="intersection")
```

How many nodes are in the intersection?

Now we will make the union of the intersection network, which contains the disease scores, and the experimental network. Make sure that the new merged network has the same number of nodes and edges as ‘String Network’, and that some nodes have a disease score.
```{r}
mergeNetworks(c("Merged Network", "String Network"), "Union")
```
Now, we can change the visualization of the merged network to be able to identify high disease score proteins. Specifically, we will change the size and color of the nodes to reflect their disease score. 

```{r}

disease.score.table <- getTableColumns('node', "stringdb::disease score")

min <- min(disease.score.table, na.rm = T)
max <- max(disease.score.table, na.rm = T)
mid <- min + (max - min)/2

node.data.values = c(min, max)
node.sizes = c(35,50)
setNodeSizeMapping("stringdb::disease score", node.data.values, node.sizes, default.size="30", style.name = "default")

color.data.values = c(min, mid, max)
node.colors <- c(brewer.pal(length(color.data.values), "Blues"))
setNodeColorMapping("stringdb::disease score", color.data.values, node.colors, style.name = "default")
```

Let's focus on only the connected nodes:

```{r}
createSubnetwork(edges='all', subnetwork.name='Union sub')
```

<center>
![](https://cytoscape.github.io/cytoscape-automation/for-scripters/R/notebooks/data/img/stringApp_union.png){width=80%}
</center>

<!-- # Exercise 4 -->
<!-- In this exercise, we will retrieve virus-host networks for two closely related viruses, merge them into a single network, and then will retrieve the functional enrichment for the host proteins in this network. -->
