silico.biotoul.fr
 

M1 BBS Graphes TP Librairies R

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (ggtree)
m (ggtree)
(13 intermediate revisions not shown)
Line 5: Line 5:
=== ggtree ===
=== ggtree ===
-
[https://bioconductor.org/packages/release/bioc/html/ggtree.html 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 <tt>plot</tt> de R, qu'il vous est conseillé d'expérimenter par vous-même.
+
[https://bioconductor.org/packages/release/bioc/html/ggtree.html ggtree] est une librairie R/Bioconductor pour la '''visualisation''' et l''''annotation''' d'arbres phylogénétiques. Elle est basée sur <tt>ggplot2</tt>, une librairie pour faire des graphiques plus élaborés qu'avec la fonction <tt>plot</tt> de R, qu'il vous est conseillé d'expérimenter par vous-même.
-
 
+
-
Installation de la librairie
+
-
<source lang='rsplus'>
+
-
source("https://bioconductor.org/biocLite.R")
+
-
biocLite("ggtree")
+
-
</source>
+
-
 
+
[[Image:Tree_toy_example.png|thumb|Toy example. Télécharger au format [[Media:Tree_toy_example.svg|SVG]] ou [[Media:Tree_toy_example.nw|Newick]].]]
[[Image:Tree_toy_example.png|thumb|Toy example. Télécharger au format [[Media:Tree_toy_example.svg|SVG]] ou [[Media:Tree_toy_example.nw|Newick]].]]
Line 19: Line 12:
  ((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;
  ((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 :''' <tt>read.tree</tt> au lien de <tt>read.newick</tt> (ape plus ancienne)
+
'''Remarque :''' <tt>read.tree</tt> au lieu de <tt>read.newick</tt> (ape plus ancienne)
 +
Chargement de l'arbre
 +
<source lang="rsplus">
 +
t = read.newick('http://silico.biotoul.fr/site/images/f/f0/Tree_toy_example.nw')
 +
t
 +
</source>
Premier plot avec les paramètres par défaut
Premier plot avec les paramètres par défaut
-
```{r}
+
<source lang="rsplus">
ggtree(t)
ggtree(t)
-
```
+
</source>
N'affiche que l'arbre, sans les étiquettes ni autre.
N'affiche que l'arbre, sans les étiquettes ni autre.
Line 38: Line 36:
ggtree(t)+ geom_tiplab(size=3) + geom_text2(aes(subset=!isTip, label=nodelabels), hjust=-.5, size=3)  
ggtree(t)+ geom_tiplab(size=3) + geom_text2(aes(subset=!isTip, label=nodelabels), hjust=-.5, size=3)  
</source>
</source>
 +
'''Remarque :''' Il doit bien y avoir une méthode plus indiquée mais je ne l'ai pas trouvée en si peu de temps.
-
Affichage sans prendre en compte les longueurs de branche
+
Affichage sans prendre en compte les longueurs de branche (cladogramme)
<source lang="rsplus">
<source lang="rsplus">
ggtree(t, branch.length="none")
ggtree(t, branch.length="none")
Line 101: Line 100:
   geom_hilight(node=eukaryota, fill="pink", alpha=.1) +  
   geom_hilight(node=eukaryota, fill="pink", alpha=.1) +  
   geom_hilight(node=archaea, fill="green", alpha=.1)  
   geom_hilight(node=archaea, fill="green", alpha=.1)  
 +
</source>
 +
 +
Vous devriez obtenir la figure suivante :
 +
 +
[[Image:iTOL.annotated.png|800px]]
 +
 +
==== Annexe ====
 +
 +
Installation de la librairie (nécessaire le 12/12/2017 et ne fonctionne que sous windows cette année)
 +
<source lang='rsplus'>
 +
source("https://bioconductor.org/biocLite.R")
 +
biocLite("ggtree")
</source>
</source>
Line 107: Line 118:
=== igraph ===
=== igraph ===
 +
La librairie [http://igraph.sourceforge.net 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 :
 +
<source lang='rsplus'>
 +
library(igraph)
 +
</source>
 +
 +
Pour charger un graphe (différents format possibles : pajek, newick, ...) :
 +
<source lang='rsplus'>
 +
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)
 +
</source>
 +
 +
'''Consulter''' l'aide de la fonction (<tt>?read.graph</tt>) pour voir les autres formats supportés.
 +
 +
 +
Pour l'afficher, il faut au préalable en effectuer le dessin (layout) :
 +
<source lang='rsplus'>
 +
# 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)
 +
</source>
 +
 +
'''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 : <tt>V(g)</tt>
 +
* la liste des arêtes : <tt>E(g)</tt>
 +
 +
 +
 +
* Quel est l'ordre du graphe ? Combien a-t-il d'arêtes ?
 +
 +
 +
On peut assigner des étiquettes aux sommets : <tt>V(g)$name = vector_of_labels</tt>
 +
 +
Charger les étiquettes des sommets :
 +
* [[Media:Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.cod]]
 +
 +
<source lang='rsplus'>
 +
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)
 +
</source>
 +
 +
Les paramètres de la fonction plot de igraph sont décrits dans la [http://igraph.org/r/doc/plot.igraph.html documentation].
 +
 +
Il est possible de stocker de l'information sur le graphe, les sommets et/ou les arêtes avec les [http://igraph.sourceforge.net/doc/R/attributes.html fonctions dédiées] :
 +
<source lang="rsplus">
 +
vertex_attr(g, name="name")
 +
vertex_attr(g, name="name", index=10)
 +
</source>
 +
 +
 +
* Quel est le '''diamètre''' du graphe ?
 +
* Lister les '''points d'articulation''' (<tt>articulation.points</tt>)
 +
* Longueur moyenne des plus courts chemins (sans valuation) (<tt>average.path.length</tt>)
 +
 +
* Afficher sa '''représentation canonique''' (<tt>canonical.permutation</tt>) et la matrice d'adjacence correspondante
 +
* Lister les cliques maximales (<tt>maximal.cliques</tt>)
 +
* Lister les '''composantes connexes''' (<tt>clusters</tt>) '''avant et après suppression du point d'articulation'''
 +
 +
* Obtenir le line graph (<tt>line.graph</tt>)
 +
* Arbre couvrant de poids minimum (<tt>minimum.spanning.tree</tt>)
 +
 +
Fichier dressing au format ncol : [[Media:M1BBS_Graphe_dressing.ncol|dressing.ncol]], [[Media:M1BBS_Graphe_Bellman-Ford.ncol|Bellman-Ford.ncol]] et celui pour Floyd-Warshall [[Media:M1BBS_Graphe_Graphe.Floyd-Warshall.ncol|Floyd-Warshall.ncol]]
 +
 +
[[Image:igraph.dressing.png|thumb|dressing]]
 +
 +
[[Image:igraph.Bellman-Ford.png|thumb|Graphe utilisé pour l'algorithme Bellman-Ford]]
 +
 +
[[Image:igraph.Floyd-Warshall.png|thumb|Graphe utilisé pour l'algorithme Floyd-Warshall]]
 +
 +
* Effectuer le tri topologique du graphe dressing.
 +
 +
* Parcours en largeur et en profondeur, et Bellman-Ford (<tt>graph.bfs</tt>, <tt>graph.dfs</tt>, <tt>distances</tt>)
 +
 +
* 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 <tt>betweenness</tt> et  <tt>edge.betweenness</tt> pour calculer cette valeur et l'ajouter au dessin.
 +
 +
Pour la centralité des sommets : [[Image:igraph.node.betweenness.png|400px]]  Pour la centralité des arêtes : [[Image:igraph.edge.betweenness.png|400px]]
 +
 +
 +
==== 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).
 +
 +
<source lang='rsplus'>
 +
# community detection with edge-betweenness
 +
com = edge.betweenness.community(g)
 +
par(cex=.5)
 +
plot(as.dendrogram(com))
 +
</source>
 +
 +
Affichage du meilleur partitionnement
 +
<source lang='rsplus'>
 +
plot(com, g, vertex.size=6, vertex.label=NA)
 +
</source>
 +
 +
[[Image:igraph.edge.betweenness.communities.png|600px]]
 +
 +
décomposition en recherchant le partitionnement ayant meilleure modularité
 +
<source lang='rsplus'>
 +
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
 +
</source>
 +
 +
Autre : marquage de certains sommets (au hasard la 1ère communauté)
 +
<source lang='rsplus'>
 +
gr1 = groups(com)[1]
 +
plot(g, mark.groups=gr1, vertex.label=NA)
 +
</source>
 +
 +
=== librairie visNetwork ===
 +
 +
<source lang='rsplus'>
 +
library(visNetwork)
 +
</source>
 +
 +
 +
conversion du graphe précédent
 +
<source lang='rsplus'>
 +
nodes = data.frame(id = as.vector(V(g)$name))
 +
links = as.data.frame(get.edgelist(g))
 +
colnames(links) = c('from','to')
 +
</source>
 +
 +
paramètres pour l'affichage
 +
<source lang='rsplus'>
 +
links$smooth = F
 +
nodes$color.background <- c("red", "green", "blue", "orange", "yellow")[membership(com)]
 +
</source>
 +
 +
Affichage interactif
 +
<source lang='rsplus'>
 +
visNetwork(nodes, links, width="1024",height="768", main="graph with visNetwork")
 +
</source>
 +
 +
 +
Pour aller plus loin : http://kateto.net/network-visualization
-
=== visNetwork ===
 
<!--[[File:R.libs.TP.Rmd]]-->
<!--[[File:R.libs.TP.Rmd]]-->

Revision as of 09:54, 3 December 2018

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.

Toy example. Télécharger au format SVG ou Newick.

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 lieu de read.newick (ape plus ancienne)

Chargement de l'arbre

t = read.newick('http://silico.biotoul.fr/site/images/f/f0/Tree_toy_example.nw')
t

Premier plot avec les paramètres par défaut

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)

Remarque : Il doit bien y avoir une méthode plus indiquée mais je ne l'ai pas trouvée en si peu de temps.

Affichage sans prendre en compte les longueurs de branche (cladogramme)

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)

Vous devriez obtenir la figure suivante :

Annexe

Installation de la librairie (nécessaire le 12/12/2017 et ne fonctionne que sous windows cette année)

source("https://bioconductor.org/biocLite.R")
biocLite("ggtree")

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, Bellman-Ford.ncol et celui pour Floyd-Warshall Floyd-Warshall.ncol

dressing
Graphe utilisé pour l'algorithme Bellman-Ford
Graphe utilisé pour l'algorithme Floyd-Warshall
  • Effectuer le tri topologique du graphe dressing.
  • Parcours en largeur et en profondeur, et Bellman-Ford (graph.bfs, graph.dfs, 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.

Pour la centralité des sommets : Pour la centralité des arêtes :


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)]

Affichage interactif

visNetwork(nodes, links, width="1024",height="768", main="graph with visNetwork")


Pour aller plus loin : http://kateto.net/network-visualization