diff --git a/R_amplicon_analysis b/R_amplicon_analysis
index 04f8cceb5d0fbda429e11c1e179f0204f9f7a65b..736cbd06d0a3abaec074b8deda11f1f1bacd6d37 100644
--- a/R_amplicon_analysis
+++ b/R_amplicon_analysis
@@ -27,8 +27,12 @@ library("caTools")
 library("tcltk")
 library("vegan")
 library("ggplot2")
+library("plyr")
 
 
+setwd("//pommard/agro-users$/evanschaik/Mes documents/Solca/Amplicon/16s")
+a=read.table("betad.txt", sep = "\t", header = T,stringsAsFactors = F)
+b=read.table("essai_treat.txt", header = T,stringsAsFactors = F, sep="\t")
 
 
 
@@ -91,48 +95,48 @@ add.alpha <- function(col, alpha=0.3){
 ###chi square test with only 2 variables
 
 chi_sam<-function(b)
+  #########################################
 #########################################
-#########################################
-{
-t<-apply(b,1,sum)
-selec<-which(t>0)
-b<-b[selec,]
-t<-apply(b,2,sum)
-selec<-which(t>0)
-b<-b[,selec]
-y<-chisq.test(b)
-print(y)
-print(y$expected)
-print(y$observed)
-print (y$p.value)
-if(y$p.value<0.05)
-{
-if(dim(b)[2]>2)
-{
-z<-combn(dim(b)[2],2)
-for(i in 1:dim(z)[2])
-{
-r<-b[,c(z[1,i],z[2,i])]
-s<-apply(r,1,sum)
-selec<-which(s>0)
-r<-r[selec,]
-x<-chisq.test(r)
-#print(paste(colnames(b)[z[1,i]],colnames(b)[z[2,i]],sep=" vs "))
-#print(x$expected)
-#print(x$observed)
-if(x$p.value<(0.05/dim(z)[2]))
-{
-print(paste(colnames(b)[z[1,i]],colnames(b)[z[2,i]],sep=" vs "))
-print(x$p.value)
-}
-}
-print(paste("corrected p-value ",(0.05/dim(z)[2]),sep=""))
-}
-else
 {
-print (y$p.value)
-}
-}
+  t<-apply(b,1,sum)
+  selec<-which(t>0)
+  b<-b[selec,]
+  t<-apply(b,2,sum)
+  selec<-which(t>0)
+  b<-b[,selec]
+  y<-chisq.test(b)
+  print(y)
+  print(y$expected)
+  print(y$observed)
+  print (y$p.value)
+  if(y$p.value<0.05)
+  {
+    if(dim(b)[2]>2)
+    {
+      z<-combn(dim(b)[2],2)
+      for(i in 1:dim(z)[2])
+      {
+        r<-b[,c(z[1,i],z[2,i])]
+        s<-apply(r,1,sum)
+        selec<-which(s>0)
+        r<-r[selec,]
+        x<-chisq.test(r)
+        #print(paste(colnames(b)[z[1,i]],colnames(b)[z[2,i]],sep=" vs "))
+        #print(x$expected)
+        #print(x$observed)
+        if(x$p.value<(0.05/dim(z)[2]))
+        {
+          print(paste(colnames(b)[z[1,i]],colnames(b)[z[2,i]],sep=" vs "))
+          print(x$p.value)
+        }
+      }
+      print(paste("corrected p-value ",(0.05/dim(z)[2]),sep=""))
+    }
+    else
+    {
+      print (y$p.value)
+    }
+  }
 }
 
 
@@ -154,17 +158,20 @@ ensemble<-function(a,b,name,threshold)
   sum_a<-apply(a,2,sum)
   if(length(which(sum_a==0))>0)
   {
-  a<-a[,-which(sum_a==0)] 
+    a<-a[,-which(sum_a==0)] 
   }
   #heatmap
   if(dim(a)[2]<threshold)
   {
     aprime<-t(a)
-  } else {
+  }
+  else
+  {
     aprime<-t(Ten_max(t(a),threshold))
   }
+  
   write.table(aprime, paste(Sys.Date(),name,"Ten_max",threshold,".txt",sep="_"), row.names=T, sep="\t",quote=F)
-
+  
   logaprime <- log2(aprime)
   logaprime[logaprime=="-Inf"] <- c(0)
   bprime <- data.frame(unclass(b))
@@ -205,12 +212,7 @@ ensemble<-function(a,b,name,threshold)
   
   #dev.off()
   
-  if(dim(a)[2]<threshold)
-  {
-    aprime<-t(Ten_max(t(a),dim(a)[2]))
-  } else {
-    aprime<-t(Ten_max(t(a),threshold))
-  }
+
   v<-apply(aprime,1,sum)
   w<-apply(a,1,sum)
   pourseq<-mean(v/w)*100
@@ -229,28 +231,36 @@ ensemble<-function(a,b,name,threshold)
     z=1
     y=1
     pdf(file=paste("boxplot",Sys.Date(),name,colnames(b)[x],y,".pdf",sep="_"),width=15,height=15)
-    par(mfrow=c(3,2))
+    par(mfrow=c(3,2), mar=c(7,4,6,4)+2)
     for(m in 1:dim(aprime)[2])
     {
-      test <- kruskal.test(aprime[,m]~bprime[,x])
-      if(test$p.value<0.05)
-      {
-      if(z%%6==0)
-      {
-        boxplot(aprime[,m]~bprime[,x],main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3)  
-        dev.off()
-        y=y+1
-        z=1
-        pdf(file=paste("boxplot",Sys.Date(),name,colnames(b)[x],y,".pdf",sep="_"),width=30,height=15)
-        par(mfrow=c(3,2))
-      }
+      condition <- length(which(!is.na(unique(bprime[,x]))==1))
+      if(condition[1]>1)
+        {
+        test <- kruskal.test(aprime[,m]~bprime[,x])
+        if(test$p.value<0.05)
+          {
+          if(z%%6==0)
+            {
+            boxplot(aprime[,m]~bprime[,x],main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3,las=2)  
+            dev.off()
+            y=y+1
+            z=1
+            pdf(file=paste("boxplot",Sys.Date(),name,colnames(b)[x],y,".pdf",sep="_"),width=30,height=15)
+            par(mfrow=c(3,2), mar=c(7,4,6,4)+2)
+            }
+          else
+            {
+            boxplot(aprime[,m]~bprime[,x],main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3,las=2)  
+            z <- z+1
+            }
+          }
+        }
       else
-      {
-        boxplot(aprime[,m]~bprime[,x],main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3)  
-        z <- z+1
-      }
+        {
+        print("kruskal not allowed")
+        }
       }
-    }
     dev.off()
     
     
@@ -269,16 +279,16 @@ ensemble<-function(a,b,name,threshold)
     d <- c()
     for(i in 1:dim(c)[1]){d[i] <- paste(c[i,c(1)],collapse="_")}
     rownames(temp) <- d
+    if(dim(t(temp))[2]>1)
+    {
     chi_sam(t(temp))
-    
-    
     pdf(file=paste("BGA",Sys.Date(),name,colnames(bprime)[x],".pdf",sep="_"),width=15,height=15)
     {
       bga1 <- bca(x = datas, fac = as.factor(bprime[,x]), scannf = FALSE, nf = 5)
       
       s.class(datas$li,
               fac = as.factor(bprime[,x]),  # colorer par groupes
-              col = c("black","blue","red","green","orange","yellow","grey","darkblue"))
+              col = c("black","blue","red","darkgreen","orange","purple","grey","darkblue"))
     }
     dev.off()
     
@@ -316,7 +326,7 @@ ensemble<-function(a,b,name,threshold)
     
     print(graph)
     dev.off()
-    
+    }
   }
   
   
