---
title: "Yersinia pestis"
author: "Yves Quentin"
date: "2016"
output: html_document
---

```{r global_options, include=FALSE}
library(knitr)
knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/', echo=FALSE, warning=FALSE, message=FALSE, cache=TRUE)

cachedata=TRUE
savepar <- par()
par(las=1)
color <- "coral"

```

#Introduction

Nous allons r?aliser une analyse de g?nomique comparative qui s'inspire d'un travail publi? par [Snipen et Liland (2015)](http://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0517-0). 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](https://cran.r-project.org/web/packages/micropan/index.html)). Ils ont mis ? disposition une ?tude de cas [casestudy. pdf](http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/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](http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/Ypestis). Les prot?ines en format FASTA sont dans le r?pertoire (**data/proteins**). Les comparaisons entre pairs 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. 

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 oeuvre

Nous allons utiliser les librairies :
*[micropan](https://cran.r-project.org/web/packages/micropan/index.html)
*[iGraph](http://igraph.sourceforge.net/).
et le fichier [my_micropan.R](http://silico.biotoul.fr/enseignement/m1-bioinfo/Graphe/micropan/Ypestis/my_micropan.R) qui contient des fonctions de micropan modifi?es pour prendre en compte les pr?fixes de nos prot?ines.
*source("my_micropan.R")

```{r, echo=FALSE, cache=cachedata}
root_dir <- "/autogsi/home/quentin/micropan/Ypestis"
#root_dir <- "C:/Users/quentin/Documents/Enseignement/CoursYQ/Graph/micropan/Ypestis"

setwd(root_dir)
library(micropan)
source("my_micropan.R")

options (stringsAsFactors = FALSE)
```

##Tableau de donn?es
Lecture du tableau d?crivant les souches utilis?es dans l'analyse.

```{r, echo=TRUE, cache=cachedata}
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.

```{r, echo=FALSE, cache=cachedata}
genome_table <- read.table(file=genome_table_file , 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).  

```{r, echo=FALSE, cache=cachedata}
in_files = c();
for (i in 1:nrow(genome_table) ) {
	in_files[i] <- file.path("data/proteins", genome_table$File[i])
	nb_entries <- length(grep('>', readLines(in_files[i])))
	genome_table[i, 5] <- nb_entries
	#cat ("Strain", genome_table$File[i], ':', nb_entries, "entries\n")
}
colnames(genome_table)[5] <- c('Proteome')
kable(genome_table)
```

  Tracer la distribution de la taille des prot?omes (**boxplot()**)

```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
boxplot(genome_table$Proteome, col=color, xlab="Proteome size")
```

##Comparaisons des pairs de prot?omes

Pour chaque pair 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. 

```{r, eval=FALSE, echo=FALSE, cache=cachedata}
out_folder <- "blast"
myblastAllAll(in_files, out_folder, e.value=1,job=1)
```

##Inspection des r?sultats du BlastP All/All

La commande **readBlastTable()** de *micropan* peut ?tre utilis?e pour relire les fichiers en 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*). 

```{r, echo=FALSE, cache=cachedata}
nb_hitsAB <- length(readLines('blast/YpesA_vs_YpesB.txt'))
nb_hitsAB / genome_table[genome_table$File == 'YpesA', 5]

nb_hitsBA <- length(readLines('blast/YpesB_vs_YpesA.txt'))
nb_hitsBA / genome_table[genome_table$File == 'YpesB', 5]
```

  Tracez la distribution du nombre de hits par prot?ine pour le fichier *YpesA_vs_YpesB.txt*.

```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
dataAB <- readBlastTable('blast/YpesA_vs_YpesB.txt')
head(dataAB)
xmax <- ceiling(max(table(dataAB$Query))/50)*50
hist(table(dataAB$Query), nclass = 100, col="tan4", xlim=c(0,xmax), xlab ="nombre de hits", main="Distribution du nombre de hits par prot?ine")
```

  Est-ce que toutes les prot?ines de A ont un hit dans B?
  
```{r, echo=FALSE, cache=cachedata}
genome_table[genome_table$File == 'YpesA', 5] - length(unique(dataAB$Query))
```

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

  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.
  
Premi?re ?tape: extraire la liste des noms des prot?ines.
Seconde ?tape: extraire les annotations Pfam des prot?ines de cette liste.
Troisi?me ?tape: 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()**.

```{r eval=TRUE, cache=cachedata, echo=FALSE}
pfamA <- readHmmer('pfam/YpesA_vs_Pfam-A.hmm.txt')
head(pfamA)
ll <- tail(sort(table(dataAB[,1])), n=200)
pfam_annot <- lapply(names(ll), function(x) pfamA[pfamA$Query == x,])
descriptions <- sort(table(unlist(lapply(pfam_annot, function(x) x$Description[1]))))
```

  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)*)
  
```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
globalpar <- par()
par <- par(mar=c(5.1,20.1,4.1,2.1))
barplot(descriptions, horiz=T, cex.names=0.8, col=color, main="Pfam domains", xlab="Frequence", xlim=c(0,100))
par <- par(mar=c(5.1,4.1,4.1,2.1))
```

##Distances entre les pairs de g?nes

Le package *micropan* a une fonction (**bDist()**) qui permet de calculer les distances entre pairs 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 pairs de prot?ines pour toutes les comparaisons (*blast_files*).

Le *data.frame* est "sauv?e" dans un fichier avec la fonction **save()**.

```{r, echo=FALSE, cache=cachedata}
threshold <- 1e-05
blast_distances10_5_file <- "res/blast_distances_10_5"
if ( file.exists(blast_distances10_5_file)) {
  cat("skip ditances", blast_distances10_5_file, "exists\n")
  load(blast_distances10_5_file)
} else {
  blast_files <- file.path("blast", dir("blast"))
  blast_distances10_5 <- mybDist(blast_files, e.value = threshold)
  save(blast_distances10_5, file = blast_distances10_5_file)
}
head(blast_distances10_5)
```

  Tracez la distribution des distances (**hist()**). 
  
```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center', fig.height=2, fig.width=5}
hist(blast_distances10_5[,3], breaks=100, col=color, xlab = "normalized distances", ylab = "Number of pairs", main="Distribution of distances between gene pairs")
```

  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?
  
```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
tab <- table(blast_distances10_5[,1])
hist(tab, nclass = 1000, col=color)
partial <- hist(tab[tab <= 20], breaks=c(0:20), plot =F)
barplot(partial$count, names.arg=partial$breaks[-1], las=1)

ll <- (tab[tab == 1])
pfam_annot <- lapply(names(ll), function(x) pfamA[pfamA$Query == x,])
lapply(pfam_annot, function(x) x$Description[1])

descriptions <- table(unlist(lapply(pfam_annot, function(x) x$Description[1])), useNA="always")

globalpar <- par()
par <- par(mar=c(5.1,20.1,4.1,2.1))
barplot(descriptions, horiz=T, cex.names=0.8, col=color, main="Pfam domains", xlab="Frequence", xlim=c(0,100))
par <- par(mar=c(5.1,4.1,4.1,2.1))

```

```{r eval=FALSE, cache=cachedata, echo=FALSE}
### en commentaires
ll <- head(sort(table(blast_distances10_5[,1])), n=200)
ll <- tail(sort(table(blast_distances10_5[,1])), n=200)
pfam_annot <- lapply(names(ll), function(x) pfamA[pfamA$Query == x,])
lapply(pfam_annot, function(x) x$Description[1])
#table(unlist(lapply(pfam_annot, function(x) x$Description[1])))
```

##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!). 
  
```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
blast_distances10_5$Similarity <- 1-blast_distances10_5$Distance
head(blast_distances10_5[,c(1,2,4)])

hist(blast_distances10_5[,4], breaks=100, col=color, xlab = "normalized similarities", ylab = "Number of pairs", main="Distribution of similarities between gene pairs")
```

  Ecrire les pairs de similarit?s dans un fichier sous la forme d'une table (**write.table(**).  

```{r eval=TRUE, cache=cachedata, echo=FALSE}
similarity_file10_5 <- "res/blast_similarity10_5.txt"
if ( file.exists(blast_distances10_5_file)) {
  cat("Skip", similarity_file10_5, "is found\n")
} else {
  write.table(blast_distances10_5[,c(1,2,4)], file=similarity_file10_5, append=FALSE, row.names=FALSE, col.names=FALSE, quote=FALSE)
}
```

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

```{r eval=TRUE, cache=cachedata, echo=FALSE}
graph_file <- similarity_file10_5
if ( file.exists(graph_file)) {
  graph <- read.graph(graph_file, format = "ncol", directed = FALSE)
  cat("graphe avec", vcount(graph), "noeuds et", ecount(graph), "ar?tes\n");
} else {
  message <- paste("Stop, le fichier", graph_file, "n'existe pas!\n", sep=" ")
  stop(message)
}
```

  V?rifiez graphiquement les poids des ar?tes (*E(graph)$weight*).
  
```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
hist(E(graph)$weight, breaks=100, col="tan4", las=1, xlab = "edge weights", ylab = "Number of pairs", main="Distribution of edge weights between gene pairs")
```

  Appliquez la fonction **simplify()** sur le graphe.
  Affichez le nombre de vertices et d'edges et comparez avec le graphe initial.
  
```{r eval=TRUE, cache=cachedata, echo=FALSE}
gr <- simplify(graph)
cat("graphe simplifi?", vcount(gr), "noeuds et", ecount(gr), "ar?tes\n");
```

##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 cc.

```{r eval=TRUE, cache=cachedata, echo=FALSE}
cc <- decompose.graph(gr, mode="strong", min.vertices = 1);
cat(length(cc), "composantes connexes dans le graphe\n");
```

Le r?sultat cc est la liste des CC du graphe (cc[[i]] est la i?me CC)

##Propri?t?es des composantes connexes

  Calculez le nombre de nouds 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()**). 

