update on figures number, and labels of solutions in the plot, to be...

update on figures number, and labels of solutions in the plot, to be consistent with final manusctipt

statistical-analysis.Rmd

@@ -23,6 +23,25 @@ library(gridExtra)
 set.seed(13)
 
 geom_Mean = function(x){prod(x^(1/length(x)))}
+
+get.custom.labels=function(name){
+  my.labels = list("D"="D",
+                   "DNASHIELD"="DNA/RNA shield",
+                   "ETOH"="95% Ethanol",
+                   "FECALSWAB"="Fecal swab",
+                   "MTM"="MTM",
+                   "NORGEN"="Norgen",
+                   "OMNIGEN"="Omnigene",
+                   "RNALATER"="RNA later",
+                   "S"="S",
+                   "STRATEC"="Stratec",
+                   "TRISEDTA"="Tris EDTA",
+                   "WHATMAN"="FTA card")
+  return(unlist(ifelse(name %in% names(my.labels),
+                       my.labels[name],
+                       name)))
+  
+}
 ```
 
 # Présentation de l'étude
@@ -83,7 +102,8 @@ plot.ratio = ggplot(data=ratios.melt, aes(x=solution, y=ratio)) +
   theme_bw() + 
   theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                       hjust = 1))+
-  geom_rect(data=normal.range, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), fill = "green", alpha=0.2, inherit.aes=FALSE)
+  geom_rect(data=normal.range, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), fill = "green", alpha=0.2, inherit.aes=FALSE)+
+  scale_x_discrete(labels = get.custom.labels)
 plot.ratio
 
 concentrations.melt =  melt(t(concentrations))
@@ -93,7 +113,8 @@ plot.conc = ggplot(data=concentrations.melt, aes(x=solution, y=concentrations))
   geom_jitter(shape=16, position = position_jitter(0.2), col="gray10") + 
   theme_bw() + 
   theme(axis.text.x = element_text(angle = 45, vjust = 1, 
-                                      hjust = 1))
+                                      hjust = 1))+
+  scale_x_discrete(labels = get.custom.labels)
 
 plot.conc
 
@@ -200,8 +221,6 @@ ps
 ## Filtre sur les échantillons
 On retire les échantillons qui ont moins de 30 000 reads.
 ```{r, fig.height = 10, fig.width = 10, fig.align = "center"}
-
-
 ps.filt = prune_samples(sample_sums(ps)>=30000, ps)
 
 # Si on veut retirer les outliers
@@ -209,16 +228,11 @@ col=rep("black", nsamples(ps.filt))
 col[sample_data(ps.filt)$Patient %in% c("SS01", "SS02", "SS03", "SS07", "SS11", "SS13", "SS12", "SS14")] = "gray70"
 col[sample_names(ps.filt) %in% c("SS01Sc1-3", "SS02Sc1-3", "SS04Dc3-3")] = "red"
 
-plot(ape::as.phylo(hclust(phyloseq::distance(rarefy_even_depth(ps.filt, rngseed = T), method = "bray"))), type="radial", cex = 0.5, tip.color = col)
-
-ps.filt = prune_samples(!sample_names(ps.filt) %in% c("SS01Sc1-3", "SS02Sc1-3", "SS04Dc3-3"), ps.filt)
-```
-
-
-```{r include=FALSE}
 pdf(file = "figures-article/figureS5.pdf", width = 9, height = 9, pointsize = 10)
 plot(ape::as.phylo(hclust(phyloseq::distance(rarefy_even_depth(ps.filt, rngseed = T), method = "bray"))), type="radial", cex = 0.5, tip.color = col)
 dev.off()
+
+ps.filt = prune_samples(!sample_names(ps.filt) %in% c("SS01Sc1-3", "SS02Sc1-3", "SS04Dc3-3"), ps.filt)
 ```
 ## Filtre sur les OTUs
 On ne garde que les OTUs présents à une proportion suppérieure à $10^{-4}$ dans au moins 5% des échantillons. 
@@ -267,11 +281,31 @@ Il s'agit de donner une première impréssion sur les données. Nous utiliserons
 L'indice de Bray-Curtis est une mesure de distance très utilisée en écologie, qui se base des données de comptage (nombre de reads associés à chaque OTUs, il faut donc utiliser les données raréfiées). Il renseigne sur les différences quantitatives. 
 
 ```{r, fig.height=15, fig.width=20}
+
+sample_data(ps.rarefied)$type = "Stabilized"
+sample_data(ps.rarefied)$type[sample_data(ps.rarefied)$Solution_agg == "Dc"] = "Dc replicates"
+sample_data(ps.rarefied)$type[sample_data(ps.rarefied)$Solution_agg == "Sc"] = "Sc replicates"
+sample_data(ps.rarefied)$type[sample_data(ps.rarefied)$Solution_agg == "REF"] = "Reference"
+sample_data(ps.rarefied)$label = as.character(sample_data(ps.rarefied)$Solution_agg)
+sample_data(ps.rarefied)$label[sample_data(ps.rarefied)$Solution_agg == "Dc"] = ""
+sample_data(ps.rarefied)$label[sample_data(ps.rarefied)$Solution_agg == "Sc"] = ""
+sample_data(ps.rarefied)$label[sample_data(ps.rarefied)$Solution_agg == "REF"] = ""
+sample_data(ps.rarefied)$label = get.custom.labels(sample_data(ps.rarefied)$label)
+
 dist.bc = phyloseq::distance(ps.rarefied, method = "bray")
 PCoA.bc  = ordinate(ps.rarefied, "PCoA", distance=dist.bc)
