BioNAR: Biological Network Analysis in R

Colin Mclean, Anatoly Sorokin, T. I. Simpson, J. Douglas Armstrong, Oksana Sorokina

Introduction

Proteomic studies typically generate a massive list of proteins being identified within a specific tissue, compartment or cell type, often accompanied by additional qualitative and/or quantitative information. Conversion of these data into meaningful biological insight requires processing in several stages to identify possible structural and/or functional dependencies.

One of the most popular ways of representing proteomic data is a protein-protein interaction network, which allows to study its topology and how it correlates with functional annotation mapped onto the network.

Many existing packages support different steps of the network building and analysis process, but few packages combine network analysis methodology into a single coherent pipeline.

We designed BioNAR to support a range of network analysis functionality, complementing existing R packages and filling the methodological gaps necessary to interrogate biomedical networks with respect to functional and disease domains. For that purpose, we do not implement network reconstruction directly (unless for synaptic networks), as other tools such as Cytoscape and Network Analyst do this already. Rathher, we provide a detailed topologically-based network analysis package, enabling the researcher to load networks generated from the lab’s own meta-data, thus making the tool as widely applicable and flexible as possible. We also provide a synaptic proteome network of our own for validation.

Overview of capabilities

The BioNAR’s pipeline starts with importing the graph of interest (typically built from nodes/proteins and edges/PPI interactions), and annotating its vertices with available metadata [annotate_vertex].

This is followed by the analysis of general network properties, such as estimating a network’s “scale-free” property. For this we used the R “PoweRlaw” package (version 0.50.0) (Gillespie, 2015) [FitDegree] and a network entropy analysis (Teschendorff et al, 2014) [getEntopy].

The package allows estimation of the main network vertex centrality measures: degree, betweenness centrality, clustering coefficient, semilocal centrality, mean shortest path, page rank, and standard deviation of the shortest path. Values for centrality values measures can be added as vertex attributes [calcCentrality] or returned as an R matrix [getCentralityMatrix], depending on user’s preferences. Any other vertex meta-data, which can be represented in matrix form, can also be stored as a vertex attribute.

To compare observed networks vertex centrality values against those of equivalently sized but randomised graphs, we support three varying randomisation models including G(n,p) Erdos-Renyi model , Barabasi-Albert model, and new random graph from a given graph by randomly adding/removing edges [getRandomGraphCentrality].

Additionally, to allow comparison of networks with different structures, we implemented normalized modularity measure (Parter et al., 2007, Takemoto, 2012, Takemoto, 2013, Takemoto and Borjigin, 2011)[normModularity].

BioNAR then proceeds to analyse the network’s community structure, which can be performed via nine different clustering algorithms, simultaneously [calcAllClustering] or with a chosen algorithm [calcClustering], community membership being stored as a vertex attribute. In situations where the network is dense and clusters are large and barely tractable, reclustering can be applied [calcReclusterMatrix]. The obtained community structure can be visualized with [layoutByCluster], and communiities further tested for robustness [getRobustness] by comparing against randomised networks. As a result, a consensus matrix can be estimated [makeConsensusMatrix], which is needed for the next step -identifying the “influential” or “bridging” proteins (Nepusz et al., 2008).

For this, we enabled a function for calculating the “bridgeness” [getBridgeness] metric, which takes into account the probability a vertex to belongs to different communities [getBridgeness], such that a vertex can be ranked under assumption the higher its community membership the more influence it has to the network topology and signaling (Han et al., 2004). “Bridgeness” can be plotted against any other centrality measure, e.g. semi-local centrality (plot is implemented), to enable useful indication of vertex (protein) importance within the network topology.

To provide a perspective of the molecular signature of multiple diseases (or biological functions) and how they might interact or overap at the network level, we implemented a disease-disease overlap analysis by measuring the mean shortest distance for each disease (⟨d⟩), using the shortest distance between each gene-disease association (GDA) to its next nearest GDA neighbor (Menche et al., 2015) [calcDiseasePairs], [runPermDisease] to be used for obtaining significance values. In the case example of the presynaptic network we found a set of neurological disorders to overlap with a high significance, e.g. AD-PD, SCH-BP, ASD-ID, and, indeed, their comorbidity was confirmed by the literature. Note that while developed for disease-disease correlation, the analysis can be performed using any in-house vertex meta-data, including biological function terms.

To study the distribution of the specific annotation(s) over a clustered graph (typically disease of biologial process/function), we enabled overrepresentation analysis [clusterORA], which helps identify the network communities enriched for specific function or disease, or any other annotation.

The case study illustrates the package functionality for the protein-protein interaction network for the presynaptic compartment of the synapse generated from Synaptic proteome database (Sorokina et al., 2021) with Synaptome.db package. The network has 1073 vertices and 3346 edges, step by step analysis is shown below.

Build the network

BioNAR allows building a network from a data frame, where the rows correspond to the edges of the graph; alternatively for our synaptic proteome exemple, a list of vertices (genes) is needed, for which the information will be retrieved from SynaptomeDB package.

Build a network from a given data frame

The command listed below builds a graph from provided data frame, simplifies the graph (removing multiple edges and loops) and return its MCC (maximum connected component)

file <- system.file("extdata", "PPI_Presynaptic.csv", package = "BioNAR")
tbl <- read.csv(file, sep="\t")
head(tbl)
#>   entA  entB We
#> 1  273  1759  1
#> 2  273  1785  1
#> 3  273 26052  1
#> 4  273   824  1
#> 5  273   161  1
#> 6  273  1020  1
gg <- buildNetwork(tbl)
summary(gg)
#> IGRAPH 76471fe UN-- 1780 6620 -- 
#> + attr: name (v/c)

Build a network given a node list extracted from SynaptomeDB

The command provides proxy for Synaptome.db package, enabling building the protein-protein interaction (PPI) network for entire synaptic compartment (pre-, or post-synaptic)