```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
cc_property <- cbind(unlist(lapply(cc, vcount)), unlist(lapply(cc, ecount)))
colnames(cc_property) <- c("vertex", "edges")

vmax <- 500 
subset <- cc_property[cc_property[,1] <= vmax, ]
rank <- which(subset[, 2] == max(subset[, 2]))
xlimit <- subset[rank, 1]
plot(cc_property[cc_property[,1] <= vmax, ], las=1, xlim=c(0, vmax), pch=20, col="coral")
curve(x*(x-1)/2, from = 0, to = xlimit, n = 101, col="purple", add=T)
```

```{r eval=FALSE, cache=cachedata, echo=FALSE, fig.align='center'}
# en commentaire
clique <- unlist(lapply(cc, function(y) {
	x <- vcount(y)
	x*(x-1)/2}
	)
)

cc_property2 <- cbind(unlist(lapply(cc, vcount)), unlist(lapply(cc, ecount)), clique)

ecart <- cc_property2[,3] - cc_property2[,2]
hist(ecart[ecart <= 100], breaks=100, col="tan4", las=1)
length(ecart[ecart == 0])
# 3016
clique.number(gr)
#1158
```

  Repr?sentez graphiquement une clique (cc[[5]] par exemple).
  
```{r eval=FALSE, cache=cachedata, echo=TRUE, fig.align='center'}
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?

```{r eval=FALSE, cache=cachedata, echo=TRUE}
count_paralogs(V(mycc)$name)
```

##Distribution de la paralogie

  Repr?sentez la distribution de la paralogie dans les CC.
  
```{r eval=TRUE, cache=cachedata, echo=TRUE, fig.align='center'}
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")
	}
}
partial <- hist(pdm$Pcount, nclass=10, plot =F)
barplot(partial$count, names.arg=partial$breaks[-1], xlab="Pcount")
```

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

```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
plot(pdm[,c(1,2)], las=1, pch=20, col=color)
```

  Quelle est le taux minimum de paralogie en fonction du nombre de prot?ines par CC?
  Tracez la courbe correspondante.
  
```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
# x-length(strain_list))/x
vmax <- max(pdm$Vcount)
plot(pdm[,c(1,2)], las=1, pch=20, col=color)
curve((x-nrow(genome_table))/x, from=nrow(genome_table)+1, to=vmax, n=101, col="purple", add=T)
```

  Que pouvez-vous dire des CC localis?es sur la courbe?

```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
hist(pdm$Vcount%%nrow(genome_table), breaks=c(0:11), col=color, xlab="Vcount%%nrow(genome_table)", xlim=c(0, nrow(genome_table)))
```

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

```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
vmax = 3000
pdm <- data.frame(Vcount=numeric(), Pcount=numeric(), Density=numeric(), Modularity=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 > 12) {
		j <- j+1
		list_name <- V(CC)$name
		fastgreedy <- cluster_fast_greedy(CC, modularity=TRUE)
		
  	pdm[j, 1] <- vcount(cc[[i]])
		pdm[j, 2] <- count_paralogs(list_name)
		pdm[j, 3] <- graph.density(cc[[i]],loops=FALSE)
		pdm[j, 4] <- max(fastgreedy$modularity)
		#cat("CC", cc_num, " ", vcount(CC), " ", pdm[j, 2], " ", pdm[j, 4], "\n")
		
		if ( pdm[j, 4] > 2 ) {
			plot(fastgreedy, CC, vertex.label=NA, vertex.size=3, edge.color="red", edge.width=E(CC)$weight)
		}
	}
}
head(pdm)
hist(pdm$Modularity, nclass=10, col=color)
```

  Tracez la modularit? des CC en fonction du taux de paralogie.
  
```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
plot(pdm[,c(2,4)], las=1, pch=20, col=color)
```

  Tracez la droite de r?gression *Modularity* en fonction de *Pcount*.

```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center', fig.height=2, fig.width=5}
cor <- lm(pdm$Modularity ~ pdm$Pcount)
summary(cor)
plot(pdm[,c(2,4)], las=1, pch=20, xlim=c(0,1), col=color)
#abline(coef(cor), xlim=c(0,1))
curve(cor$coefficients[2]*x +cor$coefficients[1], from=0, to=1, col="purple", add=T)
```

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

```{r eval=TRUE, cache=cachedata, echo=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")
		}
	}
}
```

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

```{r eval=TRUE, cache=cachedata, echo=TRUE}
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(m?thode)*.

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

