---
title: "R Notebook - Prise en main igraph"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_notebook:
    theme: paper
    highlight: tango
---

```{r echo=F}
blue = '<font style="color: #000099;">'
black = '<font style="color: #000000;">'
red = '<font style="color: #990000;">'
```


`r blue`

Prise en main de la librairie R - `igraph`
========================================

La librairie [igraph](https://igraph.org) met à disposition tout un ensemble de fonctions pour le traitement et la visualisation de graphes. Elle est codée en C. Nous allons utiliser aujourd'hui son interfaçage avec R. Pour la charger :

```{r}
library(igraph)
```

Différentes manières de l'installer :

-   `install.packages` ou avec rstudio
-   `conda` : voir <http://silico.biotoul.fr/p/Conda>

Pour charger un graphe (différents format possibles : pajek, newick, ...) :

```{r}
g = read.graph("http://silico.biotoul.fr/site/images/9/9f/Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.tgr", directed=FALSE)
g
```

**Consulter** l'aide de la fonction (`?read_graph`) pour voir les autres formats supportés.

`r black`

Pour l'afficher, `plot(g)` suffit, les paramètres par défauts seront utilisés et un algorithme de dessin sera sélectionné :
```{r}
plot(g)
```

`r blue`

Pour l'afficher, il faut au préalable en effectuer le dessin (layout) : Soit en une ligne en passant la fonction de dessin :

```{r}
plot(g, layout=layout.fruchterman.reingold)
```

soit en sauvegardant ce layout dans une variable :

```{r}
lfr = layout.fruchterman.reingold(g)
plot(g, layout=lfr, vertex.size=3, vertex.label=NA)
```

**Remarque :** le fait de sauvegarder la disposition des sommets, (i) évite de la recalculer pour chaque plot, (ii) permet d'avoir la même disposition à chaque plot (certains algorithmes, tel que FR, ne donnent pas la même disposition à chaque fois).

**Consulter** l'aide des fonction `plot.igraph` et `layout.fructhterman.reingold` pour voir les options ainsi que les autres algorithmes de dessin disponibles.

Vous trouverez la documentation de la librairie sur le site dédié. Pour celle de l'interface R en ligne : <https://igraph.org/r/>

-   Pour obtenir la liste des sommets : `V(g)`
-   la liste des arêtes : `E(g)`

`r black`
```{r}
V(g)
```


```{r}
E(g)
```


`r red` Quel est l'ordre du graphe ? Combien a-t-il d'arêtes ?

`r black`
ordre

```{r}
length(V(g))
```

nb d'arêtes

```{r}
length(E(g))
```

`r blue`


On peut assigner des étiquettes aux sommets : `V(g)$name = vector_of_labels`

Charger les étiquettes des sommets :

-   [Cleandb\_Luca\_1\_S\_1\_1\_65\_Iso\_Tr\_1-CC1.cod](http://silico.biotoul.fr/site/images/6/61/Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.cod)

`r black`
En décomposant, appeçu du contenu du fichier :
```{r}
head(read.table("http://silico.biotoul.fr/site/images/6/61/Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.cod"))
```
La 2ème colonne contient le nom des sommets. Il s'agit donc de l'assigner aux sommets du graphe que l'on a déjà chargé : 


```{r}
V(g)$name=as.character(read.table("http://silico.biotoul.fr/site/images/6/61/Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.cod")[,2])
plot(g, layout=lfr, vertex.size=3, vertex.label=V(g)$name)
```

`r blue`

Les paramètres de la fonction plot de igraph sont décrits dans la [documentation](https://igraph.org/r/doc/plot.igraph.html).

On peut redéfinir les paramètres par défaut :

```{r}
igraph.options(vertex.label.cex=.8) # font size
igraph.options(vertex.label.family='sans')
igraph.options(vertex.size=3) 
igraph.options(vertex.color=NA) 
```


Il est possible de stocker de l'information sur le graphe, les sommets et/ou les arêtes avec les [fonctions dédiées](https://igraph.org/r/doc/graph_attr_names.html) :

```{r}
vertex_attr(g, name="name")
```

```{r}
vertex_attr(g, name="name", index=10)
```

`r red` Quel est le **diamètre** du graphe ? (cf. `diameter` ou `distances`)

```{r}
diameter(g)
```

`r black` la même chose avec `distances`

```{r}
distances(g)[1:6, 1:6]
```


```{r}
max(distances(g))
```

`r red` Lister les **points d'articulation** (`articulation.points`)

```{r}
articulation.points(g)
```

Longueur moyenne des plus courts chemins (sans valuation) (`average.path.length`)

```{r}
average.path.length(g)
```

`r black` Pour obtenir la même chose à partir de distances (il faut ajouter la diagonale (=0) correspondant à la moyenne)

```{r}
mean( distances(g) + diag(42)*1.702671)
```

`r red` Afficher sa '''représentation canonique''' (`canonical.permutation`) et la matrice d'adjacence correspondante

`r black`
La permutation

```{r}
canonical.permutation(g)
```

La matrice d'adjacence correspondant à cette permutation

```{r}
get.adjacency(permute.vertices(g, canonical.permutation(g)$labeling))
```

`r red` Lister les **composantes connexes** (`clusters`) avant et après suppression du point d'articulation

```{r}
clusters(g)
```

`r black`
Une seule, bien sûr.

Suppression du point d'articulation

```{r}
g1 = delete_vertices(g, articulation_points(g))
plot(g1)
```

ou aussi

```{r}
plot(g - articulation_points(g))
```

Composantes connexes

```{r}
clusters(g1)
```

`r red` Obtenir le line graph (`line.graph`)

```{r}
lg = line.graph((g))
plot(lg, vertex.label=NA)
```

`r black`
Il semble ressembler à g mais en plus dense.

`r blue`

Fichier dressing au format ncol : [dressing.ncol](http://silico.biotoul.fr/site/images/0/02/M1BBS_Graphe_dressing.ncol), [Bellman-Ford.ncol](http://silico.biotoul.fr/site/images/b/b0/M1BBS_Graphe_Bellman-Ford.ncol) et celui pour Floyd-Warshall [Floyd-Warshall.ncol](http://silico.biotoul.fr/site/images/d/dd/M1BBS_Graphe_Graphe.Floyd-Warshall.ncol)

![dressing](http://silico.biotoul.fr/site/images/7/73/Igraph.dressing.png)

![Graphe utilisé pour l'algorithme Bellman-Ford](http://silico.biotoul.fr/site/images/a/ae/Igraph.Bellman-Ford.png)

![Graphe utilisé pour l'algorithme Floyd-Warshall](http://silico.biotoul.fr/site/images/d/df/Igraph.Floyd-Warshall.png)

`r red` Effectuer le tri topologique du graphe dressing.

`r black`
le graphe

```{r}
dressing = read.graph('http://silico.biotoul.fr/site/images/0/02/M1BBS_Graphe_dressing.ncol', format='ncol', directed=T)
dressing
```

tri topologique (la réponse à la question s'il est sans circuit est plus bas)

```{r}
topological.sort(dressing)
```

`r red`
Parcours en largeur et en profondeur, et Bellman-Ford (`graph.bfs`, `graph.dfs`, `distances`)

`r black`
Parcours en largeur

```{r}
p=bfs(dressing, root='sous-vetements', rank=TRUE, dist=TRUE, unreachable=FALSE)
p
```

Affichage

```{r}
plot(dressing, vertex.label = paste(V(dressing)$name,':', p$dist) , vertex.size=6)
```

Parcours en profondeur

```{r}
dfs(dressing, root='sous-vetements', unreachable=FALSE, father=T, dist=T)
```

`unreachable` est TRUE par défaut et correspond à refaire les appels récursifs depuis les sommets non atteints à partir du sommet source fourni.

Sans circuit/cycle ?

```{r}
is_dag(dressing)
```

Bellman-Ford

```{r}
bf = read.graph('http://silico.biotoul.fr/site/images/b/b0/M1BBS_Graphe_Bellman-Ford.ncol', format='ncol', directed=T)
bf.layout=layout.circle(bf)
bf
```

Poids présents ?

```{r}
E(bf)$weight
```

Affichage

```{r}
plot(bf, layout=bf.layout, edge.label=E(bf)$weight, vertex.size=15)
```

remarque : pas pratique, on ne sait pas trop à quel arc correspond quel poids...

```{r}
plot(bf, layout=bf.layout, edge.label=E(bf)$weight, vertex.size=15, edge.color=E(g), edge.label.color=E(g))
```

On a encore des arcs qui se superposent
```{r}
E(bf)$curved = 0.2
E(bf)$color = E(bf)
V(bf)$color = V(g)
plot(bf, layout=bf.layout, edge.label=E(bf)$weight, vertex.size=15)
```



Application de BF

```{r}
d=distances(bf, v='A', mode='out', algo='bellman-ford', weights = E(bf)$weight)
d
```

Affichage

```{r}
plot(bf, layout=bf.layout, edge.label=E(bf)$weight, vertex.label = paste(V(bf)$name,':', d), vertex.size=35 )
```

Pour le graphe FW, se sera plutôt l'algorithme de Johnson qui sera utilisé

```{r}
fw = read.graph('http://silico.biotoul.fr/site/images/d/dd/M1BBS_Graphe_Graphe.Floyd-Warshall.ncol', format='ncol', directed=T)
fw.layout=layout.circle(bf)
fw
```

Poids

```{r}
E(fw)$weight
```

Affichage

```{r}
plot(fw, layout=bf.layout, edge.label=E(fw)$weight, vertex.size=15)
```

Plus courtes distances

```{r}
shortest.paths(fw, mode='out')
```

Autre exemple à partir d'une matrice d'adjacence ou d'un data.frame
===================================================================

Création d'un graphe à partir d'une **matrice d'adjacence** :

```{r}
M = matrix(c(0,0,0,0, 1,0,1,0, 0,0,1,0, 4,2,0,0), ncol=4)
M
```

Affichage

```{r}
g=graph.adjacency(M, mode='directed', weighted=TRUE)
plot(g, edge.width=E(g)$weight)
```

Librairie très utilisée maintenant pour la manipulation de données :

```{r}
library(tidyverse)
```

Création à partir d'un data.frame (celui du 1er TP)

```{r}
links.detailed = read_delim("511145.protein.links.detailed.v11.5.txt", delim=" ")
links.detailed
```

Comme lors du 1er TP, ne gardons que les liens de coexpression et experimental \>=800 :

```{r}
links.filtered = links.detailed %>% 
  filter(protein1<protein2 & (coexpression>=800 | experimental>=800)) %>%
  select(protein1, protein2, coexpression, experimental)
links.filtered
```

Utilisons des noms de sommets plus lisibles :

```{r}
proteins = read_delim("511145.protein.info.v11.5.txt", delim="\t")
proteins
```

```{r}
g = graph_from_data_frame(d = links.filtered, directed = F, vertices = proteins)
V(g)$name = V(g)$preferred_name
```

**Exploration**

Accès aux sommets/arêtes : `V(g)` et `E(g)`.

```{r}
length(V(g))
```

Arêtes

```{r}
length(E(g))
```

Beaucoup de sommets isolés (suite au filtre score \>=800)

```{r eval=F}
plot(g, vertex.label=NA, vertex.size=2)
```

Ne gardons que la plus grande CC

```{r}
CC = clusters(g)
which ( CC$csize == max(CC$csize) )
```

Extraction du sous-graphe induit

```{r}
v5 = which(CC$membership == 5)
g5 = induced_subgraph(g, v5)
fr = layout.fruchterman.reingold(g5)
plot(g5, layout=fr, vertex.label=NA)
```

Diamètre

```{r}
diameter(g5)
```

C'est le plus long des plus courts chemins

```{r}
shortest_distances = distances(g5)
max(shortest_distances)
```

#### la betweenness d'un sommet ou d'une arête est le nombre de plus courts chemins passant par le sommet ou l'arête. Utiliser les fonctions `betweenness` et `edge.betweenness` pour calculer cette valeur et l'ajouter au dessin.

Centralité des sommets

```{r}
node.betweenness = betweenness(g5)
# normalization min-max
minmax = function(vec, newmin=0, newmax=1) newmin + (newmax - newmin) * (vec - min(vec)) / (max(vec) - min(vec))
nb = minmax(node.betweenness, newmin=2, newmax=15)
plot(g5, layout=fr, vertex.size = nb, vertex.label=NA)
```

Centralité des arêtes

```{r}
V(g5)$size = nb
eb = edge_betweenness(g5)
E(g5)$width = edge.width=minmax(eb, newmin=1, newmax=10)
plot(g5, layout=fr,  vertex.label=NA)
```

#### Arbre couvrant de poids minimum (`minimum.spanning.tree`)

```{r}
mst = minimum.spanning.tree(g5)
mst
```

Affichage

```{r}
plot(mst, vertex.label=NA)
```

Partitionnement et détection de communautéss
============================================

Pour partitionner le graphe en communautés, différentes méthodes sont disponibles. Vous allez utilisez la betweenness des arêtes pour effectuer le partitionnement du graphe. Cette méthode sélectionne les arêtes dont la betweenness est la plus importantes afin de former des communautés (clusters).

```{r}
com = edge.betweenness.community(g5)
par(cex=.5)
com
```

Contenu

```{r}
names(com)
```

```{r}
modularity(com)
```

Affichage du meilleur partitionnement

```{r}
plot(com, g5, layout=fr, vertex.size=6, vertex.label=NA, main=paste0("Q = ", round(modularity(com), 3)))
```

Layout
======

L'algorithme de Fruchterman-Reingold donne un résultats assez stisfaisant généralement sur les réseaux biologiques. Il modélise assez simplement un modèle physique (cf. cours). Nous allons comparer les performances sur une autre structure de grapĥe : un treillis.

```{r}
lat=graph.lattice(length=10,dim=2)
# layout
lat.FR=layout.fruchterman.reingold(lat)
lat.MDS=layout.mds(lat)
# tracé
plot(lat, layout=lat.FR, vertex.size=3, vertex.label=NA)
```

MDS

```{r}
plot(lat, layout=lat.MDS, vertex.size=3, vertex.label=NA)
```

Sur un petit graphe, les résultats sont assez comparables. Essayons sur un graphe un peu plus grand

```{r}
lat=graph.lattice(length=30,dim=2)
# layout
lat.FR=layout.fruchterman.reingold(lat)
lat.MDS=layout.mds(lat)
# tracé
plot(lat, layout=lat.FR, vertex.size=3, vertex.label=NA)
```

→ on obtient une illustration de Fruchterman-Reingold qui se retrouve "coincé"

MDS

```{r}
plot(lat, layout=lat.MDS, vertex.size=3, vertex.label=NA)
```

FR se retrouve bloqué dans une configuration stable non optimale.

Une implémentation plus efficace de MDS consiste à ne considérer qu'un petit nombre de sommets "pivots" : High Dimensional Embedding (HDE).

Principes :

-   calculer les plus courts chemins entre chaque paire de sommets (utiliser la fonction `shortest.paths(g)`)

-   déterminer les sommets pivots

    -   le premier est pris au hasard
    -   les suivants sont sélectionnés l'un après l'autre en maximisant à chaque fois leur distance aux pivots déjà sélectionnés

-   effectuer une ACP sur les coordonnées des sommets du graphes sur chacun des axes pivots

-   effectuer la projection par rapport aux 2 ou 3 composantes principales

```{r}
pivot_selection = function(g, n.pivots) {
	g.sp = shortest.paths(g)
	g.pivots = c( 1 )
	g.others = 2:length(V(g))
	while ( length(g.pivots) < n.pivots ) {
		p.new.i = 1 # index in g.others
		p.new.best.i = g.others[1] # intialize with index of first non pivot
		p.new.best.dist = 0
		# iterate over non pivots (g.others)
		while ( p.new.i <= length(g.others) ) {
			p.new = g.others[ p.new.i ]
			p.new.dist = 0
			for (p.old in g.pivots) { # add shortest path to each selected pivots
				p.new.dist = p.new.dist + g.sp[ p.old, p.new ]
			}
			# keep it if it maximizes distance
			if (p.new.dist > p.new.best.dist) {
				p.new.best.dist = p.new.dist
				p.new.best.i = p.new.i
			}
			p.new.i = p.new.i + 1
		}
		# add best pivot to selected pivots
		g.pivots = c(g.pivots, g.others[p.new.best.i])
		# remove best pivot from non pivots
		g.others = g.others[ -p.new.best.i ]
	}
	g.pivots
}

pca = function(m, normalization=FALSE) {
  m = as.matrix(m)
  if (normalization) {
    m = apply(m, 2, function(x) x-mean(x))
    var.n = function(x) sum((x-mean(x))**2)/length(x)
    m = apply(m, 2, function(x) x/sqrt(var.n(x)) )
  }
  n = nrow(m)
  COV = 1/n * t(m) %*% m
  E=eigen(COV)
  E$values/sum(E$values)
  m.proj = m %*% E$vectors
  list(
    variances = E$values**2,
    eigen.vectors = E$vectors,
    coordinates = m.proj
  )
}

layout.hde = function(g, n.pivots=50, ndim=2) {
  # adjust number of pivots
  if (length(V(g)) <= n.pivots) {
    n.pivots = round( length(V(g))/4 )
  }
  g.pivots = pivot_selection(g, n.pivots)
  # PCA
  # coordinates in g.pivots dimensions
  g.sp = shortest.paths(g)
  U = g.sp[, g.pivots]
  pca.res = pca(U, normalization = T)
  pca.res$coordinates[, 1:ndim]
}
```

Utilisation

```{r}
lat.HDE=layout.hde(lat, n.pivots=10)
plot(lat, layout=lat.HDE, vertex.size=3, vertex.label=NA)
```

# Pour aller plus loin dans la visualisation

Notamment, par exemple, pour obtenir un rendu interactif, je vous conseille le tuto https://kateto.net/network-visualization

```{r}
nodes = tibble(id=as.numeric(V(bf)), label=V(bf)$name)
nodes
```


```{r}
foo = as_long_data_frame(bf)
foo
```


```{r}
edges = tibble(from=foo$from, to=foo$to, weight=foo$weight)
edges
```

Nécessite la librairie visNetwork
```{r}
library(visNetwork)
visNetwork(nodes, edges, width="100%", height="400px")
```

Customisation
```{r}
vis.nodes <- nodes
vis.links <- edges

vis.nodes$shape  <- "dot"  
vis.nodes$shadow <- TRUE # Nodes will drop shadow
vis.nodes$borderWidth <- 2 # Node border width

vis.nodes$color.border <- "black"
vis.nodes$color.highlight.background <- "orange"
vis.nodes$color.highlight.border <- "darkred"


vis.links$width <- 1+abs(edges$weight)/2 # line width
vis.links$color <- "gray"    # line color  
vis.links$arrows <- "to" # arrows: 'from', 'to', or 'middle'
vis.links$smooth <- FALSE    # should the edges be curved?
vis.links$shadow <- FALSE    # edge shadow


visnet <- visNetwork(vis.nodes, vis.links)
visnet
```

