Making sense of gene and proteins lists with functional enrichment analysis

14 STRINGdb

STRING is a database of known and predicted protein-protein interactions, including both direct (physical) and indirect (functional) associations, drawn from experimental repositories, computational prediction methods, and text-mined publications. Each interaction is given a combined confidence score integrating all of the underlying evidence. STRING v12.0 covers over 59 million proteins from more than 12,500 organisms.

Earlier in this workshop you used the STRING website directly to explore both a protein-protein interaction network and functional enrichment for a gene list. The STRINGdb R/Bioconductor package gives scripted, reproducible access to the same underlying database and, being an R package, integrates directly with your existing gene lists, expression values, and enrichment workflows. STRINGdb uses R’s Reference Classes system, so most functionality is accessed as methods on a STRINGdb object (string_db$method_name(...)) rather than as standalone functions, and represents the network internally using igraph.

This chapter is based on the STRINGdb package vignette, reworked to use our own Pezzini et al. 2017 SH-SY5Y neuronal differentiation dataset in place of the vignette’s example data, so that you can compare the R output directly against what you saw on the STRING website earlier.

 

In this chapter we will:

  1. Initialise a STRINGdb object for human (NCBI taxonomy ID 9606)
  2. Load the differential expression results and map gene symbols to STRING identifiers
  3. Visualise the protein-protein interaction network
  4. Colour the network by up/down-regulation with the payload mechanism
  5. Run functional enrichment across GO, KEGG, Reactome and other databases
  6. Set a custom statistical background
  7. Retrieve functional annotations for the gene set
  8. Cluster the network and visualise individual clusters
  9. Query specific proteins of interest: interactions, interaction partners, paralogs and homologs

All of the data used in this chapter is provided in the data/ folder of this workshop’s materials, so there’s nothing extra to download - just read on and run the code chunks below in order (in RStudio, open this same .Rmd file directly to run and edit the code yourself). Note that several STRINGdb methods (network images, enrichment figures, the payload colouring) query the STRING web server directly, so an internet connection is required.

 

14.1 Load R packages

14.2 1. Initialise the STRINGdb object

To begin, you should know the NCBI taxonomy identifier of your organism of interest (e.g. 9606 for human, 10090 for mouse). If you don’t know it, search NCBI Taxonomy or the STRING organisms list.

In the constructor below we also set:

  • score_threshold: the minimum combined confidence score (0-1000) for an interaction to be loaded (400 is STRING’s own “medium confidence” default, and matches what you used on the website)
  • network_type: "full" for the full functional network (as opposed to "physical", which restricts to physical complexes only)
  • link_data: "detailed", so that individual evidence scores (e.g. experimental, database, text-mining) are loaded alongside the combined score

14.3 2. Explore the available methods

Since STRINGdb uses Reference Classes, all its functionality is exposed as methods on the string_db object. Two commands are useful for exploring what’s available:

For every method used below, you can call STRINGdb$help("method_name") for the full list of parameters and defaults.

14.4 3. Load the differential expression data

We reuse the same differential expression results behind the DE_198set gene list you used for ORA in g:Profiler and STRING earlier: genes from the Pezzini et al. 2017 SH-SY5Y differentiation experiment, each with a log2 fold change and an FDR-adjusted p-value.

The earlier DE_198set gene list was obtained by filtering this table for FDR < 0.0001 and an absolute log2 fold change ≥ 1 (i.e. at least a 2-fold change). We recreate that same filtered set here, keeping the pvalue (we’ll use FDR as our p-value, as in the STRING web activity) and logFC columns, so we can later colour the network by direction of change:

14.5 4. Map gene identifiers to STRING IDs

As a first step, we map our gene identifiers to STRING database identifiers using the map method. Here we map from HUGO gene symbols, but map also supports other common identifiers (e.g. Entrez Gene ID, Ensembl protein, RefSeq transcript, etc).

The map method adds a STRING_id column to the data frame passed as its first argument.

