Analisis de redes biolóligicas

Introducción

Hasta ahora lo que hemos visto es el comportamiento de un gen en especifico, miendo su actividad como una relación directa a la cantidad de expresión relativa de este (qRT-PCR). Despues hemos visto como a partir de una secuenciación masiva como los microarreglos podemos determinar los genes que estan siendo reprimidos o sobre expresados. Finalmente usamos los análisis de enriquecimiento para dejar de un lado esta vista reduccionista del sistema (por genes individuales) y de una forma más holística vimos como el análisis de enriquecimiento nos podía ayudar para entender estos cambios de expresión traducidos en procesos biológicos más complejos. Por otro lado, el análisis de redes regulatorias trata de atacar este problema desde un punto de vista o estrategia totalmente distinta. En primer lugar, debemos entender para este capítulo el objetivo de nuestro estudio, que son las redes regulatorias. Entenderemos a los organismos como conjuntos de redes regulatorias, que están en constante cambio y regulación (inter e intra-red). Uno de los puntos claves que veremos en este capítulo es que, a veces, un grupo de genes menos representados en los analisis de enriquecimiento, puede jugar un papel trascendental en la organización o la integridad de esta red (los llamados broker). Lo segundo que debemos entender, es que esta joven rama de la biología está aún desarrollando las herramientas para estos análisis. Una de estas herramientas es la teoría de grafos que la toma de la matemática y la trata de ajustar a un entorno biológico que nos permita caracterizar de forma más complejas o integrativa dichas redes biológicas. La transferencia aún no es perfecta pero es un tema en vertiginoso ascenso. Actualemente, existen varias herramientas que nos permiten trabajar con estas redes en R. Una de ellas es el paquete igraph. Que es una librería tambien disponible para el entorno de python (como varias de las librerías actualmente).

Elementos de los grafos

Nodos o vertices

Estan vas a ser los puntos a los cuales se van a unir/anclar nuestras relaciones (que luego llamaremos links o edges). Los llamo puntos porque es un concepto más teórico que físico o real. Estas entidades pueden tomar forma de: genes, marcas epigenéticas, metabolitos, mRNA, microRNA, regiones de la corteza cerebral, celulas únicas, genomas de organismos, etc. Como el lector podrá entender, este es una amplia y diversa gama de opciones posibles.

Redes multicapa

Estas redes implican dos o mas capas de redes regulatorias. Por ejemplo, en una podemos tener las redes regulatorias de la co-expresión de genes y en la otra capa las redes regulatorias de las marcas epigenéticas para estos genes. Cada una tendría unidades de medida distintas (si ponderamos los links) y los vertices son distintos (en la primera son los genes y en la segunda son las marcas epigenéticas). Aqui ganamos un nivel mas de interacción, “inter-capa”, pero también un nivel mas de complejidad.

Concepto básicos de redes

  • Centrality: Es un valor cuantitativo que indica cuan importante es un vertice dentro de una red.
  • Degree: Es el número de links/edges que llegan (in-degree) o salen del vertice (out-degree). \[[total]_{degree}=[out]_{degree}+[in]_{degree}\]
  • hubs: Son los vertices con inusual número de conexiones (“high degree”)
  • geodesic distance: Es el menor número de conexiones posibles que unen dos puntos (small-world effect, six degrees of separation).
  • community or clusters: Son subredes que pueden ser identificadas dentro de la red general.
  • connected network: Cuando todos los vertices de la red estan conectados, por al menos un solo edge.

Creación de una red simple

Descargando paquetería

Como mencionamos antes vamos a trabajar con el paquete igraph

install.packages("igraph", dependencies = T)

Funciones básicas de igraph

Creación de un gráfico vacío

library(igraph)
g <- igraph::make_empty_graph()
plot(g)
title("Grafico VACIO")

Creación de un grafico con “n” nodos definidos e interacciones

Vamos a crear un gráfico con 13 nodos, con uniones no ponderadas y no direccionadas entre los nodos 3 y 5, y entre los nodos 12 y 6.

g <- igraph::make_graph(edges = c(3,5,12,6),
                        n = 13,
                        directed = FALSE)
g
## IGRAPH a32aafc U--- 13 2 -- 
## + edges from a32aafc:
## [1] 3-- 5 6--12
plot(g)

Otra forma que tenemos para crear una red es introduccioendo directamente las interacciones. Por ejemplo, si queremos una red 10 nodos, con interacciones No Direccionadas entre los nodos 3 y 5, y los nodos 6 y 9.

g <- igraph::make_graph(~1,2,3,4,5,6,7,8,9,10,3--5,6--9)
plot(g)

