R Notebook - Prise en main igraph

Prise en main de la librairie R - igraph

La librairie igraph 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 :

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, …) :

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

IGRAPH 1898aee U--- 42 408 -- 
+ edges from 1898aee:
  [1] 1-- 2 1-- 3 1-- 4 1-- 5 1-- 6 1-- 7 1-- 8 1-- 9 1--10 1--11 1--12 1--13 1--14 1--15 1--16 1--17 1--18 1--19 1--20 1--21
 [21] 1--22 1--23 1--24 1--25 1--26 1--27 2-- 6 2-- 7 2--10 2--16 2--17 2--18 2--21 2--22 2--23 2--26 3-- 4 3-- 8 3-- 9 3--11
 [41] 3--12 3--14 3--15 3--16 3--19 3--20 3--24 3--25 3--27 3--28 3--29 3--30 3--31 3--32 3--33 3--34 3--35 3--36 3--37 3--38
 [61] 3--39 4-- 8 4-- 9 4--11 4--12 4--14 4--15 4--16 4--19 4--20 4--24 4--25 4--27 4--28 4--29 4--30 4--31 4--32 4--33 4--34
 [81] 4--35 4--36 4--37 4--38 4--39 5-- 6 5-- 7 5--13 5--16 5--17 5--18 5--21 5--22 5--26 6-- 7 6--10 6--13 6--14 6--16 6--17
[101] 6--18 6--21 6--22 6--23 6--26 7--10 7--13 7--14 7--16 7--17 7--18 7--21 7--22 7--23 7--26 8-- 9 8--11 8--12 8--14 8--15
[121] 8--19 8--20 8--24 8--25 8--27 8--28 8--29 8--30 8--31 8--32 8--33 8--34 8--35 8--36 8--37 8--38 8--39 9--11 9--12 9--14
[141] 9--15 9--19 9--20 9--24 9--25 9--27 9--28 9--29 9--30 9--31 9--32 9--33 9--34 9--35 9--36 9--37 9--38 9--39
+ ... omitted several edges

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

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

plot(g)

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 :

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)

V(g)

+ 42/42 vertices, from afe8fa2:
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
[38] 38 39 40 41 42

E(g)

+ 408/408 edges from afe8fa2:
  [1]  1-- 2  1-- 3  1-- 4  1-- 5  1-- 6  1-- 7  1-- 8  1-- 9  1--10  1--11  1--12  1--13  1--14  1--15  1--16
 [16]  1--17  1--18  1--19  1--20  1--21  1--22  1--23  1--24  1--25  1--26  1--27  2-- 6  2-- 7  2--10  2--16
 [31]  2--17  2--18  2--21  2--22  2--23  2--26  3-- 4  3-- 8  3-- 9  3--11  3--12  3--14  3--15  3--16  3--19
 [46]  3--20  3--24  3--25  3--27  3--28  3--29  3--30  3--31  3--32  3--33  3--34  3--35  3--36  3--37  3--38
 [61]  3--39  4-- 8  4-- 9  4--11  4--12  4--14  4--15  4--16  4--19  4--20  4--24  4--25  4--27  4--28  4--29
 [76]  4--30  4--31  4--32  4--33  4--34  4--35  4--36  4--37  4--38  4--39  5-- 6  5-- 7  5--13  5--16  5--17
 [91]  5--18  5--21  5--22  5--26  6-- 7  6--10  6--13  6--14  6--16  6--17  6--18  6--21  6--22  6--23  6--26
[106]  7--10  7--13  7--14  7--16  7--17  7--18  7--21  7--22  7--23  7--26  8-- 9  8--11  8--12  8--14  8--15
[121]  8--19  8--20  8--24  8--25  8--27  8--28  8--29  8--30  8--31  8--32  8--33  8--34  8--35  8--36  8--37
[136]  8--38  8--39  9--11  9--12  9--14  9--15  9--19  9--20  9--24  9--25  9--27  9--28  9--29  9--30  9--31
+ ... omitted several edges

