M1 MABS Graphes TP Dessin et Introduction iGraph
From silico.biotoul.fr
m |
m |
||
Line 1: | Line 1: | ||
+ | = Utilisation de la bibliothèque iGraph sous R = | ||
+ | |||
+ | Pour cela nous allons travailler avec le graphe ci-contre dans lequel les sommets correspondent à des gènes de différents organismes procaryotes et les liens correspondent à une relation d'isorthologie inférée entre les gènes. | ||
+ | |||
+ | [[Image:Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.FR.png|thumb|Graphe des isorthologues des protéines affines de systèmes ABC de la sous-famille 1 (import de sucres) dessiné par l'algorithme de Fruchterman-Reingold]] | ||
+ | |||
+ | Fichiers au format ''edgelist'' : | ||
+ | * [[Media:Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.gr]] (comporte les étiquettes des sommets) | ||
+ | * [[Media:Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.tgr]] (le même mais les sommets sont numérotés) | ||
+ | |||
+ | == Premier pas == | ||
+ | |||
+ | La librairie [http://igraph.sourceforge.net iGraph] met à disposition tout un ensemble de fonction 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("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.grid) | ||
+ | |||
+ | # 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.sourceforge.net/doc/R/00Index.html | ||
+ | |||
+ | * Pour obtenir la liste des sommets : <tt>V(g)</tt> | ||
+ | * la liste des arêtes : <tt>E(g)</tt> | ||
+ | |||
+ | On peut assigner des étiquettes aux sommets : <tt>V(g)$name = vector_of_labels</tt> | ||
+ | |||
+ | '''Remarque :''' la librairie étant codée en C puis interfacée avec R, cela créé au moins un petit inconvénient : les indices en C vont de 0 à n-1 alors que en R ils vont de 1 à n. | ||
+ | |||
+ | == Partitionnement == | ||
+ | Pour partitionner le graphe en communautés, différentes méthodes sont là aussi disponibles. Pour la méthode spinglass : | ||
+ | <source lang='rsplus'> | ||
+ | g.com = spinglass.community(g)$membership+1 | ||
+ | |||
+ | color = c('red', 'blue', 'yellow','black','pink','purple','orange','green','white'); | ||
+ | palette(color); | ||
+ | if ( length(unique(g.com)>length(color)) ) { | ||
+ | palette(rainbow(length(unique(g.com)))) | ||
+ | } | ||
+ | |||
+ | plot(g, layout=layout.fruchterman.reingold.grid, vertex.color=g.com, vertex.size=3, vertex.label=NA) | ||
+ | </source> | ||
+ | |||
+ | == High Dimensional Embedding == | ||
+ | |||
+ | A présent, vous allez devoir implémenter l'algorithme de dessin HDE vu en cours puisqu'il n'est pas disponible dans iGraph. | ||
+ | |||
+ | On rappelle son principe : | ||
+ | * calculer les plus courts chemins entre chaque paire de sommets (utiliser la fonction <tt>shortest.paths(g)</tt>) | ||
+ | * 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 | ||
+ | |||
+ | '''Comparer''' le résultat obtenu par rapport à Fruchterman-Reingold sur le graphe précédent ainsi que sur une grille régulière de 30x50 générée avec la fonction suivante : | ||
+ | <source lang='rsplus'> | ||
+ | # GENERATE REGULAR GRAPH HAVING A GRID STRUCTURE | ||
+ | graph.grid=function(rows = 10, cols = 20) { | ||
+ | N=rows*cols | ||
+ | M=matrix(data=0, nrow=N, ncol=N) | ||
+ | i=1 | ||
+ | while (i<=N) { | ||
+ | j=i+1 | ||
+ | while (j<=N) { | ||
+ | ri=(i-1) %/% cols | ||
+ | ci=(i-1) %% cols | ||
+ | rj=(j-1) %/% cols | ||
+ | cj=(j-1) %% cols | ||
+ | # adjacency on same row | ||
+ | if ( abs(ci-cj)==1 && ri-rj==0) { | ||
+ | M[i,j] = 1 | ||
+ | M[j,i] = 1 | ||
+ | } | ||
+ | # adjacency on same col | ||
+ | if ( abs(ri-rj)==1 && ci-cj==0) { | ||
+ | M[i,j] = 1 | ||
+ | M[j,i] = 1 | ||
+ | } | ||
+ | j=j+1 | ||
+ | } | ||
+ | i=i+1 | ||
+ | } | ||
+ | graph.adjacency(M, 'undirected') | ||
+ | } | ||
+ | </source> | ||
+ | |||
+ | |||
+ | |||
= Dessin de graphes = | = Dessin de graphes = | ||
[[Image:Sample_Tree_with_Branchlengths.png|thumb|Un arbre phylogénétique récupéré sur le site PATRIC. [[Media:Sample_Tree_with_Branchlengths.nw|Télécharger]] au format Newick.]] | [[Image:Sample_Tree_with_Branchlengths.png|thumb|Un arbre phylogénétique récupéré sur le site PATRIC. [[Media:Sample_Tree_with_Branchlengths.nw|Télécharger]] au format Newick.]] | ||
Line 112: | Line 218: | ||
* Pour les coordonnées sur l'axe des X : | * Pour les coordonnées sur l'axe des X : | ||
** Calculer la profondeur ou distance de chaque sommet à la racine (normalisée entre 0 et 1) | ** Calculer la profondeur ou distance de chaque sommet à la racine (normalisée entre 0 et 1) | ||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
Revision as of 14:38, 1 March 2013
Contents |
Utilisation de la bibliothèque iGraph sous R
Pour cela nous allons travailler avec le graphe ci-contre dans lequel les sommets correspondent à des gènes de différents organismes procaryotes et les liens correspondent à une relation d'isorthologie inférée entre les gènes.
Fichiers au format edgelist :
- Media:Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.gr (comporte les étiquettes des sommets)
- Media:Cleandb_Luca_1_S_1_1_65_Iso_Tr_1-CC1.tgr (le même mais les sommets sont numérotés)
Premier pas
La librairie iGraph met à disposition tout un ensemble de fonction 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("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.grid) # 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.sourceforge.net/doc/R/00Index.html
- Pour obtenir la liste des sommets : V(g)
- la liste des arêtes : E(g)
On peut assigner des étiquettes aux sommets : V(g)$name = vector_of_labels
Remarque : la librairie étant codée en C puis interfacée avec R, cela créé au moins un petit inconvénient : les indices en C vont de 0 à n-1 alors que en R ils vont de 1 à n.
Partitionnement
Pour partitionner le graphe en communautés, différentes méthodes sont là aussi disponibles. Pour la méthode spinglass :
g.com = spinglass.community(g)$membership+1 color = c('red', 'blue', 'yellow','black','pink','purple','orange','green','white'); palette(color); if ( length(unique(g.com)>length(color)) ) { palette(rainbow(length(unique(g.com)))) } plot(g, layout=layout.fruchterman.reingold.grid, vertex.color=g.com, vertex.size=3, vertex.label=NA)
High Dimensional Embedding
A présent, vous allez devoir implémenter l'algorithme de dessin HDE vu en cours puisqu'il n'est pas disponible dans iGraph.
On rappelle son principe :
- 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
Comparer le résultat obtenu par rapport à Fruchterman-Reingold sur le graphe précédent ainsi que sur une grille régulière de 30x50 générée avec la fonction suivante :
# GENERATE REGULAR GRAPH HAVING A GRID STRUCTURE graph.grid=function(rows = 10, cols = 20) { N=rows*cols M=matrix(data=0, nrow=N, ncol=N) i=1 while (i<=N) { j=i+1 while (j<=N) { ri=(i-1) %/% cols ci=(i-1) %% cols rj=(j-1) %/% cols cj=(j-1) %% cols # adjacency on same row if ( abs(ci-cj)==1 && ri-rj==0) { M[i,j] = 1 M[j,i] = 1 } # adjacency on same col if ( abs(ri-rj)==1 && ci-cj==0) { M[i,j] = 1 M[j,i] = 1 } j=j+1 } i=i+1 } graph.adjacency(M, 'undirected') }
Dessin de graphes
Dans cette partie, nous allons voir comment lire/écrire un arbre au format Newick. Puis, comment afficher l'arborescence avec différents algorithmes :
- textuelle comme par exemple une arborescence de fichiers
- graphique
Pour cela, nous allons utiliser les arbres ci-contre.
Chargement d'un arbre au format Newick
La première étape consiste donc à pouvoir charger l'arbre au format Newick dans un script python.
Pour cela, il semble opportun de se créer une petite bibliothèque pour manipuler les arbres : Tree.py
def createTree(): return { 'root': None } def createTreeNode(): return { 'parent': None, 'children': [], 'label': '', 'distance': 1 } def isLeaf(node): return len(node['children'])==0 def firstChild(node): return node['children'][0] def lastChild(node): return node['children'][ len(node['children'])-1 ] # ROUGH Newick PARSING: ## no syntax checking, ## unsecure, ## absolutely not fault tolerant/robust, ## not guaranteed to parse any newick compliant file!!! def parseNewick(text): T = createTree() i=0 stack=[] current=createTreeNode() # root while i < len(text): if text[i]==' ': # skip spaces i+=1 if text[i]==';': # finished ! i+=1 elif text[i]=='(': # new child N=createTreeNode() current['children'].append(N) N['parent'] = current stack.insert(0,current) current = N i+=1 elif text[i]==',': # end of label or branch length... next child incoming N=createTreeNode() parent = stack[0] # parent is at the top of the stack parent['children'].append(N) N['parent'] = parent current = N i+=1 elif text[i]==')': # new child complete current=stack.pop(0) i+=1 elif text[i]==':': # incoming branch length # parse number nbText = '' i+=1 while text[i] in ['0','1','2','3','4','5','6','7','8','9','0','.','E','e','-']: nbText+=text[i] i+=1 current['distance'] = float(nbText) else: # maybe incoming node label label='' while not text[i] in [' ',')',':',',',';']: label+=text[i] i+=1 current['label'] = label return current def loadNewick(filename): with open(filename) as f: text = f.readline().rstrip() T = parseNewick(text) return T
Affichage au format texte
Écrire une procédure pour afficher un arbre sous forme d'arborescence de fichier. Le résultat pour toyExample devrait ressembler à ce qui suit :
|-0.2- A | |-0.4- F | | |-0.2- B | | |-0.4- J | | |-0.2- E | | | |-0.2- C | | | |-0.4- D | | | |-0.2- H | | | |-0.4- I | |-0.2- G | | |-0.2- K | | |-0.4- L
Écrire une procédure pour exporter l'arbre au format Newick.
Dessin de l'arborescence sous forme non circulaire
Pour cela, nous allons utiliser une astuce qui consiste à générer du texte au format SVG (du xml) pour visualiser le résultat avec Firefox ou un autre logiciel capable d'afficher ce format. Il s'agit donc de placer les étiquettes sur un canevas ainsi que de tracer des lignes correspondant aux branches.
Pour vous faire une idée du format SVG, consulter le contenu du fichier SVG généré pour toyExample.
Principe :
- Calculer la hauteur et le nombre de feuilles. Ceci va servir à normaliser entre 0 et 1 les coordonnées (x,y) de chaque sommet et ainsi on pourra afficher sur un canevas de taille (largeur,hauteur) en multipliant : x*largeur et y*hauteur pour avoir les coordonnées dans cette taille d'image.
- Pour les coordonnées sur l'axe des Y :
- répartir les feuilles uniformément (entre 0 et 1)
- à l'aide d'une procédure récursive, une fois les coordonnées Y des descendants connues, la coordonnées Y du sommet est entre le premier et le dernier descendant.
- Pour les coordonnées sur l'axe des X :
- Calculer la profondeur ou distance de chaque sommet à la racine (normalisée entre 0 et 1)