Como vemos, acabamos de crear un gráfico no direccionado. Luego del caracter especial ~ comenzamos a enumerar los nodos y las interacciones. No necesita ser en este orden, puede ser uno y luego el otro, o mezclados. Para definir una relacion no direccionada usamos los caracteres “–”, mientras que para una realción direccionada usamos “-+”

g2 <- igraph::make_graph(~1,2,3,4,5,6,7,8,9,6-+9,5-+3)
plot(g2)

Ahora vemos que hay puntas de flechas apuntando a 9 y 3. Si comparamos los dos objetos

g
## IGRAPH 52ce254 UN-- 10 2 -- 
## + attr: name (v/c)
## + edges from 52ce254 (vertex names):
## [1] 3--5 6--9
g2
## IGRAPH 138471c DN-- 9 2 -- 
## + attr: name (v/c)
## + edges from 138471c (vertex names):
## [1] 5->3 6->9

En la primera linea del OUTPUT, vemos que ambos son objetos IGRAPH, pero el primero es Unidireccionado con Nombres, por lo que empieza con la sigla UN. Mientras que el segundo, es Direccionado y Nombrado por lo que empieza con la sigla DN. Luego podemos ver dos números, el primero define el número de vertices y el segundo el número de conexiones o links. Debajo de los + edges, en cada objeto se define los links, nuevamente se usa “–” parea los links no direccionados, y “-+” para los direccionados, donde “+” indica hacia donde apunta la flecha.

Añadiendo vertices y conexiones

Para agregar un vertice podemos usar la función add_vertice(). Podemos usar dos formas.

g      # INPUT
## IGRAPH 52ce254 UN-- 10 2 -- 
## + attr: name (v/c)
## + edges from 52ce254 (vertex names):
## [1] 3--5 6--9
g <- igraph::add_vertices(graph = g,     # grafica de origen
                          nv = 3)        # numero de vertices para agregar
g      # revisamos el resultado
## IGRAPH 2e586af UN-- 13 2 -- 
## + attr: name (v/c)
## + edges from 2e586af (vertex names):
## [1] 3--5 6--9

Para agregar conexiones podemos usar la función add_edges

g <- igraph::add_edges(graph = g,            # grafica de origen
                       edges = c(1,3,9,4)) # conexiones a agregar 
g
## IGRAPH 7ae5f62 UN-- 13 4 -- 
## + attr: name (v/c)
## + edges from 7ae5f62 (vertex names):
## [1] 3--5 6--9 1--3 4--9

Tambien podemos usar el operador “+” para agregar edges

g <- g + igraph::edges(c(3,4,7,2))
g
## IGRAPH c47eff9 UN-- 13 6 -- 
## + attr: name (v/c)
## + edges from c47eff9 (vertex names):
## [1] 3--5 6--9 1--3 4--9 3--4 2--7

Tambien podemos usar la paquetería magrittr para concatenar varias acciones en una sola linea de comando. Por ejemplo, supongamos que deseamos agregar uniones entre los vertices 5 y 7, luego agregar los vertices 14,15 y 16 , para finalmente agregar las conexciones entre los vertices 14,15 y los vertices 7 y 16.

g <- g %>% add_edges(edges = c(5,7)) %>% add_vertices(3) %>% add_edges(edges = c(14,15 , 7,16))
g
## IGRAPH 28e4a1c UN-- 16 9 -- 
## + attr: name (v/c)
## + edges from 28e4a1c (vertex names):
## [1] 3--5 6--9 1--3 4--9 3--4 2--7 5--7
## + ... omitted several edges

Cambio de nombre de los vertices

Como vemos algunos de los vertices que creamos posterior a la generación del objeto igraph quedaron sin nombrar. Aunque no aparezca un ID por cada punto, R (internamente) conserva un identificador para cada punto. Si nosotros deseamos cambiar o asignar el nombre a los vertices podemos hacerlo con la funcion V()$name

V(g) # verificamos que no hay nombres en los vertices
## + 16/16 vertices, named, from 28e4a1c:
##  [1] 1    2    3    4    5    6    7    8    9    10   <NA> <NA> <NA> <NA> <NA>
## [16] <NA>
length(g) # numero de vertices
## [1] 16
V(g)$name <- LETTERS[1:(length(g))]  # asignamos los nombres
V(g)
## + 16/16 vertices, named, from 28e4a1c:
##  [1] A B C D E F G H I J K L M N O P

Ahora, si volvemos a graficar veremos que todos los vertices tienen nombre

plot(g)

Ahora podemos usar los nombres de los vertices para definir las conexiones.

g <- igraph::add_edges(graph = g,edges = c("L","M", "K", "H", "M","K", "J", "H"))
plot(g)

