silico.biotoul.fr
 

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

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Partitionnement en communautés)
m (Partitionnement en communautés)
Line 350: Line 350:
Layout <- layout.fruchterman.reingold(CC)  
Layout <- layout.fruchterman.reingold(CC)  
titre <- paste(cc_num, "fastgreedy", round(Modularity5[j, 1],2), Size5[j, 1], sep=" ")
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, edge.color=edgecolor, edge.width=E(CC)$weight, main=titre)
+
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")
cat ("Press [enter] to continue or quit")
line <- readline()
line <- readline()

Revision as of 10:31, 11 December 2018

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.
  • Identifiez la protéine qui a le plus d'alignements et analysez les alignements obtenus.
  • Est-ce que toutes les protéines de A ont un alignement dans B?

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"
 
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_cci 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()**).
  • 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?
  • Combien de CC de ce type sont observées dans le jeu de données?
  • Représentez graphiquement une clique (iso_cc5 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")
  • 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)
{
  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.
  • Ajoutez la droite de régression *Modularity* en fonction de *Pcount*.
  • 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 peut être obtenu avec *length(méthode)*.

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

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

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

Nous allons fixer les paramètres suivant:

  `my_gamma <- 1`
  `my_spins <- 25`
  `stop_temp <- 0.1`
  `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`

Boucle pour spinglass:

  my_gamma <- 1
  my_spins <- 25
  stop_temp <- 0.1
  start_temp <- 1
  cool_fact <- 0.95
  update.rule <- "config"
  iter <- 10
 
  # 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)


par(mar=c(1,1,1,1))
layout(matrix(1:4, 2, 2))
vmin = 24
vmax = 50
pdm5 <- data.frame(Vcount=numeric(), Pcount=numeric())
Modularity5 <- data.frame(fastgreedy=numeric(), walktrap=numeric(), edge=numeric(), spinglass=numeric())
Size5 <- data.frame(fastgreedy=numeric(), walktrap=numeric(), edge=numeric(), spinglass=numeric())
 
my_gamma <- 1
my_spins <- 25
stop_temp <- 0.1
start_temp <- 1
cool_fact <- 0.95
update.rule <- "config"
iter <- 10
 
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 <- cluster_fast_greedy(CC, modularity=TRUE)
		fastgreedy <- fastgreedy.community(CC, modularity=TRUE)
		if ( max(fastgreedy$modularity) > 0.3 ) {
	  	j <- j+1
      #walktrap   <- cluster_walktrap(CC, modularity=TRUE)
      #edge       <- cluster_edge_betweenness(CC, modularity=TRUE)
      walktrap   <- walktrap.community(CC, modularity=TRUE)
      edge       <- edge.betweenness.community(CC, modularity=TRUE)
 
      Qiter <- array(0, c(iter))
      spinglass_iter <- list()  
      for ( k in 1:iter )
      {
        #spinglass <- cluster_spinglass(CC, , spins=my_spins,
        spinglass <- spinglass.community(CC, , 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]])
 
    	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)
  		Modularity5[j, 2] <- round(max(walktrap$modularity), 3)
  		Modularity5[j, 3] <- round(max(edge$modularity), 3)
    	Modularity5[j, 4] <- round(max(Qiter), 3)
  	  Size5[j, 1] <- length(fastgreedy)
  	  Size5[j, 2] <- length(walktrap)
	    Size5[j, 3] <- length(edge)
      Size5[j, 4] <- length(unique(membership))
 
  		cat("CC", cc_num, vcount(CC), Modularity5[j, 1], Modularity5[j, 2], Modularity5[j, 3], Modularity5[j, 4], "\n")
		  if ( max(fastgreedy$modularity) > 0.2 ) {
      	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, edge.color="white", edge.width=E(CC)$weight, main=titre)
 
      	titre <- paste("walktrap",  round(Modularity5[j, 2],2), Size5[j, 2], sep=" ")
      	plot(walktrap, CC, layout=Layout, vertex.label=NA, vertex.size=3, edge.color="white", edge.width=E(CC)$weight, main=titre)
 
      	titre <- paste("edge",  round(Modularity5[j, 3],2), Size5[j, 3], sep=" ")
      	plot(edge, CC, layout=Layout, vertex.label=NA, vertex.size=3, edge.color="white", edge.width=E(CC)$weight, main=titre)
 
      	titre <- paste("spinglass",  round(Modularity5[j, 4],2), Size5[j, 4], sep=" ")
      	plot(spinglass_iter[[step]], CC, layout=Layout, vertex.label=NA, vertex.size=3, edge.color="white", 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)

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

