Geneplast: evolutionary rooting and plasticity analysis of orthologous groups.
In geneplast: Evolutionary and plasticity analysis of orthologous groups

Overview

Geneplast is designed for evolutionary and plasticity analysis based on the distribuion of orthologous groups in a given species tree. It uses Shannon information theory to estimate the Evolutionary Plasticity Index (EPI) [@Dalmolin2011, @Castro2008].

Figure 1 shows a toy example to illustrate the analysis. The observed itens in Figure 1a are distributed evenly among the different species (i.e. high diversity), while Figure 1b shows the opposite case. The diversity is given by the normalized Shannon’s diversity and represents the distribution of orthologous and paralogous genes in a set of species. High diversity represents an homogeneous distribution among the evaluated species, while low diversity indicates that few species concentrate most of the observed orthologous genes.

The EPI characterizes the evolutionary history of a given orthologous group (OG). It accesses the distribution of orthologs and paralogs and is defined as,

$$EPI=1-\frac{H\alpha}{\sqrt{D\alpha}}, (1)$$

where Dα represents the OG abundance and Hα the OG diversity. Low values of Dα combined with high values for Hα indicates an orthologous group of low plasticity, that is, few OG members distributed over many species. It also indicates that the OG might have experienced few modifications (i.e. duplication and deletion episodes) during the evolution. Note that 0 ≤ Hα ≤ 1 and Dα ≥ 1. As a result, 0 ≤ EPI ≤ 1. For further information about the EPI, please see [@Dalmolin2011].

title Figure 1. Toy examples showing the distribution of orthologous and paralogous genes in a given species tree. (a) OG of low abundance Dα, high diversity (Hα) and consequently low plasticity (PI). In this hypothetical case, the OG comprises orthologous genes observed in all species, without apparent deletion or duplication episodes. (b) in this example the OG is observed in many species, but not all, with many paralogs in some of them. Green numbers represents the number of orthologous genes in each species.

geneplast also implements a new algorithm called Bridge in order to interrogate the evolutionary root of a given gene based on the distribution of its orthologs. The Bridge algorithm assesses the probability that an ortholog of a given gene is present in each last common ancestor (LCA) of a given species (in a given species tree). As a result, this approach infers the evolutionary root representing the gene emergence. The method is designed to deal with large scale queries in order to interrogate, for example, all genes annotated in a network (please refer to [@Castro2008] for a case study illustrating the advantages of using this approach).

To illustrate the rooting inference consider the evolutionary scenarios presented in Figure 2 for the same hypothetical OGs. These OGs comprise a number of orthologous genes distributed among 13 species, and the pattern of presence or absence is indicated by green and grey colours, respectively. Observe that at least one ortholog is present in all extant species in Figure 2a. To explain this common genetic trait, one possible evolutionary scenario could assume that the ortholog was present in the LCA of all species and was genetically transmitted up to the descendants. For this case, the evolutionary root might be placed at the bottom of the species tree (i.e. node g). The same reasoning can be done in Figure 2b, but with the evolutionary root placed at node d. The geneplast rooting pipeline is designed to infer the most consistent rooting scenario for the observed orthologs in a given species tree. The pipeline provides a consistency score called Dscore which estimates the stability of the inferred root, as well as an associated empirical p-value computed by permutation analysis.

title Figure 2. Possible evolutionary rooting scenarios for the same toy examples depitected in Figure 1. (a, b) Red circles indicate the evolutionary roots that best explain the observed orthologs in this species tree.

Quick start

The orthology data required to run geneplast is available in the gpdata.gs dataset. This dataset includes four objects containing information about Clusters of Orthologous Groups derived from the STRING database, release 9.1. geneplast can also be used with other sources of orthology information, provided that the input is set according to the gpdata.gs data structure (note: in order to reduce the processing time this example uses a subset of the STRING database).

library(geneplast)
data(gpdata.gs)

Evolutionary plasticity inference

