M1 MABS Graphes TP Recherche de communautés dans les graphes
From silico.biotoul.fr
Contents |
Contexte biologique
La grande disponibilité de données issues de séquençage de génomes complets rend les approches basées sur la génomique comparative particulièrement attractives pour aborder les questions relatives à l’adaptation des bactéries à leur environnement. Dans ce cadre, les systèmes codés par de grandes familles de gènes paralogues, comme les transporteurs ABC, sont ceux qui présentent le plus grand potentiel. En effet, ces systèmes sont impliqués dans l'import ou l'export de molécules diverses à travers la membrane cytoplasmique. Ils sont ubiquitaires chez les archaea et les bactéries, ils seraient donc apparus très tôt dans l’évolution. Cependant, ils peuvent avoir des distributions différentes même entre génomes apparentés suggérant des vagues successives de gains/pertes de gènes.
Les transporteurs ABC peuvent être classés en une trentaine de sous-familles aux quelles nous pouvons associer des types de substrats de nature spécifiques (sucres, acides aminées, peptides, ions, ferrichromes...). Cependant, nous disposons de très peu d'informations sur la nature exacte de la molécule transportée par un transporteur. Or cette information est nécessaire si l'on veut pouvoir étudier la contribution de ces transporteurs à l'adaptation des bactéries à leur environnement. Une approche bioinformatique pour tenter de répondre à cette question consiste à regrouper les différents systèmes en groupe de gènes orthologues. Sous l'hypothèse que ces gènes orthologues ont conservé la même fonction, et si le substrat transporté par le produit d'un des gènes d'un groupe a été identifié expérimentalement, alors on peut propager cette annotation à tous les autres membres du groupe. Une autre approche consiste a exploiter le contexte génomique. En effet, une forte association entre les gènes appartenant au même groupe d'orthologues et un autre gène peut révéler un lien fonctionnel entre le transporteur et le produit de ce gène, comme par exemple, si le gène code pour une enzyme, alors on peut penser que le substrat de cette enzyme a de forte chance d'être la molécule prise en charge par le transporteur ABC.
Cette approche nécessite comme préalable l'identification de groupes de gènes orthologues. Si cela peut paraitre assez facile à réaliser dans le cas de gènes présentant pas ou peu de duplication dans les génomes, la tache est beaucoup plus difficile à effectuer avec des familles de gènes fortement paralogues. Néanmoins, dans tous les cas, il est nécessaire d’établir de façon précise les liens de parenté entre les séquences. Cela peu se faire à l’aide de représentation arborée ou à l’aide de graphes. Les arbres fournissent une bonne modélisation des trajectoires évolutives des séquences mais sont très sensibles à la qualité des données utilisées pour les construire. A l’inverse, les graphes sont plus difficiles à manipuler et à interpréter, mais présentent l’avantage d’être beaucoup plus robustes et de pouvoir inclure un très grand nombre de données. Généralement, les gènes (ou les protéines) constituent les sommets du graphe. Ils sont reliés par des arcs pouvant traduire une relation évolutive entre les gènes. Ces liens sont qualitatifs (homologie, orthologie, paralogie). Si les données sont suffisamment dissemblables, le graphe peut se décomposer en un ensemble de composantes connexes (ensemble de sommets connectés par au moins un arc), on a alors une classification de type COG. Sinon, les groupes de gènes orthologues ne sont pas déconnectés dans le graphe, mais correspondent à des régions de plus forte densité (communautés). L’identification de groupes de gènes orthologues revient alors à l’identification de communautés dans un graphe.
Jeu de données
A titre d’exemple, nous allons analyser une sous familles de transporteurs ABC, sous famille 1, dont les membres sont impliqués dans l’import de sucre. A l’aide de la base de données ABCdb, nous avons sélectionné une liste de transporteurs de la sous famille 1. Pour chaque transporteur, nous avons extrait les protéines codées par les gènes, localisés sur le chromosome, en amont et en aval des gènes codant pour la transporteur (4 gènes de par et d’autre). Pour toutes ces protéines, nous disposons des pairs de liens d’isorthologie.
Remarques:
- Les liens d’homologies n’ont de sens qu’au niveau des gènes, mais ils ont été établis au niveau des protéines en raison de leur meilleure conservation.
- Les gènes orthologues a et b, appartenant aux génomes A et B, sont isorthologues s'il n'existe pas de paires de gènes paralogues (a, y) dans A et (x, b) dans B avec une distance phylogénétique inférieure à celle de la paire (a, b). La liste de ces pairs d’isorthologues constitue la liste d’arête d’un graphe dont les sommets sont les protéines.
- Dans la base de données ABCdb, les noms des protéines sont construits de la façon suivante:
- Première lettre en majuscule : la première lettre du genre (E pour Escherichia)
- Les trois lettres suivantes en minuscule : les trois premières lettres de l’espèce (col pour coli)
- La cinquième lettre majuscule est le code souche (A pour K12)
- Les deux chiffres suivant se réfèrent au numéro du chromosome (01 pour le chromosome 1)
- Un point est utilisé comme séparateur, il est suivit du nom de la protéine
Exemple : EcolA01.RBSA
- Les annotations des gènes/protéines sont compilées dans un fichier, dont les colonnes correspondent à :
- le nom de la protéine
- annotation COG
- le non du gène (protéine) d'ancrage, utilisé comme référence pour la région chromosomique contenant les gènes codant pour un transporteur ABC,
- le nom du système ABC (le nom d'un des partenaires),
- le domaine (NBD, MSD, SBP)
- sous famille
Fichiers
Mise en œuvre
Nous allons utiliser la librairie : igraphe
Prise en main du graphe
La fonction read.graph permet de lire la liste des arêtes à partir d’un fichier.
graph <- read.graph(graph_file, format = "ncol", directed = FALSE);
Nous allons utiliser les fonctions vcount et ecount pour compter le nombre de sommets et d’arêtes du graphe.
Avant de poursuivre, il est nécessaire d’applique la fonction simplify() sur notre graphe.
*Quelle est l’effet de cette fonction ? *Discutez cette observation en la ramenant à la façon dont a été construit le graphe.
Taux de paralogie
On s’attend à ce qu’un assemblage de gènes isorthologues ne comporte pas de gènes paralogues. Nous allons donc écrire une fonction qui calcule le taux de gènes paralogues dans une liste de gènes. Pour cela nous allons compter le nombre de gènes présent dans la liste et le comparer au nombre de souches. S’il n’y a pas de paralogue, chaque souche est représentée par un seul gènes, sinon, on peut avoir plus de un gène par souche. Pour obtenir le taux de paralogues, on divise la différence entre nombre de gènes et nombre de souches par le nombre de gènes.
Les fonctions a utiliser sont : substr() et unique(). L’expression V(graph)$name permet d’obtenir la liste des sommets du graphe.
count_paralogs <-function(list=list) { ... return(paralogs) }
Utilisez cette fonction pour calculer le taux de paralogues dans le graphe.
*Quelle peut être l’origine de ce taux élevé da paralogie ?
Représentation graphique
Nous allons observer graphiquement le graphe. Pour faciliter les comparaisons entre différentes « colorations » du graphe, nous allons fixer la disposition (le layout) du graphe à l’aide de l’expression suivante :
globalGraphLayout = layout.fruchterman.reingold(graph); plot(graph, layout=globalGraphLayout, vertex.label=NA, vertex.size=3);
*Décrire le graphe obtenu.
Annotations
Nous allons utiliser les annotations en domaines contenu dans le fichier annot_file pour colorier le graphe.
Lecture du fichier :
annot_file <- 'chemin du fichier/Annotations.txt'; annotation <- read.table(annot_file, header=TRUE, row.names=1);
Annotation en domaines
Les colonnes du fichier sont : 'protein name', 'anchor', 'cog', 'assembly', 'domain', 'subfamily'. La difficulté est que certaines protéines n’ont pas de domaines annotés (autre que partenaire d'un système ABC) dans le fichier.
domains = annotation[V(graph)$name, 4];
Pensez à associer une couleur à chaque domaine (color = rainbow(x)).
plot(graph, layout=globalGraphLayout, vertex.color=domains, vertex.label=NA, vertex.size=3); legend("topright", legend=levels(domains), cex=0.7, fill=color);
Annotation en sous familles
ssfam = annotation[V(graph)$name, 5];
* Commentez le résultat obtenu ?
Extraction des composantes connexes
Nous allons maintenant extraire les composantes connexes (CC) du graphe à l’aide de la fonction decompose.graph() en utilisant le mode = "strong".
Nous ne conserverons que les CC contenant au moins 5 sommets.
* Combien de composantes connexes avez-vous obtenu ?
Propriétés des composantes connexes
Les CC apparaissent hétérogènes au niveau de leur taille mais aussi au niveau de leur densité en arêtes.
* Quelle est le nombre d’arêtes maximum possible dans un graphe de n sommets ?
Vous pouvez établir une fonction donnant le nombre d’arêtes maximum d’un graphe en fonction du nombre de sommets. Nous allons donc étudier la distribution du nombre d’arêtes en fonction du nombre de sommets de nos CC et nous allons tracer la courbe de la fonction d’arêtes maximum.
On peut par exemple placer ces informations dans un tableau :
cc_info <- array(0, c(nb_cc, 2)); for ( i in 1:nb_cc ) { cc_info[i,1] <- vcount(cc[[i]]) cc_info[i,2] <- ecount(cc[[i]]) }
Puis calculer la courbe de la fonction:
maxx <- max(cc_info[,1]) maxy <- max(cc_info[,2]) courbe <- c(); i <- 0; j <- 0; while(j < maxy && i < maxx ) { i <- i + 1; j <- i*(i-1)/2; courbe[i] <- j; }
On pourra utiliser la fonction points() pour ajouter la courbe à la distribution des CC.
* Commentez le résultat obtenu.
Pour compléter l’analyse, nous allons ajouter le taux de paralogies de chaque CC.
paralogs <- array(0, c(nb_cc, 2)); for ( i in 1:nb_cc ) { list_name <- V(cc[[i]])$name; paralogs[i, 1] <- vcount(cc[[i]]); paralogs[i, 2] <- count_paralogs(list_name); }
Pour superposer le résultat sur le graphique précédant, nous pouvons utiliser :
par(new=T) plot(paralogs, pch=18, col='red', axes=F, xlim=xlim, ylim=zlim, ,xlab= '', ylab=''); axis(side=4) mtext('paralogs', side=4, line=2, col='red')
Il peut être nécessaire d’augmenter la place disponible à droite du graphique pour afficher la légende (fonction mar() à utiliser dans un par()).
* Quelle propriété du graphe semble liée à la distribution du taux de paralogie ?
Taux de paralogie et modularité
Nous allons étudier la distribution du taux de paralogie en fonction de la modularité des CC. Pour cela nous allons utiliser fastgreedy.community()qui permet de faire le découpage en communautés d’un graphe.
Remarque: il faut utiliser get.adjlist() pour passer des noms des protéines à la liste des sommets.
modularity <- c(); for ( i in 1:nb_cc ) { g <- graph.adjlist(get.adjlist(cc[[i]]), directed = FALSE) com <- fastgreedy.community(g, merges=TRUE, modularity=TRUE); modularity[i] <- max(com$modularity); }
* Que pouvez-vous en conclure?
Décomposition en communautés
Nous allons représenter graphiquement le résultat du découpage en communautés réalisée par fastgreedy, uniquement pour le CC qui ont un taux de paralogie supérieur à 0.
* Pouvez-vous justifier ce choix ?
Pour la suite nous allons conserver les graphes modifiés pour être compatibles avec les méthodes de communautés et les dispositions pour faciliter les comparaisons entre les méthodes (utilisation de listes).
G_list <- list(); # liste des graphes GLayout_list <- list(); # liste des layout for ( i in 1:nb_cc ) { G_list[[i]] <- graph.adjlist(get.adjlist(cc[[i]]), directed = FALSE) GLayout_list[[i]] <- layout.fruchterman.reingold(G_list[[i]]); }