```{r eval=TRUE, cache=cachedata, echo=TRUE}
layout(matrix(1:3,1,3))
par <- par(mar=c(1,1,1,1))

vmin = 24
vmax = 50
pdm <- data.frame(Vcount=numeric(), Pcount=numeric())
Modularity <- data.frame(fastgreedy=numeric(), walktrap=numeric(), edge=numeric())
Size <- 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)
  	  Size[j, 1] <- length(fastgreedy)
  	  Size[j, 2] <- length(walktrap)
	    Size[j, 3] <- length(edge)
	    
		  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()
      }
  		
  		
	    #cat("CC", cc_num, vcount(CC), Size[j, 1], Size[j, 2], Size[j, 3], "\n")
		}
	}
}
par <- par(mar=c(5.1,4.1,4.1,2.1))
layout(1)
```

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

```{r eval=TRUE, cache=cachedata, echo=TRUE}
Size_stat <- apply(Size, 1, function(x) rank(signif(x[1:3], digits=5), ties.method = "max"))

Size_tab <- apply(Size_stat, 1, function(x) { 
  his <- hist(x,breaks=(0:3), plot=F)
  his$count
})
barplot(Size_tab)
```

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:
```{r eval=FALSE, cache=cachedata, echo=TRUE}
  my_gamma <- 1
  my_spins <- 25
  stop_temp <- 0.1
  start_temp <- 1
  cool_fact <- 0.95
  update.rule <- "config"
  iter <- 10
  
  # test 
  CC <- cc[[1216]]
  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)
  	spinglass
  	cluster_fast_greedy(CC, modularity=TRUE)
    Qiter[k] <- modularity(spinglass)
    spinglass_iter[[k]] <- spinglass
  }
  step <- which.max(Qiter)
  membership <- membership(spinglass_iter[[step]])
  modularity <- round(max(Qiter), 3)
```

