M1 MABS Graphes TP Dessin et Introduction iGraph
From silico.biotoul.fr
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)
Premiers pas
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("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.sourceforge.net/doc/R/00Index.html
- 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
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.
Charger les étiquettes des sommets :
V(g)$name=as.character(read.table("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 :
get.vertex.attribute(g, name="name") get.vertex.attribute(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)
- Lister les cliques maximales (maximal.cliques)
- Lister les composantes connexes (clusters)
- Parcours en largeur et en profondeur (graph.bfs et graph.dfs)
- Obtenir le line graph (line.graph)
- Arbre couvrant de poids minimum (minimum.spanning.tree)
- Plus courts chemins : utiliser les fonctions shortest.paths et get.shortest.paths pour obtenir la longueur des plus courts chemin entre tous les sommets, ainsi que le plus court chemin entre les sommets BantA01.AAP27658.1 et HmarA01.YUFN.
- 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.
Suppression du point d'articulation, extraction du sous graphe induit et affichage avant l'analyse des composantes connexes :
vids = which( V(g) != V(g)[ articulation.points(g) ]) g1 = induced.subgraph(g, vids ) plot(g1, vertex.size=3, vertex.label=V(g)$name, vertex.label.cex=.5) clusters(g1)
Partitionnement
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).
Cette méthode ne fournit pas le meilleur partitionnement : seulement le clustering hiérarchique. Il va donc falloir déterminer où couper l'arbre pour former les clusters. Pour cela, il va falloir déterminer à quelle étape du clustering hiérarchique se situe le maximum de la mesure de qualité du partitionnement (on utilisera la modularité).
# community detection with edge-betweenness com = edge.betweenness.community(g) plot(as.dendrogram(com)) # find best partition com$best_step=0 com$best_modularity = -Inf for (i in 1:nrow(com$merges)) { ctm = community.to.membership(g, com$merges, steps=i) mod = modularity(g, ctm$membership) # faire ctm$membership+1 si crash il y a (dépend de la version igraph, membership commence à 1 ou 0) print(paste("step: ",i,", modularity:",mod)) if (mod>com$best_modularity) { com$best_step = i com$best_modularity = mod com$membership=ctm$membership } } # plot best partition plot(g, layout=g.FR, vertex.color=com$membership, vertex.size=3, vertex.label=NA, main=paste("step",com$best_step,", modularity:",com$best_modularity))
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
Sélection des pivots:
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 }
Comparer le résultat obtenu par rapport à Fruchterman-Reingold sur le graphe précédent ainsi que sur une grille régulière de 10x10 générée avec la fonction suivante :
lat=graph.lattice(length=10,dim=2) # layout lat.FR=layout.fruchterman.reingold(lat) lat.HDE=layout.hde(lat) # tracé plot(lat, layout=lat.FR, vertex.size=3, vertex.label=NA) plot(lat, layout=lat.HDE, vertex.size=3, vertex.label=NA)
Augmenter la taille de la grille (par exemple 30), et apprécier le temps pris pour le dessin par l'algorithme de FR, et par HDE.
Disposant d'une matrice de distance, il également envisageable de faire du multidimensional scaling (cf. TD de Traitement de données biologiques sur les analyses multivariées) :
# MDS Layout layout.mds = function(g, ndim=2) { g.sp = shortest.paths(g) fit <- cmdscale(g.sp,eig=TRUE, k=ndim) # k est le nombre de dimensions fit$points[,1:ndim] } lat.MDS=layout.mds(lat) plot(lat, layout=lat.MDS, vertex.size=3, vertex.label=NA)
Avec HDE et MDS, on peut aussi dessiner le graphe en 3 dimensions :
# grille lat=graph.lattice(length=10,dim=2) lat.HDE=layout.hde(lat,ndim=3) rglplot(lat, layout=lat.HDE, vertex.size=3, vertex.label=NA) # cube lat=graph.lattice(length=6,dim=3) lat.HDE=layout.hde(lat,ndim=3) rglplot(lat, layout=lat.HDE, vertex.size=3, vertex.label=NA)
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 d'ajouter les classes RootedTreeNode et RootedTree à notre bibliothèque pour manipuler les arbres :
class RootedTreeNode(Node): # TP3 # constructor def __init__(self, id, attributes = None): super(RootedTreeNode, self).__init__(id) self.attributes['parent'] = None self.attributes['children'] = [] self.attributes['distance'] = 1.0 def isLeaf(self): return len(self.attributes['children']) == 0 class RootedTree(Graph): # TP3 # constructor def __init__(self): super(Graph, self).__init__() Graph.__init__(self) self.attributes['root'] = None def parseNewick(self, text): # TP3 """ Parses the given text to create the tree. Warning, nothing is checked for correctness!! """ i=0 stack=[] self.attributes['root'] = RootedTreeNode('root') # root current = self.attributes['root'] while i < len(text): if text[i]==' ': # skip spaces i+=1 if text[i]==';': # finished ! i+=1 elif text[i]=='(': # new child N = RootedTreeNode('N') current.attributes['children'].append(N) N.attributes['parent'] = current stack.insert(0,current) current = N i+=1 elif text[i]==',': # end of label or branch length... next child incoming N = RootedTreeNode('N') parent = stack[0] # parent is at the top of the stack parent.attributes['children'].append(N) N.attributes['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.attributes['distance'] = float(nbText) else: # maybe incoming node label label='' while not text[i] in [' ',')',':',',',';']: label+=text[i] i+=1 current.id = label return self.attributes['root'] def loadNewick(self, filename): # TP3 """ Load the tree contained in the provided file. """ with open(filename) as f: text = f.readline().rstrip() T = self.parseNewick(text) return T
Affichage au format texte
Surchargez la méthode dump 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)
Pour visualiser le résultat, vous pouvez positionner les attributs rank (pour y entre 0 et 1) et depth (pour x entre 0 et 1) sur chacun des sommets et utiliser les méthodes suivantes à ajouter à la classe RootedTree :
# rendering def renderSVG(self, filename=None, width=800, height=600, circular=False): # TP3 svg = '<svg xmlns="http://www.w3.org/2000/svg" version="1.1" width="'+str(width)+'" height="'+str(height)+'">\n' if circular: svg += self.renderSVG_circular(self.attributes['root'], width, height) else: svg += self.renderSVG_rec(self.attributes['root'], width-80, height) svg += '</svg>' if filename is None: print svg else: f = open(filename, 'w') f.write(svg) f.close() def renderSVG_circular(self, node, width=800, height=600): # TP3 def svgText(x,y,text): return '<text x="'+str(x)+'" y="'+str(y)+'">'+text+"</text>\n" def svgLine(x1,y1,x2,y2): return '<line x1="'+str(x1)+'" y1="'+str(y1)+'" x2="'+str(x2)+'" y2="'+str(y2)+'" style="stroke:rgb(0,0,0);stroke-width:1"/>\n' def svgArc(x1,y1, x2,y2, rx,ry, cx,cy, angle): large=0 if angle>math.pi: large=1 return '<path d="M '+str(x1)+' '+str(y1)+ ' a '+str(rx)+' '+str(ry)+' 0 '+str(large)+' 1 '+str(x2-x1)+' '+str(y2-y1)+'" stroke="black" fill="none" stroke-width="1"/>\n' svg = '' x0 = width/2 y0 = height/2 # node label angle = node.attributes['rank'] * 2 * math.pi dx = node.attributes['depth']*math.cos(angle)*x0 dy = node.attributes['depth']*math.sin(angle)*y0 svg += svgText(x0+dx, y0+dy, node.id) # node branch dx2=0 dy2=0 if node.attributes['parent'] is not None: dx2 = node.attributes['parent'].attributes['depth']*math.cos(angle)*x0 dy2 = node.attributes['parent'].attributes['depth']*math.sin(angle)*y0 svg += svgLine(x0+dx,y0+dy, x0+dx2,y0+dy2) if not node.isLeaf(): for child in node.attributes['children']: svg += self.renderSVG_circular(child, width, height) # arc from 1st child to last child branches ## from 1st child angle1 = node.firstChild().attributes['rank']*2*math.pi dx1 = node.attributes['depth']*math.cos(angle1)*x0 dy1 = node.attributes['depth']*math.sin(angle1)*y0 ## to last child angle2 = node.lastChild().attributes['rank']*2*math.pi dx2 = node.attributes['depth']*math.cos(angle2)*x0 dy2 = node.attributes['depth']*math.sin(angle2)*y0 ## radius rx = node.attributes['depth']*x0 ry = node.attributes['depth']*y0 svg +=svgArc(x0+dx1,y0+dy1,x0+dx2,y0+dy2,rx,ry,x0,y0, angle2-angle1) return svg def renderSVG_rec(self, node, dScale=800, rScale=600, dY=20): # TP3 def svgText(x,y,text): return '<text x="'+str(x)+'" y="'+str(y)+'">'+text+"</text>\n" def svgLine(x1,y1,x2,y2): return '<line x1="'+str(x1)+'" y1="'+str(y1)+'" x2="'+str(x2)+'" y2="'+str(y2)+'" style="stroke:rgb(0,0,0);stroke-width:1"/>\n' # dScale: scale for depth (x axis) # rScale: scale for rank (y axis) # dY: shifts 20 on y for labels at y=0 (rank 0) if node.isLeaf(): # leaf return svgText(node.attributes['depth']*dScale+5, node.attributes['rank']*rScale+5+dY, node.id) else: # internal node svg = '' # draw children for child in node.attributes['children']: svg += self.renderSVG_rec(child, dScale, rScale) # branch line svg += svgLine(node.attributes['depth']*dScale, child.attributes['rank']*rScale+dY, child.attributes['depth']*dScale, child.attributes['rank']*rScale+dY) # branch length x=node.attributes['depth']*dScale + (child.attributes['depth']*dScale - node.attributes['depth']*dScale)/2 svg += svgText(x, child.attributes['rank']*rScale+dY, str(child.attributes['distance'])) # internal node label svg += svgText(node.attributes['depth']*dScale, node.attributes['rank']*rScale+5+dY, node.id) # line from 1st child to last child branches svg += svgLine(node.attributes['depth']*dScale, node.firstChild().attributes['rank']*rScale+dY, node.attributes['depth']*dScale, node.lastChild().attributes['rank']*rScale+dY) return svg
Annexes
