%\VignetteIndexEntry{STRINGdb Vignette}
%\VignetteDepends{}
%\VignetteKeywords{STRINGdb}
%\VignettePackage{STRINGdb}
\documentclass[10pt]{article}
\textwidth=6.2in
\textheight=8.5in
\oddsidemargin=0.2in
\evensidemargin=0.2in
\headheight=0in
\headsep=0in

\begin{document}
\SweaveOpts{concordance=TRUE, strip.white=false}
\title{STRINGdb Package Vignette}
\author{Andrea Franceschini}
\date{15 March 2015}
\maketitle

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

\section{INTRODUCTION}


STRING (https://www.string-db.org) is a database of known and predicted protein-protein interactions. 
The interactions include direct (physical) and indirect (functional) associations. 
The database contains information from numerous sources, including experimental repositories, computational prediction methods and public text collections. Each interaction is associated with a combined confidence score that integrates the various evidences.  We currently cover over 24 milions proteins from 5090 organisms.

As you will learn in this guide, the STRING database can be usefull to add meaning to list of genes (e.g. the best hits coming out from a screen or the most differentially expressed genes coming out from a Microarray/RNAseq experiment.)

We provide the STRINGdb R package in order to facilitate our users in accessing the STRING database from R.
In this guide we explain, with examples, most of the package's features and functionalities.\\

In the STRINGdb R package we use the new ReferenceClasses of R (search for "ReferenceClasses" in the R documentation.).  Besides we make use of the iGraph package (http://igraph.sourceforge.net) as a data structure to represent our protein-protein interaction network. \\

To begin, you should first know the NCBI taxonomy identifiers of the organism on which you have performed the experiment (e.g. 9606 for Human, 10090 for mouse). 
If you don't know that, you can search the NCBI Taxonomy (http://www.ncbi.nlm.nih.gov/taxonomy)  or start looking at our species table (that you can also use to verify that your organism is represented in the STRING database). \\
Hence, if your species is not Human (i.e. our default species), you can find it and their taxonomy identifiers on STRING webpage under the 'organisms' section (https://string-db.org/cgi/input.pl?input\_page\_active\_form=organisms) or download the full list in the download section of STRING website.

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

<<initialization>>=
library(STRINGdb)
string_db <- STRINGdb$new( version="11.5", species=9606, 
                           score_threshold=200, network_type="full", input_directory="")
@

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

As it has been shown in the above commands, you start instantiating the STRINGdb reference class.
In the constructor of the class you can also define the STRING version to be used and a threshold for the combined scores of the interactions, such that any interaction below that threshold is not loaded in the object (by default the score threshold is set to 400).\\
You can also specify the network type "functional" for full functional STRING network or "physical" for physical subnetwork, which link only the proteins which share a physical complex.\\ 
Besides, if you specify a local directory to the parameter input-directory, the database files will be downloaded into this directory and most of the methods can be used off-line. Otherwise, the database files will be saved and cached in a temporary directory that will be cleaned automatically when the R session is closed.\\\\

For a better understanding of the package two other commands can be useful:

<<help>>=
STRINGdb$methods()              # To list all the methods available.
STRINGdb$help("get_graph")      # To visualize their documentation.
@

For all the methods that we are going to explain below, you can always use the help function in order to get additional information/parameters with respect to those explained in this guide.\\\\


As an example, we use the analyzed data of a microarray study taken from GEO (Gene Expression Omnibus, GSE9008).
This study investigates the activity of Resveratrol, a natural phytoestrogen found in red wine and a variety of plants, in A549 lung cancer cells. 
Microarray gene expression profiling after 48 hours exposure to Revestarol has been performed and compared to a control composed by A549 lung cancer cells threated only with ethanol.
This data is already analyzed for differential expression using the limma package: the genes are sorted by fdr corrected pvalues and the log fold change of the differential expression is also reported in the table.

<<load_data>>=
data(diff_exp_example1)
head(diff_exp_example1)
@

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

As a first step, we map the gene names to the STRING database identifiers using the "map" method.
In this particular example, we map from gene HUGO names, but our mapping function supports several other common identifiers (e.g. Entrez GeneID, ENSEMBL proteins, RefSeq transcripts ... etc.).\\\\
The map function adds an additional column with STRING identifiers to the dataframe that is passed as first parameter.

<<map>>=
example1_mapped <- string_db$map( diff_exp_example1, "gene", removeUnmappedRows = TRUE )
@
\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

As you may have noticed, the previous command prints a warning showing the number of genes that we failed to map. In this particular example, we cannot map all the probes of the microarray that refer to position of the chromosome that are not assigned to a real gene (i.e. all the LOC genes). If we remove all these LOC genes before the mapping we obtain a much lower percentage of unmapped genes (i.e. < 6 \%). \\
If you set to FALSE the "removeUnmappedRows" parameter, than the rows which corresponds to unmapped genes are left and you can manually inspect them.\\
Finally, we extract the most significant 200 genes and we produce an image of the STRING network for those.
The image shows clearly the genes and how they are possibly functionally related.
On the top of the plot, we insert a pvalue that represents the probability that you can expect such an equal or greater number of interactions by chance.

<<echo=false>>=
options(SweaveHooks=list(fig=function()
par(mar=c(2.1, 0.1, 4.1, 2.1))))
@
\setkeys{Gin}{width=1.2\linewidth}

<<get_hits>>=
hits <- example1_mapped$STRING_id[1:200]  
@


<<plot_network, fig=TRUE>>=
string_db$plot_network( hits )
@



\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

\section{PAYLOAD MECHANISM}

This R library provides the ability to interact with the STRING payload mechanism. The payload appears as an additional colored "halo" around the bubbles.

For example, this allows to color in green the genes that are down-regulated and in red the genes that are up-regulated.
For this mechanism to work, we provide a function that posts the information on our web server.

<<add_diff_exp_color>>=
# filter by p-value and add a color column 
# (i.e. green down-regulated gened and red for up-regulated genes)
example1_mapped_pval05 <- string_db$add_diff_exp_color( subset(example1_mapped, pvalue<0.05), 
                                                            logFcColStr="logFC" )    
@

<<post_payload>>=
# post payload information to the STRING server
payload_id <- string_db$post_payload( example1_mapped_pval05$STRING_id, 
                                        colors=example1_mapped_pval05$color )
@

<<plot_halo_network, fig=TRUE>>=
# display a STRING network png with the "halo"
string_db$plot_network( hits, payload_id=payload_id )
@



\begin{Schunk}
\begin{Sinput}


\end{Sinput}
\end{Schunk}

\section{ENRICHMENT}

We provide a method to compute the enrichment in Gene Ontology (Process, Function and Component), KEGG and Reactome pathways, PubMed publications, UniProt Keywords, and PFAM/INTERPRO/SMART domains for your set of proteins all in one simple call. The enrichment itself is computed using an hypergeometric test and the FDR is calculated using Benjamini-Hochberg procedure.

<<enrichment>>=
enrichment <- string_db$get_enrichment( hits )
head(enrichment, n=20)
@

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

If you have performed your experiment on a predefined set of proteins, it is important to run the enrichment statistics using that set as a background (otherwise you would get a wrong p-value !). Hence, before to launch the method above, you may want to set the background:

<<background, eval=FALSE>>=
backgroundV <- example1_mapped$STRING_id[1:2000]   # as an example, we use the first 2000 genes                                                    
string_db$set_background(backgroundV)
@

You can also set the background when you instantiate the STRINGdb object:

<<new_background_inst, eval=FALSE>>=
string_db <- STRINGdb$new( score_threshold=200, backgroundV = backgroundV )
@

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

If you just want to know terms are assigned to your set of proteins (and not necessary enriched) you can use "get\_annotations" method. This method will output all the terms from most of the categories (the exceptions are KEGG terms due to licensing issues and PubMed due to the size of the output) that are associated with your set of proteins. 

<<enrichment>>=
annotations <- string_db$get_annotations( hits )
head(annotations, n=20)
@

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}






\section{CLUSTERING}

The iGraph package provides several clustering/community algorithms: "fastgreedy", "walktrap", "spinglass", "edge.betweenness".
We encapsulate this in an easy-to-use function that returns the clusters in a list.

<<clustering1>>=
# get clusters
clustersList <- string_db$get_clusters(example1_mapped$STRING_id[1:600])
@

<<echo=false>>=
options(SweaveHooks=list(fig=function()
par(mar=c(2.1, 0.1, 4.1, 2.1))))
@
\setkeys{Gin}{width=1.2\linewidth}
<<clustering2, fig=TRUE>>=
# plot first 4 clusters
par(mfrow=c(2,2))
for(i in seq(1:4)){
 string_db$plot_network(clustersList[[i]])
}
@

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}


\section{ADDITIONAL PROTEIN INFORMATION}

You can get a table that contains all the proteins that are present in our database of the species of interest. 
The protein table also include the preferred name, the size and a short description of each protein.

<<proteins>>=
string_proteins <- string_db$get_proteins()
@
\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

In the following section we will show how to query STRING with R on some specific proteins.
In the examples, we will use the famous tumor proteins TP53 and ATM .\\\\


First we need to get the STRING identifier of those proteins, using our mp method:
<<atmtp>>=
tp53 = string_db$mp( "tp53" )
atm = string_db$mp( "atm" )
@

The mp method (i.e. map proteins) is an alternative to our map method, to be used when you need to map only one or few proteins.\\
It takes in input a vector of protein aliases and returns a vector with the STRING identifiers of those proteins.

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}