If a warning about unmapped genes is printed, this is expected: not every symbol will resolve to a STRING identifier (for example, deprecated gene symbols or non-protein-coding genes). Setting removeUnmappedRows = FALSE instead would keep the unmapped rows so you can inspect them manually.

14.6 5. Visualise the interaction network

We extract the STRING identifiers as our hits, and plot the network. The p-value printed above the plot is the probability of observing at least this many interactions among these proteins by chance - the same statistic reported on the STRING website’s network page.

Compare this network to the one you generated on the STRING website earlier using the same DE_198set gene list - the nodes, edges, and clustering should look familiar!

If you want more control over the appearance of the network image, get_png and get_link expose additional styling options, letting you switch between evidence/confidence edge styles, hide labels or disconnected proteins, use a flat node design, and adjust label size.

14.7 6. Payload mechanism: colour by direction of change

STRINGdb supports the “payload” mechanism used on the website to add a coloured halo around network nodes - for example, colouring down-regulated genes green and up-regulated genes red. This works by posting the colour information to the STRING server, then requesting the network image again with the resulting payload_id.

14.8 7. Functional enrichment

get_enrichment computes enrichment across Gene Ontology (Process, Function, Component), KEGG, Reactome, PubMed publications, UniProt Keywords, and PFAM/InterPro/SMART domains, in a single call. It uses a hypergeometric test with Benjamini-Hochberg FDR correction - the same underlying statistics as the STRING website enrichment tab you used earlier.

Question: How closely do these results match the enrichment table you downloaded from the STRING website earlier for the same gene list?

You can also retrieve a ready-made enrichment figure directly from STRING for a chosen category:

14.9 8. Setting a custom statistical background

As discussed earlier in the workshop, using an appropriate background/universe is important for avoiding false positives in enrichment tests. By default, get_enrichment tests against the whole proteome. If your experiment only measured a defined subset of genes (as ours does), you should set that as the background instead:

Once a background is set, it is used automatically by subsequent calls to get_enrichment on the same string_db object - you do not need to create a new object.

Compare this to the enrichment table from step 7: results are often less dramatic once tested against a realistic background rather than the whole proteome, because the test is now calibrated against the genes that were actually measurable in this experiment.

14.10 9. Retrieve functional annotations

If you simply want to know which terms are assigned to your proteins, rather than which are statistically enriched, use get_annotations. This returns terms from most categories associated with your protein set (KEGG and PubMed are excluded, for licensing and output-size reasons respectively).

14.11 10. Clustering

The underlying igraph package provides several clustering/community-detection algorithms ("fastgreedy", "walktrap", "spinglass", "edge.betweenness"), which get_clusters wraps in an easy-to-use function returning a list of clusters.

We visualise the first four clusters individually:

14.12 11. Querying specific proteins of interest

CRABP2 and SEMA3A are among the most significantly differentially expressed genes in this dataset (see the DE table above). Let’s look at them specifically, using mp (“map proteins”) - a convenient alternative to map for mapping just one or a few gene symbols directly to STRING identifiers:

14.12.1 Interactions between named proteins

get_interactions returns a data frame describing the interactions among the queried proteins, with STRING identifiers (from/to), preferred names (from_name/to_name), the combined_score, and any additional evidence scores available under the current link_data setting.

14.12.2 Interaction partners

To retrieve the first-shell interaction partners of one or more proteins from the locally loaded network, use get_interaction_partners. The output has the same structure as get_interactions.

You can restrict to interactions supported by a particular type of evidence - here, only those with experimental support:

14.12.3 Paralogs and homologs

STRING stores all-against-all sequence alignments within species, so you can retrieve paralogs of a protein of interest:

STRING also stores cross-species best hits (by bitscore). get_homologs_besthits retrieves these:

… or restrict to a specific target species, e.g. mouse (taxonomy ID 10090):

14.13 12. Additional protein information

You can retrieve a table of every protein in the STRING database for your species of interest, including preferred name, protein size, and a short description: