silico.biotoul.fr
 

M1 BBS Graphes TP Recherche de communautés dans les graphes2.0

From silico.biotoul.fr

(Difference between revisions)
Jump to: navigation, search
m (Les données)
m (Propriétés des composantes connexes)
 
(36 intermediate revisions not shown)
Line 56: Line 56:
=== Chargement des fichiers ===
=== Chargement des fichiers ===
-
Vous pouvez effectuer les téléchargements des fichiers dans une Chunk bash dans Rstudio.
+
Téléchargement et décompression des fichiers dans un shell.
 +
<source lang='bash'>
 +
rootdir=$HOME/Graphe/micropan/Ypestis
 +
echo $rootdir
 +
mkdir -p $rootdir
 +
cd $rootdir
 +
pwd
-
<source lang='rsplus'>
 
if [ ! -f my_micropan.R ];
if [ ! -f my_micropan.R ];
then
then
Line 81: Line 86:
</source>
</source>
-
Nous allons utiliser les librairies : '''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.  
+
Nous allons nous placer dans le répertoire de travail.
<source lang='rsplus'>
<source lang='rsplus'>
 +
setwd(rootdir)
getwd()
getwd()
list.files()
list.files()
 +
source("my_micropan.R")
source("my_micropan.R")
</source>
</source>
 +
 +
Nous allons utiliser les librairies : '''igraph''', '''tidyverse''' 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 dans les résultats des blast.
<source lang='rsplus'>
<source lang='rsplus'>
-
library(dplyr)
+
library(tidyverse)
library(igraph)
library(igraph)
 +
library(ggvenn)
options (stringsAsFactors = FALSE)
options (stringsAsFactors = FALSE)
-
savepar <- par()
+
savepar <- par() # save Graphical Parameters
par(las=1)
par(las=1)
color <- "coral"
color <- "coral"
Line 120: Line 130:
++
++
 +
<source lang='rsplus'>
<source lang='rsplus'>
in_files <- c()
in_files <- c()
Line 132: Line 143:
kable(genome_table)
kable(genome_table)
</source>
</source>
-
 
*Tracer la distribution de la taille des protéomes (avec '''boxplot()''').
*Tracer la distribution de la taille des protéomes (avec '''boxplot()''').
++
++
 +
<source lang='rsplus'>
 +
boxplot(genome_table$Proteome, ylab='Proteome size', col='coral', las=1)
 +
</source>
 +
<!--
<source lang='rsplus'>
<source lang='rsplus'>
rainbowcols <- rainbow(6, s = 0.5)
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")
barplot(table(table(homologs$Query)), ylim=c(0, 50), col = rainbowcols, main='', xlab="nombre d'alignements", ylab="fréquence")
</source>
</source>
 +
-->
Remarque : le package '''seqinr''' a été développé pour l'exploration, l'analyse et la visualisation des séquences biologiques. Il offre des fonctions pour faciliter la manipulation de fichiers en format FASTA (ex. '''read.fasta()''').
Remarque : le package '''seqinr''' a été développé pour l'exploration, l'analyse et la visualisation des séquences biologiques. Il offre des fonctions pour faciliter la manipulation de fichiers en format FASTA (ex. '''read.fasta()''').
=Orthologs 1:1=
=Orthologs 1:1=
-
Dans un premier temps, nous allons réaliser 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.
+
Dans un premier temps, nous allons réaliser 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==
==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é.
+
La commande '''readBlastTable()''' de ''micropan'' peut être utilisée pour relire les fichiers de sortie de blast en format tabulé.
<source lang='rsplus'>
<source lang='rsplus'>
Line 162: Line 177:
*Combien d'alignement(s) attendez-vous par protéines?
*Combien d'alignement(s) attendez-vous par protéines?
-
''table(homologs$Query)'' vous donne le nombre de fois qu'une protéine est présente dans le fichier.
+
Vous pouvez regrouper les lignes par requête (''group_by''(Query)) et compter le nombre de cibles (''summarise''(Hits=n())) pour chaque requête.  
-
*Tracez la distribution du nombre d'alignement par protéine.
+
 +
++
 +
<source lang='rsplus'>
 +
nbaliprot <- homologs %>%
 +
  group_by(Query) %>%
 +
  summarise(Hits=n()) %>%
 +
  arrange(desc(Hits))
 +
</source>
 +
*Tracez la distribution du nombre d'alignement par protéine avec ''geom_histogram(binwidth=1, col=color, fill=color)''.
 +
 +
++
 +
<source lang='rsplus'>
 +
ggplot(nbaliprot, aes(x=Hits)) +
 +
  geom_histogram(binwidth=1, col=color, fill=color) +
 +
  scale_y_log10() +
 +
  theme_bw()
 +
</source>
 +
<!--
 +
ancienne version
<source lang='rsplus'>
<source lang='rsplus'>
rainbowcols <- rainbow(6, s = 0.5)
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")
barplot(table(table(homologs$Query)), ylim=c(0, 50), col = rainbowcols, main='', xlab="nombre d'alignements", ylab="fréquence")
</source>
</source>
 +
-->
 +
* Identifiez la protéine qui a le plus d'alignements et analysez les alignements obtenus.
-
*Identifiez la protéine qui a le plus d'alignements et analysez les alignements obtenus.
+
<!--
-
++
+
old version
<source lang='rsplus'>
<source lang='rsplus'>
max_hits <- which.max(table(homologs$Query))
max_hits <- which.max(table(homologs$Query))
subtab <- homologs[homologs$Query == names(max_hits), ]
subtab <- homologs[homologs$Query == names(max_hits), ]
kable(subtab %>% arrange(Query.start))
kable(subtab %>% arrange(Query.start))
 +
</source>
 +
-->
 +
* sélectionner la protéine qui a le plus de hits (Hits==max(Hits))
 +
* avec son nom, filtrer la table ''homologs'' (Query==max_hits$Query)
 +
* sélectionner les colonnes: Query, Hit, Query.start, Query.end, Hit.start et Hit.end
 +
* éditer cette table avec kable.
 +
 +
++
 +
<source lang='rsplus'>
 +
max_hits <- nbaliprot %>%
 +
  filter(Hits==max(Hits)) %>%
 +
  select(Query)
 +
max_hits$Query
 +
 +
subtab <- homologs %>%
 +
  filter(Query==max_hits$Query)
 +
head(subtab)
 +
 +
compact <- subtab %>%
 +
  select(Query, Hit, Query.start, Query.end, Hit.start, Hit.end) %>%
 +
  arrange(Query.start)
 +
 +
kable(compact)
</source>
</source>
*Est-ce que toutes les protéines de A ont un alignement dans B?
*Est-ce que toutes les protéines de A ont un alignement dans B?
 +
++
++
 +
<source lang='rsplus'>
 +
YpesA <- genome_table %>%
 +
  filter(File == "YpesA")
 +
 +
YpesA.Protein.Hits <- nrow(homologs %>% distinct(Query))
 +
 +
cat("Sur les ", YpesA$Proteome, " proteines de YpesA, seules ", YpesA.Protein.Hits, " ont au moins un alignement avec YpesB.", sep="")
 +
</source>
 +
<!--
 +
old
<source lang='rsplus'>
<source lang='rsplus'>
cat("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="")
cat("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="")
</source>
</source>
 +
-->
==Extraction du graphe des orthologues 1:1==
==Extraction du graphe des orthologues 1:1==
Line 194: Line 263:
sum <- 0
sum <- 0
 +
dir.create("res", showWarnings = TRUE)
 +
graph_file <- "res/iso_gr"
 +
if (file.exists(graph_file)) file.remove(graph_file)
 +
blast_files_list <- file.path("homologs", dir("homologs"))
blast_files_list <- file.path("homologs", dir("homologs"))
for( i in 1:length(blast_files_list)) {
for( i in 1:length(blast_files_list)) {
Line 199: Line 272:
   cat("read blast file", blast_file, "...\n")
   cat("read blast file", blast_file, "...\n")
   blast_res <- read.table(blast_file, sep = "\t", quote = "", header = F, col.names = columnames, stringsAsFactors = F)
   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)]
+
   blast_iso <- blast_res %>%
 +
    filter(Homology.type == 'Iso') %>%
 +
    select(Query, Hit)
   write.table(as.data.frame(blast_iso), file=graph_file, append=TRUE, row.names=FALSE, col.names=FALSE, quote=FALSE)
   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)
   sum <- sum + nrow(blast_iso)
Line 218: Line 293:
Affichez le nombre de sommets et d’arêtes (''vcount()'' et ''ecount()'').  
Affichez le nombre de sommets et d’arêtes (''vcount()'' et ''ecount()'').  
-
++
+
* Le graphe est composé de `r vcount(gr)` sommets et de `r ecount(gr)` arêtes.
-
Le graphe est composé de `r vcount(gr)` sommets et de `r ecount(gr)` arêtes.
+
Appliquez la fonction ''simplify()'' sur le graphe.
Appliquez la fonction ''simplify()'' sur le graphe.
Line 230: Line 304:
Affichez le nombre de sommets et d’arêtes du graphe simplifié.
Affichez le nombre de sommets et d’arêtes du graphe simplifié.
-
++
+
* Le graphe simplifié est composé de `r vcount(gr_iso)` sommets et de `r ecount(gr_iso)` arêtes.
-
Le graphe simplifié est composé de `r vcount(gr_iso)` sommets et de `r ecount(gr_iso)` arêtes.
+
*Comparez le nombre de sommets et d’arêtes des deux graphes et expliquez les différences observées.
*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()?
*Quel-est le but de la fonction simplify()?
 +
==Composantes connexes==
==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 allons maintenant extraire les composantes connexes (CC) du graphe à l'aide de la fonction ''decompose.graph()'' en utilisant le ''mode="strong"''.
Line 245: Line 319:
La liste des objets graphe sera conservée dans la variable ''iso_cc''.
La liste des objets graphe sera conservée dans la variable ''iso_cc''.
-
''iso_cc'' est la liste des CC du graphe.  
+
''iso_cc'' est la '''liste''' des CC du graphe.  
-
<source lang='rsplus'>iso_cc[[i]]</source> est la ième CC.
+
<source lang='rsplus'>
 +
iso_cc[[i]]
 +
</source> est la ième CC.
==Propriétés des composantes connexes==
==Propriétés des composantes connexes==
Line 252: Line 328:
Calculez le nombre de nœuds et le nombre d'arêtes de chaque CC et représentez graphiquement le résultat.  
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 pour calculer le nombre de sommets et le nombre d'arêtes de chaque CC.
+
Vous pouvez utiliser ''unlist(lapply())'' ou ''sapply()'' sur la liste de CC pour calculer le nombre de sommets et le nombre d'arêtes de chaque CC.
++
++
Line 260: Line 336:
vmax <- max(iso_cc_property[,1])
vmax <- max(iso_cc_property[,1])
 +
iso_cc_property <- cbind(sapply(iso_cc, vcount), sapply(iso_cc, ecount))
 +
colnames(iso_cc_property) <- c("vertex", "edges")
 +
</source>
 +
 +
Tracer le nombre d'arêtes en fonction du nombre de sommets.
 +
 +
<source lang='rsplus'>
 +
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")
plot(iso_cc_property[iso_cc_property[,1] <= vmax, ], las=1, xlim=c(0, vmax), pch=20, col="coral")
 +
</source>
 +
ou
 +
<source lang='rsplus'>
 +
iso_cc_property.df <- as.data.frame(iso_cc_property)
 +
Vermax <- max(iso_cc_property.df$vertex)
 +
Edgmax <- max(iso_cc_property.df$edges)
 +
 +
ggplot(iso_cc_property.df, aes(x=vertex, y=edges)) +
 +
  geom_point(
 +
    aes(color = color, fill = color),
 +
    size = 1.5, alpha = 0.5,
 +
    shape = 21) +
 +
    xlim(1, Vermax) +
 +
    ylim(0, Edgmax) +
 +
  theme_bw()
</source>
</source>
Line 266: Line 365:
*Ajoutez la courbe correspondant au nombre d'arêtes maximum d'un graphe en fonction du nombre de sommets (fonction ''curve()'').
*Ajoutez la courbe correspondant au nombre d'arêtes maximum d'un graphe en fonction du nombre de sommets (fonction ''curve()'').
-
++
 
<source lang='rsplus'>
<source lang='rsplus'>
xlimit <- 22
xlimit <- 22
curve(<fonction(x)>, from = 0, to = xlimit, n = 200, col="purple", add=T)
curve(<fonction(x)>, from = 0, to = xlimit, n = 200, col="purple", add=T)
 +
</source>
 +
ou
 +
<source lang='rsplus'>
 +
ggplot(iso_cc_property.df, aes(x=vertex, y=edges)) +
 +
  geom_point(
 +
    aes(color = color, fill = color),
 +
    size = 1.5, alpha = 0.5,
 +
    shape = 21) +
 +
    xlim(1, Vermax) +
 +
    ylim(0, Edgmax) +
 +
  geom_function(
 +
    fun = function(x) x*(x-1)/2,
 +
    show.legend = NA,
 +
    colour = "purple"
 +
    ) +
 +
  theme_bw()
</source>
</source>
Line 277: Line 391:
++
++
 +
<source lang='rsplus'>
 +
clique.12 <- iso_cc_property.df %>%
 +
  filter(vertex==12 & edges==66 )
 +
 +
cat("Il y a ", nrow(clique.12), " CC composées de 12 sommets et 66 arêtes (", round(nrow(clique.12)/length(iso_cc)*100, digits=2), "% des CC)", sep="")
 +
</source>
 +
<!--
<source lang='rsplus'>
<source lang='rsplus'>
table(iso_cc_property[,1], iso_cc_property[,2])
table(iso_cc_property[,1], iso_cc_property[,2])
Line 285: Line 406:
maxcol <- colnames(cross)[maxidx[2]]
maxcol <- colnames(cross)[maxidx[2]]
#cat("Il y a ", max(cross), " CC composées de ", maxrow, " sommets et ", maxcol, " arêtes (", max(cross)/length(iso_cc)*100, "% des CC)", sep="")
#cat("Il y a ", max(cross), " CC composées de ", maxrow, " sommets et ", maxcol, " arêtes (", max(cross)/length(iso_cc)*100, "% des CC)", sep="")
 +
</source>
 +
-->
 +
 +
*Distribution du nombre de sommets dans les CC
 +
<source lang='rsplus'>
 +
OG <- iso_cc_property.df %>%
 +
  group_by(vertex, edges) %>%
 +
  summarise(Count=n()/nrow(iso_cc_property.df)*100) %>%
 +
  arrange(desc(Count))
 +
 +
ggplot(OG, aes(x=vertex, y=Count)) +
 +
  geom_bar(stat = "identity", col=color, fill=color) +
 +
  xlab("Hits") +
 +
  ylab("fréquence") +
 +
  theme_bw()
</source>
</source>
Line 341: Line 477:
   return(paralogs)
   return(paralogs)
}
}
 +
 +
count_paralogs_CC <- function(CC)
 +
{
 +
  lst <- V(CC)$name
 +
  all_strains <- substr(lst,1,5)
 +
  strain_list <- unique(all_strains)
 +
  paralogs <- (length(all_strains) - 
 +
              length(strain_list))/length(all_strains)
 +
  return(paralogs)
 +
}
 +
</source>
</source>
Line 364: Line 511:
}
}
hist(pdm1$Pcount, nclass=10)
hist(pdm1$Pcount, nclass=10)
 +
</source>
 +
Une version plus compacte
 +
<source lang='rsplus'>
 +
CC.df <-
 +
  data.frame(Num=(1:length(iso_cc)), Vcount=sapply(iso_cc, vcount, USE.NAMES = TRUE), Ecount=sapply(iso_cc, ecount), Pcount=sapply(iso_cc, count_paralogs_CC))
 +
row.names(CC.df) = 1:nrow(CC.df)
 +
 +
head(CC.df)
 +
 +
nbclass <- 10
 +
rainbowcols <- rainbow(nbclass, s = 0.5)
 +
hist(CC.df$Pcount, nclass=nbclass, ylim=c(0,50), xlab="taux de paralogie", main='', col=rainbowcols)
</source>
</source>
Line 370: Line 529:
++
++
<source lang='rsplus'>
<source lang='rsplus'>
-
vmax <- max(pdm1$Vcount)
+
Vmax <- max(CC.df$Vcount)
-
plot(pdm1[,c(1,2)], las=1, pch=20, col=color, ylim=c(0,1))
+
Pmax <- max(CC.df$Pcount)
 +
ggplot(CC.df, aes(x=Vcount, y=Pcount)) +
 +
  geom_point(
 +
    aes(color = color, fill = color),
 +
    size = 1.5, alpha = 0.5,  
 +
    shape = 21) +
 +
    xlim(1, Vmax) +
 +
    ylim(0, Pmax) +
 +
  theme_bw()
</source>
</source>
*Quelle est le taux minimum de paralogie en fonction du nombre de protéines par CC?
*Quelle est le taux minimum de paralogie en fonction du nombre de protéines par CC?
*Tracez la courbe correspondante.
*Tracez la courbe correspondante.
 +
<source lang='rsplus'>
 +
  geom_function(
 +
    colour = "purple",
 +
    fun = function(x) (x-12)/x,
 +
    show.legend = NA,
 +
    ) +
 +
</source>
-
++
 
Jusqu'à 12 protéines présentes (le nombre de souches) le taux minimum de paralogie est de 0. Au-delà, il est minimum si toutes les souches sont présentes au moins une fois.  
Jusqu'à 12 protéines présentes (le nombre de souches) le taux minimum de paralogie est de 0. Au-delà, il est minimum si toutes les souches sont présentes au moins une fois.  
 +
*Sélectionner les CC avec un taux de paralogie supérieur à 0.5
 +
++
 +
<source lang='rsplus'>
 +
kable(CC.df %>%
 +
        filter(Pcount > 0.5) %>%
 +
        arrange(desc(Pcount)) %>%
 +
        mutate(Pcount=round(Pcount,3)),
 +
      caption="CC avec un taux de paralogie supérieur à 0.5")
 +
</source>
 +
<!--
<source lang='rsplus'>
<source lang='rsplus'>
vmax <- max(pdm1$Vcount)
vmax <- max(pdm1$Vcount)
Line 385: Line 568:
curve((x-nrow(genome_table))/x, from=nrow(genome_table), to=vmax, n=200, col="purple", add=T)
curve((x-nrow(genome_table))/x, from=nrow(genome_table), to=vmax, n=200, col="purple", add=T)
</source>
</source>
 +
-->
==Modularité des CC==
==Modularité des CC==
Line 391: Line 575:
++
++
 +
<source lang='rsplus'>
 +
fastgreedy.Q <- function(CC)
 +
{
 +
fastgreedy <- cluster_fast_greedy(CC, modularity=TRUE)
 +
return(round(max(fastgreedy$modularity), 3))
 +
}
 +
 +
CC.df$fastgreedy <- sapply(iso_cc, fastgreedy.Q)
 +
kable(head(CC.df))
 +
</source>
 +
<!--
 +
old
<source lang='rsplus'>
<source lang='rsplus'>
pdm2 <- data.frame(Vcount=numeric(), Pcount=numeric(), Density=numeric(), Modularity=numeric())
pdm2 <- data.frame(Vcount=numeric(), Pcount=numeric(), Density=numeric(), Modularity=numeric())
Line 411: Line 607:
kable(head(pdm2))
kable(head(pdm2))
</source>
</source>
-
 
+
-->
*Tracer la distribution de Q.
*Tracer la distribution de Q.
++
++
<source lang='rsplus'>
<source lang='rsplus'>
-
hist(pdm2$Modularity, nclass=10, col=color, ylim=c(0, 50))
+
ggplot(CC.df, aes(x=fastgreedy)) +
 +
  geom_histogram(binwidth=0.1, col="grey", fill=color) +
 +
  xlab("Modularity") +
 +
  ylab("frequence") +
 +
  xlim(0, 1) +
 +
  ylim(0, 30) +
 +
  theme_bw()
</source>
</source>
Line 422: Line 624:
++
++
 +
<source lang='rsplus'>
 +
p <- ggplot(CC.df, aes(x=Pcount, y=fastgreedy)) +
 +
  geom_point(
 +
    aes(color = color, fill = color),
 +
    size = 1.5, alpha = 0.5,
 +
    shape = 21) +
 +
  geom_smooth(method=lm , color="purple", fill="grey", se=FALSE) +
 +
  theme_bw()
 +
plot(p)
 +
</source>
 +
<!--
 +
old
<source lang='rsplus'>
<source lang='rsplus'>
plot(pdm2[,c(2,4)], las=1, pch=20, col=color)
plot(pdm2[,c(2,4)], las=1, pch=20, col=color)
</source>
</source>
-
 
+
-->
*Ajoutez la droite de régression *Modularity* en fonction de *Pcount*.
*Ajoutez la droite de régression *Modularity* en fonction de *Pcount*.
++
++
 +
<source lang='rsplus'>
 +
p <- p +  geom_smooth(method=lm , color="purple", fill="grey", se=FALSE) +
 +
plot(p)
 +
 +
</source>
 +
<!--
 +
old
<source lang='rsplus'>
<source lang='rsplus'>
cor <- lm(pdm2$Modularity ~ pdm2$Pcount)
cor <- lm(pdm2$Modularity ~ pdm2$Pcount)
Line 435: Line 656:
curve(cor$coefficients[2]*x +cor$coefficients[1], from=0.05, to=1, col="purple", add=T)
curve(cor$coefficients[2]*x +cor$coefficients[1], from=0.05, to=1, col="purple", add=T)
</source>
</source>
-
 
+
-->
*Commentez le résultat obtenu.
*Commentez le résultat obtenu.
Line 445: Line 666:
*`edge <-  cluster_edge_betweenness(CC, modularity=TRUE)`
*`edge <-  cluster_edge_betweenness(CC, modularity=TRUE)`
*`louvain    <- cluster_louvain(CC)`
*`louvain    <- cluster_louvain(CC)`
-
 
+
<!--
Le nombre de communautés trouvée dans la CC j peut être obtenu avec *length(méthode)*.
Le nombre de communautés trouvée dans la CC j peut être obtenu avec *length(méthode)*.
*`size[j, 1] <- length(fastgreedy)`
*`size[j, 1] <- length(fastgreedy)`
Line 457: Line 678:
*`Modularity5[j, 3] <- round(max(edge$modularity), 3)`
*`Modularity5[j, 3] <- round(max(edge$modularity), 3)`
*`Modularity5[j, 4] <- round(max(edge$modularity), 4)`
*`Modularity5[j, 4] <- round(max(edge$modularity), 4)`
 +
-->
 +
 +
Nous avons déja écrit notre fonction pour fastgreedy. Nous allons ajouter les fonctions poure les autres méthodes.
-
===Un exemple avec fastgreedy===
 
-
Un exemple simple:
 
<source lang='rsplus'>
<source lang='rsplus'>
-
stop = 100
+
walktrap.Q <- function(CC)
-
for (cc_num in 2:stop) {
+
{
-
  CC <- iso_cc[[cc_num]]
+
walktrap <- cluster_walktrap(CC, modularity=TRUE)
-
  fastgreedy <- fastgreedy.community(CC, modularity=TRUE)
+
return(round(max(walktrap$modularity), 3))
-
  Q = round(max(fastgreedy$modularity), 3)
+
-
  cat("CC", cc_num, " modularité", Q, "\n")
+
}
}
-
</source>
 
-
Un code un peu plus complexe:
+
betweenness.Q <- function(CC)
-
<source lang='rsplus'>
+
{
-
vmin <- 12
+
edge <- cluster_edge_betweenness(CC, modularity=TRUE)
-
vmax <- max(sapply(iso_cc, vcount, simplify=TRUE))
+
return(round(max(edge$modularity), 3))
-
edgecolor <- "red"
+
}
-
pdm5 <- data.frame(Vcount=numeric(), Pcount=numeric())
+
louvain.Q <- function(CC)
-
Modularity5 <- data.frame(fastgreedy=numeric())
+
{
-
Size5 <- data.frame(fastgreedy=numeric())
+
louvain <- cluster_louvain(CC)
-
 
+
return(round(modularity(louvain), 3))
-
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 ) {
+
-
 
+
-
fastgreedy <- fastgreedy.community(CC, modularity=TRUE)
+
-
Q = round(max(fastgreedy$modularity), 3)
+
-
if ( Q > 0.45 ) {
+
-
j <- j+1
+
-
list_name <- V(CC)$name     
+
-
pdm5[j, 1] <- cc_num
+
-
pdm5[j, 2] <- nb_vertex
+
-
pdm5[j, 3] <- count_paralogs(list_name)
+
-
 
+
-
Modularity5[j, 1] <- Q
+
-
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("user interruption")
+
-
}
+
-
}
+
}
}
</source>
</source>
-
Réinitialiser ''layout''
+
<!--
-
<source lang='rsplus'>
+
Nous allons extraire un sous ensemble de CC pour réaliser les tests (Q>=0.5).
-
par <- par(mar=c(5.1,4.1,4.1,2.1))
+
-
layout(1)
+
-
</source>
+
-
Avec un code couleur pour les souches:
 
-
<!--
 
<source lang='rsplus'>
<source lang='rsplus'>
-
vmin <- 12
+
Q.thres = 0.5
-
vmax <- max(sapply(iso_cc, vcount, simplify=TRUE))
+
CC.df.light <- CC.df[CC.df$fastgreedy >= Q.thres, ]
-
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)
+
-
Q = round(max(fastgreedy$modularity), 3)
+
-
if ( Q > 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] <- Q
+
-
Size5[j, 1] <- length(fastgreedy)
+
-
+
-
cat("CC", cc_num, vcount(CC), Modularity5[j, 1], Modularity5[j, 2], Modularity5[j, 3], Modularity5[j, 4], "\n")
+
-
+
-
local_gene_list <- substr(V(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))
+
-
Layout <- layout.fruchterman.reingold(CC)
+
kable(CC.df.light)
-
titre <- paste(cc_num, "fastgreedy", round(Modularity5[j, 1],2), Size5[j, 1], sep=" ")
+
-
plot(CC, layout=Layout,
+
-
    vertex.size=15, vertex.shape="sphere", vertex.color=tocolor,
+
-
    vertex.label.cex=0.8, vertex.label=local_gene_list, vertex.label.dist=0,
+
-
    edge.color=edgecolor, edge.width=E(CC)$weight, main=titre)
+
-
+
-
cat ("Press [enter] to continue or quit")
+
-
line <- readline()
+
-
if ( line == 'quit') stop("user interruption")
+
-
}
+
-
}
+
-
}
+
</source>
</source>
-->
-->
 +
Pour chaque CC trier les méthodes en fonction des modularités obtenues. Représentez graphiquement les résultats.
-
===Ajouter les autres méthodes===
+
Rang des modularités obtenues par les méthodes pour chaque CC
-
Le principe :
+
<source lang='rsplus'>
<source lang='rsplus'>
-
stop = length(iso_cc)
+
Q.thres = 0.1
-
for (cc_num in 2:stop) {
+
CC.df.positif <- CC.df[CC.df$fastgreedy >= Q.thres, ]
-
  CC <- iso_cc[[cc_num]]
+
Modu_stat <- apply(CC.df.positif, 1, function(x) rank(signif(x[5:8], digits=5), ties.method = "max"))
-
nb_vertex <- vcount(CC)
+
-
  fastgreedy <- fastgreedy.community(CC, modularity=TRUE)
+
-
  Q = round(max(fastgreedy$modularity), 2)
+
-
  if ( Q > 0.45 ) {
+
-
    list_name <- V(CC)$name
+
-
    P <- round(count_paralogs(list_name), 2)
+
-
 
+
-
    walktrap <- walktrap.community(CC, modularity=TRUE)
+
-
    W <- round(max(walktrap$modularity), 2)
+
-
    edge <- edge.betweenness.community(CC, modularity=TRUE)
+
-
    E <- round(max(edge$modularity), 2)
+
-
    louvain <- cluster_louvain(CC)
+
-
    L <- round(max(louvain$modularity), 2)
+
-
    cat("CC", cc_num, nb_vertex, P, Q, W, E, L, "\n")
+
-
  }
+
-
}
+
-
</source>
+
-
l'application:
+
-
++
+
# fréquence des différents rangs
-
<source lang='rsplus'>
+
Modu_tab <- apply(Modu_stat, 1, function(x) {
-
par(mar=c(1,1,1,1))
+
  his <- hist(x,breaks=(0:4), plot=F)
-
layout(matrix(1:4, 1, 4))
+
   his$count
-
vmin <- 12
+
})
-
vmax <- max(sapply(iso_cc, vcount, simplify=TRUE))
+
rainbowcols <- rainbow(4, s = 0.5)
-
edgecolor <- "red"
+
barplot(Modu_tab, col=rainbowcols, ylab="frequency")
-
pdm5 <- data.frame(CC=numeric(), Vcount=numeric(), Pcount=numeric())
+
-
Modularity5 <- data.frame(fastgreedy=numeric(), walktrap=numeric(), edge=numeric(), louvain=numeric())
+
-
Size5 <- data.frame(fastgreedy=numeric(), walktrap=numeric(), edge=numeric(), louvain=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 <- cluster_fast_greedy(CC, modularity=TRUE)
+
-
fastgreedy <- cluster_fast_greedy(CC, modularity=TRUE)
+
-
if ( max(fastgreedy$modularity) > 0.45 ) {
+
-
j <- j+1
+
-
walktrap  <- cluster_walktrap(CC, modularity=TRUE)
+
-
edge      <- cluster_edge_betweenness(CC, modularity=TRUE)
+
-
louvain    <- cluster_louvain(CC)
+
-
 
+
-
     
+
-
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(louvain$modularity), 3)
+
-
 
+
-
  if ( min(Modularity5[j,]) < max(Modularity5[j,])) {
+
-
  Size5[j, 1] <- length(fastgreedy)
+
-
  Size5[j, 2] <- length(walktrap)
+
-
  Size5[j, 3] <- length(edge)
+
-
  Size5[j, 4] <- length(louvain)
+
-
 
+
-
 
+
-
 
+
-
   
+
-
    #cat("CC", cc_num, vcount(CC), Modularity5[j, 1], Modularity5[j, 2], Modularity5[j, 3], "\n")
+
-
 
+
-
  Layout = layout.fruchterman.reingold(CC)
+
-
  titre <- paste(cc_num, "fastgreedy", round(Modularity5[j, 1],3), 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)
+
-
    
+
-
  titre <- paste("walktrap",  round(Modularity5[j, 2],3), Size5[j, 2], sep=" ")
+
-
  plot(walktrap, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, edge.width=E(CC)$weight, main=titre)
+
-
 
+
-
  titre <- paste("edge",  round(Modularity5[j, 3],3), Size5[j, 3], sep=" ")
+
-
  plot(edge, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, edge.width=E(CC)$weight, main=titre)
+
-
 
+
-
  titre <- paste("louvain",  round(Modularity5[j, 4],3), Size5[j, 3], sep=" ")
+
-
  plot(louvain, 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)
+
</source>
</source>
-
 
===spinglass===  
===spinglass===  
La méthode spinglass repose sur un plus grand nombre de paramètres:
La méthode spinglass repose sur un plus grand nombre de paramètres:
Line 710: Line 767:
*Diminuez la valeur de ''stop.temp'' (0.4 par exemple). Que pouvez-vous observer?
*Diminuez la valeur de ''stop.temp'' (0.4 par exemple). Que pouvez-vous observer?
 +
*construction de la fonction
 +
<source lang='rsplus'>
 +
spinglass.Q <- function(CC, gamma=1, spins=25, start_temp=1, stop_temp=0.8, cool_fact=0.95, update_rule="config", iter=100)
 +
{
 +
  Qiter <- array(0, c(iter))
 +
  spinglass_iter <- list()
 +
  for ( k in 1:iter )  {
 +
  spinglass <- cluster_spinglass(CC, spins=spins,
 +
  start.temp=start_temp, stop.temp=stop_temp,
 +
  cool.fact=cool_fact, update.rule=update_rule,
 +
  gamma=gamma)
 +
  Qiter[k] <- modularity(spinglass)
 +
  spinglass_iter[[k]] <- spinglass
 +
  }
 +
  return(round(max(Qiter), 3))
 +
}
 +
spinglass.Q(iso_cc[[1584]], iter=10, stop_temp=0.1)
 +
</source>
====Assemblage des cinq méthodes====
====Assemblage des cinq méthodes====
-
Nous allons fixer les paramètres suivant:
+
Nous allons réitérer spinglass 10 fois pour chaque CC. L'itération maximisant la modularité sera conservée.   
-
<pre>
+
-
  `my_gamma <- 1`
+
-
  `my_spins <- 25`
+
-
  `stop_temp <- 0.8`
+
-
  `start_temp <- 1`
+
-
  `cool_fact <- 0.95`
+
-
  `update.rule <- "config"`
+
-
</pre>
+
-
et nous allons réitérer spinglass 10 fois pour chaque CC. L'itération maximisant la modularité sera conservée.   
+
<pre>
<pre>
   `iter <- 10`
   `iter <- 10`
Line 726: Line 792:
*Assemblez les cinq méthodes dans un seul bloc.
*Assemblez les cinq méthodes dans un seul bloc.
-
Principe :
 
<source lang='rsplus'>
<source lang='rsplus'>
-
my_gamma <- 1
+
iso_cc[CC.df.positif$Num]
-
my_spins <- 25
+
CC.df.positif$spinglass <- sapply(iso_cc[CC.df.positif$Num], function(x) spinglass.Q(x, iter=10, stop_temp=0.1))
-
stop_temp <- 0.1
+
kable(head(CC.df))
-
start_temp <- 1
+
 
-
cool_fact <- 0.95
+
kable(CC.df.positif, row.names = FALSE)
-
update.rule <- "config"
+
-
iter <- 10
+
-
stop = length(iso_cc)
+
Modu_stat.spin <- apply(CC.df.positif, 1, function(x) rank(signif(x[5:9], digits=5), ties.method = "max"))
-
for (cc_num in 2:stop) {
+
Modu_tab.spin <- apply(Modu_stat.spin, 1, function(x) {
-
  CC <- iso_cc[[cc_num]]
+
  his <- hist(x, breaks=(0:5), plot=F)
-
  nb_vertex <- vcount(CC)
+
  his$count
-
  fastgreedy <- fastgreedy.community(CC, modularity=TRUE)
+
})
-
  Q = round(max(fastgreedy$modularity), 2)
+
barplot(Modu_tab.spin, col=c("blue", "green", "yellow", "orange", "red"), ylab="frequency")
-
  if ( Q > 0.45 ) {
+
</source>
-
    list_name <- V(CC)$name
+
 
-
    P <- round(count_paralogs(list_name), 2)
+
* Tracer les graphes des CC
-
 
+
 
-
    walktrap <- walktrap.community(CC, modularity=TRUE)
+
<source lang='rsplus'>
-
    W <- round(max(walktrap$modularity), 2)
+
edgecolor <- "red"
-
    edge <- edge.betweenness.community(CC, modularity=TRUE)
+
 
-
    E <- round(max(edge$modularity), 2)
+
for (i in CC.df.light$Num) {
-
    louvain <- cluster_louvain(CC)
+
  cat("num", i, "\n")
-
    L <- round(max(louvain$modularity), 2)
+
  CC <- iso_cc[[i]]
-
+
fastgreedy.comm <- cluster_fast_greedy(CC, modularity=TRUE)
-
    Qiter <- array(0, c(iter))
+
plot(fastgreedy.comm, CC, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, , main=i)
-
    spinglass_iter <- list() 
+
-
    for ( k in 1:iter ) {
+
-
      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
+
-
    }
+
-
    S <- round(max(Qiter), 2)
+
-
    if ( min(qiter) < max(qiter)) {
+
-
      cat("CC", cc_num, nb_vertex, P, Q, W, E, L, S, "\n")
+
-
    }
+
-
  }
+
}
}
</source>
</source>
-
l'ensemble
 
-
++
+
Un embellissement
 +
 
<source lang='rsplus'>
<source lang='rsplus'>
-
par(mar=c(1,1,1,1))
+
for (i in CC.df.light$Num) {
-
layout(matrix(1:5, 1, 5))
+
  cat("num", i, "\n")
-
vmin <- 12
+
  CC <- iso_cc[[i]]
-
vmax <- max(sapply(iso_cc, vcount, simplify=TRUE))
+
local_gene_list <- substr(V(CC)$name, 5, 5)
-
pdm5 <- data.frame(cc_num=numeric(), Vcount=numeric(), Pcount=numeric())
+
local_strain_list <- unique(local_gene_list)
-
Modularity5 <- data.frame(fastgreedy=numeric(), walktrap=numeric(), edge=numeric(), louvain=numeric(), spinglass=numeric())
+
tocolor <- unlist(lapply(local_gene_list, function(x,y) grep(x, y), local_strain_list))
-
Size5 <- data.frame(fastgreedy=numeric(), walktrap=numeric(), edge=numeric(), louvain=numeric(), spinglass=numeric())
+
   
   
-
my_gamma <- 1
+
plot(CC,
-
my_spins <- 25
+
    vertex.size=15,  
-
stop_temp <- 0.1
+
    vertex.shape="sphere",  
-
start_temp <- 1
+
    vertex.color=tocolor,  
-
cool_fact <- 0.95
+
    vertex.label.cex=0.8,  
-
update.rule <- "config"
+
    vertex.label=local_gene_list,  
-
iter <- 10
+
    vertex.label.dist=0,  
-
+
    edge.color=edgecolor,  
-
j <- 0
+
    main=i)
-
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)
+
-
if ( max(fastgreedy$modularity) > 0.3 ) {
+
-
j <- j+1
+
-
walktrap  <- cluster_walktrap(CC, modularity=TRUE)
+
-
edge      <- cluster_edge_betweenness(CC, modularity=TRUE)
+
-
louvain    <- cluster_louvain(CC)
+
-
+
-
Qiter <- array(0, c(iter))
+
-
spinglass_iter <- list() 
+
-
for ( k in 1:iter ) {
+
-
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)
+
-
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(louvain$modularity), 3)
+
-
Modularity5[j, 5] <- round(max(Qiter), 3)
+
-
+
-
Size5[j, 1] <- length(fastgreedy)
+
-
  Size5[j, 2] <- length(walktrap)
+
-
Size5[j, 3] <- length(edge)
+
-
Size5[j, 4] <- length(edge)
+
-
Size5[j, 5] <- length(unique(membership))
+
-
+
-
#cat("CC", cc_num, vcount(CC), Modularity5[j, 1], Modularity5[j, 2], Modularity5[j, 3], Modularity5[j, 4], "\n")
+
-
      #if ( Modularity5[j, 5] > min(Modularity5[j,])) {
+
-
  Layout = layout.fruchterman.reingold(CC)
+
-
  titre <- paste(cc_num, "fastgreedy", round(Modularity5[j, 1],3), Size5[j, 1], sep=" ")
+
-
  plot(fastgreedy, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color='grey', edge.width=E(CC)$weight, main=titre)
+
-
 
+
-
  titre <- paste("walktrap",  round(Modularity5[j, 2],3), Size5[j, 2], sep=" ")
+
-
  plot(walktrap, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color='grey', edge.width=E(CC)$weight, main=titre)
+
-
 
+
-
  titre <- paste("edge",  round(Modularity5[j, 3],3), Size5[j, 3], sep=" ")
+
-
  plot(edge, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color='grey', edge.width=E(CC)$weight, main=titre)
+
-
 
+
-
  titre <- paste("louvain",  round(Modularity5[j, 4],3), Size5[j, 3], sep=" ")
+
-
  plot(louvain, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color='grey', edge.width=E(CC)$weight, main=titre)
+
-
+
-
  titre <- paste("spinglass",  round(Modularity5[j, 4],3), Size5[j, 4], sep=" ")
+
-
  plot(spinglass_iter[[step]], CC, layout=Layout, vertex.label=NA, vertex.size=3, edge.color='grey', 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)
 
</source>
</source>
-
Pour chaque CC triez les méthodes en fonction des modularités obtenues.
+
=== Ajouter les autres méthodes à ce code ===
-
Représentez graphiquement les résultats.
+
-
++
 
<source lang='rsplus'>
<source lang='rsplus'>
-
Modu_stat <- apply(Modularity5, 1, function(x) rank(signif(x[1:5], digits=5), ties.method = "max"))
+
par(mar=c(1,1,1,1))
 +
layout(matrix(1:4, 1, 4))
-
# fréquence des différents rangs
+
edgecolor <- "red"
-
Modu_tab <- apply(Modu_stat, 1, function(x) {
+
-
  his <- hist(x, breaks=(0:5), plot=F)
+
-
  his$count
+
-
})
+
-
rainbowcols <- rainbow(5, s = 0.2)
+
-
barplot(Modu_tab, col=rainbowcols, ylab="frequency", ylim=c(0,25))
+
-
</source>
+
-
Pour chaque CC triez les méthodes en fonction du nombre de communautés obtenues.  
+
for (i in CC.df.light$Num) {
-
Représentez graphiquement les résultats.
+
  cat("num", i, "\n")
 +
  CC <- iso_cc[[i]]
 +
fastgreedy.com <- cluster_fast_greedy(CC, modularity=TRUE)
 +
Layout = layout.fruchterman.reingold(CC)
 +
fastgreedy.Q <- CC.df.light[CC.df.light$Num==i, ]$fastgreedy
 +
titre <- paste(cc_num, "fastgreedy", fastgreedy.Q, sep=" ")
 +
plot(fastgreedy.com, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, main=titre)
-
++
+
walktrap.com <- cluster_walktrap(CC, modularity=TRUE)
-
<source lang='rsplus'>
+
walktrap.Q <- CC.df.light[CC.df.light$Num==i, ]$walktrap
-
Size_stat <- apply(Size5, 1, function(x) rank(signif(x[1:4], digits=5), ties.method = "max"))
+
titre <- paste("walktrap", walktrap.Q, sep=" ")
 +
plot(walktrap.com, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, main=titre)
-
Size_tab <- apply(Size_stat, 1, function(x) {
+
betweenness.com <- cluster_edge_betweenness(CC, modularity=TRUE)
-
  his <- hist(x,breaks=(0:4), plot=F)
+
betweenness.Q <- CC.df.light[CC.df.light$Num==i, ]$betweenness
-
  his$count
+
titre <- paste("edge", betweenness.Q, sep=" ")
-
})
+
plot(betweenness.com, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, main=titre)
-
barplot(Size_tab)
+
     
 +
louvain.com <-  cluster_louvain(CC)
 +
louvain.Q <- CC.df.light[CC.df.light$Num==i, ]$louvain
 +
titre <- paste("louvain", louvain.Q, sep=" ")
 +
plot(louvain.com, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, main=titre)
 +
 
 +
cat ("Press [enter] to continue or quit")
 +
line <- readline()
 +
if ( line == 'quit') stop()
 +
}
</source>
</source>
-
*Commentez les résultats.
+
=== Intersection des communautées ===
-
 
+
-
Pour obtenir les résultats sur toutes les CC, relancez la boucle en supprimant le test:
+
<source lang='rsplus'>
<source lang='rsplus'>
-
      cat ("Press [enter] to continue or quit")
+
for (i in CC.df.light$Num) {
-
      line <- readline()
+
  cat("num", i, "\n")
-
      if ( line == 'quit') stop()
+
  CC <- iso_cc[[i]]
 +
fastgreedy.com <- cluster_fast_greedy(CC, modularity=TRUE)
 +
  fastgreedy.pairs <- unlist(lapply(fastgreedy.com[], function(x) levels(interaction(x,x))))
 +
 
 +
  betweenness.com <-  cluster_edge_betweenness(CC, modularity=TRUE)
 +
  betweenness.pairs <- unlist(lapply(betweenness.com[], function(x) levels(interaction(x,x))))
 +
 +
  walktrap.com <- cluster_walktrap(CC, modularity=TRUE)
 +
  walktrap.pairs <- unlist(lapply(walktrap.com[], function(x) levels(interaction(x,x))))
 +
 
 +
  louvain.com <-  cluster_louvain(CC)
 +
  louvain.pairs <- unlist(lapply(louvain.com[], function(x) levels(interaction(x,x))))
 +
 
 +
  p <- ggvenn(list(fastgreedy=fastgreedy.pairs, betweenness=betweenness.pairs, walktrap=walktrap.pairs, louvain=louvain.pairs),  show_percentage = FALSE, fill_alpha = 0.3, stroke_size = 0.1, set_name_size = 5)
 +
  plot(p)
 +
 
 +
cat ("Press [enter] to continue or quit")
 +
line <- readline()
 +
  if ( line == 'quit') stop()
 +
}
</source>
</source>
-
 
-
*Tracez les deux plots et commentez les résultats obtenus avec les CC de taille supérieur à 12 arêtes!
 
=End=
=End=

Current revision as of 14:17, 29 November 2022

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 (Tettelin et al., 2005). 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 (Kim et al., 2020)! 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

Nous allons travailler avec la bactérie Yersinia pestis, une bactérie à Gram négatif du genre Yersinia. Elle est responsable de la peste. Elle a été découverte en 1894 par Alexandre Yersin, un bactériologiste franco-suisse travaillant pour l'Institut Pasteur, durant une épidémie de peste à Hong Kong, en même temps que Kitasato Shibasaburō mais séparément. Kitasato tout d'abord la baptisa Pasteurella pestis en l'honneur de Pasteur. Ce n'est que plus tard qu'elle prit son nom actuel, en hommage à Yersin.

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.

Dans un premier temps nous allons télécharger my_micropan.R et les archives data et homologs.

ou sous linux:

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

Vous pouvez effectuer les téléchargements des fichiers dans une Chunk bash dans Rstudio. 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 la souche 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

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)

Définir le répertoire de travail

knitr::opts_knit$set(root.dir = "/home/user/TP5/")

Chargement des fichiers

Téléchargement et décompression des fichiers dans un shell.

rootdir=$HOME/Graphe/micropan/Ypestis
echo $rootdir
mkdir -p $rootdir
cd $rootdir
pwd
 
if [ ! -f my_micropan.R ];
then
  wget http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/Ypestis/my_micropan.R -nc
else
  echo skip my_micropan.R
fi
if [ ! -d data ];
then
  wget http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/Ypestis/data.tar.gz -nc
  tar -xvf data.tar.gz
else
  echo skip data
fi
if [ ! -d homologs ];
then
  wget http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/Ypestis/homologs.tar.gz -nc
  tar -xvf homologs.tar.gz
else
  echo skip homologs
fi

Nous allons nous placer dans le répertoire de travail.

setwd(rootdir)
getwd()
list.files()
 
source("my_micropan.R")

Nous allons utiliser les librairies : igraph, tidyverse 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 dans les résultats des blast.

library(tidyverse)
library(igraph)
library(ggvenn)
 
options (stringsAsFactors = FALSE)
savepar <- par() # save Graphical Parameters
par(las=1)
color <- "coral"

Tableau de données

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

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).
genome_table <- read.table(file=genome_table_file, sep = "\t" , header= TRUE)
head(genome_table)

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

++

boxplot(genome_table$Proteome, ylab='Proteome size', col='coral', las=1)

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

Orthologs 1:1

Dans un premier temps, nous allons réaliser 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?

Vous pouvez regrouper les lignes par requête (group_by(Query)) et compter le nombre de cibles (summarise(Hits=n())) pour chaque requête.

++

nbaliprot <- homologs %>% 
  group_by(Query) %>%
  summarise(Hits=n()) %>% 
  arrange(desc(Hits))
  • Tracez la distribution du nombre d'alignement par protéine avec geom_histogram(binwidth=1, col=color, fill=color).

++

ggplot(nbaliprot, aes(x=Hits)) + 
  geom_histogram(binwidth=1, col=color, fill=color) + 
  scale_y_log10() +
  theme_bw()
  • Identifiez la protéine qui a le plus d'alignements et analysez les alignements obtenus.
  • sélectionner la protéine qui a le plus de hits (Hits==max(Hits))
  • avec son nom, filtrer la table homologs (Query==max_hits$Query)
  • sélectionner les colonnes: Query, Hit, Query.start, Query.end, Hit.start et Hit.end
  • éditer cette table avec kable.

++

max_hits <- nbaliprot %>% 
  filter(Hits==max(Hits)) %>%
  select(Query)
max_hits$Query
 
subtab <- homologs %>% 
  filter(Query==max_hits$Query)
head(subtab)
 
compact <- subtab %>% 
  select(Query, Hit, Query.start, Query.end, Hit.start, Hit.end) %>% 
  arrange(Query.start)
 
kable(compact)
  • Est-ce que toutes les protéines de A ont un alignement dans B?

++

YpesA <- genome_table %>% 
  filter(File == "YpesA")
 
YpesA.Protein.Hits <- nrow(homologs %>% distinct(Query))
 
cat("Sur les ", YpesA$Proteome, " proteines de YpesA, seules ", YpesA.Protein.Hits, " 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
dir.create("res", showWarnings = TRUE)
graph_file <- "res/iso_gr"
if (file.exists(graph_file)) file.remove(graph_file)
 
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 %>% 
    filter(Homology.type == 'Iso') %>%
    select(Query, Hit)
  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 votre bloc de commandes.

Traitement du graphe

Relire ce fichier comme un graphe (fonction read.graph() de igraph, avec format="ncol").

++

gr <- read.graph(file=graph_file, format="ncol")

Affichez le nombre de sommets et d’arêtes (vcount() et ecount()).

  • Le graphe est composé de `r vcount(gr)` sommets et de `r ecount(gr)` arêtes.

Appliquez la fonction simplify() sur le graphe.

++

gr_iso <- simplify(gr)

Affichez le nombre de sommets et d’arêtes du graphe simplifié.

  • Le graphe simplifié est composé de `r vcount(gr_iso)` sommets et de `r ecount(gr_iso)` arêtes.
  • 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".

++

iso_cc <- decompose.graph(gr_iso, 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é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 unlist(lapply()) ou sapply() sur la liste de CC pour calculer le nombre de sommets et le nombre d'arêtes de chaque 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])
 
iso_cc_property <- cbind(sapply(iso_cc, vcount), sapply(iso_cc, ecount))
colnames(iso_cc_property) <- c("vertex", "edges")

Tracer le nombre d'arêtes en fonction du nombre de sommets.

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

ou

iso_cc_property.df <- as.data.frame(iso_cc_property)
Vermax <- max(iso_cc_property.df$vertex)
Edgmax <- max(iso_cc_property.df$edges)
 
ggplot(iso_cc_property.df, aes(x=vertex, y=edges)) +
  geom_point(
    aes(color = color, fill = color),
    size = 1.5, alpha = 0.5, 
    shape = 21) +
    xlim(1, Vermax) +
    ylim(0, Edgmax) +
  theme_bw()
  • 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(<fonction(x)>, from = 0, to = xlimit, n = 200, col="purple", add=T)

ou

ggplot(iso_cc_property.df, aes(x=vertex, y=edges)) +
  geom_point(
    aes(color = color, fill = color),
    size = 1.5, alpha = 0.5, 
    shape = 21) +
    xlim(1, Vermax) +
    ylim(0, Edgmax) +
  geom_function(
    fun = function(x) x*(x-1)/2,
    show.legend = NA,
    colour = "purple"
    ) +
  theme_bw()
  • A quoi correspondent 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?

++

clique.12 <- iso_cc_property.df %>% 
  filter(vertex==12 & edges==66 )
 
cat("Il y a ", nrow(clique.12), " CC composées de 12 sommets et 66 arêtes (", round(nrow(clique.12)/length(iso_cc)*100, digits=2), "% des CC)", sep="")
  • Distribution du nombre de sommets dans les CC
OG <- iso_cc_property.df %>% 
  group_by(vertex, edges) %>%
  summarise(Count=n()/nrow(iso_cc_property.df)*100) %>% 
  arrange(desc(Count))
 
ggplot(OG, aes(x=vertex, y=Count)) + 
  geom_bar(stat = "identity", col=color, fill=color) + 
  xlab("Hits") +
  ylab("fréquence") + 
  theme_bw()
  • 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")

Un embellissement!

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 G présent dans la liste,
  en extraire le nombre de souches S. 

S'il n'y a pas de paralogue, chaque souche est représentée par un seul gène (G=S), sinon (G>S), on a plus d'un gène par souche.

Pour estimer le taux de paralogues, on divise la différence entre nombre de gènes et nombre de souches (G-S) par le nombre de gènes G.

Les fonctions à utiliser sont :

  **substr()** pour extraire le nom de souche à partir du nom de gène 
  **unique()** pour obtenir le nombre de souche. 

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)
}
 
count_paralogs_CC <- function(CC)
{
  lst <- V(CC)$name
  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)

Une version plus compacte

CC.df <- 
  data.frame(Num=(1:length(iso_cc)), Vcount=sapply(iso_cc, vcount, USE.NAMES = TRUE), Ecount=sapply(iso_cc, ecount), Pcount=sapply(iso_cc, count_paralogs_CC))
row.names(CC.df) = 1:nrow(CC.df)
 
head(CC.df)
 
nbclass <- 10
rainbowcols <- rainbow(nbclass, s = 0.5)
hist(CC.df$Pcount, nclass=nbclass, ylim=c(0,50), xlab="taux de paralogie", main='', col=rainbowcols)

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

++

Vmax <- max(CC.df$Vcount)
Pmax <- max(CC.df$Pcount)
ggplot(CC.df, aes(x=Vcount, y=Pcount)) +
  geom_point(
    aes(color = color, fill = color),
    size = 1.5, alpha = 0.5, 
    shape = 21) +
    xlim(1, Vmax) +
    ylim(0, Pmax) +
  theme_bw()
  • Quelle est le taux minimum de paralogie en fonction du nombre de protéines par CC?
  • Tracez la courbe correspondante.
  geom_function(
    colour = "purple",
    fun = function(x) (x-12)/x,
    show.legend = NA,
    ) +

Jusqu'à 12 protéines présentes (le nombre de souches) le taux minimum de paralogie est de 0. Au-delà, il est minimum si toutes les souches sont présentes au moins une fois.

  • Sélectionner les CC avec un taux de paralogie supérieur à 0.5

++

kable(CC.df %>% 
        filter(Pcount > 0.5) %>% 
        arrange(desc(Pcount)) %>% 
        mutate(Pcount=round(Pcount,3)), 
      caption="CC avec un taux de paralogie supérieur à 0.5")

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

++

fastgreedy.Q <- function(CC)
{
		fastgreedy <- cluster_fast_greedy(CC, modularity=TRUE)
		return(round(max(fastgreedy$modularity), 3))
}
 
CC.df$fastgreedy <- sapply(iso_cc, fastgreedy.Q)
kable(head(CC.df))
  • Tracer la distribution de Q.

++

ggplot(CC.df, aes(x=fastgreedy)) + 
  geom_histogram(binwidth=0.1, col="grey", fill=color) + 
  xlab("Modularity") +
  ylab("frequence") + 
  xlim(0, 1) +
  ylim(0, 30) +
  theme_bw()
  • Tracez la modularité des CC en fonction du taux de paralogie.

++

p <- ggplot(CC.df, aes(x=Pcount, y=fastgreedy)) +
  geom_point(
    aes(color = color, fill = color),
    size = 1.5, alpha = 0.5, 
    shape = 21) +
  geom_smooth(method=lm , color="purple", fill="grey", se=FALSE) +
  theme_bw()
plot(p)
  • Ajoutez la droite de régression *Modularity* en fonction de *Pcount*.

++

p <- p +  geom_smooth(method=lm , color="purple", fill="grey", se=FALSE) +
plot(p)
  • 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, edge betweenness et louvain. Ces méthodes sont utilisées de façon similaire (vous pouvez ajouter la méthode cluster_leiden si elle disponible dans votre version d'igraph):

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

Nous avons déja écrit notre fonction pour fastgreedy. Nous allons ajouter les fonctions poure les autres méthodes.

walktrap.Q <- function(CC)
{
		walktrap <- cluster_walktrap(CC, modularity=TRUE)
		return(round(max(walktrap$modularity), 3))
}
 
betweenness.Q <- function(CC)
{
		edge <- cluster_edge_betweenness(CC, modularity=TRUE)
		return(round(max(edge$modularity), 3))
}
 
louvain.Q <- function(CC)
{
		louvain <-  cluster_louvain(CC)
		return(round(modularity(louvain), 3))
}

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

Rang des modularités obtenues par les méthodes pour chaque CC

Q.thres = 0.1
CC.df.positif <- CC.df[CC.df$fastgreedy >= Q.thres, ]
Modu_stat <- apply(CC.df.positif, 1, function(x) rank(signif(x[5:8], digits=5), ties.method = "max"))
 
# fréquence des différents rangs
Modu_tab <- apply(Modu_stat, 1, function(x) { 
  his <- hist(x,breaks=(0:4), plot=F)
  his$count
})
rainbowcols <- rainbow(4, s = 0.5)
barplot(Modu_tab, col=rainbowcols, ylab="frequency")

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)

Distribution de Q en fonction des itérations :

plot(Qiter, ylim=c(0,1), xlab="iteration", ylab="Modularity")
hist(Qiter, col=color, xlim=c(0,1))
  • Diminuez la valeur de stop.temp (0.4 par exemple). Que pouvez-vous observer?
  • construction de la fonction
spinglass.Q <- function(CC, gamma=1, spins=25, start_temp=1, stop_temp=0.8, cool_fact=0.95, update_rule="config", iter=100)
{
  Qiter <- array(0, c(iter))
  spinglass_iter <- list() 
  for ( k in 1:iter )  {
   	spinglass <- cluster_spinglass(CC, spins=spins,
   	start.temp=start_temp, stop.temp=stop_temp,
   	cool.fact=cool_fact, update.rule=update_rule, 
   	gamma=gamma)
   	Qiter[k] <- modularity(spinglass)
   	spinglass_iter[[k]] <- spinglass
  }
  return(round(max(Qiter), 3))
}
spinglass.Q(iso_cc[[1584]], iter=10, stop_temp=0.1)

Assemblage des cinq méthodes

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

  `iter <- 10`
  • Assemblez les cinq méthodes dans un seul bloc.
iso_cc[CC.df.positif$Num]
CC.df.positif$spinglass <- sapply(iso_cc[CC.df.positif$Num], function(x) spinglass.Q(x, iter=10, stop_temp=0.1))
kable(head(CC.df))
 
kable(CC.df.positif, row.names = FALSE)
 
Modu_stat.spin <- apply(CC.df.positif, 1, function(x) rank(signif(x[5:9], digits=5), ties.method = "max"))
Modu_tab.spin <- apply(Modu_stat.spin, 1, function(x) { 
  his <- hist(x, breaks=(0:5), plot=F)
  his$count
})
barplot(Modu_tab.spin, col=c("blue", "green", "yellow", "orange", "red"), ylab="frequency")
* Tracer les graphes des CC
edgecolor <- "red"
 
for (i in CC.df.light$Num) {
  cat("num", i, "\n")
  CC <- iso_cc[[i]]
	fastgreedy.comm <- cluster_fast_greedy(CC, modularity=TRUE)
	plot(fastgreedy.comm, CC, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, , main=i)
}

Un embellissement

for (i in CC.df.light$Num) {
  cat("num", i, "\n")
  CC <- iso_cc[[i]]
			local_gene_list <- substr(V(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))
 
			plot(CC,
			     vertex.size=15, 
			     vertex.shape="sphere", 
			     vertex.color=tocolor, 
			     vertex.label.cex=0.8, 
			     vertex.label=local_gene_list, 
			     vertex.label.dist=0, 
			     edge.color=edgecolor, 
			     main=i)
}

Ajouter les autres méthodes à ce code

par(mar=c(1,1,1,1))
layout(matrix(1:4, 1, 4))
 
edgecolor <- "red"
 
for (i in CC.df.light$Num) {
  cat("num", i, "\n")
  CC <- iso_cc[[i]]
	fastgreedy.com <- cluster_fast_greedy(CC, modularity=TRUE)
	Layout = layout.fruchterman.reingold(CC) 
	fastgreedy.Q <- CC.df.light[CC.df.light$Num==i, ]$fastgreedy
	titre <- paste(cc_num, "fastgreedy", fastgreedy.Q, sep=" ")
	plot(fastgreedy.com, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, main=titre)
 
	walktrap.com <- cluster_walktrap(CC, modularity=TRUE)
	walktrap.Q <- CC.df.light[CC.df.light$Num==i, ]$walktrap
	titre <- paste("walktrap", walktrap.Q, sep=" ")
	plot(walktrap.com, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, main=titre)
 
	betweenness.com <-  cluster_edge_betweenness(CC, modularity=TRUE)
	betweenness.Q <- CC.df.light[CC.df.light$Num==i, ]$betweenness
	titre <- paste("edge", betweenness.Q, sep=" ")
	plot(betweenness.com, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, main=titre)
 
	louvain.com <-  cluster_louvain(CC)
	louvain.Q <- CC.df.light[CC.df.light$Num==i, ]$louvain
	titre <- paste("louvain", louvain.Q, sep=" ")
	plot(louvain.com, CC, layout=Layout, vertex.label=NA, vertex.size=3, vertex.shape="sphere", edge.color=edgecolor, main=titre)
 
	cat ("Press [enter] to continue or quit")
	line <- readline()
	if ( line == 'quit') stop()
}

Intersection des communautées

for (i in CC.df.light$Num) {
  cat("num", i, "\n")
  CC <- iso_cc[[i]]
	fastgreedy.com <- cluster_fast_greedy(CC, modularity=TRUE)
  fastgreedy.pairs <- unlist(lapply(fastgreedy.com[], function(x) levels(interaction(x,x))))
 
  betweenness.com <-  cluster_edge_betweenness(CC, modularity=TRUE)
  betweenness.pairs <- unlist(lapply(betweenness.com[], function(x) levels(interaction(x,x))))
 
  walktrap.com <- cluster_walktrap(CC, modularity=TRUE)
  walktrap.pairs <- unlist(lapply(walktrap.com[], function(x) levels(interaction(x,x))))
 
  louvain.com <-  cluster_louvain(CC)
  louvain.pairs <- unlist(lapply(louvain.com[], function(x) levels(interaction(x,x))))
 
  p <- ggvenn(list(fastgreedy=fastgreedy.pairs, betweenness=betweenness.pairs, walktrap=walktrap.pairs, louvain=louvain.pairs),  show_percentage = FALSE, fill_alpha = 0.3, stroke_size = 0.1, set_name_size = 5)
  plot(p)
 
	cat ("Press [enter] to continue or quit")
	line <- readline()
  if ( line == 'quit') stop()
}

End