# Let's get the ID for presynaptic compartment
cid<-match('Presynaptic',getCompartments()$Name) 
cid
#> [1] 2
# Now we need to collect all the gene IDs for presinaptic  compartment
t<-getAllGenes4Compartment(cid) 
dim(t)
#> [1] 2929    8
head(t)
#> # A tibble: 6 × 8
#>   GeneID MGI       HumanEntrez MouseEntrez RatEntrez HumanName MouseName RatName
#>    <int> <chr>           <int>       <int>     <int> <chr>     <chr>     <chr>  
#> 1      1 MGI:1277…        1742       13385     29495 DLG4      Dlg4      Dlg4   
#> 2      2 MGI:88256         815       12322     25400 CAMK2A    Camk2a    Camk2a 
#> 3      3 MGI:96568        9118      226180     24503 INA       Ina       Ina    
#> 4      4 MGI:98388        6711       20742    305614 SPTBN1    Sptbn1    Sptbn1 
#> 5      5 MGI:88257         816       12323     24245 CAMK2B    Camk2b    Camk2b 
#> 6      6 MGI:1344…        1740       23859     64053 DLG2      Dlg2      Dlg2
# finally, build the graph using respecctive gene EntrezIDs as node IDs
gg<-buildFromSynaptomeByEntrez(t$HumanEntrez) 
summary(gg)
#> IGRAPH 2030fb9 UN-- 2275 9190 -- 
#> + attr: name (v/c)

Use a predifined network

Any predefined network stored as a graph file (e.g. .gml, .graphml) can be loaded for further analysis using Igraph’s functionality.

file <- system.file("extdata", "PPI_Presynaptic.gml", package = "BioNAR")
gg1 <- igraph::read.graph(file,format="gml")
summary(gg1)
#> IGRAPH 6cf5a71 UN-- 2169 8673 -- 
#> + attr: id (v/n), name (v/c), GeneName (v/c), UniProt (v/c), louvain
#> | (v/n), TopOntoOVG (v/c), TopOntoOVGHDOID (v/c)

Annotate the nodes with node attributes