Cargando archivos

Un archivo que defina una red realmente es muy sensillo. Solo necesitamos una tabla con dos columnas, la primera son los vertices de partida y la segunda columna los vertices de llegada. Este archivo puede estar en formato TXT o CSV

n0 <- read.csv("/Users/alfredocardenasrivera/Project/bioinformatics_2022/data/NetworkAnalysis/Arabidopsis_Mendosa.sif", sep = "\t", header = F)
head(n0)
##     V1 V2   V3
## 1 EMF1  R  AP1
## 2 EMF1  R  LFY
## 3 EMF1  A TFL1
## 4  AP1  R   AG
## 5  AP1  A  LFY
## 6  LFY  A  AP1

Este archivo muestra tres columnas. Para construir nuestra red solo vamos a necesitar la primera y la tercera columna, ya que en estas es donde se define el vertice de partida y el vertice de llegada, respectivamente. Para construir la red, vamos a usar la función graph.edgelist()

n1 <- igraph::graph.edgelist(el = as.matrix(n0[,c(1,3)]),   # damos los datos
                             directed = T)                  # indicamos direccion
plot(n1)

Tambien podemos cargar la información de una red que se encuentra almacenada en archivos con extensiones específicas como “ncol”

s1 <- igraph::read_graph(file = "/Users/alfredocardenasrivera/Project/bioinformatics_2022/data/NetworkAnalysis/PredictedInteractions1000000.ncol",format = "ncol")
plot(s1, vertex.label = NA, vertex.size = 2)

Determinación de número de vertices y conexiones

length(n1)           #numero de vertices
## [1] 10
V(n1)                # numero y nombre de los vertices
## + 10/10 vertices, named, from 6834815:
##  [1] EMF1 AP1  LFY  TFL1 AG   LUG  AP3  PI   UFO  SUP
E(n1)                # numero y definicion de todas las conexiones
## + 21/21 edges from 6834815 (vertex names):
##  [1] EMF1->AP1  EMF1->LFY  EMF1->TFL1 AP1 ->AG   AP1 ->LFY  LFY ->AP1 
##  [7] LFY ->TFL1 TFL1->AG   TFL1->LFY  AG  ->AP1  LUG ->AG   LFY ->AP3 
## [13] LFY ->PI   UFO ->AP3  UFO ->PI   SUP ->PI   SUP ->AP3  AP3 ->PI  
## [19] AP3 ->AP3  PI  ->AP3  PI  ->PI

Metadata de los redes

En general, los objetos igraph pueden guardar informacion relacionada (metadata) a:

  • los vertices: get.vertex.attribute()
  • las conexiones: get.edge.attribute()
  • la red: list.graph.attributes()
igraph::get.vertex.attribute(graph = n1)
## $name
##  [1] "EMF1" "AP1"  "LFY"  "TFL1" "AG"   "LUG"  "AP3"  "PI"   "UFO"  "SUP"
igraph::get.edge.attribute(graph = n1)
## list()
igraph::get.graph.attribute(graph = n1)
## named list()

En el archivo original, la segunda columna (que obviamos al momento de crear el objeto igraph), contiene la inforamción cualitativa de las conexciones. Siendo “R” para represión y “A” para activación. Para guardar esta información en nuestro objeto igraph, podemos agregar esa información como un slot en la lista que define la metadata para las conexiones. Para esto, usaremos el operador $ para crear el nombre del slot y llenar la información en un solo paso.

E(n1)$type <- n0[,2]
igraph::get.edge.attribute(graph = n1)
## $type
##  [1] "R" "R" "A" "R" "A" "A" "R" "R" "R" "R" "R" "A" "A" "A" "A" "R" "R" "A" "A"
## [20] "A" "A"

Ahora vemos que en el objeto lista que define la metadata para las conexiones de la red (edge.attributes) aparece un slot llamado type que es un vector de string y contiene la información de las conexiones.

Ahora vamos a agregar un segundo vector a los metadatos de las conexiones para definir el color que queremos que tenga cuando se imprima la gráfica. Primero llenaremos todo el vector con el color rojo y luego le diremos que cambie a verde en los lugares donde el tipo de conexión sea de activacion, es decir “A”. Para ello usaremos el siguiente código E(n1)$type == "A".

