silico.biotoul.fr
 

M1 BBS Graphes TP Recherche de communautés dans les graphes

From silico.biotoul.fr

Jump to: navigation, search

Contents

Introduction

Nous allons réaliser une analyse de génomique comparative qui s’inspire d’un travail publié par Snipen et Liland (2015). Le but est d’identifier le core et le pan génome d’un ensemble de souches bactériennes. Le core génome est l’ensemble des gènes communs à toutes les souches. Le pan génome est l’ensemble des gènes codés par au moins une des souches. Dans ce contexte un gène est en réalité une famille de gènes ou un cluster de gènes similaires. Idéalement, une famille de gènes serait une collection de gènes orthologues. Néanmoins, comme une analyse pan génomique ne produit pas toujours des groupes de gènes orthologues, il est préférable de parler de clusters de gènes. Bien que ce concept de clusters de gènes joue un rôle central dans les analyses de génomiques comparatives, il n’y a pas de consensus sur la façon de les calculer! Dans ce TP, nous allons explorer une approche basée sur la détection de communautés dans un graphe.

Snipen et Liland ont écrit un package R qui permet de répondre à cette question (micropan). Ils ont mis à disposition une étude de cas (casestudy.pdf) afin de guider les utilisateurs de ce package.

Les données

12 souches de Yersinia pestis seront analysées. Les données sont disponible sur le serveur silico dans le répertoire enseignement/m1-bioinfo/Graphe/micropan/Ypestis. Les protéines en format FASTA sont dans le répertoire data/proteins. Les comparaisons entre paires de génomes sont réalisées à l’aide de Blast+ et une annotation des domaines Pfam des protéines est réalisée avec HMMER3. Les fichiers pré-calculés sont disponibles dans les répertoires blast et pfam. Ils sont également disponibles sous forme compressée.

wget http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/Ypestis/my_micropan.R .
wget http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/Ypestis/data.tar.gz .
wget http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/Ypestis/homologs.tar.gz .

Les noms des protéines sont construits de la façon suivante:

  • Première lettre en majuscule : la première lettre du genre (Y pour Yersinia)
  • Les trois lettres suivantes en minuscule : les trois premières lettres de l’espèce (pes pour pestis)
  • La cinquième lettre majuscule est le code souche (A pour KIM)
  • 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

Mise en œuvre

Nous allons utiliser les librairies : micropan et iGraph et le fichier my_micropan.R qui contient des fonctions de micropan modifiées pour prendre en compte les préfixes de nos protéines.

my_micropan.R

source("my_micropan.R")

Pour la création d'un rapport en format HTML, utilisez la librairie knitr pour prédéfinir les formats des figures et les valeurs par défaut d'autres paramètres.

library(knitr)
knitr::opts_chunk$set(fig.align='center', fig.width=6, fig.height=4, fig.path='Figs/', echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE)
 
savepar <- par()
par(las=1)
color <- "coral"

Tableau de données

Lecture du tableau décrivant les souches utilisées dans l’analyse.

root_dir <- "<local directory>"
setwd(root_dir)
library(igraph)
options (stringsAsFactors = FALSE)
 
genome_table_file <- "data/Ypestis.csv"
  • Utilisez la fonction read.table() pour charger le fichier dans la variable genome_table. Le fichier est tabulé avec un entête (sep = "\t" , header= TRUE).

Les fichiers contenant les séquences des protéines en format FASTA sont dans le sous répertoire data/proteins. Les noms des fichiers sont dans la première colonne de genome_table.

  • Ajoutez une colonne avec le nombre de protéines dans chaque souche.

Lire chaque fichier FASTA avec la commande readLines(). Le nombre de protéines est égal au nombre de ‘>’ présents dans le fichier (utilisez grep() par exemple).

in_files = c();
for (i in 1:nrow(genome_table) ) {
	# file list
	in_files[i] <- file.path("data/proteins", genome_table$File[i])
	# number of proteins
	nb_entries <- length(grep('>', readLines(in_files[i])))
	genome_table[i, 5] <- nb_entries
}
colnames(genome_table)[5] <- c('Proteome')
kable(genome_table)
  • Tracer la distribution de la taille des protéomes (avec boxplot()).

