M1 BBS Graphes TP Librairies R
From silico.biotoul.fr
Contents |
Arbres
Il existe plusieurs librairies disponibles : ape, genoPlotR, ggtree, ... Celle choisie pour ce TP est ggtree.
ggtree
ggtree est une librairie R/Bioconductor pour la visualisation et l'annotation d'arbres phylogénétiques. Elle est basée sur ggplot2, une librairie pour faire des graphiques plus élaborés qu'avec la fonction plot de R, qu'il vous est conseillé d'expérimenter par vous-même.
Installation de la librairie
source("https://bioconductor.org/biocLite.R") biocLite("ggtree")
Au format Newick :
((B:0.2,J:0.4,(C:0.2,D:0.4,H:0.2,I:0.4)E:0.2)F:0.4,(K:0.2,L:0.4)G:0.2)A:0.2;
Remarque : read.tree au lien de read.newick (ape plus ancienne)
Premier plot avec les paramètres par défaut
```{r}
ggtree(t)
```
N'affiche que l'arbre, sans les étiquettes ni autre.
ggtree(t)+ geom_tiplab(size=3)
Et les étiquettes des noeuds internes :
nodelabels = c(t$tip.label, t$node.label) ggtree(t)+ geom_tiplab(size=3) + geom_text2(aes(subset=!isTip, label=nodelabels), hjust=-.5, size=3)
Affichage sans prendre en compte les longueurs de branche
ggtree(t, branch.length="none")
Différents algorithmes de dessin
ggtree(t) + ggtitle("default: rectangular") ggtree(t, layout="slanted") + ggtitle("slanted") ggtree(t, layout="circular") + ggtitle("circular")
Avec ou sans racine
ggtree(t, layout="unrooted") + ggtitle("unrooted layout")
Avec l'échelle pour les longueurs de branche
ggtree(t) + geom_treescale()
La même chose avec une thème différent et en stockant le plot dans une variable
p = ggtree(t) + theme_tree2() + ggtitle("rectangular tree with branch lengths and scale") p p + geom_tiplab(size=3, color="orange")
Ajout d'annotations
Pour cette partie, allez sur https://itol.embl.de/itol.cgi pour télécharger l'arbre au format Newick.
t = read.newick('iTOL.export.newick.txt') ggtree(t, layout="circular") + ggtitle("iTOL") + geom_tiplab(aes(angle=angle))+ geom_treescale()
annotations :
t$tip.label[1:10] t$node.label[1:10]
Recherche des n° de sommets correspondants aux ancêtres des bactéries, eucaryotes, et archées :
all=c(t$tip.label, t$node.label) which(all=='Bacteria') bacteria = which(all=='Bacteria')[1] which(all=='Eukaryota') eukaryota = which(all=='Eukaryota')[1] which(all=='Archaea') archaea = which(all=='Archaea')[1]
plot avec l'annotation
ggtree(t, layout="circular") + ggtitle("iTOL") + geom_tiplab(aes(angle=angle)) + geom_hilight(node=bacteria, fill="steelblue", alpha=.1) + geom_hilight(node=eukaryota, fill="pink", alpha=.1) + geom_hilight(node=archaea, fill="green", alpha=.1)
Graphes
igraph
La librairie iGraph met à disposition tout un ensemble de fonctions pour le traitement et la visualisation de graphes. Nous allons utiliser aujourd'hui son interfaçage avec R. Pour la charger :
library(igraph)
Pour charger un graphe (différents format possibles : pajek, newick, ...) :
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)
Consulter l'aide de la fonction (?read.graph) pour voir les autres formats supportés.
Pour l'afficher, il faut au préalable en effectuer le dessin (layout) :
# soit en une ligne en passant la fonction de dessin : plot(g, layout=layout.fruchterman.reingold) # soit en sauvegardant ce layout dans une variable : g.FR = layout.fruchterman.reingold(g) plot(g, layout=g.FR, vertex.size=3, vertex.label=NA)
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 : http://igraph.org/r/
- Pour obtenir la liste des sommets : V(g)
- la liste des arêtes : E(g)
- Quel est l'ordre du graphe ? Combien a-t-il d'arêtes ?
On peut assigner des étiquettes aux sommets : V(g)$name = vector_of_labels
Charger les étiquettes des sommets :
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=g.FR, vertex.size=3, vertex.label=V(g)$name)
Les paramètres de la fonction plot de igraph sont décrits dans la documentation.
Il est possible de stocker de l'information sur le graphe, les sommets et/ou les arêtes avec les fonctions dédiées :
vertex_attr(g, name="name") vertex_attr(g, name="name", index=10)
- Quel est le diamètre du graphe ?
- Lister les points d'articulation (articulation.points)
- Longueur moyenne des plus courts chemins (sans valuation) (average.path.length)
- Afficher sa représentation canonique (canonical.permutation) et la matrice d'adjacence correspondante
- Lister les cliques maximales (maximal.cliques)
- Lister les composantes connexes (clusters) avant et après suppression du point d'articulation
- Obtenir le line graph (line.graph)
- Arbre couvrant de poids minimum (minimum.spanning.tree)
Fichier dressing au format ncol : dressing.ncol et celui pour Floyd-Warshall Floyd-Warshall.ncol
- Parcours en largeur et en profondeur, et Belllman-Ford (graph.bfs, graph.dfs, distances)
- la betweenness d'une arête est le nombre de plus courts chemins passant par cette arête. Utiliser la fonction edge.betweenness pour calculer cette valeur et l'ajouter au dessin.
Partitionnement de graphe
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).
# community detection with edge-betweenness com = edge.betweenness.community(g) par(cex=.5) plot(as.dendrogram(com))
Affichage du meilleur partitionnement
plot(com, g, vertex.size=6, vertex.label=NA)
décomposition en recherchant le partitionnement ayant meilleure modularité
com$best_step=0 com$best_modularity = -Inf for (i in 1:nrow(com$merges)) { #igraph old version: ctm = community.to.membership(g, com$merges, steps=i) #igraph old version: mod = modularity(g, ctm$membership+1) # faire ctm$membership+1 si crash il y a (dépend de la version igraph, membership commence à 1 ou 0) ctm = cut_at(com, steps=i) mod = modularity(g, ctm) print(paste("step: ",i,", modularity:",mod)) if (mod>com$best_modularity) { com$best_step = i com$best_modularity = mod #igraph old version: com$membership=ctm$membership com$membership=ctm } } com$best_step ; com$best_modularity
Autre : marquage de certains sommets (au hasard la 1ère communauté)
gr1 = groups(com)[1] plot(g, mark.groups=gr1, vertex.label=NA)
librairie visNetwork
library(visNetwork)
conversion du graphe précédent
nodes = data.frame(id = as.vector(V(g)$name)) links = as.data.frame(get.edgelist(g)) colnames(links) = c('from','to')
paramètres pour l'affichage
links$smooth = F nodes$color.background <- c("red", "green", "blue", "orange", "yellow")[membership(com)]
visNetwork(nodes, links, width="1024",height="768", main="graph with visNetwork")
Pour aller plus loin : http://kateto.net/network-visualization