if ( nb_vertex > vmin && nb_vertex <= vmax ) {

et

      	cat ("Press [enter] to continue or quit")
      	line <- readline()
      	if ( line == 'quit') stop()
  • Tracez les deux plots et commentez!

BlastP scores

Comparaisons des paires de protéomes

Pour chaque paire de génomes (A, B), chaque protéine du génome A est comparée à toutes les protéines du génome B et réciproquement à l’aide de BlastP. Les alignements avec une E value < 1 sont retenus. Les résultats, sous forme tabulée, sont dans le sous répertoire blast. Les noms des fichiers sont de la forme YpesA_vs_YpesB.txt.

  • Pourquoi utiliser les protéines à la place des gènes?

La fonction blastAllAll() de micropan permet de réaliser les différents BlastP.

in_files = c();
for (i in 1:nrow(genome_table) ) {
	in_files[i] <- file.path("data/proteins", genome_table$File[i])
}
out_folder <- "blast"
myblastAllAll(in_files, out_folder, e.value=1, job=1)

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

Nous allons fixer les paramètres suivant:

  `my_gamma <- 1`
  `my_spins <- 25`
  `stop_temp <- 0.1`
  `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`

Boucle pour spinglass:

  my_gamma <- 1
  my_spins <- 25
  stop_temp <- 0.1
  start_temp <- 1
  cool_fact <- 0.95
  update.rule <- "config"
  iter <- 10
 
  # 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)

Inspection des résultats du BlastP All/All

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

readBlastTable(‘blast/YpesA_vs_YpesB.txt)
  • Calculez le nombre moyen de hits par protéines pour les deux comparaisons (YpesA_vs_YpesB et YpesB_vs_YpesA).
  • Tracez la distribution du nombre de hits par protéine pour le fichier YpesA_vs_YpesB.txt.
  • Est-ce que toutes les protéines de A ont un hit dans B?

Annotation des protéines avec les profils de PfamA

Le package micropan permet de réaliser l’annotation des protéines avec la base de données Pfam et le package HMMER. Les résultats sont dans le sous répertoire pfam. Les noms des fichiers sont du type YpesA_vs_Pfam-A.hmm.txt. Le fichier de résultat peut être chargé dans R avec la fonction readHmmer() de micropan.

  • Chargez l’annotation Pfam des protéines de YpesA.

La fonction table() permet de calculer le nombre de hits par protéines (dataAB[,1]).

  • Triez les résultats du BlastP (dataAB) par nombre de hits croissant (sort()).
  • Recherchez les annotations Pfam pour les 200 dernières protéines (tail()).
  • Calculer la fréquence de Description Pfam dans ces données.
  1. extraire la liste des noms des protéines (dataAB[,1]).
  2. extraire les annotations Pfam des protéines de cette liste (pfamA[pfamA$Query == x,]).
  3. calculer la fréquence des Description Pfam.

Pour cela, vous pouvez utiliser la fonction lapply(). Elle s’applique à une liste ou un vecteur x et elle retourne une liste de même longueur que x après avoir appliqué une fonction à chaque élément. La fonction unlist() peut être utile pour “dé-lister”" de le résultat de lapply().

  • Tracez graphiquement le résultat (barplot()).

Pour une meilleur représentation, vous pouvez basculer l’histogramme (horiz=T). Il peut être nécessaire de redéfinir les marges du graphe (mar=c(5.1,20.1,4.1,2.1))

Distances entre les paires de gènes

Le package micropan a une fonction (bDist()) qui permet de calculer les distances entre paires de gènes à partir des résultats des BlatsP. Se reporter au paragraphe BLAST distances de Microbial pangenomics in R - A case study.

Un filtre sur les E values peut être utilisé. (evalue < 1e-05). Cette distance est normalisée entre 0 et 1. La fonction file.path() renvoie la liste des fichiers présents dans un répertoire (ici blast). La fonction mybDist(blast_files, e.value = threshold) renvoie un data.frame avec les distances calculées entre paires de protéines pour toutes les comparaisons (blast_files).

Le data.frame peut être “sauvée” dans un fichier avec la fonction save() et rechargé avec la fonction load().

  • Tracez la distribution des distances (hist()).
  • Tracez le nombre de hits pour chaque protéine (hist()).

Limitez la distribution à 20 ou moins hits.

  • Quelle-est la particularité des protéines qui n’ont qu’un seul hit?

Transformation des distances en similarités

Les auteurs de micropan proposent d'utiliser une méthode de classification automatique pour découper les données en familles de gènes (fonction bClust()). Nous allons tester une méthode alternative qui repose sur le partitionnement de graphes.

  • Transformez les distances en similarités
  • Tracez la distribution des similarités (vérification!).
  • Ecrire les paires de similarités dans un fichier sous la forme d’une table (write.table()).

Graphe

Relire ce fichier comme un graphe (fonction read.graph() de igraph, avec format=“ncol”). Affichez le nombre de vertices et d’edges (vcount() et ecount()).

  • Vérifiez graphiquement les poids des arêtes (E(graph)$weight).
  • Appliquez la fonction simplify() sur le graphe.
  • Affichez le nombre de vertices et d’edges et comparez avec le graphe initial.

Composantes connexes

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

Le résultat cc est la liste des CC du graphe (cci est la ième CC).

Propriétés 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).
  • 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()).
  • Représentez graphiquement une clique (cc5 par exemple).
