---
title: "TD phyloseq 05-12-2023"
output: html_document
date: "2023-12-04"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
install.packages("phyloseq")
install.packages("magrittr")
install.packages("dendextend")
install.packages("vegan")
library(phyloseq)
library(tidyverse)
```


```{r}
ps = readRDS("/Users/djelika/Documents/2019-11 Sci Rep Villette Scarcity Paper/Metadata/2021-05 Sci Rep Villette et al 16S refinement/2021-05 VILLETTE et al 16S Refinement ps.rds")

ps
```


```{r}

human = subset_samples(ps, Platform=="BioFidal")

Donor <- ps 

Donor= Donor %>% transform_sample_counts(function(x)x/sum(x))

sample_names(Donor)= paste(Donor@sam_data$Filename, 
                           Donor@sam_data$Mechanical_Lyse, 
                           Donor@sam_data$PCR, 
                           Donor@sam_data$Concentration,sep="_")
 
sample_names(human)= paste(human@sam_data$Filename, 
                           human@sam_data$Mechanical_Lyse, 
                           human@sam_data$PCR, 
                           human@sam_data$Concentration,sep="_")

```

# Représenter la diversité alpha

```{r}
richness = as.data.frame(estimate_richness(human, measures = c("Shannon","Observed"))) 
head(richness) 
```

#Mesurer l'alpha-diversité pour différente concentration de microbiote
```{r}
# Rajouter des variables importantes pour la suite de l'analyse

richness$SampleOrigin<- human@sam_data$SampleOrigin #paste sample info 
richness$PCR= human@sam_data$PCR
richness$Concentration = human@sam_data$Concentration

ggplot(richness, aes(x=PCR, y=Observed, fill=Concentration))+
  geom_bar(stat='identity', position = "dodge")+
  facet_grid(SampleOrigin~PCR, scales = "free_x", space = "free")+
  ylab("Richness observed")+
  theme_bw()+
  theme(strip.text.y =element_text(size = 6, face = "bold"), strip.background = element_blank(), 
        legend.text = element_text(size = 10, face="bold"), 
        axis.text.x = element_text(size = 10, face = "bold"),
        legend.title = element_text(size = 10, face="bold"), 
        axis.title.x = element_blank(), 
        axis.title.y=element_text(size = 10, face="bold"), 
        axis.text.y = element_text(size = 10), 
        panel.grid = element_blank())
```

#Composition microbienne au niveau Phylum
```{r}
Phylum <- tax_glom(human, "Phylum")

plot_bar(Phylum, x = "Concentration", fill = "Phylum")+
  facet_grid(SampleOrigin~PCR, space = "free", scales = "free")+
  labs(x="Microbes quantity", y="Relative abundance")+
  theme_bw()+
  theme(strip.background = element_blank(), 
        axis.title.x = element_text(size=10, face="bold"), 
        strip.text = element_text(size = 6, face = "bold"), 
        axis.text.x = element_text(size = 10, face = "bold"), 
        axis.text.y = element_text(size = 10), 
        axis.title.y = element_text(size = 10, face="bold"), 
        legend.title = element_blank(), 
        legend.text = element_text(size = 10, face="bold.italic"), 
        legend.spacing.x = unit(0.5,"cm"), 
        panel.grid = element_blank())
```


#Composition microbienne au niveau Class
```{r}
Class<- tax_glom(human, "Class")

plot_bar(Class, x = "Concentration", fill = "Class")+
  facet_grid(SampleOrigin~PCR, space = "free", scales = "free")+
  labs(x="Microbes quantity", y="Relative abundance")+
  theme_bw()+
  theme_bw()+
  theme(strip.background = element_blank(), 
        axis.title.x = element_text(size=10, face="bold"), 
        strip.text = element_text(size = 6, face = "bold"), 
        axis.text.x = element_text(size = 10, face = "bold"), 
        axis.text.y = element_text(size = 10), 
        axis.title.y = element_text(size = 10, face="bold"), 
        legend.title = element_blank(), 
        legend.text = element_text(size = 10, face="bold.italic"), 
        legend.spacing.x = unit(0.5,"cm"), 
        panel.grid = element_blank())

```


```{r}
Genus<- tax_glom(human, "Genus")