As soon as the graph is loaded it can be annotated with any relevant annotations, such as protein names [annotateGeneNames], functionality [annotateGOont], disease associations [annotateTopOntoOVG], or any customized annotation set [annotate_vertex {BioNAR}]. We also provide two functional annotations for synaptic graphs based on published synaptic functional studies ([annotateSCHanno], and [annotatePresynaptic].

Gene name

Adding gene names to vertices.

gg<-annotateGeneNames(gg)
summary(gg)
#> IGRAPH 2030fb9 UN-- 2275 9190 -- 
#> + attr: name (v/c), GeneName (v/c)
head(V(gg))
#> + 6/2275 vertices, named, from 2030fb9:
#> [1] 273   6455  1759  1785  26052 1837
head(V(gg)$GeneName)
#> [1] "AMPH"   "SH3GL1" "DNM1"   "DNM2"   "DNM3"   "DTNA"

Diseases

Adding diseases associations to genes linked to Human Disease Ontology (HDO) terms extracted from the package (topOnto.HDO.db)[https://github.com/hxin/topOnto.HDO.db].

afile<-system.file("extdata", "flatfile_human_gene2HDO.csv", package = "BioNAR")
dis <- read.table(afile,sep="\t",skip=1,header=FALSE,strip.white=TRUE,quote="")
gg <- annotateTopOntoOVG(gg, dis)
summary(gg)
#> IGRAPH 2030fb9 UN-- 2275 9190 -- 
#> + attr: name (v/c), GeneName (v/c), TopOnto_OVG (v/c),
#> | TopOnto_OVG_HDO_ID (v/c)

Presynaptic functional annotation

Adding the presynaptic genes functional annotation derived from Boyken at al. (2013) doi:10.1016/j.neuron.2013.02.027.

sfile<-system.file("extdata", "PresynAn.csv", package = "BioNAR")
pres <- read.csv(sfile,skip=1,header=FALSE,strip.white=TRUE,quote="")
sgg <- annotatePresynaptic(gg, pres)
summary(sgg)

Functional annotation with Gene Ontology (GO)

GO annotation is specifically supported with the function [annotateGOont]:


ggGO <- annotateGOont(gg)
#however, functionality from GO: BP, MF,CC can be added

sfile<-system.file("extdata", "flatfile.go.BP.csv", package = "BioNAR")
goBP <- read.table(sfile,sep="\t",skip=1,header=FALSE,strip.white=TRUE,quote="")
sgg <- annotateGoBP(gg, goBP)
summary(sgg)
#> IGRAPH 2030fb9 UN-- 2275 9190 -- 
#> + attr: name (v/c), GeneName (v/c), TopOnto_OVG (v/c),
#> | TopOnto_OVG_HDO_ID (v/c), GO_BP_ID (v/c), GO_BP (v/c)

sfile<-system.file("extdata", "flatfile.go.MF.csv", package = "BioNAR")
goMF <- read.table(sfile,sep="\t",skip=1,header=FALSE,strip.white=TRUE,quote="")
sgg <- annotateGoMF(gg, goMF)
summary(sgg)
#> IGRAPH 2030fb9 UN-- 2275 9190 -- 
#> + attr: name (v/c), GeneName (v/c), TopOnto_OVG (v/c),
#> | TopOnto_OVG_HDO_ID (v/c), GO_MF_ID (v/c), GO_MF (v/c)

sfile<-system.file("extdata", "flatfile.go.CC.csv", package = "BioNAR")
goCC <- read.table(sfile,sep="\t",skip=1,header=FALSE,strip.white=TRUE,quote="")
sgg <- annotateGoCC(gg, goCC)
summary(sgg)
#> IGRAPH 2030fb9 UN-- 2275 9190 -- 
#> + attr: name (v/c), GeneName (v/c), TopOnto_OVG (v/c),
#> | TopOnto_OVG_HDO_ID (v/c), GO_CC_ID (v/c), GO_CC (v/c)

Estimate vertex centrality measures

Estimate centrality measures with values added as vertex attributes.

BioNAR supports centrality measures as following: * DEG - degree, * BET - betweenness, * CC - clustering coefficient, * SL - semilocal centrality, * mnSP - mean shortest path, * PR - page rank, * sdSP - standard deviation of the shortest path. These are saved as vertex atrtributes.

gg <- calcCentrality(gg)
summary(gg)
#> IGRAPH 2030fb9 UN-- 2275 9190 -- 
#> + attr: name (v/c), GeneName (v/c), TopOnto_OVG (v/c),
#> | TopOnto_OVG_HDO_ID (v/c), DEG (v/n), BET (v/n), CC (v/n), SL (v/n),
#> | mnSP (v/n), PR (v/n), sdSP (v/c)

Get vertex centralities as a matrix.

Instead of saving entrality centrality values on the graph, e.g. to provide different names for the vertex centrality attributes, they can be obtained in a matrix form:

mc <- getCentralityMatrix(gg)
head(mc)
#>      ID DEG        BET         CC     SL  mnSP           PR  sdSP
#> 1   273  14  1437.7544 0.08791209  77701 3.343 0.0006682096 0.666
#> 2  6455  24  6933.8371 0.07608696 283341 2.948 0.0010598575 0.666
#> 3  1759  16  1879.9038 0.21666667 235112 3.060 0.0007062358 0.667
#> 4  1785  32 16999.3049 0.06250000 358555 2.869 0.0014391585 0.678
#> 5 26052   7   624.3073 0.14285714  90813 3.438 0.0003420742 0.715
#> 6  1837   6   672.2350 0.13333333 176022 3.153 0.0003303946 0.706

Get the centrality measures for random graph

Sometimes one needs to compare the graph properties of the the properties of an the observed network to randomised networks of a similar size. The BioNAR command below provides three ways of generating randomization, randomised networks given an observed network including: G(n,p) Erdos-Renyi model, Barabasi-Albert model and new random graph from a given graph by randomly adding/removing edges.

{r}
ggrm <- getRandomGraphCentrality(sgg, type = c("cgnp"))
head(ggrm)

Power law fit

To examine a network’s underlying structure (i.e. not random), one can test a network’s degree distribution for evidence of scale-free structure and compare this against randomised network models. For this we used the R “PoweRlaw” package (version 0.50.0) (Gillespie, 2015). For the case study, i.e.  our presynaptic PPI network, we found evidence for disassortative mixing (Newman, 2002), i.e. a preference for high-degree genes to attach to low-degree gene(presynaptic: -0.16).

pFit <- fitDegree( as.vector(igraph::degree(graph=gg)),threads=1,Nsim=5,
                    plot=TRUE,WIDTH=2480, HEIGHT=2480)

Get entropy rate

Evidence for scale-free structure can also be tested by performing a perturbation analysis of each of the network’s vertices. In this analysis each protein is being perturbed through over-expression (red) and under-expression (green), with the global graph entropy rate (SR) after each proteins perturbation being plotted against the log of the proteins degree, as shown at the plot below. In our case study of the presynaptic PPI network we observe a bi-modal response, between gene over-expression and degree, and opposing bi-phasic response relative to over/under-expression between global entropy rate and degree. This type of bi-modal, bi-phasic behaviour has been observed only in networks with scale-free or approximate scale-free topology (Teschendorff et al, 2014).

ent <- getEntropyRate(gg)
ent
#> $maxSr
#> [1] 4.064787
#> 
#> $SRo
#> [1] 2.755673
SRprime <- getEntropy(gg, maxSr = NULL)
head(SRprime)
#>   ENTREZ.ID GENE.NAME DEGREE        UP      DOWN
#> 1       273      AMPH     14 0.6758968 0.6779169
#> 2      6455    SH3GL1     24 0.6764513 0.6776185
#> 3      1759      DNM1     16 0.6770435 0.6776884
#> 4      1785      DNM2     32 0.6772362 0.6773413
#> 5     26052      DNM3      7 0.6770322 0.6778923
#> 6      1837      DTNA      6 0.6769965 0.6779295
plotEntropy(SRprime, subTIT = "Entropy", SRo = ent$SRo, maxSr = ent$maxSr)

Get modularity. Normalised modularity.

Normalised modularity (Qm) allows the comparison of networks with varying structure. Qm based on the previous studies by Parter et al., 2007, Takemoto, 2012, Takemoto, 2013, Takemoto and Borjigin, 2011, which was defined as:

\[Qm = \frac{Q_{real}-Q_{rand}}{Q_{max}-Q_{rand}}\]

Where \(Q_{real}\) is the network modularity of a real-world signaling network and, \(Q_{rand}\) is the average network modularity value obtained from 10,000 randomized networks constructed from its real-world network. \(Q_{max}\) was estimated as: \[Q_{max}=1 − \frac{1}{M}\], where \(M\) is the number of modules in the real network.

Randomized networks were generated from a real-world network using the edge-rewiring algorithm (Maslov and Sneppen, 2002).

nm<-normModularity(gg,alg='louvain')
nm
#> [1] 0.01920175

Clustering

Clustering, or community detection, in networks has been well studied in the field of statistical physics with particular attention to methods developed for social science networks. The underlying assumption(s) of what makes a community in social science, translates remarkably well to what we think of as a community (sub-complex, module or cluster) in PPI networks. The possible algorithms of choice implemented in BioNAR are: * “lec”(‘Leading-Eigenvector, Newman, 2006), * “wt”(Walktrap, Pons & Latapy, 2006), * “fc”(Fast-Greedy Community’ algorithm, Clauset et al., 2004), * “infomap” (InfoMAP, Rosvall et al., 2007; Rosvall et al., 2010), * “louvain” (Louvain, Blondel et al., 2008), * “sgG1”, “sgG2”, “sgG5”(SpinGlass, Reichardt & Bornholdt). For each algorithm of interest the community membership can be obtained with'calcMembership command. All algorithm implementations, apart from Spectral were performed using the publicly available R package igraph (Csardi & Nepusz, 2006) (R version 3.4.2, igraph version 1.1.2). Parameters used in the fc, lec, sg, wt and lourvain algorithms were chosen as to maximise the measure Modularity (Newman & Girvan, 2004); infomap seeks the optimal community structure in the data by maximising the objective function called the Minimum Description Length (Rissanen, 1978; Grwald et al., 2005)

# choose one algorithm from the list
alg = "louvain"
mem <- calcMembership(gg, alg)
pander(head(mem))
names membership
273 1
6455 1
1759 1
1785 1
26052 1
1837 2

Due to internal random initialisation consecutive invocation of the same algorithm could produce slightly different community structures:

mem2 <- calcMembership(gg, alg)
idx<-match(mem$names,mem2$names)
idnx<-which(mem$membership!=mem2$membership[idx])
pander(head(cbind(mem[idnx,],mem2[idx[idnx],])))
  names membership names membership
11 22808 5 22808 2
12 58513 1 58513 3
15 29993 4 29993 1
17 22895 6 22895 2
18 5094 4 5094 5
20 3303 7 3303 6

To avoid inconsistency in downstream analysis we provide two additional functions calcClustering and calcAllClustering that use calcMembership to calculate community memberships and store them within the graph vertices attributes named after the algorithm. They also calculate modularity values and store them as graph vertex attributes named after the clustering algorithm. The difference between calcClustering and calcAllClustering is that calcAllClusteringallows to calculate memberships for all clustering algorithms simultaneously (may take time), and store them as graph vertices attributes, while calcClusteringcommand will work for a specific algorithm.

gg <- calcClustering(gg, alg)
summary(gg)
#> IGRAPH 2030fb9 UN-- 2275 9190 -- 
#> + attr: louvain (g/n), name (v/c), GeneName (v/c), TopOnto_OVG (v/c),
#> | TopOnto_OVG_HDO_ID (v/c), DEG (v/n), BET (v/n), CC (v/n), SL (v/n),
#> | mnSP (v/n), PR (v/n), sdSP (v/c), louvain (v/c)

Comminity membership data could be obtained from the graph vertex attribute:

mem.df<-data.frame(names=V(gg)$name,membership=as.numeric(V(gg)$louvain))

To compare different clustering algorithms,a summary matrix can be calculated with the following properties: 1. maximum Modularity obtained (mod), 2. number of detected communities (C), 3. the number of communities with size (Cn1) equal to 1, 4. the number of communities >= 100 (Cn100), 5. the fraction of edges lying between communities (mu), 6. the size of the smallest community (Min. C), 7. the size of the largest community (Max. C), 8. the average ( Mean C), median (Median C), 9. first quartile (1st Qu. C), and 10. third quartile (3rd Qu. C) of the community size.

ggc <- calcAllClustering(gg)
m<-clusteringSummary(ggc,att=c('lec','wt','fc',
                                'infomap','louvain',
                                'sgG1','sgG2','sgG5'))
pander(m)
Table continues below
  mod C Cn1 Cn100 mu Min. C 1st Qu. C
lec 0.2029 2 0 2 0.2163 552 818.2
wt 0.3602 231 116 4 0.4874 1 1
fc 0.4382 26 0 5 0.3682 3 3.25
infomap 0.4179 159 0 2 0.5449 2 5
louvain 0.4656 15 0 9 0.4313 22 77.5
sgG1 0.4755 34 0 7 0.4017 2 4
sgG2 0.4596 44 0 6 0.4948 3 16.25
sgG5 0.392 97 0 0 0.5909 2 13
  Median C Mean C 3rd Qu. C Max. C
lec 1084 1084 1351 1617
wt 1 9.39 3 594
fc 8 83.42 31 475
infomap 9 13.64 13 181
louvain 152 144.6 205.5 290
sgG1 7.5 63.79 73.75 361
sgG2 44 49.3 69.25 163
sgG5 20 22.36 28 82

It is often useful to be able to visualize clusters of the graph. The simplest way to do this is to color each cluster uniquely and plot the graph:

palette <- distinctColorPalette(max(as.numeric(mem.df$membership)))
plot(gg,vertex.size=3,layout=layout_nicely,
        vertex.label=NA,
        vertex.color=palette[as.numeric(mem.df$membership)],
        edge.color='grey95')
legend('topright',legend=names(table(mem.df$membership)),
        col=palette,pch=19,ncol = 2)

On the plot we can see some distinctive clusters but most verices are indistinguishable within the central part of the plot. So we could layout graph clusterwise:

lay<-layoutByCluster(gg,mem.df,layout = layout_nicely)
plot(gg,vertex.size=3,layout=lay,
        vertex.label=NA,
        vertex.color=palette[as.numeric(mem.df$membership)],
        edge.color='grey95')
legend('topright',legend=names(table(mem.df$membership)),
        col=palette,pch=19,ncol = 2)

It is also possible to visualize the interaction between communities:

idx<-match(V(gg)$name,mem.df$names)
cgg<-getCommunityGraph(gg,mem.df$membership[idx])
D0 = unname(degree(cgg))
plot(cgg, vertex.size=sqrt(V(cgg)$size), vertex.cex = 0.8,
vertex.color=round(log(D0))+1,layout=layout_with_kk,margin=0)

Reclustering

Reclustering a clustered graph using the same, or different, clustering algorithm:

remem<-calcReclusterMatrix(gg,mem.df,alg,10)
head(remem)
#>   names membership recluster
#> 1   273          1         1
#> 2  6455          1         1
#> 3  1759          1         1
#> 4  1785          1         1
#> 5 26052          2         8
#> 6  1837          3        18

And we can apply second order clustering layout:

lay<-layoutByRecluster(gg,remem,layout_nicely)
plot(gg,vertex.size=3,layout=lay,
        vertex.label=NA,
        vertex.color=palette[as.numeric(mem.df$membership)],
        edge.color='grey95')
legend('topright',legend=names(table(mem.df$membership)),
        col=palette,pch=19,ncol = 2)

Consensus matrix

To assess the robustness of obtained clusters, a randomization study can be performed, which applies the same clustering algorithm to N perturbed networks. The clustering results are returned as a consensus matrix where each matrix elements is assigned the frequency with which a pair of nodes vertices is found in the same cluster.

Where ‘alg’ gives the name of the clustering algorithm, ‘type’ the sampling scheme (1 sample edges, and 2 sample verices) used, ‘mask’ the percentage of edges or vertices to mask, and ‘reclust’ whether reclustering should be performed on the community set found, ‘Cnmin’ minimum cluster size and ‘Cnmax’ the maximum cluster size above which reclustering will be preformed (if reClust=TRUE).

#Build consensus matrix for louvain clustering
conmat <- makeConsensusMatrix(gg, N=5,
                                alg = alg, type = 2, 
                                mask = 10,reclust = FALSE, 
                                Cnmax = 10)

For the sake of timing we use only five randomisation rounds, for the real analysis you should use at least 500.

##Consensus matrix value distribution

Consensus matrix values can be visualised in the following way:

steps <- 100
Fn  <- ecdf(conmat[lower.tri(conmat)])
X<-seq(0,1,length.out=steps+1)
cdf<-Fn(X)
dt<-data.frame(cons=X,cdf=cdf)
ggplot(dt,aes(x=cons,y=cdf))+geom_line()+
        theme(            
        axis.title.x=element_text(face="bold",size=rel(2.5)),
        axis.title.y=element_text(face="bold",size=rel(2.5)),
        legend.title=element_text(face="bold",size=rel(1.5)),
        legend.text=element_text(face="bold",size=rel(1.5)),
        legend.key=element_blank())+
    theme(panel.grid.major = element_line(colour="grey40",size=0.2),
            panel.grid.minor = element_line(colour="grey40",size=0.1),
            panel.background = element_rect(fill="white"),
            panel.border = element_rect(linetype="solid",fill=NA))

Cluster robustness

Cluster robustness assesses the robustness of obtained clusters and can help evaluate the “goodness” of a chosen clustering algorithm.

clrob<-getRobustness(gg, alg = alg, conmat)
pander(clrob)
alg C Cn Crob CrobScaled
louvain 1 270 0.09614 0.2368
louvain 2 159 0.09543 0.2218
louvain 3 360 0.0862 0.02562
louvain 4 80 0.09042 0.1153
louvain 5 87 0.08965 0.099
louvain 6 119 0.08499 0
louvain 7 190 0.08913 0.088
louvain 8 163 0.08882 0.08142
louvain 9 127 0.08583 0.0179
louvain 10 188 0.08656 0.03328
louvain 11 72 0.09041 0.1151
louvain 12 205 0.0898 0.1022
louvain 13 162 0.08969 0.09978
louvain 14 74 0.09704 0.2561
louvain 15 19 0.1321 1

Bridgeness

Bridging proteins are known to interact with many neighbours simultaneously, organise function inside the communities they belong to, but also affect/influence other communities in the network (Nepusz et al., 2008). Bridgeness can be estimated from the consensus clustering matrix estimated above and vertex degree to calculate the vertex’s community membership, i.e.  the probability a specific vertex belongs to every community obtained by a given clustering algorithm.

The Bridgeness measure lies between 0 - implying a vertex clearly belongs in a single community, and 1 - implying a vertex forms a ‘global bridge’ across every community with the same strength.

br<-getBridgeness(gg,alg = alg,conmat)
pander(head(br))
ID GENE.NAME BRIDGENESS.louvain
273 AMPH 0.2312
6455 SH3GL1 0.5261
1759 DNM1 0
1785 DNM2 0.2995
26052 DNM3 0.2314
1837 DTNA 0.166
gg<-calcBridgeness(gg,alg = alg,conmat)
vertex_attr_names(gg)
#>  [1] "name"               "GeneName"           "TopOnto_OVG"       
#>  [4] "TopOnto_OVG_HDO_ID" "DEG"                "BET"               
#>  [7] "CC"                 "SL"                 "mnSP"              
#> [10] "PR"                 "sdSP"               "louvain"           
#> [13] "BRIDGENESS.louvain"

Bridgeness plot

Semi-local centrality measure (Chen et al., 2011) also lies between 0 and 1 indicating whether protein is important globally or locally. By plotting Bridgeness against semi-local centrality we can categorises the influence each protein found in our network has on the overall network structure: * Region 1, proteins having a ‘global’ rather than ‘local’ influence in the network (also been called bottle-neck bridges, connector or kinless hubs (0<Sl<0.5; 0.5<Br<1). * Region 2, proteins having ‘global’ and ‘local’ influence (0.5<Sl<1, 0.5<Br<1). * Region 3, proteins centred within the community they belong to, but also communicating with a few other specific communities (0<Sl<0.5; 0.1<Br<0.5). * Region 4, proteins with ‘local’ impact , primarily within one or two communities (local or party hubs, 0.5<Sl<1, 0<Br<0.5).

g<-plotBridgeness(gg,alg = alg,
               VIPs=c('8495','22999','8927','8573',
                      '26059','8497','27445','8499'),
               Xatt='SL',
               Xlab = "Semilocal Centrality (SL)",
               Ylab = "Bridgeness (B)",
               bsize = 3,
               spsize =7,
               MainDivSize = 0.8,
               xmin = 0,
               xmax = 1,
               ymin = 0,
               ymax = 1,
               baseColor="royalblue2",
               SPColor="royalblue2")
g

#Interactive view of bridgeness plot

library(plotly)
g<-plotBridgeness(gg,alg = alg,
               VIPs=c('8495','22999','8927','8573',
                      '26059','8497','27445','8499'),
               Xatt='SL',
               Xlab = "Semilocal Centrality (SL)",
               Ylab = "Bridgeness (B)",
               bsize = 1,
               spsize =2,
               MainDivSize = 0.8,
               xmin = 0,
               xmax = 1,
               ymin = 0,
               ymax = 1,
               baseColor="royalblue2",
               SPColor="royalblue2")
ggplotly(g)

Disease/annotation pairs

Given that Disease associated genes are connected within the graph, the common question is to check whether the networks for two different diseases are overlapping, which may indicate the common molecular mechanisms. Same is valid for any pair of annotations, e.g. one would ask if two different biological functions are related.

To address this question, we utilise the algorithm from Menche et al, which estimates the minimal shortest paths between two distinct annotations and compares this to the randomly annotated graph.

Below example shows the estimation of disease separation for the following diseases: DOID:10652 (Alzheimer’s_disease), (bipolar_disorder), DOID:12849 (autistic_disorder), DOID:1826 (epilepsy). Command calcDiseasePairs quickly estimates the two annotation separation on the graph and compares it with one randomly reannotated graph. This could be used for initial guess of the relationships between the annotations.

To assess the significance of the obtained separation values the command runPermDisease should be used, where the user can specify the number of permutations. TExecuting this command, which will take time depending of the number of permutations, returns table with of p-values, adjusted p-values, q-values and Bonferroni test for each disease-disease pair.

p <- calcDiseasePairs(
    gg,
    name = "TopOnto_OVG_HDO_ID",
    diseases = c("DOID:10652","DOID:3312"),
    permute = "r"
)
pander(p$disease_separation)
  DOID:10652 DOID:3312
DOID:10652 0 -0.04423
DOID:3312 NA 0

r <- runPermDisease(
    gg,
    name = "TopOnto_OVG_HDO_ID",
    diseases = c("DOID:10652","DOID:3312"),
    Nperm = 100,
    alpha = c(0.05, 0.01, 0.001)
)

pander(r$Disease_overlap_sig)
Table continues below
HDO.ID.A N.A HDO.ID.B N.B sAB Separated Overlap zScore
DOID:10652 310 DOID:10652 310 0 . . 0
DOID:10652 310 DOID:3312 153 -0.2763 . YES -3.086
DOID:3312 153 DOID:3312 153 0 . . 0
pvalue Separation/Overlap.than.chance Bonferroni p.adjusted q-value
1 larger . 1 1
0.002025 Smaller ** 0.01114 0.006076
1 larger . 1 1

Cluster overrepresentation

To identify the clusters with overrepresented function or disease we introduced the function which calculates the overrepesentation (enrichment for specified annotation). Based on R package fgsea.

ora <- clusterORA(gg, alg, name = 'TopOnto_OVG_HDO_ID', vid = "name",alpha = 0.1, col = COLLAPSE)

Session Info

Platform

Packages

ondiskversion loadedversion date source
AnnotationDbi 1.62.2 1.62.2 2023-08-10 Bioconductor
AnnotationHub 3.8.0 3.8.0 2023-08-10 Bioconductor
apcluster 1.4.10 1.4.10 2022-05-31 CRAN (R 4.3.1)
backports 1.4.1 1.4.1 2021-12-13 CRAN (R 4.3.1)
base64enc 0.1.3 0.1-3 2015-07-28 CRAN (R 4.3.1)
Biobase 2.60.0 2.60.0 2023-08-10 Bioconductor
BiocFileCache 2.8.0 2.8.0 2023-08-10 Bioconductor
BiocGenerics 0.46.0 0.46.0 2023-08-10 Bioconductor
BiocManager 1.30.22 1.30.22 2023-08-08 CRAN (R 4.3.1)
BiocParallel 1.34.2 1.34.2 2023-08-10 Bioconductor
BiocVersion 3.17.1 3.17.1 2023-08-10 Bioconductor
BioNAR 1.2.5 1.2.5 2023-08-10 Bioconductor
Biostrings 2.68.1 2.68.1 2023-08-10 Bioconductor
bit 4.0.5 4.0.5 2022-11-15 CRAN (R 4.3.1)
bit64 4.0.5 4.0.5 2020-08-30 CRAN (R 4.3.1)
bitops 1.0.7 1.0-7 2021-04-24 CRAN (R 4.3.1)
blob 1.2.4 1.2.4 2023-03-17 CRAN (R 4.3.1)
bslib 0.5.0 0.5.0 2023-06-09 CRAN (R 4.3.1)
cachem 1.0.8 1.0.8 2023-05-01 CRAN (R 4.3.1)
callr 3.7.3 3.7.3 2022-11-02 CRAN (R 4.3.1)
checkmate 2.2.0 2.2.0 2023-04-27 CRAN (R 4.3.1)
cli 3.6.1 3.6.1 2023-03-23 CRAN (R 4.3.1)
cluster 2.1.4 2.1.4 2022-08-22 CRAN (R 4.3.1)
clusterCons 1.2 1.2 2022-02-22 CRAN (R 4.3.1)
codetools 0.2.19 0.2-19 2023-02-01 CRAN (R 4.3.1)
colorspace 2.1.0 2.1-0 2023-01-23 CRAN (R 4.3.1)
cowplot 1.1.1 1.1.1 2020-12-30 CRAN (R 4.3.1)
crayon 1.5.2 1.5.2 2022-09-29 CRAN (R 4.3.1)
crosstalk 1.2.0 1.2.0 2021-11-04 CRAN (R 4.3.1)
curl 5.0.1 5.0.1 2023-06-07 CRAN (R 4.3.1)
data.table 1.14.8 1.14.8 2023-02-17 CRAN (R 4.3.1)
DBI 1.1.3 1.1.3 2022-06-18 CRAN (R 4.3.1)
dbplyr 2.3.3 2.3.3 2023-07-07 CRAN (R 4.3.1)
devtools 2.4.5 2.4.5 2022-10-11 CRAN (R 4.3.1)
digest 0.6.33 0.6.33 2023-07-07 CRAN (R 4.3.1)
doParallel 1.0.17 1.0.17 2022-02-07 CRAN (R 4.3.1)
dplyr 1.1.2 1.1.2 2023-04-20 CRAN (R 4.3.1)
dynamicTreeCut 1.63.1 1.63-1 2016-03-11 CRAN (R 4.3.1)
ellipsis 0.3.2 0.3.2 2021-04-29 CRAN (R 4.3.1)
evaluate 0.21 0.21 2023-05-05 CRAN (R 4.3.1)
fansi 1.0.4 1.0.4 2023-01-22 CRAN (R 4.3.1)
farver 2.1.1 2.1.1 2022-07-06 CRAN (R 4.3.1)
fastcluster 1.2.3 1.2.3 2021-05-24 CRAN (R 4.3.1)
fastmap 1.1.1 1.1.1 2023-02-24 CRAN (R 4.3.1)
fastmatch 1.1.3 1.1-3 2021-07-23 CRAN (R 4.3.1)
fgsea 1.26.0 1.26.0 2023-08-10 Bioconductor
filelock 1.0.2 1.0.2 2018-10-05 CRAN (R 4.3.1)
foreach 1.5.2 1.5.2 2022-02-02 CRAN (R 4.3.1)
foreign 0.8.84 0.8-84 2022-12-06 CRAN (R 4.3.1)
Formula 1.2.5 1.2-5 2023-02-24 CRAN (R 4.3.1)
fs 1.6.3 1.6.3 2023-07-20 CRAN (R 4.3.1)
generics 0.1.3 0.1.3 2022-07-05 CRAN (R 4.3.1)
GenomeInfoDb 1.36.1 1.36.1 2023-08-10 Bioconductor
GenomeInfoDbData 1.2.10 1.2.10 2023-06-20 Bioconductor
ggplot2 3.4.2 3.4.2 2023-04-03 CRAN (R 4.3.1)
ggrepel 0.9.3 0.9.3 2023-02-03 CRAN (R 4.3.1)
glue 1.6.2 1.6.2 2022-02-24 CRAN (R 4.3.1)
GO.db 3.17.0 3.17.0 2023-06-20 Bioconductor
graph 1.78.0 1.78.0 2023-08-10 Bioconductor
gridExtra 2.3 2.3 2017-09-09 CRAN (R 4.3.1)
gtable 0.3.3 0.3.3 2023-03-21 CRAN (R 4.3.1)
highr 0.10 0.10 2022-12-22 CRAN (R 4.3.1)
Hmisc 5.1.0 5.1-0 2023-05-08 CRAN (R 4.3.1)
htmlTable 2.4.1 2.4.1 2022-07-07 CRAN (R 4.3.1)
htmltools 0.5.5 0.5.5 2023-03-23 CRAN (R 4.3.1)
htmlwidgets 1.6.2 1.6.2 2023-03-17 CRAN (R 4.3.1)
httpuv 1.6.11 1.6.11 2023-05-11 CRAN (R 4.3.1)
httr 1.4.6 1.4.6 2023-05-08 CRAN (R 4.3.1)
igraph 1.5.1 1.5.1 2023-08-10 CRAN (R 4.3.1)
impute 1.74.1 1.74.1 2023-08-10 Bioconductor
interactiveDisplayBase 1.38.0 1.38.0 2023-08-10 Bioconductor
IRanges 2.34.1 2.34.1 2023-08-10 Bioconductor
iterators 1.0.14 1.0.14 2022-02-05 CRAN (R 4.3.1)
jquerylib 0.1.4 0.1.4 2021-04-26 CRAN (R 4.3.1)
jsonlite 1.8.7 1.8.7 2023-06-29 CRAN (R 4.3.1)
KEGGREST 1.40.0 1.40.0 2023-08-10 Bioconductor
knitr 1.43 1.43 2023-05-25 CRAN (R 4.3.1)
labeling 0.4.2 0.4.2 2020-10-20 CRAN (R 4.3.1)
later 1.3.1 1.3.1 2023-05-02 CRAN (R 4.3.1)
latex2exp 0.9.6 0.9.6 2022-11-28 CRAN (R 4.3.1)
lattice 0.21.8 0.21-8 2023-04-05 CRAN (R 4.3.1)
lazyeval 0.2.2 0.2.2 2019-03-15 CRAN (R 4.3.1)
lifecycle 1.0.3 1.0.3 2022-10-07 CRAN (R 4.3.1)
magrittr 2.0.3 2.0.3 2022-03-30 CRAN (R 4.3.1)
Matrix 1.6.0 1.6-0 2023-07-08 CRAN (R 4.3.1)
matrixStats 1.0.0 1.0.0 2023-06-02 CRAN (R 4.3.1)
memoise 2.0.1 2.0.1 2021-11-26 CRAN (R 4.3.1)
mime 0.12 0.12 2021-09-28 CRAN (R 4.3.1)
miniUI 0.1.1.1 0.1.1.1 2018-05-18 CRAN (R 4.3.1)
minpack.lm 1.2.3 1.2-3 2023-01-26 CRAN (R 4.3.1)
munsell 0.5.0 0.5.0 2018-06-12 CRAN (R 4.3.1)
nnet 7.3.19 7.3-19 2023-05-03 CRAN (R 4.3.1)
org.Hs.eg.db 3.17.0 3.17.0 2023-06-20 Bioconductor
pander 0.6.5 0.6.5 2022-03-18 CRAN (R 4.3.1)
pillar 1.9.0 1.9.0 2023-03-22 CRAN (R 4.3.1)
pkgbuild 1.4.2 1.4.2 2023-06-26 CRAN (R 4.3.1)
pkgconfig 2.0.3 2.0.3 2019-09-22 CRAN (R 4.3.1)
pkgload 1.3.2.1 1.3.2.1 2023-07-08 CRAN (R 4.3.1)
plotly 4.10.2 4.10.2 2023-06-03 CRAN (R 4.3.1)
png 0.1.8 0.1-8 2022-11-29 CRAN (R 4.3.1)
poweRlaw 0.70.6 0.70.6 2020-04-25 CRAN (R 4.3.1)
pracma 2.4.2 2.4.2 2022-09-22 CRAN (R 4.3.1)
preprocessCore 1.62.1 1.62.1 2023-08-10 Bioconductor
prettyunits 1.1.1 1.1.1 2020-01-24 CRAN (R 4.3.1)
processx 3.8.2 3.8.2 2023-06-30 CRAN (R 4.3.1)
profvis 0.3.8 0.3.8 2023-05-02 CRAN (R 4.3.1)
promises 1.2.1 1.2.1 2023-08-10 CRAN (R 4.3.1)
ps 1.7.5 1.7.5 2023-04-18 CRAN (R 4.3.1)
purrr 1.0.2 1.0.2 2023-08-10 CRAN (R 4.3.1)
R6 2.5.1 2.5.1 2021-08-19 CRAN (R 4.3.1)
randomcoloR 1.1.0.1 1.1.0.1 2019-11-24 CRAN (R 4.3.1)
rappdirs 0.3.3 0.3.3 2021-01-31 CRAN (R 4.3.1)
rbibutils 2.2.14 2.2.14 2023-08-07 CRAN (R 4.3.1)
RColorBrewer 1.1.3 1.1-3 2022-04-03 CRAN (R 4.3.1)
Rcpp 1.0.11 1.0.11 2023-07-06 CRAN (R 4.3.1)
RCurl 1.98.1.12 1.98-1.12 2023-03-27 CRAN (R 4.3.1)
Rdpack 2.4 2.4 2022-07-20 CRAN (R 4.3.1)
remotes 2.4.2.1 2.4.2.1 2023-07-18 CRAN (R 4.3.1)
rlang 1.1.1 1.1.1 2023-04-28 CRAN (R 4.3.1)
rmarkdown 2.23 2.23 2023-07-01 CRAN (R 4.3.1)
rpart 4.1.19 4.1.19 2022-10-21 CRAN (R 4.3.1)
RSpectra 0.16.1 0.16-1 2022-04-24 CRAN (R 4.3.1)
rSpectral 1.0.0.10 1.0.0.10 2023-01-18 CRAN (R 4.3.1)
RSQLite 2.3.1 2.3.1 2023-04-03 CRAN (R 4.3.1)
rstudioapi 0.15.0 0.15.0 2023-07-07 CRAN (R 4.3.1)
Rtsne 0.16 0.16 2022-04-17 CRAN (R 4.3.1)
S4Vectors 0.38.1 0.38.1 2023-08-10 Bioconductor
sass 0.4.7 0.4.7 2023-07-15 CRAN (R 4.3.1)
scales 1.2.1 1.2.1 2022-08-20 CRAN (R 4.3.1)
sessioninfo 1.2.2 1.2.2 2021-12-06 CRAN (R 4.3.1)
shiny 1.7.4.1 1.7.4.1 2023-07-06 CRAN (R 4.3.1)
stringi 1.7.12 1.7.12 2023-01-11 CRAN (R 4.3.1)
stringr 1.5.0 1.5.0 2022-12-02 CRAN (R 4.3.1)
survival 3.5.5 3.5-5 2023-03-12 CRAN (R 4.3.1)
synaptome.data 0.99.6 0.99.6 2023-08-09 Bioconductor
synaptome.db 0.99.12 0.99.12 2023-08-09 Bioconductor
tibble 3.2.1 3.2.1 2023-03-20 CRAN (R 4.3.1)
tidyr 1.3.0 1.3.0 2023-01-24 CRAN (R 4.3.1)
tidyselect 1.2.0 1.2.0 2022-10-10 CRAN (R 4.3.1)
urlchecker 1.0.1 1.0.1 2021-11-30 CRAN (R 4.3.1)
usethis 2.2.2 2.2.2 2023-07-06 CRAN (R 4.3.1)
utf8 1.2.3 1.2.3 2023-01-31 CRAN (R 4.3.1)
V8 4.3.3 4.3.3 2023-07-18 CRAN (R 4.3.1)
vctrs 0.6.3 0.6.3 2023-06-14 CRAN (R 4.3.1)
viridis 0.6.4 0.6.4 2023-07-22 CRAN (R 4.3.1)
viridisLite 0.4.2 0.4.2 2023-05-02 CRAN (R 4.3.1)
WGCNA 1.72.1 1.72-1 2023-01-18 CRAN (R 4.3.1)
withr 2.5.0 2.5.0 2022-03-03 CRAN (R 4.3.1)
xfun 0.40 0.40 2023-08-09 CRAN (R 4.3.1)
xtable 1.8.4 1.8-4 2019-04-21 CRAN (R 4.3.1)
XVector 0.40.0 0.40.0 2023-08-10 Bioconductor
yaml 2.3.7 2.3.7 2023-01-23 CRAN (R 4.3.1)
zlibbioc 1.46.0 1.46.0 2023-08-10 Bioconductor