igraph::E(n1)$color <- rep("red", length = length(E(n1)))    # definimos todos como rojo
igraph::get.edge.attribute(graph = n1)
## $type
##  [1] "R" "R" "A" "R" "A" "A" "R" "R" "R" "R" "R" "A" "A" "A" "A" "R" "R" "A" "A"
## [20] "A" "A"
## 
## $color
##  [1] "red" "red" "red" "red" "red" "red" "red" "red" "red" "red" "red" "red"
## [13] "red" "red" "red" "red" "red" "red" "red" "red" "red"
E(n1)$color[E(n1)$type == "A"] <- "green"
igraph::get.edge.attribute(graph = n1)
## $type
##  [1] "R" "R" "A" "R" "A" "A" "R" "R" "R" "R" "R" "A" "A" "A" "A" "R" "R" "A" "A"
## [20] "A" "A"
## 
## $color
##  [1] "red"   "red"   "green" "red"   "green" "green" "red"   "red"   "red"  
## [10] "red"   "red"   "green" "green" "green" "green" "red"   "red"   "green"
## [19] "green" "green" "green"

La función plot() reconoce automaticamente la columna color y asigna los colores a las conexiones

plot(n1)

Asignación de pesos

Los pesos en las conexiones nos sirven para ponderar la fuerza o la direccionalidad de las conexiones. Para este ejemplo, usaremos valores aleatoreos con la función sample(), que la asignaremos a la metadata de las conexiones.

E(n1)$weight <- sample(x = 1:4, 
                       size = length(E(n1)),
                       replace = T)
E(n1)$weight
##  [1] 4 4 4 4 3 4 2 2 3 2 1 1 4 1 3 2 4 2 4 4 3

Ahora para graficarlo, podemos representar esta ponderación como grosor de las conexiones o como una etiqueta de estas

plot(n1, edge.width = E(n1)$weight)

plot(n1, edge.label = E(n1)$weight)

o ambos.

plot(n1, edge.width = E(n1)$weight, edge.label = E(n1)$weight)

Generación de redes aleatorias

ERDOS RENYI MODEL

A veces es necesario generar redes aleatorias para probar nuestros pipelines o hacer algunos cálculos con nuestra red objetivo. Para ello, este paquete cuenta con la función erdos.renyi.game() que genera dichas redes. Esta función necesita al menos dos parametros. El parámetro n que define el número de vertices, y el parámetro p.or.m que indica la probabilidad que exista un enlace entre los vertices “p” y “m”. Opcionalmente podemos indicar si la red es direccionada con el parámetro directed, y si deseamos que haya loops en la red con el parámetro loops

r1 <- igraph::erdos.renyi.game(n = 200,      # numero de vertices
                               p.or.m = 0.01,# probabilidad que se de el enlace 
                               directed = F, # si es direccionada o no
                               loops = F)   # si presenta loops
plot(r1, vertex.label = NA, vertex.size = 2)

WATTS STROGATZ MODEL

Si bien es cierto el modelo de “Erdos-Renyi” es bastante facil de construir, tiene una gran desventaja: asume una misma probabilidad para todos los enlaces, cuando esto no sucede en la realidad. Una de las caracteristicas de las redes es que estas se subdibiden en comunidades o vecindarios. Teniendo esto en mente, Ud podrá imaginar o suponer que la frecuencia con la que ocurran las conexiones entre vertices de un mismo vencidario debe ser mayor que con otros vertices fuera de ese vencidario. Este principio es recogido por el modelo de Watts Strogatz

r2 <- igraph::watts.strogatz.game(dim = 1, # tamaño del lattice inicial 
                                  size = 200, # tamaño de los lattice a lo largo de cada dimensión
                                  nei = 1, # numero de vecindarios
                                  p = 0.05,
                                  loops = F,
                                  multiple = F)
plot(r2,vertex.label = NA, vertex.size = 2)

La ventaja de este modelo es que el principio de “mundo pequeño” está presente y podemos encontrar algunos hubs. Pero no podemos controlar el número de nodos (al menos en una forma sencilla)

BARABASI MODEL

El tercer modelo, el de BARABASI, simula mejor las redes biológicas.

r3 <- igraph::barabasi.game(n = 200,power = 1,m = 1, directed = F)
plot(r3, vertex.label = NA, vertex.size = 2)

Calculos básicos de las redes

Camino mas corto: Geodesic distance

Este es un concepto de teoria de grafico que nos indica cual el menor numero de conexiones que se necesita para conectar un nodo a otro. Ese conocimiento por si mismo no es tan intuitivo, pero si ahora recuerda que el broker es un jugador importante en las redes, ahora podrá entender como encontrar un broker en la red. A partir del nodo por el cual pasan mas “caminos cortos”. Para encontrar el camino corto usaremos la función get.shortest.paths().

sp1 <- get.shortest.paths(n1,from="LUG",to="AP3")
sp1
## $vpath
## $vpath[[1]]
## + 5/10 vertices, named, from 6834815:
## [1] LUG AG  AP1 LFY AP3
## 
## 
## $epath
## NULL
## 
## $predecessors
## NULL
## 
## $inbound_edges
## NULL

