Redes de Interacción

Instalar Bioconductor

Bioconductor (http://www.bioconductor.org) es un proyecto para crear y administrar paquetes de R para resolver problemas bioinformáticos. Incluye paquetes que implementan algoritmos y otros que funcioan como repositorios de datos.

Para este taller vamos a utilizar un paquete, simpIntLists, que contiene información sobre datos de interacción proteína-proteína para varos organsmos. A nosotros nos interesan los datos de Arabidopsis.

Para instalarlo primero tener que ejecutar unos scripts de inicialización, después cargar los paquetes base y en tercer lugar simpIntList:

source("http://bioconductor.org/biocLite.R")## Bioconductor version 3.2 (BiocInstaller 1.20.3), ?biocLite for help## A new version of Bioconductor is available after installing the most## recent version of R; see http://bioconductor.org/install

biocLite()

biocLite("simpIntLists")## BioC_mirror: https://bioconductor.org## Using Bioconductor 3.2 (BiocInstaller 1.20.3), R 3.2.5 (2016-04-14).## Installing package(s) 'simpIntLists'## installing the source package 'simpIntLists'## ## The downloaded source packages are in## 'C:\Users\Marcelo\AppData\Local\Temp\RtmpyqfAAO\downloaded_packages'library(simpIntLists)

También necesitamos el paquete igraph para procesar grafos

# La primera vez:# install.packages("igraph")library(igraph)

## ## Attaching package: 'igraph'

## The following objects are masked from 'package:stats': ## ## decompose, spectrum

## The following object is masked from 'package:base': ## ## union

Asignamos a una variable la lista con los datos de interacción proteína - proteína de Arabidopsis thaliana:

data(ArabidopsisBioGRIDInteractionOfficial)arab.int <-ArabidopsisBioGRIDInteractionOfficialhead(arab.int)

## [[1]] ## [[1]]$name ## [1] "BRCA2(IV)" ## ## [[1]]$interactors ## [1] "ATRAD51" "DMC1" "RAD51" "ATDSS1(I)" "ATDSS1(V)" ## ## ## [[2]] ## [[2]]$name ## [1] "ATRAD51" ## ## [[2]]$interactors ## [1] "BRCA2(IV)" "ATRAD54" "RAD54" "BRCA2B" "DMC1" "RAD51" ## ## ## [[3]] ## [[3]]$name ## [1] "DMC1" ## ## [[3]]$interactors ## [1] "BRCA2(IV)" "BRCA2B" "ATRAD51" ## ## ## [[4]] ## [[4]]$name ## [1] "TOC1" ## ## [[4]]$interactors ## [1] "PIF4" "PIL6" "TOC1" "APRR9" "ZTL" "LKP2" "FKF1" "APRR3" ## [9] "PIF3" "PIL1" "PIL2" "PIL5" "LHY" "CCA1" "GI" "APRR5" ## [17] "TIC" "ABI3" ## ## ## [[5]] ## [[5]]$name ## [1] "PIF4" ## ## [[5]]$interactors ## [1] "TOC1" "HRB1" "PIF7" "PIF3" "RGA1" "PIL6" "PHYB" "PHYA" "HFR1" "RGL2" ## ## ## [[6]] ## [[6]]$name ## [1] "PIL6" ## ## [[6]]$interactors ## [1] "TOC1" "PIF3" "PHYB" "PIF4" "PHYA" "HFR1"

length(arab.int)

## [1] 2109

Esta elemento de esta lista registra una proteína y con cuál o cuáles se detectó interacción, incluyendo interacciones consigo misma.

A partir de estos datos tenemos que construir una “lista de aristas”. Esto es, en lugar de una lista, necesitamos una matriz o un dataframe con dos columnas. Cada fila de esta estructura de datos registra una interacción. No importa que proteína va a la izquierda y cual a la derecha:

DMC1 BRCA2(IV) DMC1 BRCA2B DMC1 ATRAD51 TOC1 PIF4 TOC1 PIL6 TOC1 TOC1 TOC1 APRR9 TOC1 ZTL

Podemos ver dos formas diferentes de generar esta matriz y comparar los tiempos de ejecución. Podemos ver cuanto más lento es el loop for():

system.time({ df <- data.frame(A=NULL, B=NULL) for(i in 1:length(arab.int)){ df <- rbind(df, cbind(arab.int[[i]]$name, arab.int[[i]]$interactors)) }})

## user system elapsed ## 5.86 0.00 5.93

head(df)

## V1 V2 ## 1 BRCA2(IV) ATRAD51 ## 2 BRCA2(IV) DMC1 ## 3 BRCA2(IV) RAD51 ## 4 BRCA2(IV) ATDSS1(I) ## 5 BRCA2(IV) ATDSS1(V) ## 6 ATRAD51 BRCA2(IV)

df es un dataframe con factores.

system.time({ df.2 <- do.call("rbind", sapply(arab.int, function(x) cbind(x$name, x$interactors) ))})

## user system elapsed ## 0 0 0

head(df.2)

## [,1] [,2] ## [1,] "BRCA2(IV)" "ATRAD51" ## [2,] "BRCA2(IV)" "DMC1" ## [3,] "BRCA2(IV)" "RAD51" ## [4,] "BRCA2(IV)" "ATDSS1(I)" ## [5,] "BRCA2(IV)" "ATDSS1(V)" ## [6,] "ATRAD51" "BRCA2(IV)"

df.2 es una matriz de caracteres.

Con la matriz df.2 ya podemos construir el grafo:

arab.gr <- graph_from_edgelist(df.2, directed = F)arab.gr

## IGRAPH UN-- 2109 8100 -- ## + attr: name (v/c) ## + edges (vertex names): ## [1] BRCA2(IV)--ATRAD51 BRCA2(IV)--DMC1 BRCA2(IV)--RAD51 ## [4] BRCA2(IV)--ATDSS1(I) BRCA2(IV)--ATDSS1(V) BRCA2(IV)--ATRAD51 ## [7] ATRAD51 --ATRAD54 ATRAD51 --RAD54 ATRAD51 --BRCA2B ## [10] ATRAD51 --DMC1 ATRAD51 --RAD51 BRCA2(IV)--DMC1 ## [13] DMC1 --BRCA2B ATRAD51 --DMC1 TOC1 --PIF4 ## [16] TOC1 --PIL6 TOC1 --TOC1 TOC1 --APRR9 ## [19] TOC1 --ZTL TOC1 --LKP2 TOC1 --FKF1 ## [22] TOC1 --APRR3 TOC1 --PIF3 TOC1 --PIL1 ## + ... omitted several edges

Y ahora un gráfico:

plot(arab.gr)

Lo podemos mejorar:

plot(arab.gr, vertex.size=3, vertex.label=NA)

Este grafo presenta muchos loops. En términos biológicos son proteínas que interactúan consigo mismas. También podría ser un artefacto de la técnicas. Podemos hacer un grafo sin estos loops:

arab.gr.2 <- simplify(arab.gr)arab.gr.2

## IGRAPH UN-- 2109 3946 -- ## + attr: name (v/c) ## + edges (vertex names): ## [1] BRCA2(IV)--ATRAD51 BRCA2(IV)--DMC1 BRCA2(IV)--RAD51 ## [4] BRCA2(IV)--ATDSS1(I) BRCA2(IV)--ATDSS1(V) ATRAD51 --DMC1 ## [7] ATRAD51 --RAD51 ATRAD51 --ATRAD54 ATRAD51 --RAD54 ## [10] ATRAD51 --BRCA2B DMC1 --BRCA2B RAD51 --ATRAD54 ## [13] ATDSS1(I)--BRCA2B ATDSS1(I)--EMB2719 ATDSS1(I)--At1g75990 ## [16] ATDSS1(V)--BRCA2B ATDSS1(V)--EER5 ATDSS1(V)--EMB2719 ## [19] ATDSS1(V)--At1g75990 ATDSS1(V)--CEN2 TOC1 --PIF4 ## [22] TOC1 --PIL6 TOC1 --APRR9 TOC1 --ZTL ## + ... omitted several edges

plot(arab.gr.2, vertex.size=3, vertex.label=NA)

¿Cuáles son las diez proteínas con mayor número de interacciones?

head( sort(degree(arab.gr.2), decreasing = T), 10)

## CDC2 SKP1 CKS1 CKS2 CDKB1;1 CAK4 ATCUL1 KRP2 AHK2 ## 94 59 51 49 46 43 41 41 40 ## CYCH;1 ## 40

Podríamos ver en NCBI qué se conoce de estas proteínas de A. thaliana.

Vemos que en arab.gr.2 quedan una cuantas proteínas aisladas, o sea, con grado cero. Vamos analizar esto en más detalle. Además algunos nodos (proteínas) aparecen muy conectados y otros mucho menos. Veamos esto.

grados.ordenados <- sort(degree(arab.gr.2), decreasing=T )

¿Cuantos nodos hay de grado cero? ¿Y de grado uno?

table(grados.ordenados == 0)

## ## FALSE TRUE ## 2081 28

table(grados.ordenados == 1)

## ## FALSE TRUE ## 1138 971

Y ahora hagamos un gráfico de la distribución de todos los grados:

plot(grados.ordenados, type="h", xlab = "nodos ordenados por grado", ylab= "grado")

Otra cuestión que nos puede interesar es determinar si hay componentes. Es decir, subgrafos dentro del grafo que están separados del resto.

gr.2.components <- components(arab.gr.2)names(gr.2.components)

## [1] "membership" "csize" "no"

gr.2.components$no

## [1] 162

Hay 162 componentes. Algunos son grandes, con muchas proteínas, y otros formados por una sola.

Podemos ver que hay 28 componentes con una sola proteína, 89 con dos, etc. Esta situación es frecuente. Generalmente nos interesa el componente más grande, que se suele llamar componente gigante:

table(gr.2.components$csize)

## ## 1 2 3 4 5 6 8 9 10 11 12 15 1679 ## 28 89 19 8 2 7 2 1 2 1 1 1 1

¿Cuál es el indice del componente de mayor tamaño?

which.max(gr.2.components$csize)

## [1] 1

Y ahora extraemos los nodos de ese componente:

nodos.gigante <- gr.2.components$membership[ gr.2.components$membership == 1 ]nodos.gigante <- names(nodos.gigante)

Y con esto armamos el subgrafo inducido de este componente:

gigante.gr <- induced_subgraph(arab.gr.2, nodos.gigante)plot(gigante.gr, vertex.size=4, vertex.label=NA)

Otra cosa que nos podría interesar es el subgrafo inducido a partir de una proteína particular, por ejemplo TGA1:

TGA1.gr <- induced_subgraph( arab.gr.2, subcomponent(arab.gr.2, "TGA1") )plot(TGA1.gr, vertex.size=1, vertex.color="white")

El análisis de estos grafos se puede profundizar mucho más, por ejemplo buscando nodos de acuerdo a diferentes criterios de centralidad, buscar clusters dentro del componente gigante, ver la relación entre nodos y otra información que se disponga, por ejemplo, categoría funcional, etc., etc.