Using the following method, you can see the proteins that interact with one or more of your proteins:
<<neighbors, eval=FALSE>>=
string_db$get_neighbors( c(tp53, atm) )
@

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

It is also possible to retrieve the interactions that connect certain input proteins between each other.\\
Using the "get\_interactions" method we can clearly see that TP53 and ATM interact with each other with a good evidence/score.

<<interactions>>=
string_db$get_interactions( c(tp53, atm) )
@

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}


STRING provides a way to get homologous proteins: in our database we store ALL-AGAINST-ALL alignments within all 5090 organisms. You can retrive all of the paralogs of the protein using "get\_paralogs" method.  

<<paralogs, eval=FALSE>>=
# Get all homologs of TP53 in human.
string_db$get_paralogs(tp53)
@
\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

STRING also stores best hits (as measured by bitscore) between the proteins from different species. "get\_homologs\_besthits" lets you retrieve these homologs.

<<Closest homologs from other species, eval=FALSE>>=
# get the best hits of the following protein in all the STRING species
string_db$get_homologs_besthits(tp53)
@
\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}

... or you can specify the species of interest (i.e. all the blast hits): 

<<homologs_besthits in target species, eval=FALSE>>=
# get the homologs of the following two proteins in the mouse (i.e. species_id=10090)
string_db$get_homologs_besthits(c(tp53, atm), target_species_id=10090, bitscore_threshold=60)
@

\begin{Schunk}
\begin{Sinput}

\end{Sinput}
\end{Schunk}


\section{CITATION}

Please cite:\\\\

Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, Doncheva NT, Legeay M, Fang T, Bork P, Jensen LJ, von Mering C. 'The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets.' Nucleic Acids Res. 2021 Jan 8;49(D1):D605-12

\end{document}