The first step is to create an OGP object by running the gplast.preprocess function. This example uses 121 eukaryotic species from the STRING database and all OGs mapped to the genome stabilty gene network [@Castro2008]. Next, the gplast function perform the plasticity analysis and the gplast.get returns the results:

1 - Create an object of class OGP.

ogp <- gplast.preprocess(cogdata=cogdata, sspids=sspids, cogids=cogids, verbose=FALSE)

2 - Run the gplast function.

ogp <- gplast(ogp, verbose=FALSE)

3 - Get results.

res <- gplast.get(ogp,what="results")
head(res)

The results are returned in a 3-column data.frame with OG ids (cogids) identified in row.names. Columns are named as abundance, diversity, and plasticity.

The metric abundance simply indicates the ratio of orthologs and paralogs by species. For example, KOG0011 cromprises 201 genes distributed in 116 eukaryotic species, with a resulting abundance of 1.7328. Abundance of 1 indicates an one-to-one orthology relationship, while high abundance denotes many duplication episodes on the OG’s evolutionary history. Diversity is obtained applying normalized Shannon entropy on orthologous distribution and Plasticity is obtained by EPI index, as described equation (1).

Evolutionary rooting inference

The rooting analysis starts with an OGR object by running the groot.preprocess function. This example uses all OGs mapped to the genome stability gene network using H. sapiens as reference species [@Castro2008] and is set to perform 100 permutations for demonstration purposes (for a full analysis, please set Permutations≥1000). Next, the groot function performs the rooting analysis and the results are retrieved by groot.get, which returns a data.frame listing the root of each OG evaluated by the groot method. The pipeline also returns the inconsistency score, which estimates the stability of the rooting analysis, as well as the associated empirical p-value. Additionally, the groot.plot function allows the visualization of the inferred root for a given OG (e.g. Figure 3) and the LCAs for the reference species (Figure 4).

1 - Create an object of class OGR.

ogr <- groot.preprocess(cogdata=cogdata, phyloTree=phyloTree, spid="9606", cogids=cogids, verbose=FALSE)

2 - Run the groot function.

set.seed(1)
ogr <- groot(ogr, nPermutations=100, verbose=FALSE)

3 - Get results.

res <- groot.get(ogr,what="results")
head(res)

4 - Check the inferred root of a given OG

groot.plot(ogr,whichOG="NOG40170")

5 - Visualization of the LCAs for the reference species in the analysis (i.e. H. sapiens)

groot.plot(ogr,plot.lcas = TRUE)

title Figure 3. Inferred evolutionary rooting scenario for NOG40170. Monophyletic groups are ordered to show all branches of the tree below the queried species in the analysis.

title Figure 4. Visualization of the LCAs for the reference species in the analysis.

Case studies

High-throughput rooting inference

This example shows how to assess all OGs annotated for H. sapiens.

1 - Load orthogy data from the geneplast.data.string.v91 package.

# source("https://bioconductor.org/biocLite.R")
# biocLite("geneplast.data.string.v91")
library(geneplast.data.string.v91)
data(gpdata_string_v91)

2 - Create an object of class 'OGR' for a reference 'spid'.

ogr <- groot.preprocess(cogdata=cogdata, phyloTree=phyloTree, spid="9606")

3 - Run the groot function and infer the evolutionary roots. Note: this step should take a long processing time due to the large number of OGs in the input data (also, nPermutations argument is set to 100 for demonstration purpose only).

ogr <- groot(ogr, nPermutations=100, verbose=TRUE)

Map rooting information on PPI networks

This example aims to show the evolutionary root of a protein-protein interaction (PPI) network, mapping the appearance of each gene in a given species tree. The next steps show how to tranfer evolutionary rooting information from geneplast to a graph model. Note: to make this work the gene annotation available from the input PPI network needs to match the annotation available from the geneplast data (in this case, ENTREZ gene IDs are used to match the datasets).

1 - Load a PPI network and required packages. The igraph object called 'ppi.gs' provides PPI information for apoptosis and genome-stability genes [@Castro2008].