Quel est l’ordre du graphe ? Combien a-t-il d’arêtes ?

ordre

length(V(g))

[1] 42

nb d’arêtes

length(E(g))

[1] 408

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

En décomposant, appeçu du contenu du fichier :

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é :

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)

Les paramètres de la fonction plot de igraph sont décrits dans la documentation.

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

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 :

vertex_attr(g, name="name")

 [1] "AfulA01.AF0890"     "AperA01.APE0061"    "BantA01.AAP27658.1" "BcerA01.AAP10714.1" "BjapA01.BLR6158"   
 [6] "BmelA02.BMEII0702"  "CkorA01.ACB07774.1" "CtepA01.BMPA"       "DradA01.DR2070"     "HbutA01.ABM80170.1"
[11] "HwalA01.RBSB"       "LinnA01.TCSA"       "MlotA01.MLL0501"    "NspeA01.ALR5361"    "PabyA01.TMPC-LIKE" 
[16] "PaerA01.AAG03536.1" "ParsA01.ABP50747.1" "PcalA01.ABO08469.1" "PfurA01.PF1695"     "PhorA01.PH1714"    
[21] "PislA01.ABL88151.1" "PyaeA01.PAE3414"    "SmarA01.ABN69244.1" "SmutA01.SMU.1121C"  "TmarA01.AAD35196.1"
[26] "TneuA01.ACB40241.1" "TonnA01.ACJ16081.1" "AperA01.APE2592"    "BsubA01.YUFN"       "CkorA01.ACB07069.1"
[31] "HbutA01.ABM81242.1" "MlotA01.MLR3148"    "OiheA01.OB3388"     "ParsA01.ABP51144.1" "PcalA01.ABO09261.1"
[36] "PislA01.ABL87654.1" "PyaeA01.PAE3121"    "SmarA01.ABN69796.1" "TneuA01.ACB40566.1" "HmarA01.YUFN"      
[41] "HsalA01.RBSB"       "HspeA01.VNG0903C"  

vertex_attr(g, name="name", index=10)

[1] "HbutA01.ABM80170.1"

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

diameter(g)

[1] 3

la même chose avec distances

distances(g)[1:6, 1:6]

                   AfulA01.AF0890 AperA01.APE0061 BantA01.AAP27658.1 BcerA01.AAP10714.1 BjapA01.BLR6158
AfulA01.AF0890                  0               1                  1                  1               1
AperA01.APE0061                 1               0                  2                  2               2
BantA01.AAP27658.1              1               2                  0                  1               2
BcerA01.AAP10714.1              1               2                  1                  0               2
BjapA01.BLR6158                 1               2                  2                  2               0
BmelA02.BMEII0702               1               1                  2                  2               1
                   BmelA02.BMEII0702
AfulA01.AF0890                     1
AperA01.APE0061                    1
BantA01.AAP27658.1                 2
BcerA01.AAP10714.1                 2
BjapA01.BLR6158                    1
BmelA02.BMEII0702                  0

max(distances(g))

[1] 3

Lister les points d’articulation (articulation.points)

articulation.points(g)

+ 1/42 vertex, named, from afe8fa2:
[1] HwalA01.RBSB

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

average.path.length(g)

[1] 1.702671

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

mean( distances(g) + diag(42)*1.702671)

[1] 1.702671

Afficher sa ’‘’représentation canonique’’’ (canonical.permutation) et la matrice d’adjacence correspondante

La permutation

canonical.permutation(g)

$labeling
 [1] 41  6 40 39  4 15 14 37 36  5 42 35  7 38 34 16  8 13 33 32 12 11  9 31 30 10 29 28 24 23 27 26 22 25 21 20 19
[38] 18 17  3  2  1

$info
$info$nof_nodes
[1] 118

$info$nof_leaf_nodes
[1] 3

$info$nof_bad_nodes
[1] 0

$info$nof_canupdates
[1] 1

$info$max_level
[1] 27

$info$group_size
[1] "6067901693952000"

La matrice d’adjacence correspondant à cette permutation