plot_bar(Genus, x = "Concentration", fill = "Genus")+
  facet_grid(SampleOrigin~PCR, space = "free", scales = "free")+
  labs(x="Microbes quantity", y="Relative abundance")+
  theme_bw()+
  theme_bw()+
  theme(strip.background = element_blank(), 
        axis.title.x = element_text(size=10, face="bold"), 
        strip.text = element_text(size = 6, face = "bold"), 
        axis.text.x = element_text(size = 10, face = "bold"), 
        axis.text.y = element_text(size = 10), 
        axis.title.y = element_text(size = 10, face="bold"), 
        legend.title = element_blank(), 
        legend.text = element_text(size = 10, face="bold.italic"), 
        legend.spacing.x = unit(0.5,"cm"), 
        panel.grid = element_blank())

Genus
```

#Dsitance entre échantillon

```{r}
ord <- ordinate(human, method = "PCoA", distance = "bray")
dist <- phyloseq::distance(human, method = "bray")  

plot_ordination(human, ord, color = "SampleOrigin", shape="PCR",label = "Concentration")+
  geom_point(size=4, alpha=0.5)+ 
  theme_bw()+
  theme(legend.text = element_text(size = 7, face = "bold"),
        axis.title = element_text(size = 10, face = "bold"), 
        legend.spacing.x = unit(0.3, "cm"))
```

#Clusterisation
##Clusteriser
```{r}
ord <- ordinate(human, method = "PCoA", distance = "bray")
dist <- phyloseq::distance(human, method = "bray")  
plot_scree(ord)
```

##2 clusters
```{r}
clust <- hclust(dist, method = "ward.D2")
clust$labels = with(sample_data(human), paste(SampleOrigin, PCR, Concentration))
sample_data(human)$cluster= as.factor(cutree(clust, k=2))

tmp2 = ord$vectors[,1:2] 
tmp2 = as.data.frame(merge(tmp2, sample_data(human), by="row.names"))

tmp2 %>%
  ggplot(aes(x=Axis.1, y=Axis.2, color=SampleOrigin, group=cluster))+ 
  stat_ellipse()+
    geom_point(alpha=0.7,size=3)+
    labs(x="Axis.1", y="Axis.2", color="Cluster")+
    theme_bw(base_line_size = 0)+
      theme(legend.position = "right", 
        text = element_text(size = 20, face = "bold"))


```
##Clusteriser
```{r}

clust <- hclust(dist, method = "ward.D2")
clust$labels = with(sample_data(human), paste(SampleOrigin, PCR, Concentration))
sample_data(human)$cluster= as.factor(cutree(clust, k=3))


tmp3 = ord$vectors[,1:2] 
tmp3 = as.data.frame(merge(tmp3, sample_data(human), by="row.names"))

tmp3 %>%
  ggplot(aes(x=Axis.1, y=Axis.2, color=SampleOrigin, group=cluster))+ 
  stat_ellipse()+
    geom_point(alpha=0.7,size=3)+
    labs(x="Axis.1", y="Axis.2", color="Cluster")+
    theme_bw(base_line_size = 0)+
      theme(legend.position = "right", 
        text = element_text(size = 20, face = "bold"))


```

#Comparer la composition phylogénétique
## Préparer les annotations
```{r}
Heat<- tax_glom(ps, "Genus") 
Heat<- prune_taxa(names(sort(taxa_sums(Heat),TRUE)[1:50]), Heat) 

x= t(as.matrix(Heat@otu_table@.Data)) 

dend= hclust(distance(ps, method="bray"), "ward.D2") 
dend= hclust(vegdist(ps@otu_table, method="bray"), "ward.D2") 

dend= dendextend:: click_rotate(as.dendrogram(dend)) 

group= HeatmapAnnotation(Donor=Heat@sam_data$SampleName,
            PCR=Heat@sam_data$PCR , 
            col=list(Donor= c("HD1"="#114477", "HD2"="#AA4488"),
            PCR= c("Nested"="#44AAAA" , "Standard"="#DDDD77")),
            show_annotation_name = T, 
            annotation_name_gp =  gpar(fontsize=15, fontface="bold")) 

name= HeatmapAnnotation(Concentration= anno_text(Heat@sam_data$Concentration, 
                            gp = gpar(fontsize=14, fontface="bold"))) 
```

##Heatmap
```{r}

Heatmap(as.matrix(scale(x)), 
        top_annotation = group, 
        bottom_annotation = name,
        cluster_rows = F, 
        cluster_columns = as.dendrogram(dend), 
        row_labels = Heat@tax_table[,"Genus"],
        row_names_gp = gpar(fontsize=15, fontface="bold.italic"),
        show_column_names = F, 
        col = c("white","#88CCAA","#771122"), 
        column_dend_height = unit(3, "cm"), 
        name = "Abundance", 
        column_split = 3, 
        column_gap = unit(5, "mm"), 
        column_title = NULL)
```