mycc <- cc[[5]]
local_gene_list <- substr(V(mycc)$name, 1, 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(mycc, vertex.shape="sphere", vertex.size=15, vertex.color=tocolor, vertex.label.cex=0.8, vertex.label.dist=0.8, edge.color="red", edge.width=E(mycc)$weight)

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(list=list) 
{ 
   ... 
   return(paralogs) 
}
  • Quel est le taux de paralogie dans la CC 5?

Distribution de la paralogie

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

vmin = nrow(genome_table)
vmax = 3000
pdm <- data.frame(Vcount=numeric(), Pcount=numeric())
j <- 0
for (i in 2:length(cc))
{
    cc_num    <- i
    CC        <- cc[[cc_num]]
    nb_vertex <- vcount(CC)
    if ( nb_vertex <= vmax && nb_vertex > vmin) {
        j <- j+1
        list_name <- V(CC)$name
 
        pdm[j, 1] <- vcount(cc[[i]])
        pdm[j, 2] <- round(count_paralogs(list_name), 3)
        cat("CC", cc_num, vcount(CC), pdm[j, 2], "\n")
    }
}
  • Tracez le taux de paralogie des CC en fonction du nombre de vertices.
  • Quelle est le taux minimum de paralogie en fonction du nombre de protéines par CC?
  • Tracez la courbe correspondante.
  • Que pouvez-vous dire des CC localisées sur la courbe?
hist(pdm$Vcount%%nrow(genome_table), breaks=c(0:11), col=color, xlab="Vcount%%nrow(genome_table)", xlim=c(0, nrow(genome_table)))
  • Tracez le taux de paralogie des CC en fonction du nombre de vertices.

Modularité des CC

Nous pouvons calculer la modularité d'un graphe avec la fonction cluster_fast_greedy() d'igraph. La modularité du graphe est donnée par :

max(fastgreedy$modularity)
  • Tracez la modularité des CC en fonction du taux de paralogie.
  • Tracez la droite de régression Modularity en fonction de Pcount.
cor <- lm(pdm$Modularity ~ pdm$Pcount)
summary(cor)

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)
vmin = 24
vmax = 50
pdm <- data.frame(Vcount=numeric(), Pcount=numeric())
Modularity <- data.frame(fastgreedy=numeric(), walktrap=numeric(), edge=numeric())
 