@@ -329,41 +339,38 @@ ensemble<-function(a,b,name,threshold)
   {
     temp1 <- comb[1,x]
     temp2 <- comb[2,x]
+    print(paste(colnames(b)[temp1],colnames(b)[temp2],sep="_"))
     {
-      if(dim(a)[2]<threshold)
-      {
-        aprime<-t(Ten_max(t(a),dim(a)[2]))
-      }
-      else
-      {
-        aprime<-t(Ten_max(t(a),threshold))
-      }
       grouping <- as.factor(paste(b[,temp1],b[,temp2], sep="_"))
       z=1
       y=1
       pdf(file=paste("boxplot",Sys.Date(),name,colnames(b)[temp1],colnames(b)[temp2],y,".pdf",sep="_"),width=15,height=15)
-      par(mfrow=c(3,2))
+      par(mfrow=c(3,2), mar=c(7,4,6,4)+2)
       for(m in 1:dim(aprime)[2])
       {
-      if(test$p.value<0.05)
-       {
-      
-        test <- kruskal.test(aprime[,m]~grouping)
-        if(z%%6==0)
+        condition <- length(which(!is.na(unique(grouping))==1))
+        if(condition[1]>1)
         {
-          boxplot(aprime[,m]~grouping,main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3)  
-          dev.off()
-          y=y+1
-          z=1
-          pdf(file=paste("boxplot",Sys.Date(),name,colnames(b)[temp1],colnames(b)[temp2],y,".pdf",sep="_"),width=30,height=15)
-          par(mfrow=c(3,2))
-        }
-        else
+        test<-kruskal.test(aprime[,m]~grouping)
+        if(test$p.value<0.05)
         {
-          boxplot(aprime[,m]~grouping,main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3)  
-          z <- z+1
+          
+          if(z%%6==0)
+          {
+            boxplot(aprime[,m]~grouping,main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3,las=2)  
+            dev.off()
+            y=y+1
+            z=1
+            pdf(file=paste("boxplot",Sys.Date(),name,colnames(b)[temp1],colnames(b)[temp2],y,".pdf",sep="_"),width=30,height=15)
+            par(mfrow=c(3,2), mar=c(7,4,6,4)+2)
+          }
+          else
+          {
+            boxplot(aprime[,m]~grouping,main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3,las=2)  
+            z <- z+1
+          }
+        }
         }
-	}
       }
       dev.off()
       
@@ -373,8 +380,8 @@ ensemble<-function(a,b,name,threshold)
         
         s.class(datas$li,
                 fac = grouping,  
-                col = c("black","blue","red","green","orange","yellow","grey","darkblue"))
-      }
+                col = c("black","blue","red","darkgreen","orange","purple","grey","darkblue"))
+                }
       dev.off()
       
       c<-aggregate(a, list(a=b[,temp1],b=b[,temp2]), FUN=sum)