-plot_ordination(ps.rarefied, PCoA.bc, color = "Solution_agg") + 
-  facet_wrap(~Patient, scales = "free" , ncol = 5) + 
-  geom_text_repel(aes(label = Solution), color = 'black', size = 3) 
+PCoA.bc.plot = plot_ordination(ps.rarefied, PCoA.bc, color = "type") + 
+  geom_point(size=2.5)+
+  scale_colour_manual(values = c("red", "blue", "green", "black")) + 
+  facet_wrap(~Patient, scales = "free" , ncol = 3) + 
+  geom_text_repel(aes(label =label), color = 'black', size = 3) 
+
+PCoA.bc.plot
+
+ggsave(filename = "figures-article/figureS5.pdf",
+       plot = PCoA.bc.plot,
+       width = 15,
+       height=20)
 ```
 
 # Effet de la dilution
@@ -481,7 +515,9 @@ analyse_distance = function(phyloseq, method, title){
                                       hjust = 1))+
     geom_text(aes(label =signif), cex = 5) +
     coord_fixed() + 
-    ggtitle("B")
+    ggtitle("B")+
+  scale_x_discrete(labels = get.custom.labels)+
+  scale_y_discrete(labels = get.custom.labels)
    
   # we sort the levels of Solution in the order of median performance (ggplot will use this order)
   classement = order(apply(perf.matrix,2, FUN = median, na.rm = T))
@@ -528,7 +564,8 @@ analyse_distance = function(phyloseq, method, title){
                  min(dist.inter.patient),min(dist.inter.aliquot)),
          1.1*max(max(perf.matrix.melt$dist2ref, na.rm=T),
                  max(dist.inter.patient),max(dist.inter.aliquot)))+
-    ggtitle("A")
+    ggtitle("A")+
+  scale_x_discrete(labels = get.custom.labels)
   # we sort back the levels of Solution in the alphabetical order (to have "D" first and used as reference in stat models)
   perf.matrix.melt$Solution = factor(perf.matrix.melt$Solution, 
                                      levels = sort(levels(perf.matrix.melt$Solution)))
@@ -650,6 +687,8 @@ wilcox.phyla.siginif[wilcox.phyla.fdr<0.05]='*'
 wilcox.phyla.siginif[wilcox.phyla.fdr<0.01]='**'
 wilcox.phyla.siginif[wilcox.phyla.fdr<0.001]='***'
 
+colnames(phyla.matrix.logVsRef) = get.custom.labels(colnames(phyla.matrix.logVsRef))
+
 pheatmap(phyla.matrix.logVsRef, 
          display_numbers = wilcox.phyla.siginif, 
          na_col = "grey50", 
@@ -661,7 +700,9 @@ plot_bar(select_top_taxa(transform_sample_counts(ps.phyla.merged,
          fill="p") +
   geom_bar(aes(color=p, fill=p), stat="identity") + 
   theme_bw()+ 
-  theme(axis.text.x = element_text(angle = 90, hjust = 1)) 
+  theme(axis.text.x = element_text(angle = 90, hjust = 1))  +
+  scale_x_discrete(labels = get.custom.labels)
+
 ```
 
 ```{r include=FALSE}
@@ -727,6 +768,8 @@ wilcox.genus.siginif[wilcox.genus.fdr<0.05]='*'
 wilcox.genus.siginif[wilcox.genus.fdr<0.01]='**'
 wilcox.genus.siginif[wilcox.genus.fdr<0.001]='***'
 
+colnames(genus.matrix.logVsRef) = get.custom.labels(colnames(genus.matrix.logVsRef))
+
 pheatmap(genus.matrix.logVsRef, 
          display_numbers = wilcox.genus.siginif, 
          na_col = "grey50", 
@@ -740,7 +783,7 @@ pheatmap(genus.matrix.logVsRef,
 
 
 ```{r include=FALSE}
-pdf(file = "figures-article/figureS6.pdf", width = 9, height = 9, pointsize = 10)
+pdf(file = "figures-article/figureS7.pdf", width = 9, height = 9, pointsize = 10)
 pheatmap(genus.matrix.logVsRef, 
          display_numbers = wilcox.genus.siginif, 
          na_col = "grey50", 
@@ -761,7 +804,8 @@ plot_bar(select_top_taxa(transform_sample_counts(ps.genus.merged,
          fill="g") +
   geom_bar(aes(color=g, fill=g), stat="identity") +
   theme_bw() + 
-  theme(axis.text.x = element_text(angle = 90, hjust = 1))
+  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
+  scale_x_discrete(labels = get.custom.labels)
 
 
 rownames(genus.matrix.logVsRef) = old_rownames