get.adjacency(permute.vertices(g, canonical.permutation(g)$labeling))

42 x 42 sparse Matrix of class "dgCMatrix"

   [[ suppressing 42 column names ‘HspeA01.VNG0903C’, ‘HsalA01.RBSB’, ‘HmarA01.YUFN’ ... ]]
   [[ suppressing 42 column names ‘HspeA01.VNG0903C’, ‘HsalA01.RBSB’, ‘HmarA01.YUFN’ ... ]]

                                                                                                      
HspeA01.VNG0903C   . 1 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
HsalA01.RBSB       1 . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
HmarA01.YUFN       1 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
BjapA01.BLR6158    . . . . . . 1 1 . 1 1 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . . . . . . 1 .
HbutA01.ABM80170.1 . . . . . 1 . 1 1 1 1 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . . . . . . 1 .
AperA01.APE0061    . . . . 1 . . 1 1 1 1 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . . . . . . 1 .
MlotA01.MLL0501    . . . 1 . . . 1 1 1 1 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . . . . . . 1 .
ParsA01.ABP50747.1 . . . 1 1 1 1 . 1 1 1 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . . . . . . 1 .
SmarA01.ABN69244.1 . . . . 1 1 1 1 . 1 1 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . . . 1 . . 1 .
TneuA01.ACB40241.1 . . . 1 1 1 1 1 1 . 1 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . . . 1 . . 1 .
PyaeA01.PAE3414    . . . 1 1 1 1 1 1 1 . 1 1 1 1 1 . . . . . . . . . . . . . . . . . . . . . 1 . . 1 .
PislA01.ABL88151.1 . . . 1 1 1 1 1 1 1 1 . 1 1 1 1 . . . . . . . . . . . . . . . . . . . . . 1 . . 1 .

 ..............................
 ........suppressing 19 rows in show(); maybe adjust 'options(max.print= *, width = *)'
 ..............................

   [[ suppressing 42 column names ‘HspeA01.VNG0903C’, ‘HsalA01.RBSB’, ‘HmarA01.YUFN’ ... ]]

                                                                                                      
PhorA01.PH1714     . . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1 1 1 1 1 1 1 1 1
PfurA01.PF1695     . . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1 1 1 1 1 1 1 1
PabyA01.TMPC-LIKE  . . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1 1 1 1 1 1 1
LinnA01.TCSA       . . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1 1 1 1 1 1
DradA01.DR2070     . . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1 1 1 1 1
CtepA01.BMPA       . . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1 1 1 1
NspeA01.ALR5361    . . . . . . . . 1 1 1 1 1 1 1 1 . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1 1 1
BcerA01.AAP10714.1 . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1 1
BantA01.AAP27658.1 . . . . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . 1 1
AfulA01.AF0890     . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 . 1
HwalA01.RBSB       1 1 1 . . . . . . . . . . . . 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 .

Lister les composantes connexes (clusters) avant et après suppression du point d’articulation

clusters(g)

$membership
    AfulA01.AF0890    AperA01.APE0061 BantA01.AAP27658.1 BcerA01.AAP10714.1    BjapA01.BLR6158  BmelA02.BMEII0702 
                 1                  1                  1                  1                  1                  1 
CkorA01.ACB07774.1       CtepA01.BMPA     DradA01.DR2070 HbutA01.ABM80170.1       HwalA01.RBSB       LinnA01.TCSA 
                 1                  1                  1                  1                  1                  1 
   MlotA01.MLL0501    NspeA01.ALR5361  PabyA01.TMPC-LIKE PaerA01.AAG03536.1 ParsA01.ABP50747.1 PcalA01.ABO08469.1 
                 1                  1                  1                  1                  1                  1 
    PfurA01.PF1695     PhorA01.PH1714 PislA01.ABL88151.1    PyaeA01.PAE3414 SmarA01.ABN69244.1  SmutA01.SMU.1121C 
                 1                  1                  1                  1                  1                  1 