Remarque : le package seqinr a été développé pour l'exploration, l'analyse et la visualistation des séquences biologiques. Il offre des fonctions pour faciliter la manipulation de fichers en format FASTA (ex. read.fasta()).

Orthologs 1:1

Dans un premier temps, nous allons réalisée la recherche de groupes de gènes orthologues avec des fichiers BlastP¨tabulés dans lesquels les liens d'orthologies calculés ont été ajoutés.

Inspection des liens d'orthologies

La commande **readBlastTable()** de *micropan* peut être utilisée pour relire les fichiers de sortie de blast en format tabulé.

blast.file <- 'homologs/YpesA-YpesB.par'
 
columnames <- c("Query", "Hit", "Percent.identity", "Alignment.length", "Mismatches", "Gap.openings", "Query.start", "Query.end", "Hit.start", "Hit.end", "E.value", "Bit.score", "Homology.type", "Query.length", "Hit.length", "Query.score", "Hit.score","Paire.score")
homologs <- read.table(blast.file, sep = "\t", quote = "", header = F, col.names = columnames, stringsAsFactors = F)
names(homologs)

Trois colonnes ont été ajoutées au fichier de résultats de blastp: Homology.type, Query.length, Hit.length, Query.score, Hit.score, Paire.score

La première colonne donne le type d'homologie : Iso (orthologs 1:1), Ort (orthologs), BeH (Best Hit), IsP (Best BBH paralogs) and BeP (Best paralogs). Les autres colonnes renferment les longueurs des séquences et les scores normalisés par rapport à ces longueurs.

  • Combien d'alignement(s) attendez-vous par protéines (table(homologs$Query))?
  • Tracez la distribution du nombre d'alignement par protéine.
rainbowcols <- rainbow(6, s = 0.5)
barplot(table(table(homologs$Query)), ylim=c(0, 50), col = rainbowcols, main='', xlab="nombre d'alignements", ylab="fréquence")
  • Identifiez la protéine qui a le plus d'alignements et analysez les alignements obtenus.
max_hits <- which.max(table(homologs$Query))
subtab <- homologs[homologs$Query == names(max_hits), ]
kable(subtab %>% arrange(Query.start))
  • Est-ce que toutes les protéines de A ont un alignement dans B?
cat("Non, car sur les ", genome_table[genome_table$File == "YpesA", 5], " proteines de YpesA, seules ", length(unique(homologs$Query)), " ont au moins un alignement avec YpesB.", sep="")

Extraction du graphe des orthologues 1:1

Nous allons retenir que les liens du type Homology.type == 'Iso' qui sont des orthologues 1:1.

Écriture du graphe dans le fichier "res/iso_gr".

dir.create("res", showWarnings = TRUE)
graph_file <- "res/iso_gr"
if (file.exists(graph_file)) file.remove(graph_file)
 
sum <- 0
blast_files_list <- file.path("homologs", dir("homologs"))
for( i in 1:length(blast_files_list)) {
  blast_file <- blast_files_list[i]
  cat("read blast file", blast_file, "...\n")
  blast_res <- read.table(blast_file, sep = "\t", quote = "", header = F, col.names = columnames, stringsAsFactors = F)
  blast_iso <- blast_res[blast_res$Homology.type == 'Iso', c(1,2)]
  write.table(as.data.frame(blast_iso), file=graph_file, append=TRUE, row.names=FALSE, col.names=FALSE, quote=FALSE)
  sum <- sum + nrow(blast_iso)
}
cat(sum, " rows in files\n")
  • Décrivez les différentes étapes réalisées dans ce bloc de commandes.

Traitement du graphe

Relire ce fichier comme un graphe (fonction **read.graph()** de *igraph*, avec *format="ncol"*). Affichez le nombre de sommets et d’arêtes (**vcount()** et **ecount()**).

Appliquez la fonction **simplify()** sur le graphe. Affichez le nombre de sommets et d’arêtes du graphe simplifié.

  • Comparez le nombre de sommets et d’arêtes des deux graphes et expliquez les différences observées.
  • Quel-est le but de la fonction simplify()?