```{r eval=FALSE, cache=cachedata, echo=TRUE}
par(mar=c(1,1,1,1))
layout(matrix(1:3, 1, 3))
vmin = 24
vmax = 50
pdm <- data.frame(Vcount=numeric(), Pcount=numeric())
Modularity <- data.frame(fastgreedy=numeric(), walktrap=numeric(), edge=numeric(), spinglass=numeric())
Size <- 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(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)
  		
      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]])
      
    	pdm[j, 1] <- cc_num
    	pdm[j, 2] <- vcount(cc[[i]])
  		pdm[j, 3] <- count_paralogs(list_name)
  		
  		Modularity[j, 1] <- round(max(fastgreedy$modularity), 3)
  		Modularity[j, 2] <- round(max(walktrap$modularity), 3)
  		Modularity[j, 3] <- round(max(edge$modularity), 3)
    	Modularity[j, 4] <- round(max(Qiter), 3)
  	  Size[j, 1] <- length(fastgreedy)
  	  Size[j, 2] <- length(walktrap)
	    Size[j, 3] <- length(edge)
      Size[j, 4] <- length(unique(membership))
  		
  		cat("CC", cc_num, " ", vcount(CC), " ", Modularity[j, 1], " ", Modularity[j, 2], " ", Modularity[j, 3], " ", Modularity[j, 4], "\n")
		}
	}
}
par <- par(mar=c(5.1,4.1,4.1,2.1))
layout(1)
```

  Pour chaque CC triez les methodes en fonction des modularit?s obtenues. 
  Repr?sentez graphiquement les r?sultats.

