example code for PPI presynaptic network

[cmclean5]

reset

rm(list=ls())

library(igraph)

install.packages("pastecs")

library(pastecs)

read-in graph

gg = read.graph("PPI_Presynaptic.gml",format="gml")

No of nodes

N = length(V(gg))

Adjacency Matrix

adj = get.adjacency(gg)

Degree Matrix

D = diag(degree(gg))

Laplacian

Measures information flow through graph

L = D - as.matrix(adj)

Eigenvalues of L

The eigenvalues of the Laplacian are related to the spectrum of the graph, providing information about the graph's connectivity.

In the context of random walks, the eigenvalues of the Laplacian influence the rates of convergence of the random walk to the stationary

distribution.

Random walks on graphs represented by the normalized Laplacian tend to have more uniform behavior across nodes, reducing the bias

introduced by high-degree nodes. This can result in improved convergence properties and mixing times.

Eigen = eigen(L)

Entropy Measure given t=tau for graph Laplacian's eigenvalues

S.tau <- function(x, t, n){ mu = exp(-1xt) mu = mu/sum(mu) (-1/log(n))sum(ifelse(mu>0, mulog(mu),0)) }

tau sequence

tau = seq(0,50,10^-2)

Entropy

S = sapply(tau, function(j) S.tau(x=Eigen$values, t=j, n=N))

Specific Heat

C = -diff(S)/diff(log(tau)) tau1st = tau[-1]

1st derivative

plot(xy.coords(x=tau1st,y=C))

Find turning points of C

tp1st = turnpoints(C) if( tp1st$nturns > 0 ){ Ctau = tau1st[tp1st$tppos] }

The rate in change of Specific Heat curve??

C2 = diff(C) tau2nd = tau1st[-1]

2nd derivative

plot(xy.coords(x=tau2nd,y=C2))

tp2nd = turnpoints(C2)

if( tp2nd$nturns > 0 ){ C2tau = tau2nd[tp2nd$tppos] }

Starting Point, i.e. all nodes disconnected

S0=as.matrix(expm(-0*L)) S0=S0/sum(diag(S0))

C: The communicability matrix

C is used to quantify the efficiency of information flow between pairs of nodes in a graph.

It is derived from the exponential of the adjacency or Laplacian matrix and is related to random walks.

Approximation to C given by: C = exp(-tauL) ~ 1 - tauL

rho ~ C / Tr(C)

So Approx for St:

Eye = diag(rep(1,N))

St = Eye - lambdaC * L

St = Stationary/sum(diag(St))

Calculate rho,

St=as.matrix(expm(-Ctau[1]*L))

St=as.matrix(expm(-C2tau[1]*L)) St=St/sum(diag(St))

test... not sure what this implies yet?

rho = St - S0

rho = St

Nodes should be connected/grouped together if information flow between them as large as the minimum information at nodes

rho2=matrix(0,N,N) colnames(rho2) = V(gg)$name rownames(rho2) = V(gg)$name for( i in 1:N ){ ii = rho[i,i] for( j in 1:N ){ jj = rho[j,j] tmp = rho[i,j]/min(ii,jj) if( tmp -1 >= 0 ){ rho2[i,j]=1 } }}

build network from rho2... remember to add node names to columns/row of rho2 first

cc = graph_from_adjacency_matrix(rho2, mode=ifelse(isSymmetric(rho2),"undirected","directed"))

node groups/clusters/complexes/communities

complexes = igraph::decompose(cc)

function clusters gives a bit more information than decompose

clusters = igraph::clusters(cc)