@@ -384,9 +391,9 @@ ensemble<-function(a,b,name,threshold)
       d <- c()
       for(i in 1:dim(temp3)[1]){d[i] <- paste(temp3[i,c(1:2)],collapse="_")}
       rownames(c) <- d
-      chi_sam(t(c))
-      
-      
+      if(dim(t(c))[2]>1)
+      {
+        chi_sam(t(c))
       #NMDS
       
       
@@ -419,11 +426,11 @@ ensemble<-function(a,b,name,threshold)
       
     }
     dev.off()
+    }
   }
   
   
   
-  
   #########################################################
   #en foncion de 3 parametres
   #########################################################
@@ -433,37 +440,35 @@ ensemble<-function(a,b,name,threshold)
     temp1 <- comb[1,x]
     temp2 <- comb[2,x]
     temp3 <- comb[3,x]
-    if(dim(a)[2]<20)
-    {
-      aprime<-t(Ten_max(t(a),dim(a)[2]))
-    } else {
-      aprime<-t(Ten_max(t(a),threshold))
-    }
     grouping <- as.factor(paste(b[,temp1],b[,temp2],b[,temp3], sep="_"))
     
     z=1
     y=1
     pdf(file=paste("boxplot",Sys.Date(),name,colnames(b)[temp1],colnames(b)[temp2],colnames(b)[temp3],y,".pdf",sep="_"),width=15,height=15)
-    par(mfrow=c(3,2))
+    par(mfrow=c(3,2), mar=c(7,4,6,4)+2)
     for(m in 1:dim(aprime)[2])
     {
-      test <- kruskal.test(aprime[,m]~grouping)
+      condition <- length(which(!is.na(unique(bprime[,x]))==1))
+      if(condition[1]>1)
+      {
+        test <- kruskal.test(aprime[,m]~grouping)
       if(test$p.value<0.05)
       {
-      	if(z%%6==0)
-      	{
-        	boxplot(aprime[,m]~grouping,main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3)  
-        	dev.off()
-        	y=y+1
-        	z=1
-        	pdf(file=paste("boxplot",Sys.Date(),name,colnames(b)[temp1],colnames(b)[temp2],colnames(b)[temp3],y,".pdf",sep="_"),width=30,height=15)
-        	par(mfrow=c(3,2))
-      	}
-      	else
-      	{
-        	boxplot(aprime[,m]~grouping,main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.4)  
-        	z <- z+1
-      	}
+        if(z%%6==0)
+        {
+          boxplot(aprime[,m]~grouping,main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.3,las=2)  
+          dev.off()
+          y=y+1
+          z=1
+          pdf(file=paste("boxplot",Sys.Date(),name,colnames(b)[temp1],colnames(b)[temp2],colnames(b)[temp3],y,".pdf",sep="_"),width=30,height=15)
+          par(mfrow=c(3,2), mar=c(7,4,6,4)+2)
+        }
+        else
+        {
+          boxplot(aprime[,m]~grouping,main=paste(colnames(aprime)[m],"_kruskall_test=",formatC(test$p.value, format = "e", digits = 2),sep=""), cex.names=0.4, las=2)  
+          z <- z+1
+        }
+      }
       }
     }
     dev.off()
@@ -474,7 +479,7 @@ ensemble<-function(a,b,name,threshold)
       
       s.class(datas$li,
               fac = as.factor(grouping),  
-              col = c("black","blue","red","green","orange","yellow","grey","darkblue"))
+              col = c("black","blue","red","darkgreen","orange","purple","grey","darkblue"))
     }
     dev.off()
     
@@ -485,9 +490,9 @@ ensemble<-function(a,b,name,threshold)
     d <- c()
     for(i in 1:dim(temp4)[1]){d[i] <- paste(temp4[i,c(1:3)],collapse="_")}
     rownames(c) <- d
-    chi_sam(t(c))
-    
-    
+    if(dim(t(c))[2]>1)
+    {
+      chi_sam(t(c))
     #NMDS
     
     
@@ -519,6 +524,6 @@ ensemble<-function(a,b,name,threshold)
     print(graph)
     
     dev.off()
+    }
   }
-  
 }