Como vemos en este ejemplo para ir del nodo LUG al nodo AP3, necesitamos pasar por los nodos AG, AP1 y LFY. Como las conexiones tienen direccionalidad, entonces no es lo mismo ir de A a B que de B a A. Lo veremos con el ejemplo anterior pero en sentido inverso

sp2 <- get.shortest.paths(n1,from="AP3",to="LUG")
## Warning in get.shortest.paths(n1, from = "AP3", to = "LUG"): At
## core/paths/dijkstra.c:446 : Couldn't reach some vertices.
sp2
## $vpath
## $vpath[[1]]
## + 0/10 vertices, named, from 6834815:
## 
## 
## $epath
## NULL
## 
## $predecessors
## NULL
## 
## $inbound_edges
## NULL

Alternativamente podemos calcular cual es la distancia entre un nodo y el resto de nodos de la red

sp3 <- igraph::shortest.paths(graph = n1, v = "LUG")
sp3
##     EMF1 AP1 LFY TFL1 AG LUG AP3 PI UFO SUP
## LUG    7   3   5    3  1   0   6  8   7  10

O podemos ver una matriz con todos los caminos mas cortos en una red

sp4 <- igraph::shortest.paths(graph = n1)
sp4
##      EMF1 AP1 LFY TFL1 AG LUG AP3 PI UFO SUP
## EMF1    0   4   4    4  6   7   5  7   6   9
## AP1     4   0   3    4  2   3   4  6   5   8
## LFY     4   3   0    2  4   5   1  3   2   5
## TFL1    4   4   2    0  2   3   3  5   4   7
## AG      6   2   4    2  0   1   5  7   6   9
## LUG     7   3   5    3  1   0   6  8   7  10
## AP3     5   4   1    3  5   6   0  2   1   4
## PI      7   6   3    5  7   8   2  0   3   2
## UFO     6   5   2    4  6   7   1  3   0   5
## SUP     9   8   5    7  9  10   4  2   5   0

Con esta herramienta, podemos comenzar a caracterizar nuestra red. Por ejemplo podemos ver la distribución de las frecuencias relativas a lo largo de nuestras redes

par(mfrow = c(2,2))
plot(density(shortest.paths(graph = n1)), main = "Arabidopsis")
plot(density(shortest.paths(graph = r1)), main = "Edos Renyi")
plot(density(shortest.paths(graph = r2)), main = "Watts Strogatz")
plot(density(shortest.paths(graph = r3)), main = "Barabasi")

Otro calculo que nos podría ayudar a describir nuestra red, es el promedio de caminos cortos

igraph::average.path.length(graph = n1)
## [1] 4.289474
igraph::average.path.length(graph = r1)
## [1] 6.896615
igraph::average.path.length(graph = r2)
## [1] 32.72201
igraph::average.path.length(graph = r3)
## [1] 8.247337

Grados

Un parámetro importante es la cantidad de conexiones que recibe, envia o simplementa conecta un nodo con el resto de la red. En primer lugar, estamos viendo el lado opuesto del broker. En este caso, tratamos de ver cual es el nodo mas relacionado con el resto de la red. En otras palabras, si fuera una red social, queremos identificar a la persona más popular. Para ello vamos a usar la función degree().

igraph::degree(graph = n1,v = "AP3")
## AP3 
##   7
igraph::degree(graph = n1,v = "LUG")
## LUG 
##   1

En este ejemplo, vemos que AP3 esta mucho mas conectado (tiene mas edges/links/conexiones) con la red que LUG, que solo cuenta con 1. Si trabajamos con una red direccionada, podemos identificar a la persona que trata de relacionarse más con la red (OUT), o cual es la persona más consultada en la red (IN).

igraph::degree(graph = n1,v = "AP3", mode = "in")
## AP3 
##   5
igraph::degree(graph = n1,v = "AP3", mode = "out")
## AP3 
##   2

Finalmente, podemos perfilar la distribución de las conexiones con la función degree.distribution

par(mfrow=c(2,2))
plot(degree.distribution(graph = n1), type="l", ylab = "", xlab = "Grado", main = "Arabidpsis")
plot(degree.distribution(graph = r1), type="l", ylab = "", xlab = "Grado", main = "Erdos-Renyi")
plot(degree.distribution(graph = r2), type="l", ylab = "", xlab = "Grado", main = "Watts-Strogatz")
plot(degree.distribution(graph = r3), type="l", ylab = "", xlab = "Grado", main = "Barbasi")

Conectividad