TmarA01.AAD35196.1 TneuA01.ACB40241.1 TonnA01.ACJ16081.1    AperA01.APE2592       BsubA01.YUFN CkorA01.ACB07069.1 
                 1                  1                  1                  1                  1                  1 
HbutA01.ABM81242.1    MlotA01.MLR3148     OiheA01.OB3388 ParsA01.ABP51144.1 PcalA01.ABO09261.1 PislA01.ABL87654.1 
                 1                  1                  1                  1                  1                  1 
   PyaeA01.PAE3121 SmarA01.ABN69796.1 TneuA01.ACB40566.1       HmarA01.YUFN       HsalA01.RBSB   HspeA01.VNG0903C 
                 1                  1                  1                  1                  1                  1 

$csize
[1] 42

$no
[1] 1

Une seule, bien sûr.

Suppression du point d’articulation

g1 = delete_vertices(g, articulation_points(g))
plot(g1)

ou aussi

plot(g - articulation_points(g))

Composantes connexes

clusters(g1)

$membership
    AfulA01.AF0890    AperA01.APE0061 BantA01.AAP27658.1 BcerA01.AAP10714.1    BjapA01.BLR6158  BmelA02.BMEII0702 
                 1                  1                  1                  1                  1                  1 
CkorA01.ACB07774.1       CtepA01.BMPA     DradA01.DR2070 HbutA01.ABM80170.1       LinnA01.TCSA    MlotA01.MLL0501 
                 1                  1                  1                  1                  1                  1 
   NspeA01.ALR5361  PabyA01.TMPC-LIKE PaerA01.AAG03536.1 ParsA01.ABP50747.1 PcalA01.ABO08469.1     PfurA01.PF1695 
                 1                  1                  1                  1                  1                  1 
    PhorA01.PH1714 PislA01.ABL88151.1    PyaeA01.PAE3414 SmarA01.ABN69244.1  SmutA01.SMU.1121C TmarA01.AAD35196.1 
                 1                  1                  1                  1                  1                  1 
TneuA01.ACB40241.1 TonnA01.ACJ16081.1    AperA01.APE2592       BsubA01.YUFN CkorA01.ACB07069.1 HbutA01.ABM81242.1 
                 1                  1                  1                  1                  1                  1 
   MlotA01.MLR3148     OiheA01.OB3388 ParsA01.ABP51144.1 PcalA01.ABO09261.1 PislA01.ABL87654.1    PyaeA01.PAE3121 
                 1                  1                  1                  1                  1                  1 
SmarA01.ABN69796.1 TneuA01.ACB40566.1       HmarA01.YUFN       HsalA01.RBSB   HspeA01.VNG0903C 
                 1                  1                  2                  2                  2 

$csize
[1] 38  3

$no
[1] 2

Obtenir le line graph (line.graph)

lg = line.graph((g))
plot(lg, vertex.label=NA)

Il semble ressembler à g mais en plus dense.

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.

le graphe

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

IGRAPH b34321a DN-- 8 8 -- 
+ attr: name (v/c)
+ edges from b34321a (vertex names):
[1] chaussettes   ->chaussures pantalon      ->chaussures chemise       ->ceinture   cravate       ->veste     
[5] ceinture      ->veste      pantalon      ->ceinture   sous-vetements->chaussures sous-vetements->pantalon  

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

topological.sort(dressing)

+ 8/8 vertices, named, from b34321a:
[1] chaussettes    chemise        cravate        sous-vetements pantalon       chaussures     ceinture      
[8] veste         

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

Parcours en largeur

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

$root
[1] 8

$neimode
[1] "out"

$order
+ 8/8 vertices, named, from b34321a:
[1] sous-vetements chaussures     pantalon       ceinture       veste          <NA>           <NA>          
[8] <NA>          

$rank
   chaussettes     chaussures       pantalon        chemise       ceinture        cravate          veste 
           NaN              2              3            NaN              4            NaN              5 
sous-vetements 
             1 

$father
NULL

$pred
NULL

$succ
NULL

$dist
   chaussettes     chaussures       pantalon        chemise       ceinture        cravate          veste 
           NaN              1              1            NaN              2            NaN              3 