Composantes connexes

Nous allons maintenant extraire les composantes connexes (CC) du graphe à l'aide de la fonction **decompose.graph()** en utilisant le *mode="strong"*.

La liste des objets graphe sera conservée dans la variable iso_cc.

iso_cc est la liste des CC du graphe.

iso_cc[[i]]
est la ième CC.

Propriétées des composantes connexes

Calculez le nombre de nœuds et le nombre d'arêtes de chaque CC et représentez graphiquement le résultat (vous pouvez utiliser **lapply()** sur la liste de CC).

iso_cc_property <- cbind(unlist(lapply(iso_cc, vcount)), unlist(lapply(iso_cc, ecount)))
colnames(iso_cc_property) <- c("vertex", "edges")
vmax <- max(iso_cc_property[,1])
 
plot(iso_cc_property[iso_cc_property[,1] <= vmax, ], las=1, xlim=c(0, vmax), pch=20, col="coral")
  • Quel est le nombre d'arêtes maximum d'un graphe en fonction du nombre de sommets?
  • Ajoutez la courbe correspondant au nombre d'arêtes maximum d'un graphe en fonction du nombre de sommets (fonction **curve()**).
xlimit <- 22
curve(x*(x-1)/2, from = 0, to = xlimit, n = 200, col="purple", add=T)
  • A quoi correspond les zones du graphe délimitées par cette courbe et quelle est la propriété des CC localisées sur la courbe?.
  • Quel est le nombre maximum de sommets pour qu'une CC soit une clique dans notre jeu de données?
  • Combien de CC de ce type sont observées dans le jeu de données?
table(iso_cc_property[,1], iso_cc_property[,2])
cross <- table(iso_cc_property[,1], iso_cc_property[,2])
cat("Il y a ", max(cross), " CC composées de 12 sommets et 66 arêtes (", max(cross)/length(iso_cc)*100, "% des CC)", sep="")


  • Représentez graphiquement une clique.
iso_cc[[5]]
par exemple.
myiso_cc <- iso_cc[[5]]
plot.igraph(myiso_cc, vertex.shape="sphere", vertex.size=15, vertex.label.cex=0.8, vertex.label.dist=0, edge.color="red")


myiso_cc <- iso_cc[[5]]
local_gene_list <- substr(V(myiso_cc)$name, 5, 5)
local_strain_list <- unique(local_gene_list)
tocolor <- unlist(lapply(local_gene_list, function(x,y) grep(x, y), local_strain_list))
mycolor = rainbow(length(local_strain_list));
palette(mycolor);
plot.igraph(myiso_cc, vertex.shape="sphere", vertex.size=15, vertex.color=tocolor, vertex.label.cex=0.8, vertex.label=local_gene_list, vertex.label.dist=0, vertex.label.color="white", edge.color="red", edge.width=E(myiso_cc)$weight)
  • Décrivez les différentes étapes réalisées dans ce bloc de commandes.
  • Combien de sommets et d'arêtes possède la CC 5?

Taux de paralogie

Nous allons é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ène, sinon, on peut avoir plus d'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 à utiliser sont : **substr()** et **unique()**. L'expression `V(graph)$name` permet d'obtenir la liste des sommets du graphe.

count_paralogs <- function(lst)
{
  ...
  return(paralogs)
}
count_paralogs <- function(lst)
{
  all_strains <- substr(lst,1,5)
  strain_list <- unique(all_strains)
  paralogs <- (length(all_strains) - length(strain_list))/length(all_strains)
  return(paralogs)
}

Distribution de la paralogie

Représentez la distribution de la paralogie dans les CC.

pdm1 <- data.frame(Vcount=numeric(), Pcount=numeric())
j <- 0
for (i in 1:length(iso_cc))
{
	cc_num    <- i
	CC        <- iso_cc[[cc_num]]
	nb_vertex <- vcount(CC)
 
	j <- j+1
	list_name <- V(CC)$name
 
  	pdm1[j, 1] <- vcount(CC)
	pdm1[j, 2] <- round(count_paralogs(list_name), 3)
}
hist(pdm1$Pcount, nclass=10)