Hablamos de una “connected network” cuando todos los vertices estan conectados a la red por lo menos con un edge. Veamos el ejemplo de las dos redes que tenemos hasta el momento. Primero analizemos la red del archivo importado

igraph::is.connected(n1)
## [1] TRUE
plot(n1)

Como vemos esta red esta conectada TRUE, y vemos que todos los nodos estan conectados por lo menos con un enlace al resto de la red. Por otro lado, la red g

igraph::is.connected(g)
## [1] FALSE
plot(g)

Esta red no está conectada. Existen vertices que no se unen al resto de la red.

Cluster

Para saber cuantas subredes no conectadas hay en nuestra red podemos usar la función clusters. Podemos tambien determinar cuales estan fuertemente conectado o debilmente conectados, usando parámetro mode con las opciones strong o weak, respectivamente.

igraph::clusters(r2, mode = "strong")
## $membership
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $csize
## [1] 200
## 
## $no
## [1] 1
igraph::clusters(r2, mode = "weak")
## $membership
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $csize
## [1] 200
## 
## $no
## [1] 1

Esta información la podemos guardar como metadata de los vertices, asignandolo como un slot mas

V(r2)$cluster <- igraph::clusters(graph = r2,mode = "strong")$membership
V(r2)$cluster
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Ahora podemos seleccionar solo aquellos que tengan

r2.strong <- igraph::delete.vertices(graph = r2,
                                     v = V(r2)$cluster < 3)
set.seed(1)
plot(r2.strong, vertex.size = 2, vertex.label = NA)

set.seed(1)
plot(r2, vertex.size = 2, vertex.label = NA)

Eficiencia

La eficiencia de una red es un parámetro que nos ayuda a evaluar con cuanta eficacia se comunica o pasa la información a traves de la red. En otras palabras, lo que queremos evaluar es cuan distante, en promedio, estan separados dos nodos. Si la distancia promedio (medido en caminos mas cortos) es menor en una red que en otra, entonces decimos que la primera tiene mayor eficiencia que la otra. Esto nos habla de la fortaleza de la red. Para evaluar este parámetro de la red, podemos usar las funciones global_efficiency o local_efficiency, para evaluar la eficiencia de la red como conjunto o para evaluar la de un vertice en particular, respectivamente. Tambien podemos calcular el promedio de las eficiencias locales de todos los vertices de la red con average_local_efficiency. Finalmente, una caracteristica interesante que podemos calcular, es como cambia la eficiencia global de la red si quitamos uno a uno los vertices de la red. Esto nos permitiría evaluar la relevancia de cada vertice para la eficiencia de la red. Recuerde que el resultado de estos (como vectores numéricos) los podemos guardar en la sección de attributes del slot vertices de la red.

igraph::global_efficiency(graph = n1)
## [1] 0.1348854
igraph::local_efficiency(graph = n1,vids = "AP3")
##        AP3 
## 0.09027778
igraph::local_efficiency(graph = n1,vids = "LUG")
## LUG 
##   0
igraph::local_efficiency(graph = n1)
##       EMF1        AP1        LFY       TFL1         AG        LUG        AP3 
## 0.31111111 0.11111111 0.07500000 0.11666667 0.05714286 0.00000000 0.09027778 
##         PI        UFO        SUP 
## 0.18750000 0.37500000 0.37500000
igraph::average_local_efficiency(graph = n1)
## [1] 0.169881

Transitividad o coeficiente de conectividad

Se define como el número de triangulos que se pueden formar uniendo los vertices conectados al vertice objetivo. Lo que trata de evaluar este parámetro es la formación de comunidades o vencindarios en la red. Es decir, nos da un valor númerico para definir las sub-redes de nuestra red. Con la función transivity() podemos calcular la conectividad global o la conectividad local, entre otras.

# conectividad global de la red de Arabidpsis
igraph::transitivity(graph = n1, type = "global")
## [1] 0.4166667
# conectividad local de la red de Arabidopsis
igraph::transitivity(graph = n1, type = "local")
##  [1] 0.6666667 0.3333333 0.3000000 0.3333333 0.0000000       NaN 0.5000000
##  [8] 0.5000000 1.0000000 1.0000000

Esta información la podemos agregar a nuestro objeto igraph, asignandolo dentro del slot vertices

V(n1)$cc <- igraph::transitivity(graph = n1, type = "local")

El coeficiente de conectividad es uno de los métodos para definir sub-redes. Este tema, probablemente, es una de los más importantes en los análisis de redes. Existen al menos otros 3 algoritmos que nos permiten realizar este análisis:

  • fastgreedy.community(): rápido pero no muy exacto.
  • walktrap.community(): “caminante aleatorio”, la idea es que caminos cortos aleatorios tienden a estar en las mismas comunidades.
  • spinglass.community(): computacionalmente mas exigente.