j <- 0
for (i in 2:length(cc))
{
    cc_num    <- i
    CC        <- cc[[cc_num]]
    nb_vertex <- vcount(CC)
    if ( nb_vertex > vmin && nb_vertex <= vmax ) {
        list_name <- V(CC)$name
        fastgreedy <- cluster_fast_greedy(CC, modularity=TRUE)
        if ( max(fastgreedy$modularity) > 0.3 ) {
        j <- j+1
      walktrap   <- cluster_walktrap(CC, modularity=TRUE)
      edge       <- cluster_edge_betweenness(CC, modularity=TRUE)
 
        pdm[j, 1] <- cc_num
        pdm[j, 2] <- vcount(cc[[i]])
        pdm[j, 3] <- count_paralogs(list_name)
 
        Modularity[j, 1] <- max(fastgreedy$modularity)
        Modularity[j, 2] <- max(walktrap$modularity)
        Modularity[j, 3] <- max(edge$modularity)
 
        cat("CC", cc_num, vcount(CC), Modularity[j, 1], Modularity[j, 2], Modularity[j, 3], "\n")
        }
    }
}

Comparaison des méthodes

Modularité

  • Pour chaque CC triez les méthodes en fonction des modularités obtenues.
  • Représentez graphiquement les résultats.
Modu_stat <- apply(Modularity, 1, function(x) rank(signif(x[1:3], digits=3), ties.method = "max"))
Modu_tab <- apply(Modu_stat, 1, function(x) { 
  his <- hist(x,breaks=(0:3), plot=F)
  his$count
})
barplot(Modu_tab)

Nombre de communautés

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

size[i, 1] <- length(fastgreedy)
size[i, 2] <- length(walktrap)
size[i, 3] <- length(edge)
  • Pour chaque CC triez les méthodes en fonction du nombre de communautés obtenues.
  • Représentez graphiquement les résultats.

Vous pouvez visualiser les partitionnements obtenus avec les différences méthodes.

En entête :

par <- par(mar=c(1,1,1,1))
layout(matrix(1:3,1,3))

Dans la boucle :

      if ( max(fastgreedy$modularity) > 0.3 ) {
      	Layout = layout.fruchterman.reingold(CC) 
      	titre <- paste("fastgreedy", round(Modularity[j, 1],2), Size[j, 1], sep=" ")
      	plot(fastgreedy, CC, layout=Layout, vertex.label=NA, vertex.size=3, edge.color="white", edge.width=E(CC)$weight, main=titre)
 
      	titre <- paste("walktrap",  round(Modularity[j, 2],2), Size[j, 2], sep=" ")
      	plot(walktrap, CC, layout=Layout, vertex.label=NA, vertex.size=3, edge.color="white", edge.width=E(CC)$weight, main=titre)
 
      	titre <- paste("edge",  round(Modularity[j, 3],2), Size[j, 3], sep=" ")
      	plot(edge, CC, layout=Layout, vertex.label=NA, vertex.size=3, edge.color="white", edge.width=E(CC)$weight, main=titre)
 
      	cat ("Press [enter] to continue or quite")
      	line <- readline()
      	if ( line == 'quite') stop()
      }

A la fin :

par <- par(mar=c(5.1,4.1,4.1,2.1))
layout(1)

Application

Nous allons utiliser fastgreedy pour définir les clusters de gènes. Pour chaque communauté,

mem <- membership(fastgreedy)

renvoie la classification des protéines appartenant à la CC.

classes <- unique(mem)

donne la liste des classes trouvées et

class_list <- lapply(classes, function(x, y) names(y[y==x]), y=mem)

permet d'extraire la liste des classes associées aux protéines de la classe.

La liste gene_cluster_list permet de conserver ces listes. Le data.frame gene_classes conserve les propriétés des nouvelles classes.

Les vmin et vmax sont bornés dans la phase de test!

vmin = 12
vmax = 60
ccmax = length(cc)
 
gene_classes <- data.frame(Vcount=numeric(), Strains=numeric(), Pcount=numeric())
 