library(RedeR)
library(igraph)
library(RColorBrewer)
data(ppi.gs)

2 - Map rooting information on the igraph object.

g <- ogr2igraph(ogr, cogdata, ppi.gs, idkey = "ENTREZ")

3 - Adjust colors for rooting information.

pal <- brewer.pal(9, "RdYlBu")
color_col <- colorRampPalette(pal)(25) #set a color for each root!
g <- att.setv(g=g, from="Root", to="nodeColor", cols=color_col, na.col = "grey80", breaks = seq(1,25))

4 - Aesthetic adjusts for some graph attributes.

g <- att.setv(g = g, from = "SYMBOL", to = "nodeAlias")
E(g)$edgeColor <- "grey80"
V(g)$nodeLineColor <- "grey80"

5 - Send the igraph object to RedeR interface.

rdp <- RedPort()
calld(rdp)
resetd(rdp)
addGraph(rdp, g)
addLegend.color(rdp, colvec=g$legNodeColor$scale, size=15, labvec=g$legNodeColor$legend, title="Roots represented in Fig4")

6 - Get apoptosis and genome-stability sub-networks.

g1  <- induced_subgraph(g=g, V(g)$name[V(g)$Apoptosis==1])
g2  <- induced_subgraph(g=g, V(g)$name[V(g)$GenomeStability==1])

7 - Group apoptosis and genome-stability genes into containers.

myTheme <- list(nestFontSize=25, zoom=80, isNest=TRUE, gscale=65, theme=2)
addGraph(rdp, g1, gcoord=c(25, 50), theme = c(myTheme, nestAlias="Apoptosis"))
addGraph(rdp, g2, gcoord=c(75, 50), theme = c(myTheme, nestAlias="Genome Stability"))
relax(rdp, p1=50, p2=50, p3=50, p4=50, p5= 50, ps = TRUE)

title Figure 5. Inferred evolutionary roots of a protein-protein interaction network.

Map rooting information on regulatory networks

This example aims to show the evolutionary root of regulons [@Fletcher2013]. The idea is to map the appearance of each regulon (and the corresponding target genes) in a species tree. The next steps show how to tranfer evolutionary rooting information from geneplast to a graph model. Note: to make this work the gene annotation available from the input regulatory network needs to match the annotation available from the geneplast data (in this case, ENTREZ gene IDs are used to match the datasets).

1 - Load a TNI class object and required packages. The rtni1st object provides regulons available from the Fletcher2013b data package computed from breast cancer data [@Fletcher2013].

library(RTN)
library(Fletcher2013b)
library(RedeR)
library(igraph)
library(RColorBrewer)
data("rtni1st")

2 - Extract two regulons from rtni1st into an igraph object.

regs <- c("FOXM1","PTTG1")
g <- tni.graph(rtni1st, gtype = "rmap", tfs = regs)

3 - Map rooting information on the igraph object.

g <- ogr2igraph(ogr, cogdata, g, idkey = "ENTREZ")

4 - Adjust colors for rooting information.

pal <- brewer.pal(9, "RdYlBu")
color_col <- colorRampPalette(pal)(25) #set a color for each root!
g <- att.setv(g=g, from="Root", to="nodeColor", cols=color_col, na.col = "grey80", breaks = seq(1,25))

5 - Aesthetic adjusts for some graph attributes.

idx <- V(g)$SYMBOL %in% regs
V(g)$nodeFontSize[idx] <- 30
V(g)$nodeFontSize[!idx] <- 1
E(g)$edgeColor <- "grey80"
V(g)$nodeLineColor <- "grey80"

6 - Send the igraph object to RedeR interface.

rdp <- RedPort()
calld(rdp)
resetd(rdp)
addGraph( rdp, g, layout=NULL)
addLegend.color(rdp, colvec=g$legNodeColor$scale, size=15, labvec=g$legNodeColor$legend, title="Roots represented in Fig4")
relax(rdp, 15, 100, 20, 50, 10, 100, 10, 2, ps=TRUE)