Skip to content

GitLab

Commit 29013a53 authored by magali's avatar magali
Browse files

fichier rds et html rotenone

parent 45d0eaa3
Pipeline #44571 passed with stage
in 3 minutes and 29 seconds
Show whitespace changes
Inline Side-by-side
202111_Rotenone/report_Rotenone.Rmd
View file @ 29013a53
...@@ -91,11 +91,7 @@ git pull ...@@ -91,11 +91,7 @@ git pull
# Rappel du contexte de l'étude # Rappel du contexte de l'étude
Du microbiote d'enfants atteints de troubles du spectre autistique a été donné à des souris axéniques, afin de voir l'impact de ce transfert sur leur comportement mais aussi sur divers paramètres biochimiques Deux groupes de 10 souris ont été traitées ou non avec de la rotenone (modèle de Parkinson)
Chez les donneurs, il y avait 4 groupes composés chacun de 4 enfants
Un pool des 4 échantillons d'un même groupe a été fait pour le donner aux souris
# Metabarcoding analysis # Metabarcoding analysis
...@@ -207,7 +203,8 @@ Les adaptateurs Illumina sont présents à différents niveaux pour certains éc ...@@ -207,7 +203,8 @@ Les adaptateurs Illumina sont présents à différents niveaux pour certains éc
Ces séquences seront supprimées plus tard avec FROGS.</div> Ces séquences seront supprimées plus tard avec FROGS.</div>
<div class="alert alert-success" role="alert">En résumé, à part pour l'échantillon 2.1, les séquences sont de très bonnes qualités</div>
<div class="alert comment">En résumé, à part pour l'échantillon 2.1, les séquences sont de très bonnes qualités</div>
```{bash} ```{bash}
#Creer archive #Creer archive
...@@ -420,6 +417,7 @@ scp mmonnoye@migale.jouy.inrae.fr:~/save/analyses_16s/202111_Rotenone/data.tar.g ...@@ -420,6 +417,7 @@ scp mmonnoye@migale.jouy.inrae.fr:~/save/analyses_16s/202111_Rotenone/data.tar.g
# Metagenomic phyloseq analysis # Metagenomic phyloseq analysis
## Import packages et données ## Import packages et données
```{r load-packages, eval=TRUE, cache = FALSE} ```{r load-packages, eval=TRUE, cache = FALSE}
...@@ -444,16 +442,16 @@ library(ape) ...@@ -444,16 +442,16 @@ library(ape)
if(file.exists("frogs.data.rds")){ if(file.exists("frogs.data.rds")){
frogs.data <- readRDS("frogs.data.rds") frogs.data <- readRDS("frogs.data.rds")
}else{ }else{
frogsBiom <- "cleaned_biom-2021-10-19.biom" frogsBiom <- "cleaned_biom-2021-11-17.biom"
frogs.data.biom <- import_frogs(frogsBiom, taxMethod = "blast") frogs.data.biom <- import_frogs(frogsBiom, taxMethod = "blast")
metadata<-read.table("metadatasBALBc.txt",row.names=1, header=T) metadata<-read.table("metadatasRotenone.txt",row.names=1, header=T)
sample_data(frogs.data.biom)<-metadata sample_data(frogs.data.biom)<-metadata
sample_names(frogs.data.biom) <- get_variable(frogs.data.biom, "Name") sample_names(frogs.data.biom) <- get_variable(frogs.data.biom, "Name")
frogs.data<-frogs.data.biom frogs.data<-frogs.data.biom
treefile<- read.tree("tree.nhx") treefile<- read.tree("FROGS/tree.nhx")
phy_tree(frogs.data) <- treefile phy_tree(frogs.data) <- treefile
tax_table(frogs.data) <- gsub("\"", "", tax_table(frogs.data)) tax_table(frogs.data) <- gsub("\"", "", tax_table(frogs.data))
...@@ -478,43 +476,28 @@ qui contient les métadonnées suivantes: ...@@ -478,43 +476,28 @@ qui contient les métadonnées suivantes:
```{r metadata, eval=TRUE, echo=TRUE} ```{r metadata, eval=TRUE, echo=TRUE}
sample_variables(frogs.data) sample_variables(frogs.data)
``` ```
Composition des groupes:
* AG : enfants atteints de TSA et ayant aussi des troubles gastro-intestinaux ## Sélection des données
* S-AG: leurs frères et sœurs comme contrôles <div class="alert alert-danger" role="alert">Je supprime l'échantillon 2.1 car il possède trop peu de reads</div>
* A: Enfants atteints de TSA sans troubles gastro-intestinaux
* S-A: leurs frères et sœurs comme contrôles ```{r sample selection,fig.width=10,fig.height=4, eval=TRUE}
#Je supprime les échantillons TF51_2, T026, TF50_2 et T023_2 qui ont un nombre de séquences < 5000
frogs.data <- subset_samples(frogs.data, !(sample_names(frogs.data) %in% c("2.1")) )
Pour vérifier l'implantation du microbiote et voir les éventuelles différences entre les groupes, les fèces ont été prélévées à 3 et 6 semaines après le transfert et une partie du contenu cæcal a été gardée pour faire l'analyse 16S. frogs.data
```
* Nombre d'échantillons par groupe: * Nombre d'échantillons par groupe:
```{r variable-1, eval = TRUE, echo=TRUE} ```{r variable-1, eval = TRUE, echo=TRUE}
table(get_variable(frogs.data, "Group")) table(get_variable(frogs.data, "Group"))
``` ```
Le groupe **Vehicule** sans rotenone
* Nombre d'échantillons par statut: Le groupe **Rotenone** avec rotenone
```{r variable-2, eval = TRUE, echo=TRUE}
table(get_variable(frogs.data, "Statut"))
```
* Nombre d'échantillons par cohorte:
```{r variable-3, eval = TRUE, echo=TRUE}
table(get_variable(frogs.data, "Cohorte"))
```
## Sélection des données
```{r sample selection 1,fig.width=10,fig.height=4, eval=TRUE, echo = TRUE}
frogs.data <- subset_samples(frogs.data, Cohorte == "Cohorte_1")
frogs.data
```
## Visualisation de l'abondance par échantillon {.tabset} ## Visualisation de l'abondance par échantillon {.tabset}
...@@ -543,12 +526,15 @@ plot(p) ...@@ -543,12 +526,15 @@ plot(p)
#Get count of phyla #Get count of phyla
table<-table(phyloseq::tax_table(frogs.data)[, "Phylum"]) table<-table(phyloseq::tax_table(frogs.data)[, "Phylum"])
table table
``` ```
#### Tableaux relatives abondances Phylum #### Tableaux relatives abondances Phylum
* Relatives Abondances par échantillons: * Relatives Abondances par échantillons:
<div class="alert alert-info" role="alert">Abandance Relative = nb total des individus d'une éspèce par rapport au nb total des individus de toutes les espèces présentes</div>
```{r relative abundance-phylum, eval = TRUE} ```{r relative abundance-phylum, eval = TRUE}
#Convert to relative abundance #Convert to relative abundance
...@@ -585,20 +571,20 @@ rel_abun2<-phyloseq::otu_table(data_phylum_2) ...@@ -585,20 +571,20 @@ rel_abun2<-phyloseq::otu_table(data_phylum_2)
rel_abun2 %>% rel_abun2 %>%
DT::datatable(extensions = 'Buttons', DT::datatable(extensions = 'Buttons',
options = list(dom = 'Bfrtip', options = list(dom = 'Bfrtip',
pageLength = 20, pageLength = 10,
buttons = list('copy', 'csv', 'excel'))) buttons = list('copy', 'csv', 'excel')))
``` ```
#### Plot composition Phylum #### Plot composition Phylum
```{r 05-composition-phylum, eval = TRUE, fig.height=15, fig.width=15} ```{r 03-composition-phylum, eval = TRUE, fig.height=5, fig.width=10}
correct.order <- c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec" ) correct.order <- c("Vehicule","Rotenone")
sample_data(frogs.data)$Group <- factor(sample_data(frogs.data)$Group, sample_data(frogs.data)$Group <- factor(sample_data(frogs.data)$Group,
levels = correct.order) levels = correct.order)
p <- plot_composition(frogs.data, "Kingdom", "Bacteria", "Phylum",numberOfTaxa = 6, fill = "Phylum") p <- plot_composition(frogs.data, "Kingdom", "Bacteria", "Phylum",numberOfTaxa = 6, fill = "Phylum")
p<- p + facet_wrap(~Group, scales = "free_x", nrow = 3) + p<- p + facet_wrap(~Group, scales = "free_x", nrow = 1) +
theme(strip.text.x = element_text(size = 14, color = "black")) + theme(strip.text.x = element_text(size = 14, color = "black")) +
scale_y_continuous(label = scales::percent) scale_y_continuous(label = scales::percent)
...@@ -606,16 +592,16 @@ plot(p) ...@@ -606,16 +592,16 @@ plot(p)
``` ```
```{r 06-composition-phylum-merged, eval = TRUE, fig.height=10, fig.width=10} ```{r 04-composition-phylum-merged, eval = TRUE, fig.height=5, fig.width=6}
frogs.data.merged<-merge_samples(frogs.data,group="Group") frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged) sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)
correct.order <- c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec") correct.order <- c("Vehicule","Rotenone")
sample_data(frogs.data.merged)$Group <- factor(sample_data(frogs.data.merged)$Group, sample_data(frogs.data.merged)$Group <- factor(sample_data(frogs.data.merged)$Group,
levels = correct.order) levels = correct.order)
p <- plot_composition(frogs.data.merged, "Kingdom", "Bacteria", "Phylum",numberOfTaxa = 6, fill = "Phylum") p <- plot_composition(frogs.data.merged, "Kingdom", "Bacteria", "Phylum",numberOfTaxa = 6, fill = "Phylum")
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 3) p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
p <- p + ggtitle("Phylum Composition (12 top Phylum)")+ p <- p + ggtitle("Phylum Composition (12 top Phylum)")+
scale_y_continuous(label = scales::percent) + theme(axis.text.x = element_blank()) scale_y_continuous(label = scales::percent) + theme(axis.text.x = element_blank())
...@@ -628,39 +614,41 @@ plot(p) ...@@ -628,39 +614,41 @@ plot(p)
##### Boxplot Phylum avec stats ##### Boxplot Phylum avec stats
On va se focaliser sur les 6 Phylums et représenter leurs abondances dans les différents groupe sous forme de boxplot. On va se focaliser sur les 6 Phylums et représenter leurs abondances dans les différents groupe sous forme de boxplot.
On va utiliser **S-A_F3w** comme groupe de référence et le comparer à tout les autres avec un test de wilcoxon pour voir desquels il diffère.
Le nombre d'étoiles codent le niveau de significativité, les comparaisons non-significatives ne sont pas indiqués Le nombre d'étoiles codent le niveau de significativité, les comparaisons non-significatives ne sont pas indiqués
<div class="alert alert-info" role="alert">Le test wilcox compare S-A_F3w aux autres groupes</div> <div class="alert alert-info" role="alert">Le test wilcox compare Vehicule à Rotenone</div>
```{r 04-Phylum stats,fig.width=12,fig.height=20, eval = TRUE} ```{r 05-Phylum stats,fig.width=6,fig.height=8, eval = TRUE}
library(ggpubr) library(ggpubr)
library(cowplot)
#Select groups
frogs.data_bis <- subset_samples(frogs.data, Group %in% c("Vehicule","Rotenone"))
## Select only some families ## Select only some families
families <- c("Actinobacteriota","Bacteroidota", "Desulfobacterota","Fibrobacterota" ,"Firmicutes", "Proteobacteria") families <- c("Actinobacteriota","Bacteroidota","Deferribacterota", "Desulfobacterota" ,"Firmicutes", "Proteobacteria")
phy <- frogs.data %>% subset_taxa(Phylum %in% families) %>% tax_glom(taxrank="Phylum") #agglomerate at family level phy <- frogs.data_bis %>% subset_taxa(Phylum %in% families) %>% tax_glom(taxrank="Phylum") #agglomerate at family level
## Transform count to relative abundances ## Transform count to relative abundances
depth <- sample_sums(frogs.data)[1] depth <- sample_sums(frogs.data)[1]
plotdata<-psmelt(phy) %>% plotdata<-psmelt(phy) %>%
mutate(Abundance = Abundance / depth, mutate(Abundance = Abundance / depth,
Group = factor(Group, labels = c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec" ))) Group = factor(Group, labels = c("Vehicule","Rotenone" )))
p <- ggplot(plotdata,aes(x = Group,y=Abundance, color = Group, Group = Group)) + p <- ggplot(plotdata,aes(x = Group,y=Abundance, color = Group, Group = Group)) +
stat_boxplot(geom = "errorbar", width = 0.5) + stat_boxplot(geom = "errorbar", width = 0.5) +
geom_boxplot(outlier.alpha = 1, geom_boxplot(outlier.alpha = 1,
outlier.size = 0.8) + outlier.size = 0.8) +
facet_wrap(~Phylum, scales = "free_y", ncol = 1) + facet_wrap(~Phylum, scales = "free_y", ncol = 3) +
scale_color_manual(values = c("S-A_F3w" = "#c61104","A_F3w" = "#027d49", "S-AG_F3w" = "#015e90", "AG_F3w" = "#74037c", scale_color_manual(values = c("Vehicule" = "#227432", "Rotenone" = "#8C67A9"), guide = "none") +
"S-A_F6w" = "#fc392c","A_F6w" = "#03c372","S-AG_F6w" = "#0190dd", "AG_F6w" = "#b803c4",
"S-A_Caec" = "#ff7b73","A_Caec" = "#7afcc5", "S-AG_Caec" = "#6ecafc","AG_Caec" = "#f79efd"), guide = "none") +
# theme_classic() + ## fond blanc # theme_classic() + ## fond blanc
labs(x = NULL) + labs(x = NULL) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10)) theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10))
## Compare reference level T0 to all other levels using a wilcoxon test and adjust p-values using the holm correction. ## Compare reference level T0 to all other levels using a wilcoxon test and adjust p-values using the holm correction.
p <- p + stat_compare_means(aes(label = ..p.signif..), method = "wilcox.test", p.adjust.method = "holm", ref.group = "S-A_F3w", hide.ns = T, p <- p + stat_compare_means(aes(label = ..p.signif..), method = "wilcox.test", p.adjust.method = "holm", ref.group = "Vehicule", hide.ns = T,
label.y.npc = c(0.90), label.y.npc = c(0.90),
size = 7, size = 7,
fontface = "bold") fontface = "bold")
...@@ -678,7 +666,7 @@ data.phylum <- frogs.data %>% ...@@ -678,7 +666,7 @@ data.phylum <- frogs.data %>%
transform_sample_counts(function(x) { x / sum(x)}) %>% ## transform counts to proportions transform_sample_counts(function(x) { x / sum(x)}) %>% ## transform counts to proportions
fast_tax_glom(taxrank = "Phylum") %>% fast_tax_glom(taxrank = "Phylum") %>%
psmelt() %>% psmelt() %>%
mutate(Group = factor(Group, labels = c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec"))) mutate(Group = factor(Group, labels = c("Vehicule","Rotenone")))
``` ```
```{r test kruskal-phylum, eval=TRUE} ```{r test kruskal-phylum, eval=TRUE}
...@@ -686,63 +674,15 @@ data.test <- compare_means(Abundance ~ Group, data = data.phylum, method = "krus ...@@ -686,63 +674,15 @@ data.test <- compare_means(Abundance ~ Group, data = data.phylum, method = "krus
data.test data.test
``` ```
##### Dunn Post-Hoc test
Dunn Post-Hoc test sur le résultat Kruskall-Wallis de l'impact de `Group` sur les abondances Actinobacteriota
```{r Post-Hoc test Actinobacteriota, eval=TRUE}
filter(data.phylum, Phylum == "Actinobacteriota") %>%
dunn.test::dunn.test(x = .$Abundance, g = .$Group, method = "bh")
```
Dunn Post-Hoc test sur le résultat Kruskall-Wallis de l'impact de `Group` sur les abondances Bacteroidota
```{r Post-Hoc test Bacteroidetes, eval=TRUE}
filter(data.phylum, Phylum == "Bacteroidota") %>%
dunn.test::dunn.test(x = .$Abundance, g = .$Group, method = "bh")
```
Dunn Post-Hoc test sur le résultat Kruskall-Wallis de l'impact de `Group` sur les abondances Desulfobacterota
```{r Post-Hoc test Desulfobacterota, eval=TRUE}
filter(data.phylum, Phylum == "Desulfobacterota") %>%
dunn.test::dunn.test(x = .$Abundance, g = .$Group, method = "bh")
```
Dunn Post-Hoc test sur le résultat Kruskall-Wallis de l'impact de `Group` sur les abondances Fibrobacterota
```{r Post-Hoc test Fibrobacterota, eval=TRUE}
filter(data.phylum, Phylum == "Fibrobacterota") %>%
dunn.test::dunn.test(x = .$Abundance, g = .$Group, method = "bh")
```
Dunn Post-Hoc test sur le résultat Kruskall-Wallis de l'impact de `Group` sur les abondances Firmicutes
```{r Post-Hoc test firmicutes, eval=TRUE}
filter(data.phylum, Phylum == "Firmicutes") %>%
dunn.test::dunn.test(x = .$Abundance, g = .$Group, method = "bh")
```
Dunn Post-Hoc test sur le résultat Kruskall-Wallis de l'impact de `Group` sur les abondances Proteobacteria
```{r Post-Hoc test Proteobacteria, eval=TRUE}
filter(data.phylum, Phylum == "Proteobacteria") %>%
dunn.test::dunn.test(x = .$Abundance, g = .$Group, method = "bh")
```
### Représentation des communautés au niveau Family {.tabset} ### Représentation des communautés au niveau Family {.tabset}
* Nombre d'OTUs par Family: * Nombre d'OTUs par Family:
```{r count-family,echo=FALSE, eval=TRUE} ```{r count-family, eval=TRUE}
#Get count of phyla #Get count of phyla
table(phyloseq::tax_table(frogs.data)[, "Family"]) table(phyloseq::tax_table(frogs.data)[, "Family"])
``` ```
#### Tableaux relatives abondances Family #### Tableaux relatives abondances Family
...@@ -770,7 +710,7 @@ rel_abun_3 %>% ...@@ -770,7 +710,7 @@ rel_abun_3 %>%
* Relative Abondance par groupe: * Relative Abondance par groupe:
```{r 2.1-relative abundance-family-2,echo=FALSE, eval=TRUE} ```{r 2.1-relative abundance-family-2, eval=TRUE}
#Convert to relative abundance #Convert to relative abundance
data_rel_abund = phyloseq::transform_sample_counts(frogs.data, function(x) 100 * x/sum(x)) data_rel_abund = phyloseq::transform_sample_counts(frogs.data, function(x) 100 * x/sum(x))
...@@ -794,29 +734,29 @@ rel_abun_4 %>% ...@@ -794,29 +734,29 @@ rel_abun_4 %>%
#### Plot composition Family #### Plot composition Family
```{r 8-composition-phylum,fig.width=10,fig.height=10,echo=FALSE, eval=TRUE} ```{r 06-composition-phylum,fig.width=10,fig.height=5, eval=TRUE}
correct.order <- c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec") correct.order <- c("Vehicule","Rotenone")
sample_data(frogs.data)$Group <- factor(sample_data(frogs.data)$Group, sample_data(frogs.data)$Group <- factor(sample_data(frogs.data)$Group,
levels = correct.order) levels = correct.order)
p <- plot_composition(frogs.data, "Kingdom", "Bacteria", "Family",numberOfTaxa = 10, fill = "Family") p <- plot_composition(frogs.data, "Kingdom", "Bacteria", "Family",numberOfTaxa = 10, fill = "Family")
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 3) p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
plot(p) plot(p)
``` ```
```{r 7-composition-family merged,fig.width=8,fig.height=10,echo=FALSE, eval=TRUE} ```{r 07-composition-family merged,fig.width=6,fig.height=5, eval=TRUE}
frogs.data.merged<-merge_samples(frogs.data,group="Group") frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged) sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)
correct.order <- c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec") correct.order <- c("Vehicule","Rotenone")
sample_data(frogs.data.merged)$Group <- factor(sample_data(frogs.data.merged)$Group, sample_data(frogs.data.merged)$Group <- factor(sample_data(frogs.data.merged)$Group,
levels = correct.order) levels = correct.order)
p <- plot_composition(frogs.data.merged, "Kingdom", "Bacteria", "Family",numberOfTaxa = 10, fill = "Family") p <- plot_composition(frogs.data.merged, "Kingdom", "Bacteria", "Family",numberOfTaxa = 10, fill = "Family")
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 3) + theme(axis.text.x = element_blank()) p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1) + theme(axis.text.x = element_blank())
plot(p) plot(p)
``` ```
...@@ -827,47 +767,47 @@ plot(p) ...@@ -827,47 +767,47 @@ plot(p)
On va se focaliser sur les 10 Family et représenter leurs abondances dans les différents groupe sous forme de boxplot. On va se focaliser sur les 10 Family et représenter leurs abondances dans les différents groupe sous forme de boxplot.
On va utiliser **S-A_F3w** comme groupe de référence et le comparer à tout les autres avec un test de wilcoxon pour voir desquels il diffère.
Le nombre d'étoiles codent le niveau de significativité, les comparaisons non-significatives ne sont pas indiqués Le nombre d'étoiles codent le niveau de significativité, les comparaisons non-significatives ne sont pas indiqués
```{r 04-Family stats,fig.width=12,fig.height=22, eval = TRUE, eval=TRUE} <div class="alert alert-info" role="alert">Le test wilcox compare Vehicule à Rotenone </div>
library(ggpubr)
```{r 08-Family stats,fig.width=10,fig.height=8, eval = TRUE}
#Select groups
frogs.data_quat <- subset_samples(frogs.data, Group %in% c("Vehicule","Rotenone"))
## Select only some families ## Select only some families
families <- c("Atopobiaceae","Bacteroidaceae", "Bifidobacteriaceae", "Desulfovibrionaceae","Lachnospiraceae" , "Marinifilaceae", "Prevotellaceae", "Rikenellaceae","Ruminococcaceae", "Tannerellaceae") families <- c("Muribaculaceae","Marinifilaceae", "Rikenellaceae", "Tannerellaceae","Prevotellaceae" , "Bacteroidaceae", "Lachnospiraceae", "Ruminococcaceae","Anaerovoracaceae", "Peptostreptococcaceae")
phy <- frogs.data %>% subset_taxa(Family %in% families) %>% tax_glom(taxrank="Family") #agglomerate at family level phy <- frogs.data_quat %>% subset_taxa(Family %in% families) %>% tax_glom(taxrank="Family") #agglomerate at family level
## Transform count to relative abundances ## Transform count to relative abundances
depth <- sample_sums(frogs.data)[1] depth <- sample_sums(frogs.data_quat)[1]
plotdata<-psmelt(phy) %>% plotdata<-psmelt(phy) %>%
mutate(Abundance = Abundance / depth, mutate(Abundance = Abundance / depth,
Group = factor(Group, labels = c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec" ))) Group = factor(Group, labels = c("Vehicule","Rotenone")))
p <- ggplot(plotdata,aes(x = Group,y=Abundance, color = Group, Group = Group)) + p <- ggplot(plotdata,aes(x = Group,y=Abundance, color = Group, Group = Group)) +
stat_boxplot(geom = "errorbar", width = 0.5) + stat_boxplot(geom = "errorbar", width = 0.5) +
geom_boxplot(outlier.alpha = 1, geom_boxplot(outlier.alpha = 1,
outlier.size = 0.8) + outlier.size = 0.8) +
facet_wrap(~Family, scales = "free_y", ncol = 1) + facet_wrap(~Family, scales = "free_y", ncol = 5) +
scale_color_manual(values = c("S-A_F3w" = "#c61104","A_F3w" = "#027d49", "S-AG_F3w" = "#015e90", "AG_F3w" = "#74037c", scale_color_manual(values = c("Vehicule" = "#227432", "Rotenone" = "#8C67A9"), guide = "none") +
"S-A_F6w" = "#fc392c","A_F6w" = "#03c372","S-AG_F6w" = "#0190dd", "AG_F6w" = "#b803c4",
"S-A_Caec" = "#ff7b73","A_Caec" = "#7afcc5", "S-AG_Caec" = "#6ecafc","AG_Caec" = "#f79efd"), guide = "none") +
# theme_classic() + ## fond blanc # theme_classic() + ## fond blanc
labs(x = NULL) + labs(x = NULL) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10)) theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10))
## Compare reference level T0 to all other levels using a wilcoxon test and adjust p-values using the holm correction. ## Compare reference level T0 to all other levels using a wilcoxon test and adjust p-values using the holm correction.
p <- p + stat_compare_means(aes(label = ..p.signif..), method = "wilcox.test", p.adjust.method = "holm", ref.group = "S-A_F3w", hide.ns = T, p <- p + stat_compare_means(aes(label = ..p.signif..), method = "wilcox.test", p.adjust.method = "holm", ref.group = "Vehicule", hide.ns = T,
label.y.npc = c(0.90), label.y.npc = c(0.90),
size = 7, size = 7,
fontface = "bold") fontface = "bold")
plot(p) plot(p)
``` ```
##### Kruskal test ##### Kruskal test
Test de Kruskal (Anova non paramétrique) pour chaque Phylum pour tester si les abondances sont similaires (pour ce Phylum) entre les différents groupes. Test de Kruskal (Anova non paramétrique) pour chaque Family pour tester si les abondances sont similaires (pour la Family) entre les différents groupes.
```{r data.-family, eval=TRUE} ```{r data.-family, eval=TRUE}
# depth <- sample_sums(frogs.data)[1] # depth <- sample_sums(frogs.data)[1]
...@@ -875,7 +815,7 @@ data.phylum <- frogs.data %>% ...@@ -875,7 +815,7 @@ data.phylum <- frogs.data %>%
transform_sample_counts(function(x) { x / sum(x)}) %>% ## transform counts to proportions transform_sample_counts(function(x) { x / sum(x)}) %>% ## transform counts to proportions
fast_tax_glom(taxrank = "Family") %>% fast_tax_glom(taxrank = "Family") %>%
psmelt() %>% psmelt() %>%
mutate(Group = factor(Group, labels = c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec"))) mutate(Group = factor(Group, labels = c("Vehicule","Rotenone")))
``` ```
```{r test kruskal-family, eval=TRUE} ```{r test kruskal-family, eval=TRUE}
...@@ -889,147 +829,135 @@ Les 6 phylums d'intérêts sont les Actinobacteriota, les Bacteroidota, les Desu ...@@ -889,147 +829,135 @@ Les 6 phylums d'intérêts sont les Actinobacteriota, les Bacteroidota, les Desu
#### Actinobacteriota #### Actinobacteriota
```{r 15-composition-Actinobacteriota,echo=FALSE,fig.width=10,fig.height=10, eval=TRUE} ```{r 09-composition-Actinobacteriota,fig.width=8,fig.height=5, eval=TRUE}
correct.order <- c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec") correct.order <- c("Vehicule","Rotenone")
sample_data(frogs.data)$Group <- factor(sample_data(frogs.data)$Group, sample_data(frogs.data)$Group <- factor(sample_data(frogs.data)$Group,
levels = correct.order) levels = correct.order)
p <- plot_composition(frogs.data, "Phylum", "Actinobacteriota", "Family",numberOfTaxa = 5, fill = "Family") p <- plot_composition(frogs.data, "Phylum", "Actinobacteriota", "Family",numberOfTaxa = 5, fill = "Family")
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 3) p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
plot(p) plot(p)
``` ```
```{r 7-composition-Actinobacteriota merged,fig.width=8,fig.height=10,echo=FALSE, eval=TRUE}