Tracez le taux de paralogie des CC en fonction du nombre de vertices.

vmax <- max(pdm1$Vcount)
plot(pdm1[,c(1,2)], las=1, pch=20, col=color, ylim=c(0,1))
  • Quelle est le taux minimum de paralogie en fonction du nombre de protéines par CC?
  • Tracez la courbe correspondante.
vmax <- max(pdm1$Vcount)
plot(pdm1[,c(1,2)], las=1, pch=20, col=color, ylim=c(0,1))
curve((x-nrow(genome_table))/x, from=nrow(genome_table), to=vmax, n=200, col="purple", add=T)

Modularité des CC

Nous pouvons calculer la modularité d'un graphe avec la fonction **cluster_fast_greedy()** d'*igraph* (**fastgreedy.community()** suivant la version). La modularité du graphe est donnée par `max(fastgreedy$modularity)`.

pdm2 <- data.frame(Vcount=numeric(), Pcount=numeric(), Density=numeric(), Modularity=numeric())
j <- 0
for (i in 1:length(iso_cc))
{
	cc_num    <- i
	CC        <- iso_cc[[cc_num]]
	nb_vertex <- vcount(CC)
 
	j <- j+1
	list_name <- V(CC)$name
	fastgreedy <- fastgreedy.community(CC, modularity=TRUE)
 
  	pdm2[j, 1] <- vcount(CC)
	pdm2[j, 2] <- count_paralogs(list_name)
	pdm2[j, 3] <- graph.density(CC,loops=FALSE)
	pdm2[j, 4] <- max(fastgreedy$modularity)
}
kable(head(pdm2))
 
hist(pdm2$Modularity, nclass=10, col=color)
  • Tracez la modularité des CC en fonction du taux de paralogie.
plot(pdm2[,c(2,4)], las=1, pch=20, col=color)
  • Ajoutez la droite de régression *Modularity* en fonction de *Pcount*.
cor <- lm(pdm2$Modularity ~ pdm2$Pcount)
summary(cor)
plot(pdm2[,c(2,4)], las=1, pch=20, xlim=c(0,1), ylim=c(0,1), col=color)
#abline(coef(cor), xlim=c(0,1))
curve(cor$coefficients[2]*x +cor$coefficients[1], from=0.05, to=1, col="purple", add=T)
  • Commentez le résultat obtenu.

Partitionnement en communautés

Plusieurs méthodes sont disponibles dans igraph pour partitionner un graphe en communautés. Nous allons, dans un premier temps, comparer les méthodes fastgreedy, walktrap et edge betweenness. Ces trois méthodes sont utilisées de façon similaire:

  • `fastgreedy <- cluster_fast_greedy(CC, modularity=TRUE)`
  • `walktrap <- cluster_walktrap(CC, modularity=TRUE)`
  • `edge <- cluster_edge_betweenness(CC, modularity=TRUE)`

Le nombre de communautés trouvée par une méthode maximisant Q peut être obtenu avec *length(méthode)*.

  • `size[i, 1] <- length(fastgreedy)`
  • `size[i, 2] <- length(walktrap)`
  • `size[i, 3] <- length(edge)`

avec la modularité Q associée:

  • `Modularity5[j, 1] <- round(max(fastgreedy$modularity), 3)`
  • `Modularity5[j, 2] <- round(max(walktrap$modularity), 3)`
  • `Modularity5[j, 3] <- round(max(edge$modularity), 3)`

Un exemple avec fastgreedy

vmin <- 12
vmax <- max(sapply(iso_cc, vcount, simplify=TRUE))
edgecolor <- "red"
 
pdm5 <- data.frame(Vcount=numeric(), Pcount=numeric())
Modularity5 <- data.frame(fastgreedy=numeric())
Size5 <- data.frame(fastgreedy=numeric())
 