gene_cluster_list <- c()
 
cluster <- 1
for (i in 2:ccmax)
{
	cc_num    <- i
	CC        <- cc[[cc_num]]
	nb_vertex <- vcount(CC)
	list_name <- V(CC)$name
	if ( nb_vertex > vmin && nb_vertex <= vmax ) {
    paralogs <- round(count_paralogs(list_name), 2)
    if ( paralogs > 0 ) {
      fastgreedy <- cluster_fast_greedy(CC, modularity=TRUE)
 
  		mem <- membership(fastgreedy)
  		classes <- unique(mem)
  		class_list <- lapply(classes, function(x, y) names(y[y==x]), y=mem)
  		for (k in classes ) {
        paralogs <- round(count_paralogs(class_list[[k]]), 2)
  		  gene_cluster_list[[cluster]] <- class_list[[k]]
      	#cat(cc_num, nb_vertex, "mod", round(max(fastgreedy$modularity), 3), "size", length(fastgreedy), "class", k, "paralogs", paralogs, "cluster", cluster, "\n")
      	all_strains <- substr(class_list[[k]],1,5)
        strain_list <- unique(all_strains)
 
      	gene_classes[cluster,1] <- length(class_list[[k]])
      	gene_classes[cluster,2] <- length(strain_list)
      	gene_classes[cluster,3] <- count_paralogs(class_list[[k]])
 
  		  cluster <- cluster+1
  		}
  	} else {
    	gene_cluster_list[[cluster]] <- list_name
 
      all_strains <- substr(list_name,1,5)
      strain_list <- unique(all_strains)
      gene_classes[cluster,1] <- length(list_name)
      gene_classes[cluster,2] <- length(strain_list)
      gene_classes[cluster,3] <- count_paralogs(list_name)
 
    	#cat(cc_num, nb_vertex, "paralogs", paralogs, "cluster", cluster, "\n")
  	  cluster <- cluster + 1
  	}
	}
}
  • Tracez la distribution des taux de paralogie.
  • Tracez la distribution du nombre de souches présents dans les clusters de gènes.

Liens d'orthologies

Nous allons rependre l'analyse réalisée avec les résultats des BlastP All/All, mais avec des fichiers BlastP¨tabulés dans lesquels les orthologs ont été calculés (répertoire homologs).

*[http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/Ypestis/homologs.tar.gz homologs.tar.gz]

Les résultats, par paires de génomes, sont dans le format BlastP tabulé avec 6 colonnes en plus: le type de lien d'homologie, la taille des protéines Query et Hit, les scores normalisés des Query, Hit et des deux.

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)

Extraction du graphe des orthologues 1:1

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

  • Relire les fichiers et ré-écrire les paires de gènes orthologues 1:1 dans le fichier "res/iso_gr".

Traitement du graphe

  • Relire ce fichier comme un graphe.
  • Affichez le nombre de vertices et d'edges.
  • Appliquez la fonction **simplify()** sur le graphe.
  • Affichez le nombre de vertices et d'edges et comparez avec le graphe initial.
  • Extraire les composantes connexes (CC).
  • Calculez le nombre de nœuds et le nombre d'arêtes de chaque CC et représentez graphiquement le résultat.
  • Ajoutez la courbe correspondant au nombre d'arêtes maximum d'un graphe en fonction du nombre de sommets.
  • Représentez la distribution de la paralogie dans les CC.
  • Tracez la modularité des CC en fonction du taux de paralogie.
  • Tracez la droite de régression *Modularity* en fonction de *Pcount*.

Partitionnement en communautés

En plus des méthodes déjà vues, nous allons utiliser la méthode spinglass. Elle 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)`

Nous allons fixer les paramètres suivant:

  `my_gamma <- 1`
  `my_spins <- 25`
  `stop_temp <- 0.1`
  `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`

Boucle pour spinglass:

  my_gamma <- 1
  my_spins <- 25
  stop_temp <- 0.1
  start_temp <- 1
  cool_fact <- 0.95
  update.rule <- "config"
  iter <- 10
 
  # 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)