fg <- igraph::fastgreedy.community(graph = r3)
wt <- igraph::walktrap.community(graph = r3)
sg <- igraph::spinglass.community(graph = r3)

par(mfrow=c(1,3))
set.seed(124)
plot(fg, r3, vertex.label=NA, vertex.size=5, main="Fast Greedy")
set.seed(124)
plot(wt, r3, vertex.label=NA, vertex.size=5, main="Random Walker")
set.seed(124)
plot(sg, r3, vertex.label=NA, vertex.size=5, main="Spin Glass")

Para obetner el numero del grupo de cada uno de los vertices en los modelos podemos usar la función membership. Nuevamente podemos almacernar esta infomación en el objeto igraph original en la el slot de vertices V()

V(r3)$fg <- igraph::membership(fg)
V(r3)$wt <- igraph::membership(wt)
V(r3)$sg <- igraph::membership(sg)
igraph::get.vertex.attribute(r3)
## $fg
##   [1]  6 12  6  8 13 14 13  4  8 10 10  6  4  1  1  3  8  3 11  2  2 12  1 13  2
##  [26]  8  2 11  5  8  1 13 10  9  1  6  7 11  2  9 11  1  7  6 15 12  4  7  1 11
##  [51]  9  4  2 10  9  6  2 12  7  7 14  6 10 10  3  4  9 13  5 10 14  4  6 13 15
##  [76]  7 12  5 10  2  4  4  2  4 15  5  3  7  7  6  5  1  1  5 13 15 13  4  5 11
## [101]  5  9  7 10  8  4  5  5  2  9  4  6  5 10  9 12  8  2  4  2  8  6  5  2 15
## [126]  7  2  8 11  4  1 11 14  1 12  3  2 10 14  4  9  3  2  8  3  6  9  8  7  3
## [151]  1 12  1  1 11  5  6  5 11 14  9  3  7  2 14  1  2 15  3  7  3  3  4  3 12
## [176]  3  3  8  9  5  1  4  3 11  1 14  1 13 14  9 12  8 15  5  5  3 14  1  6  7
## 
## $wt
##   [1]  5  2  5  5  5  2  5  4  5  5  1  5  4  5 15  1 13  1  3  7  7  2  5 16  7
##  [26] 13  8  3  1  5  5 16  5  6 12  5  5  3  7  6  3 12  5  5  1  2  4 14 15  3
##  [51]  6  4  7  1  6  5  8  2  9 14  2  5  1  1 10  4  6  5  1  5  2  4  5 16  1
##  [76]  9  2  1  5  7  4  4  8  4  1  1 10  9 14  5  1  5 12  1  5  1 16  4  1  3
## [101]  1  6  5  1 13  4  1  1  8  6  4  5  1  5  6  2  5  8  4  8 13  5  1  7  1
## [126]  9  8  5  3  4 12  3  2 15  2  1  7  5  2  4  6  1  8  5  1  5  6 13  9 11
## [151]  5  2 12  5  3  1  5  1  3  2  6 11 14  7  2 12  7  1 11  5 10 10  4 11  2
## [176] 10 11 13  6  1  5  4 11  3 15  2 12  5  2  6  2  5  1  1  1 10  2 15  5  9
## 
## $sg
##   [1] 13 12  1 21  1  6  1  4 21 11 11  1 16  8  8 17 20 17 22 14 14 12 15  1 14
##  [26] 20  9 22  3 21  8  1 11 11 15  1  7 22 14 11 22 15  7  1 18 12  4 10  8 22
##  [51] 11  4 14 11  5 13  9 12  7 10  6 13 11 11 19  4  5  1  3 11  6 16  1  1 18
##  [76]  7 12  3 11 14  4  4  9 16 18  3 19  7 10 13  2  8 15  3  1 18  1  4  2 22
## [101]  3  5  7 11 20 16  3  3  9  5  4 13  2 11 11 12 21  9 16  9 20  1  3 14 18
## [126]  7  9 21 22  4 15 22  6  8 12 17 14 11  6  4  5 17  9 21 17  1  5 20  7 17
## [151]  8 12 15  8 22  3 13  2 22  6  5 17 10 14  6 15 14 18 17  7 19 19  4 17 12
## [176] 19 17 20 11  3  8  4 17 22  8  6 15  1  6 11 12 21 18  3  2 19  6  8  1  7

Graficando redes