sous-vetements 
             0 

Affichage

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

Parcours en profondeur

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

$root
[1] 7

$neimode
[1] "out"

$order
+ 8/8 vertices, named, from b34321a:
[1] sous-vetements chaussures     pantalon       ceinture       veste          <NA>           <NA>          
[8] <NA>          

$order.out
NULL

$father
+ 8/8 vertices, named, from b34321a:
[1] chaussettes    chaussures     pantalon       chemise        ceinture       cravate        veste         
[8] sous-vetements

$dist
   chaussettes     chaussures       pantalon        chemise       ceinture        cravate          veste 
           NaN              1              1            NaN              2            NaN              3 
sous-vetements 
             0 

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 ?

is_dag(dressing)

[1] TRUE

Bellman-Ford

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

IGRAPH 7d9b848 DNW- 5 10 -- 
+ attr: name (v/c), weight (e/n)
+ edges from 7d9b848 (vertex names):
 [1] A->B A->E B->C B->D B->E C->B D->A D->C E->C E->D

Poids présents ?

E(bf)$weight

 [1]  6  7  5 -4  8 -2  2  7 -3  9

Affichage

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…

On a encore des arcs qui se superposent

Application de BF

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

  A B E C  D
A 0 2 7 4 -2

Affichage

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é

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

IGRAPH 6c344b1 DNW- 5 9 -- 
+ attr: name (v/c), weight (e/n)
+ edges from 6c344b1 (vertex names):
[1] A->B A->C A->E B->D B->E C->B D->A D->C E->D

Poids

E(fw)$weight

[1]  3  8 -4  1  7  4  2 -5  6

Affichage

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

Plus courtes distances

shortest.paths(fw, mode='out')

  A  B  C  E D
A 0  1 -3 -4 2
B 3  0 -4 -1 1
C 7  4  0  3 5
E 8  5  1  0 6
D 2 -1 -5 -2 0

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

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

M = matrix(c(0,0,0,0, 1,0,1,0, 0,0,1,0, 4,2,0,0), ncol=4)
M

     [,1] [,2] [,3] [,4]
[1,]    0    1    0    4
[2,]    0    0    0    2
[3,]    0    1    1    0
[4,]    0    0    0    0

Affichage

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 :

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.6     ✓ dplyr   1.0.7
✓ tidyr   1.1.4     ✓ stringr 1.4.0
✓ readr   2.1.0     ✓ forcats 0.5.1
── Conflicts ────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::as_data_frame() masks tibble::as_data_frame(), igraph::as_data_frame()
x purrr::compose()       masks igraph::compose()
x tidyr::crossing()      masks igraph::crossing()
x dplyr::filter()        masks stats::filter()
x dplyr::groups()        masks igraph::groups()
x dplyr::lag()           masks stats::lag()
x purrr::simplify()      masks igraph::simplify()

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

links.detailed = read_delim("511145.protein.links.detailed.v11.5.txt", delim=" ")

Rows: 1083186 Columns: 10
── Column specification ────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: " "
chr (2): protein1, protein2
dbl (8): neighborhood, fusion, cooccurence, coexpression, experimental, database, textmining, combined_score

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

links.detailed

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

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 :

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

Rows: 4127 Columns: 4
── Column specification ────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): #string_protein_id, preferred_name, annotation
dbl (1): protein_size

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

proteins

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

length(V(g))

[1] 4127

Arêtes

length(E(g))

[1] 5425

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

plot(g, vertex.label=NA, vertex.size=2)

Ne gardons que la plus grande CC

CC = clusters(g)
which ( CC$csize == max(CC$csize) )

[1] 5

Extraction du sous-graphe induit

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

Diamètre

diameter(g5)

[1] 17

C’est le plus long des plus courts chemins

shortest_distances = distances(g5)
max(shortest_distances)

[1] 17

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

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

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)

mst = minimum.spanning.tree(g5)
mst