```{r eval=FALSE, cache=cachedata, echo=TRUE}
Modu_stat <- apply(Modularity, 1, function(x) rank(signif(x[1:3], digits=5), ties.method = "max"))

Modu_tab <- apply(Modu_stat, 1, function(x) { 
  his <- hist(x,breaks=(0:3), plot=F)
  his$count
})
barplot(Modu_tab)
```

  Pour chaque CC triez les methodes en fonction du nombre de communaut?es obtenues. 
  Repr?sentez graphiquement les r?sultats.

```{r eval=FALSE, cache=cachedata, echo=TRUE}
Size_stat <- apply(Size, 1, function(x) rank(signif(x[1:4], digits=5), ties.method = "max"))

Size_tab <- apply(Size_stat, 1, function(x) { 
  his <- hist(x,breaks=(0:4), plot=F)
  his$count
})
barplot(Size_tab)
```


```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center'}
vmin = 12
vmax = 60
ccmax = length(cc)
#ccmax = 20
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
  	}
	}
}
```

```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center', fig.height=2, fig.width=5}
hist(gene_classes$Vcount, col=color, nclass=10)
```

```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center', fig.height=2, fig.width=5}
hist(gene_classes$Pcount, col=color, nclass=10, xlim=c(0,1))
```
```{r eval=TRUE, cache=cachedata, echo=FALSE, fig.align='center', fig.height=2, fig.width=5}
hist(gene_classes$Strains, col=color, nclass=12, xlim=c(1,12))
```