En la programción orientada a objetos, las funciones deben crearse y definirse cuando se construye el objeto. Esto no es diferente en el paquete igraph, pero que pasa cuando usamos la función plot en un objeto igraph y obtenemos como un output un gráfico. Lo que realmente está pasando aquí, es que estamos ejecutando la función plot.igraph de los objetos igraph. Para este ejemplo vamos a trbajar con datos de interactoma de ser humano. Primero importamos la base de datos en formato .ncol con la función `rea.graph``

gw<-read.graph("/Users/alfredocardenasrivera/Project/bioinformatics_2022/data/NetworkAnalysis/PredictedInteractions1000000.ncol",
               format="ncol")
str(gw)
## Class 'igraph'  hidden list of 10
##  $ : num 277
##  $ : logi FALSE
##  $ : num [1:498] 1 3 5 7 7 8 10 8 10 11 ...
##  $ : num [1:498] 0 2 4 6 4 3 9 2 7 7 ...
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ : NULL
##  $ :List of 4
##   ..$ : num [1:3] 1 0 1
##   ..$ : Named list()
##   ..$ :List of 1
##   .. ..$ name: chr [1:277] "RFC4" "RFC3" "PSMA2" "PSMA5" ...
##   ..$ :List of 1
##   .. ..$ weight: num [1:498] 2412885 2344830 804075 804075 804075 ...
##  $ :<environment: 0x1410a72b0>
V(gw)
## + 277/277 vertices, named, from eab0bf7:
##   [1] RFC4      RFC3      PSMA2     PSMA5     LSM7      LSM2      SNRPD2   
##   [8] LSM5      PSMA3     SNRPD1    SNRPD3    LSM3      F2        F10      
##  [15] FGA       FGB       PSMA4     PSMA7     PSMA6     LSM6      LSM4     
##  [22] PSMB7     MCM3      MCM2      MCM4      SNRPG     JAK2      CSF1R    
##  [29] STAT5A    FYN       PDGFRB    PLG       PSMB1     LSM1      STAT5B   
##  [36] SNRPF     RFC5      PSMB3     ZAP70     TEC       CRKL      PDGFRA   
##  [43] FLT1      RAF1      MAP2K2    KIT       CDC6      STAT1     GRB2     
##  [50] F9        LSM8      TAGLN2    ACTG1     BRAF      RFC2      CHMP4C   
##  [57] CHMP4A    CHMP4B    CHMP6     PIK3R2    FGR       VAV1      FRK      
##  [64] INSR      CBL       PSMB2     PSMB6     PSMC2     PSMC1     WASPIP   
## + ... omitted several vertices
E(gw)
## + 498/498 edges from eab0bf7 (vertex names):
##  [1] RFC4  --RFC3   PSMA2 --PSMA5  LSM7  --LSM2   SNRPD2--LSM5   LSM7  --LSM5  
##  [6] PSMA5 --PSMA3  SNRPD1--SNRPD3 PSMA2 --PSMA3  LSM5  --SNRPD3 LSM5  --LSM3  
## [11] F2    --F10    FGA   --FGB    PSMA2 --PSMA4  PSMA5 --PSMA4  PSMA5 --PSMA7 
## [16] PSMA3 --PSMA6  LSM5  --LSM6   LSM7  --LSM4   PSMA2 --PSMA7  LSM2  --SNRPD2
## [21] LSM2  --SNRPD1 LSM2  --SNRPD3 PSMA5 --PSMB7  MCM3  --MCM2   MCM3  --MCM4  
## [26] MCM2  --MCM4   LSM2  --LSM5   SNRPD3--SNRPG  JAK2  --CSF1R  JAK2  --STAT5A
## [31] FYN   --PDGFRB F2    --PLG    PSMB7 --PSMB1  LSM4  --LSM1   JAK2  --STAT5B
## [36] LSM2  --LSM3   LSM2  --LSM4   SNRPD1--SNRPF  SNRPD2--SNRPG  PSMA7 --PSMA6 
## [41] PSMA3 --PSMA4  LSM7  --LSM3   SNRPD3--SNRPF  LSM7  --SNRPD2 RFC3  --RFC5  
## [46] SNRPD3--LSM6   LSM5  --SNRPG  RFC4  --RFC5   PSMB7 --PSMB3  JAK2  --ZAP70 
## + ... omitted several edges

Ahora removemos el grupo dos, y proyectamos la red resultante.

#separar los datos con mayor conectividad
V(gw)$cluster<-clusters(gw)$membership
g.human<-simplify(delete.vertices(gw,V(gw)$cluster!=2))


par(mfrow=c(2,2))
plot(gw, vertex.label=NA, vertex.size = 2, main = "Raw")
plot(g.human, vertex.label=NA, vertex.size = 2, main = "Clean")

plot(g.human,vertex.label=NA, 
     vertex.size=sample(1:5, vcount(g.human), replace=T))