IGRAPH d3b3f6d UN-- 688 687 -- 
+ attr: name (v/c), preferred_name (v/c), protein_size (v/n), annotation (v/c), size (v/n),
| coexpression (e/n), experimental (e/n), width (e/n)
+ edges from d3b3f6d (vertex names):
 [1] carA--carB ftsL--ftsQ ftsI--ftsQ ftsW--ftsQ ftsA--ftsZ aceE--aceF aceF--lpd  ftsI--mrcB tsf --pyrH fabZ--lpxA
[11] bamA--rcsF metQ--metN metI--metN dnaE--dnaQ yaaJ--prpC secA--phoA sbcC--sbcD phoB--phoR dnaJ--yajL aceE--yajL
[21] cyoE--cyoB cyoD--cyoB clpP--clpX dnaK--htpG glsA--ybaT carA--allB purK--purE fepA--entF fes --entC entF--entC
[31] entD--entE entF--entE entF--entB fepG--entB entF--entA entB--entH ftsQ--mrdA ileS--leuS yaaJ--gltA gltA--sucA
[41] sdhC--sucA sdhD--sucA sdhA--sucA sdhB--sucA sucA--sucB sucA--sucC sucA--sucD tolA--tolB tolB--pal  metN--mntR
[51] dps --clpS clpP--clpA clpS--clpA ahpC--trxB pal --lolA dnaK--rpsA mukF--mukB mukE--mukB degP--ompF leuS--asnS
[61] ftsZ--zapC secA--ompA degP--ompA skp --ompA pal --ompA hyaA--hyaB hyaB--hyaD hyaD--hyaE hyaD--hyaF yccJ--wrbA
+ ... omitted several edges

Affichage

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

com = edge.betweenness.community(g5)
par(cex=.5)
com

IGRAPH clustering edge betweenness, groups: 22, mod: 0.38
+ groups:
  $`1`
   [1] "yaaJ" "lpxD" "fabZ" "lpxA" "yahF" "prpC" "cyoE" "cyoD" "cyoC" "cyoB" "cyoA" "fdrA" "gltA" "sdhC" "sdhD"
  [16] "sdhA" "sdhB" "sucA" "sucB" "sucC" "sucD" "ybjT" "hcr"  "serS" "pyrD" "fabA" "pyrC" "plsX" "fabH" "fabD"
  [31] "acpP" "fabF" "dauA" "fabI" "pfo"  "paaE" "pqqL" "sodB" "yeaX" "preT" "preA" "yfaE" "nuoN" "nuoM" "nuoL"
  [46] "nuoK" "nuoJ" "nuoI" "nuoH" "nuoG" "nuoF" "nuoE" "nuoC" "nuoB" "nuoA" "ackA" "pta"  "fabB" "aroC" "eutD"
  [61] "maeB" "hyfC" "hyfD" "hyfH" "hmp"  "acpS" "hycF" "hycD" "sdhE" "tdcD" "yraR" "bioH" "yhjJ" "spoT" "glmU"
  [76] "atpB" "fre"  "sodA" "fpr"  "frdD" "frdC" "frdB" "frdA" "mscM" "ytfE" "qorB" "ylbE"
  
  $`2`
   [1] "dnaK" "dnaJ" "ftsA" "ftsZ" "yacG" "htpG" "zapC" "sulA" "minE" "minD" "minC" "yebG" "gyrA" "zipA" "rodZ"
  + ... omitted several groups/vertices

Contenu

names(com)

[1] "removed.edges"    "edge.betweenness" "merges"           "bridges"          "modularity"      
[6] "membership"       "names"            "vcount"           "algorithm"       

modularity(com)

[1] 0.3838229

Affichage du meilleur partitionnement

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.

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

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

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

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

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

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

nodes = tibble(id=as.numeric(V(bf)), label=V(bf)$name)
nodes

foo = as_long_data_frame(bf)
foo

edges = tibble(from=foo$from, to=foo$to, weight=foo$weight)
edges

Nécessite la librairie visNetwork

library(visNetwork)
visNetwork(nodes, edges, width="100%", height="400px")

Customisation

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