j <- 0
for (i in 2:length(iso_cc)) {
	cc_num    <- i
	CC        <- iso_cc[[cc_num]]
	nb_vertex <- vcount(CC)
	if ( nb_vertex > vmin && nb_vertex <= vmax ) {
		list_name <- V(CC)$name
		fastgreedy <- fastgreedy.community(CC, modularity=TRUE)
		if ( max(fastgreedy$modularity) > 0.45 ) {
			j <- j+1
 
			pdm5[j, 1] <- cc_num
			pdm5[j, 2] <- vcount(CC)
			pdm5[j, 3] <- count_paralogs(list_name)
 
			Modularity5[j, 1] <- round(max(fastgreedy$modularity), 3)
			Size5[j, 1] <- length(fastgreedy)
 
			cat("CC", cc_num, vcount(CC), Modularity5[j, 1], Modularity5[j, 2], Modularity5[j, 3], Modularity5[j, 4], "\n")
 
			Layout <- layout.fruchterman.reingold(CC) 
			titre <- paste(cc_num, "fastgreedy", round(Modularity5[j, 1],2), Size5[j, 1], sep=" ")
			plot(fastgreedy, CC, layout=Layout, 
			vertex.label=NA, vertex.size=3, vertex.shape="sphere",
			edge.color=edgecolor, edge.width=E(CC)$weight, main=titre)
 
			cat ("Press [enter] to continue or quit")
			line <- readline()
			if ( line == 'quit') stop()
		}
	}
}
par <- par(mar=c(5.1,4.1,4.1,2.1))
layout(1)

Ajouter les autres méthodes à ce code

spinglass

La méthode spinglass repose sur un plus grand nombre de paramètres:

spinglass <- cluster_spinglass(CC, spins=my_spins, start.temp=start_temp, stop.temp=stop_temp, cool.fact=cool_fact, update.rule=update.rule, gamma=my_gamma)`

Illustration avec la CC 1584

my_gamma <- 1
my_spins <- 25
stop_temp <- 0.8
start_temp <- 1
cool_fact <- 0.95
update.rule <- "config"
iter <- 100
 
# test 
CC <- iso_cc[[1584]]
Qiter <- array(0, c(iter))
spinglass_iter <- list() 
 
for ( k in 1:iter )  {
 	spinglass <- cluster_spinglass(CC, weights=E(CC)$weight+10, spins=my_spins,
 	start.temp=start_temp, stop.temp=stop_temp,
 	cool.fact=cool_fact, update.rule=update.rule, 
 	gamma=my_gamma)
 	Qiter[k] <- modularity(spinglass)
 	spinglass_iter[[k]] <- spinglass
}
step <- which.max(Qiter)
membership <- membership(spinglass_iter[[step]])
modularity <- round(max(Qiter), 3)
 
plot(Qiter, ylim=c(0,1), xlab="iteration", ylab="Modularity")
hist(Qiter, col=color, xlim=c(0,1))
  • Diminuez la valeur de stop.temp. Que pouvez-vous observer?

Assemblage des quatre méthodes

Nous allons fixer les paramètres suivant:

  `my_gamma <- 1`
  `my_spins <- 25`
  `stop_temp <- 0.8`
  `start_temp <- 1`
  `cool_fact <- 0.95`
  `update.rule <- "config"`

et nous allons réitérer spinglass 10 fois pour chaque CC. L'itération maximisant la modularité sera conservée.

  `iter <- 10`
  • Assemblez les quatre méthodes dans un seul bloc.

Pour chaque CC triez les méthodes en fonction des modularités obtenues. Représentez graphiquement les résultats.

Modu_stat <- apply(Modularity5, 1, function(x) rank(signif(x[1:4], digits=5), ties.method = "max"))
 
Modu_tab <- apply(Modu_stat, 1, function(x) { 
  his <- hist(x,breaks=(0:4), plot=F)
  his$count
})
barplot(Modu_tab)

Pour chaque CC triez les méthodes en fonction du nombre de communautés obtenues. Représentez graphiquement les résultats.


  • Commentez les résultats.

Relancez la boucle en supprimant le test

      	cat ("Press [enter] to continue or quit")
      	line <- readline()
      	if ( line == 'quit') stop()
  • Tracez les deux plots et commentez les résultats obtenus avec les CC de taille supérieur à 12 arêtes!

End