DMM_clustering.Rmd

---
title: "DMM_clustering"
author: "ruben"
date: "30/03/2022"
output: html_document
---

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

## Librairies

```{r,message=F,warning=F}

library(tidyverse)
library(data.table)
library(devtools)
library(microbiome)
library(DirichletMultinomial)
library(reshape2)
library(magrittr)
library(phyloseq)
library(dplyr)
library(tibble)
library(rstatix)
library(ggpubr)
library(ggplot2)
library(tidyverse)
#library(readxl)
library("optparse")
library(varhandle)
library(tidyr)
library("stringr")

library(gplots)

source(file = "functions.R")

#setwd("/home/ladeirru/GitHub/microbiome.actibiome")

```


## Data

```{r,message=F,warning=F}

seed = 12000

set.seed(seed)

threshold = 500

# we load this file to have taxa table specifically
load(file = "study_after_QC.rda")

# we load only the phy object
load(file = "phy_curated_data_after_QC.rda")

phy <- subset_samples(phy, body_site == "stool")

nrow(data.frame(phy@sam_data))

## diversity vegan

phy_diversity <- estimate_richness(phy, split = TRUE, measures = "Shannon")

#phy <- subset_samples(phy, age_category == "adult")

#phy <- subset_samples(phy, disease == "Control")

phy <- subset_samples(phy, antibiotics_current_use == "no")

nrow(data.frame(phy@sam_data))

#phy <- subset_samples(phy, westernized == "yes")

```

## Filtering curated/ascinar data

```{r,message=F,warning=F}

pseq <- subset_taxa(phy, Genus %in% c("Bifidobacterium"))

#pseq = phy

# 
# which_study <- readline(prompt=paste("Do you want to cluster the data on ASNICAR (Predict) or whole CURATED data without predict ? : ", "pre - PREDICT", "cur - CURATED minus PREDICT","who - WHOLE CURATED DATA SET", "6clusters - 6 clusters model","who - WHOLE CURATED DATA SET", "", sep="\n"))
# 
# while(!which_study %in% c("cur", "pre", "who", "6clusters")){
#   which_study <- readline(prompt=paste("Error, please redo : ", "pre - PREDICT", "cur -
#                                        CURATED minus PREDICT","who - WHOLE CURATED DATA SET","6clusters - 6 clusters model", "", sep="\n"))
# }

which_study = "6clusters"


if(which_study == "cur"){
  ## we exclude predict data
  
  pseq <- subset_samples(pseq, !PMID == 33432175)
  
  ## keep only adults subjects
  
  pseq <- subset_samples(pseq, age_category == "adult")
  
}



if(which_study == "pre"){
  ## we exclude predict data
  
  pseq <- subset_samples(pseq, PMID == 33432175)
  

  ## keep only adults subjects
  
  pseq <- subset_samples(pseq, age_category == "adult")
  
}






```

## File names

```{r,message=F,warning=F}

## file names for saving seeds and DMM results

dmm_name = paste0("DMM_", which_study)

seeds_list_name = paste0("seeds_", which_study)

```

## Transformation

```{r,message=F,warning=F}

# phy.comp <- microbiome::transform(phy, "compositional")
# compositional transformation, ie, relative abundance

#taxa <- core_members(phy.comp, detection = 0.1/100, prevalence = 50/100)

#pseq <- prune_taxa(taxa, phy)

# Pick the OTU count matrix
# and convert it into samples x taxa format
dat <- abundances(pseq)
count_raw <- as.matrix(t(dat))


# rowsums equal zero (sample with 0 count everywhere) need to be deleted

sample_low_reads <- names(which(rowSums(count_raw) <= 500))

# colsums equal zero (bacteria not present in any sample) need to be deleted

count <- count_raw %>%
  as.data.frame(.) %>%
  filter_all(., any_vars(. != 0))



## select ppl with 0 bifid, they are selected as a no bifid cluster

no_bifid <- count_raw %>%
  as.data.frame(.) %>%
  filter(!row.names(.) %in% rownames(count))

# colsums equal zero (bacteria not present in any sample) need to be deleted

otu_not_present <- names(which(colSums(count) == 0))

count <- count %>%
  select(-otu_not_present) %>%
  ## filtering ppl with nb reads < threshold
  filter(!row.names(.) %in% sample_low_reads) %>%
  as.matrix(.)

# transform colnames to species name

colnames(count) <- get_species(as.data.frame(tax),colnames(count))

nrow(count)

```


## Code à faire marcher la nuit

```{r,message=F,warning=F}

fit_genus_list = vector("list",5)


set.seed(seed)
seeds=sample(1:1000, 5)

#setwd("/home/ladeirru/GitHub/microbiome.actibiome")

if(!(dmm_name %in% list.files("./DMM_files/")) & !(seeds_list_name %in% list.files("./DMM_files/"))) {

  for(i in 1:5) {

    set.seed(seeds[i])

    fit_genus <- parallel::mclapply(1:30, dmn, count=count, verbose=FALSE, mc.cores=12)

    fit_genus_list[[i]] = fit_genus
    
    print(i)

  }
  
  save(fit_genus_list, file=paste0("./DMM_files/", dmm_name))
  
  save(seeds, file = paste0("./DMM_files/", seeds_list_name))

} else {
  load(paste0("./DMM_files/", dmm_name))
  load(paste0("./DMM_files/", seeds_list_name))
  }


# collect scores to find the best fit
lplc = vector("list",5)
bic <- vector("list",5)
aic <- vector("list",5)

for(i in 1:5) {

  
  lplc[[i]] <- base::sapply(fit_genus_list[[i]], DirichletMultinomial::laplace)

  aic[[i]]  <- base::sapply(fit_genus_list[[i]], DirichletMultinomial::AIC)

  bic[[i]]  <- base::sapply(fit_genus_list[[i]], DirichletMultinomial::BIC)

}

```

## Model fit

```{r,message=F,warning=F}

number_of_clusters <- as.character(length(fit_genus_list[[1]]))

## we will take the number k for each model and for each method

list_k <- vector("list",15)

for (i in 1:5) {
  
  seed_experiment <-seeds[[i]]
  
  list_k[[i]] <- c("lplc", rank(lplc[[i]]), seed_experiment)
  
  list_k[[i+5]] <- c("aic", rank(aic[[i]]), seed_experiment)
  
  list_k[[i+10]] <- c("bic", rank(bic[[i]]), seed_experiment)
  
  print(i)
  
}

test <- data.frame(t(as.data.frame(list_k)))

rownames(test) <- NULL

colnames(test) <- c("method", 1:number_of_clusters, "seed")


test_unfactored <- test %>%
  select('1':as.character(number_of_clusters)) %>%
  varhandle::unfactor(.) %>%
  cbind(method = test$method, . , seed = test$seed)
rownames(test_unfactored) <- NULL
  

## calcul des rangs moyens

mean_rank <- test_unfactored %>%
  mutate_at(-1, as.numeric) %>%
  select(-seed) %>%
  mutate_at(2:31, as.numeric) %>%
  #reshape2::melt(id.vars="method") %>%
  #dplyr::rename(k=variable) %>%
  #mutate(value=as.numeric(value)) %>%
  group_by(method) %>%
  #summarise_all(mean) %>%
  #summarise(value=mean(value)) %>%
  summarise(across(everything(), mean)) %>%
  #summarise_all(mean) %>%
  gather(., key="k", value="mean_rank", '1':as.character(number_of_clusters))


sd_rank <- test_unfactored %>%
  mutate_at(-1, as.numeric) %>%
  select(-seed) %>%
  mutate_at(2:31, as.numeric) %>%
  #reshape2::melt(id.vars="method") %>%
  #dplyr::rename(k=variable) %>%
  #mutate(value=as.numeric(value)) %>%
  group_by(method) %>%
  #summarise_all(mean) %>%
  #summarise(value=mean(value)) %>%
  summarise(across(everything(), sd)) %>%
  #summarise_all(mean) %>%
  gather(., key="k", value="sd_rank", '1':as.character(number_of_clusters))




```


## Plot des rangs moyens

```{r,message=F,warning=F}

merge(mean_rank, sd_rank, by=c("method","k")) %>%
  filter(method=="bic") %>%

  ggplot() +
    geom_line(aes(x=k, y=mean_rank, group=method, color =method)) +
  #geom_line(aes(x=k, y=mean_rank))
    #geom_line(aes(x=k, y=mean_rank + (1.96*(sd_rank/sqrt(30))), group=method, color =method, linetype=2)) +
    geom_ribbon(aes(x=k, ymin = mean_rank - (1.96*(sd_rank/sqrt(30))), ymax = mean_rank + (1.96*(sd_rank/sqrt(30))),  group=method, fill =method), alpha = 0.2) +
  #geom_ribbon(aes(x=k, ymin = mean_rank - (1.96*(sd_rank/sqrt(30))), ymax = mean_rank + (1.96*(sd_rank/sqrt(30)))), alpha = 0.2) +
    theme_classic() +
    scale_x_discrete("Number of Dirichlet component",limits = as.character(1:number_of_clusters)) +
    scale_color_brewer("criteria",type="qual") +
    scale_fill_brewer("criteria",type="qual") +
    ylab("model optimum average rank") + guides(color=FALSE,fill=FALSE) -> criteria_dmm_plot

save(criteria_dmm_plot, file="criteria_dmm_plot.rda")

criteria_dmm_plot

ggsave(filename = "criteria_dmm_plot.pdf" )


```
## Model selection

```{r,message=F,warning=F}

## k to select is the minimum for BIC

k_BIC <- mean_rank %>%
  filter(method == "bic") %>%
  filter(mean_rank == min(.$mean_rank)) %>%
  select(k) %>%
  as.numeric(.)

## we select k=6 with BIC method

model_k <- vector("list",5)

for (i in 1:5) {
  model_k[[i]] <- which(bic[[i]] == min(bic[[i]]))
}

index <- which(abs(as.data.frame(model_k) - k_BIC) == min(abs(as.data.frame(model_k) - k_BIC)) )

## the model 1 is the best fit model


if (length(index) > 0){
  ## if several models has a min corresponding to the mean rank
  ## choose a random model
  
  index <- sample(index, 1)
  
  }

fit <- fit_genus_list[[index]]

best <- fit[[which.min(unlist(bic[[index]]))]]


## save best

name_best <- paste0("./DMM_files/best_fit_DMM_", which_study)

if(!name_best %in% list.files("./DMM_files/")){
  
  save(best, file = name_best)
}

## best is the best model based on BIC criteria

pi_theta <- data.frame(mixturewt(best)) %>%
  rownames_to_column("cluster") %>%
  gather(key='method', value='value',pi:theta)

pi_theta$cluster <- as.character(pi_theta$cluster)

```

## teta comparison

# calculate a model k = 1

```{r,message=F,warning=F}

set.seed(seeds)

fit_genus_k1 <- parallel::mclapply(1:1, dmn, count=count, verbose=FALSE, mc.cores=12)

```

# plot

```{r,message=F,warning=F}

## we plot pi and theta by clusters. theta : give the idea of the homogeneity of each cluster

ggplot(pi_theta) +
  geom_line(aes(x=cluster, y=value, group=method, color =method)) +
  theme_minimal() +
  scale_x_discrete(limits = as.character(1:6)) +
  labs(color="", y = "")

```


## Check figure

```{r,message=F,warning=F}

laplace_plot <- data.frame(Fit = lplc[[index]], k = c(1:length(lplc[[index]])))
aic_plot <- data.frame(Fit = aic[[index]], k = c(1:length(aic[[index]])))
bic_plot <-data.frame(Fit = bic[[index]], k = c(1:length(bic[[index]])))

## on bind les 3 par lignes

binded <- bind_rows(list(laplace = laplace_plot, aic = aic_plot, bic = bic_plot), .id = "type_check")

## the minimum for laplace check

mins_laplace <- binded %>%
  filter(type_check == 'laplace' ) %>%
  filter(Fit == min(Fit)) %>%
  select(Fit, k)

mins_bic <- binded %>%
  filter(type_check == 'bic' ) %>%
  filter(Fit == min(Fit)) %>%
  select(Fit, k)


nb_clusters <- ggplot(binded) +
  geom_line(aes(x=k, y=Fit, group=type_check, color =type_check)) +
  theme_minimal() +
  scale_x_continuous(breaks = seq(0, max(binded$k) + 5, by = 5))  +
  
  ### we add laplace minimum
  annotate("point", x=mins_laplace$k, y=mins_laplace$Fit, color='orange', size=4) +
  geom_hline(yintercept = mins_laplace$Fit, linetype='dashed') +
  annotate("label", x=max(laplace_plot$k) + 3, y=mins_laplace$Fit, color='black', label= c("laplace minimum"), size=4) +
  annotate(geom = "text", x = mins_laplace$k, y = mins_laplace$Fit, label = paste0("k = ", mins_laplace$k), vjust = -0.5) +
  
    ### we add bic minimum
  annotate("point", x=mins_bic$k, y=mins_bic$Fit, color='orange', size=4) +
  geom_hline(yintercept = mins_bic$Fit, linetype='dashed') +
  annotate("label", x=max(laplace_plot$k) + 3, y=mins_bic$Fit, color='black', label= c("bic minimum"), size=4) +
  annotate(geom = "text", x = mins_bic$k, y = mins_bic$Fit, label = paste0("k = ", mins_bic$k), vjust = -0.5) +
  
  coord_cartesian(xlim =c(min(binded$k),max(binded$k) + 6))  +
  
  labs(x = "Number of Dirichlet Component - K", y = "Model Fit", title = "", color = "")

nb_clusters

```

## Laplace check

```{r,message=F,warning=F}

## sample component assignments

model <- apply(mixture(best), 1, which.max)

save(model, file="model_DMM")

## write the enterotype file

enterotype <- data.frame(bifidotype = model)

enterotype <- enterotype %>%
  rownames_to_column("sample_id")

## add no bifid cluster as cluster 7 to enterotype

no_bifid_ <- no_bifid %>%
  rownames_to_column("sample_id") %>%
  select(sample_id) %>%
  mutate(bifidotype = rep(max(enterotype$bifidotype) +1,nrow(.)))

# rbind no bifid to enterotype
enterotype <- enterotype %>%
  rbind(., no_bifid_)

name_file <- "enterotype_DMM"

write.csv(enterotype, file="enterotype_DMM", row.names=TRUE, col.names=TRUE, quote=FALSE)

enterotype <- fread("enterotype_DMM") %>%
  column_to_rownames("V1")

enterotype$bifidotype <- as.character(enterotype$bifidotype)

```


## Clusters distributions

```{r,message=F,warning=F}

distribution <- enterotype %>%
  select(-sample_id) %>%
  group_by(bifidotype) %>%
  table(.) %>%
  data.frame(.)

distribution$. <- as.character(distribution$.)

distr_clusters <- ggplot(distribution, aes(x = ., y = Freq, fill= .)) +  # Plot with values on top
  geom_bar(stat = "identity") +
  geom_text(aes(label = Freq, vjust=1.5),
            size=5,
            color="white") +
  labs( y = "Subject count", x = "Bifidobacterium based gut microbiome partitions") +
  theme_classic() +
  scale_fill_brewer("", type="qual", palette=2) +
  scale_x_discrete(limits = as.character(1:max(distribution$.))) +
  ylab("number of individuals") +
  guides(fill=FALSE) + annotate("text", x=7, y=700, label="(undetected\nbifid\npartition)", size=4)

distr_clusters

```

## Clusters distributions

```{r,message=F,warning=F, fig.width=14, fig.height = 4}

figure <- ggarrange(nb_clusters, distr_clusters, ncol = 2,
                              labels = c("A","B"),
                    nrow = 1)

figure

```

## Bifido relative abundance per cluster

```{r,message=F,warning=F,fig.width=15, fig.height=8}

phy_relative <- phy %>%
  prune_samples(sample_names(.) %in% enterotype$sample_id, .) %>%
  microbiome::transform(., "compositional")

tax <- data.frame(phy_relative@tax_table)

otu <- data.frame(phy_relative@otu_table)

```


## Heatmaps

```{r,message=F,warning=F,fig.width=9, fig.height=9}

matrix_int <- fitted(best, scale= TRUE) %>%
  log10(.) %>%
  round(.,2)

matrix_int <- round(matrix_int,2)

## replace col and rownames

colnames(matrix_int) <- paste("partitions ",1:ncol(matrix_int))

rownames(matrix_int) <- gsub("Bifidobacterium","B. ",rownames(matrix_int))


my_palette <- colorRampPalette(c("white", "yellow", "red"))(n = 400)

my_palette <- viridis_pal()

col_breaks = c(seq(0.81,1,length=100), # for red
seq(0,0.8,length=100),  # for yellow
seq(-1,0,length=100)) # for green

#png(filename = "heatmap_dmm.png", width = 1000, height = 1000, units = "px")

heatmap.2(matrix_int,
          Colv=FALSE,   # do not cluster partitions
  cellnote = matrix_int,  # same data set for cell labels
  main = "", # heat map title
  notecol="black",      # change font color of cell labels to black
  density.info="none",  # turns off density plot inside color legend
  trace="none",         # turns off trace lines inside the heat map
  margins = c(9,4),     # widens margins around plot
  col=my_palette,       # use on color palette defined earlier
  key.title = "", symm=F,symkey=F,symbreaks=F, scale="none",
  ## change legend scale size
  lhei=c(1,9), keysize=2, key.par = list(cex=0.5),
  key.xlab = "") -> bif_DMM_alpha_weight_df
  #key.xlab = expression(paste("log"[10],"(Value)")))           # turn off column clustering

#matrix_int %>%
  #pheatmap::pheatmap(show_colnames = FALSE, main = "eggnog annotations")


bif_DMM_alpha_weight_df$carpet %>%
  reshape2::melt() %>%
  ggplot() + geom_tile(aes(x=Var1,y=Var2,fill=10^value)) +
  geom_text(aes(x=Var1,y=Var2,label=round(10^value,4))) +
  scale_fill_viridis_c("DMM model\nalpha weights", trans="log10") +
  xlab("Bifidobacterium based gut micrbiome partitions") + 
  ylab(expression(paste(italic("Bifidobacterium"), " species"))) +
  theme_classic() +
  theme(axis.text.y = element_text(face="italic")) -> bif_DMM_alpha_weight_plot
 
bif_DMM_alpha_weight_plot

```

## Figure for manuscript



```{r fig.height=5, fig.width=14}

((criteria_dmm_plot /

distr_clusters) |

bif_DMM_alpha_weight_plot) + plot_annotation(tag_levels = "A")

```




## Bifido relative abundance per cluster

```{r,message=F,warning=F}

## select only bifid species and do sum of all bifid species abundances

bifido_names <- tax %>%
  filter(Genus == "Bifidobacterium") %>%
  rownames(.)

otu_bifid <- otu %>%
  filter(row.names(.) %in% bifido_names) %>%
  colSums(.) %>%
  data.frame(bif_rel = .) %>%
  rownames_to_column("sample_id")

sam_data <- phy_relative@sam_data %>%
  data.frame(.) %>%
  rownames_to_column("sample_id") %>%
  left_join(., otu_bifid, by="sample_id") %>%
  ## merge with cluster
  left_join(., enterotype, by="sample_id")
  

```

## Bifido five most abundant species per cluster

```{r,message=F,warning=F}

## top x number

x_nb <- "15"

colnames(matrix_int) <- 1:ncol(matrix_int)

top_5_ab <- matrix_int %>%
  data.frame(.)
colnames(top_5_ab) <- c( 1:ncol(matrix_int))

rownames(top_5_ab) = gsub("B.", "Bifidobacterium", rownames(top_5_ab))

## we select the commensal + animaliss species

species_to_select <- c("Bifidobacterium  animalis","Bifidobacterium  adolescentis","Bifidobacterium  pseudocatenulatum","Bifidobacterium  dentium","Bifidobacterium  longum","Bifidobacterium  bifidum","Bifidobacterium  catenulatum","Bifidobacterium  angulatum", "Bifidobacterium  breve")

data_sorted <- top_5_ab %>%
  rownames_to_column('species') %>%
  gather(key="cluster", value="value", 2:(ncol(matrix_int)+1)) %>%
  select(-value) %>%
  dplyr::rename('bifidotype' = 'cluster') %>%
  filter(species %in% species_to_select)

data_sorted$species <- str_squish(data_sorted$species)

  ## option where we select the top 5 abundant species per cluster
top_5 <- top_5_ab %>%
  rownames_to_column('species') %>%
  gather(key="cluster", value="value", 2:(ncol(matrix_int)+1)) %>%
  group_by(cluster) %>%
  top_n(x_nb, value) %>%
  select(-value) %>%
  dplyr::rename('bifidotype' = 'cluster')
top_5$bifidotype = gsub("cluster", "", top_5$bifidotype)
## remove multiple spaces
top_5$species <- str_squish(top_5$species)


## second way it to select all the commensal + animalis species

```

## search for the mean abundance of each species by cluster

```{r,message=F,warning=F}

otu_most <- otu

rownames(otu_most) <- get_species(tax, rownames(otu_most))

## group by cluster the otu table

otu_most <- otu_most %>%
  t(.) %>%
  data.frame(.) %>%
  select(contains("Bifidobacterium"))
colnames(otu_most) <- gsub("\\.", " ", colnames(otu_most))
## remove multiple spaces
colnames(otu_most) <- str_squish(colnames(otu_most))

# which(is.na(otu_most))

bifidotype <- enterotype

otu_most <- otu_most %>%
  rownames_to_column("sample_id") %>%
  right_join(., bifidotype, by = "sample_id") %>%
  filter(!is.na(bifidotype)) %>%
  gather(key="species", value="abundance", 2:33)
  

```

## search for the mean abundance of each species by cluster

```{r,message=F,warning=F}

data_sorted_abundance <- data_sorted %>%
  left_join(., otu_most, by = c("species", "bifidotype")) %>%
  filter(!is.na(abundance)) %>%
  select(-sample_id)

top_5_abundance <- top_5 %>%
  left_join(., otu_most, by = c("species", "bifidotype")) %>%
  filter(!is.na(abundance)) %>%
  select(-sample_id)

test <- top_5_abundance %>%
  group_by(bifidotype, species) %>%
  summarise(abundance = sum(abundance))

top_5_abundance$species <- gsub("Bifidobacterium", "B.", top_5_abundance$species)

data_sorted_abundance$species <- gsub("Bifidobacterium", "B.", data_sorted_abundance$species)

## change top_5 to lopg value, and a pseudo count of 10-6


top_5_abundance$abundance <- log10(top_5_abundance$abundance + 10^(-6))

data_sorted_abundance$abundance <- log10(data_sorted_abundance$abundance + 10^(-6))

p <- ggplot(top_5_abundance, aes(x = bifidotype,
                       y = abundance, fill = factor(species))) +
  geom_boxplot(outlier.shape = NA) +
  theme_minimal() +
  labs(fill= "Bifidobacterium species") #+
 # coord_cartesian(ylim = c(0, 0.25))
  
p

```

## ggplot

```{r,message=F,warning=F}

p <- ggplot(top_5_abundance, aes(x = species,
                       y = abundance, fill = species)) +
  geom_boxplot(outlier.shape = NA) +
  theme_minimal() +
  theme(axis.title.x = element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  labs(fill= "Bifidobacterium species",) +
  #scale_x_discrete(guide = guide_axis(n.dodge=3)) +
  facet_wrap(~bifidotype)
  
 # coord_cartesian(ylim = c(0, 0.25))
  
p


```
## ggplot

```{r,message=F,warning=F}

p <- ggplot(data_sorted_abundance, aes(x = bifidotype,
                       y = abundance, fill = bifidotype)) +
  #geom_boxplot() +
  geom_boxplot(outlier.shape = NA) +
  theme_minimal() +
  labs(fill= "Bifidotype", x="", y=expression(paste("log"[10],"(Abundance)"))) +
  facet_wrap(~species)
  
 # coord_cartesian(ylim = c(0, 0.25))
  
p

size = 3

p <- ggplot(data_sorted_abundance %>% filter(grepl("adolescentis|bifidum|catenulatum|longum|breve",species)), aes(x = bifidotype,
                       y = abundance, fill = bifidotype)) +
  geom_violin() +
  #geom_boxplot(outlier.shape = NA) +
  theme_minimal() +
  labs(fill= "Bifidotype", x="", y=expression(paste("log"[10],"(Abundance)"))) +
  stat_compare_means(label.y = 0.3, label.x ="2", size = size) +
  coord_cartesian(ylim=c(-6,1)) +
  facet_wrap(~species)
  
 # coord_cartesian(ylim = c(0, 0.25))
  
p

## facet_wrap par espèce
## et utiliser scale = free-Y
```

## tests

```{r,message=F,warning=F}

test <- top_5_abundance %>% filter(species %in% c("B. angulatum"))

```

## Heatmaps

```{r,message=F,warning=F}

sam_data$bifidotype <- as.character(sam_data$bifidotype)
  
p_relative_abundance <- ggplot(sam_data, aes(x = bifidotype,
                       y = bif_rel+10^-5, fill=bifidotype))+
  geom_boxplot(outlier.shape = NA) + scale_y_log10() +
  scale_fill_brewer(type="qual", palette=2) +
  #labs(x = "Bifidotype", y = NULL, title = "Relative abundance of Bifidobacterium",fill= "Bifidotype") + 
  #coord_cartesian(ylim=c(0,0.46)) +
  #scale_x_discrete(limits = as.character(1:max(sam_data$bifidotype))) +
  stat_compare_means(label.y = 0.45, label.x ="3") +
    theme_classic() +
  theme(legend.position = "none") +
  ylab("Bifidobacterium genus proportion") +
  xlab("Bifidobacterium based gut microbiome partitions")
  #stat_compare_means(label = "p.signif", method = "wilcox", ref.group = "5", label.y = 0.42)
  

p_relative_abundance


pairwise.wilcox.test(sam_data$bif_rel, sam_data$bifidotype) %>% 
  broom::tidy() %>%
  mutate(fdr=p.adjust(p.value, method="fdr")) %>%
  readr::write_tsv(file = "bif_rel_pairwise_wilcox.tsv")


```

## Boxplot Shannon

```{r,message=F,warning=F}


my_comparisons <- list( c("1", "2"), c("1", "3"), c("1", "4"), c("1","5"), c("1","6"), c("1","7"),
                        c("2", "3"), c("2", "4"), c("2","5"), c("2","6"), c("2","7"))

diversity <- phy_diversity %>%
  rownames_to_column("sample_id") %>%
  merge(., enterotype, by="sample_id") %>%
  select(-sample_id) %>%
  mutate(grp_healthy = case_when(
    bifidotype %in% c(1,2,6) ~ "high diversity",
    TRUE ~ "low diversity" )) 

diversity$bifidotype <- as.character(diversity$bifidotype)

p_diversity <- ggplot(diversity, aes(x = bifidotype,
                       y = Shannon)) +
  geom_boxplot(outlier.shape = NA) +
  theme_classic() +
  #stat_compare_means(label.y = 4.5, label.x ="3") +
  #labs(x = "Bifidotype", y = "", title = "Shannon diversity",color= "") + 
  #coord_cartesian(ylim=c(0.8,4.5)) +
  #stat_compare_means(label = "p.signif", method = "wilcox", ref.group = "1", label.y = 4.1)
  stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox")

p_diversity


p_diversity2 <- ggplot(diversity, aes(x = bifidotype,
                       y = Shannon, fill=bifidotype)) +
  geom_boxplot(outlier.shape = NA) +
  theme_classic() +
  stat_compare_means(label.y = 4.5, label.x ="3") +
  ylab("Alpha-diversity (Shannon)" ) +
  xlab("Bifidobacterium based gut microbiome partition") +
  scale_fill_brewer(type="qual", palette = 2) +
  guides(fill=FALSE)
  #labs(x = "Bifidotype", y = "", title = "Shannon diversity",color= "") + 
  #coord_cartesian(ylim=c(0.8,4.5)) +
  #stat_compare_means(label = "p.signif", method = "wilcox", ref.group = "1", label.y = 4.1)
  #stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "wilcox")

p_diversity2



pairwise.wilcox.test(diversity$Shannon, diversity$bifidotype) %>% 
  broom::tidy() %>%
  mutate(fdr=p.adjust(p.value, method="fdr")) %>%
  readr::write_tsv(file = "bif_shannon_pairwise_wilcox.tsv")

```

## Boxplot Shannon

```{r,message=F,warning=F, fig.width = 10, fig.height = 4}
figure <- ggarrange(p_relative_abundace, p_diversity,
                    labels = c("A","B"),
                    ncol =2,
                    nrow = 1)
figure
```


```{r,message=F,warning=F, fig.width = 6, fig.height = 8}

(p_relative_abundance/p_diversity2) + plot_annotation(tag_levels = "A")



```






## Description des clusters

```{r,message=F,warning=F, fig.height=9, fig.width = 15}

## set the geom text size to adjust graph

size_geom_text <- 2.5

distribution <- function(category, type){
  

  map <- phy@sam_data %>%
    data.frame(.) %>%
    rownames_to_column("sample_id") %>%
    merge(., enterotype, by="sample_id") %>%
    select(bifidotype,all_of(category)) %>%
    group_by_all(.) %>%
    summarise(rows = n(), .groups = 'drop') %>%
    group_by_at('bifidotype') %>%
    mutate(percent = rows/sum(rows)) %>%
    dplyr::rename(., 'category' = category)
  
  if(category=="age_category") {
  map <- map %>%
    mutate(category = factor(category, levels=c("newborn","child","schoolage","adult","senior"))) }
  
    p2 <- ggplot(map, aes(x = bifidotype,y=percent, fill = category)) +
    geom_bar(stat="identity", position=position_dodge()) +
    #geom_text(aes(label = paste0(rows)),vjust=-2, color="black",
    #        position = position_dodge(0.9),
    #        #size=size_geom_text) +
    #        size=ifelse(map$percent>0.035, size_geom_text,0)) +
    #  geom_text(aes(label = paste0(round(percent*100,0), "%")),vjust=-0.6, color="black",
    #        position = position_dodge(0.9),
    #        ## add a condition, if 
    #        #size=size_geom_text) +
    #        size=ifelse(map$percent>0.035, size_geom_text,0)) +
    #geom_text(aes(label = paste0(round(percent*100,1),"%"), vjust = 1), position = position_fill(vjust = 0.5), size =3) +
    theme_classic() +
    scale_fill_brewer(type = type) +
    labs(x = "Cluster K", y = "percent", fill = category) +
      #guides(y = "none") +
      ylim(0, 1.2)# + scale_y_continuous("prevalence per partition", label=scales::percent())
  
  p2

}

distribution('westernized', 'qual')


```

## Description des clusters

```{r,message=F,warning=F, fig.height=4, fig.width = 6}

distr_disease <- distribution('disease', "qual") +
  labs(x = "Bifidobacterium based gut microbiome partitions", y = "", fill= "disease") +
  coord_cartesian(ylim=c(0,1.02)) + scale_fill_brewer("Health status", type="qual")

distr_disease

```
## Description des clusters

```{r,message=F,warning=F, fig.height=4, fig.width = 6}

distr_age <- distribution('age_category', "qual") +
  labs(x = "Bifidobacterium based gut microbiome partitions",fill= "age category") +
  coord_cartesian(ylim=c(0,0.95)) +  scale_fill_brewer("Age category", type="div", palette = 7)


distr_age 

```
## Description des clusters

```{r,message=F,warning=F, fig.height=4, fig.width = 6}

distr_westernized <- distribution('westernized', "qual") +
  labs(x = "Bifidobacterium based gut microbiome partitions", fill= "westernized") +
  coord_cartesian(ylim=c(0,1.03)) +  scale_fill_brewer("Westernized\nlifestyle", type="qual", palette = 4)


distr_westernized

```

## descriptive metadata barplot 

```{r,message=F,warning=F, fig.width = 6, fig.height = 8}

figure <- ggarrange(distr_disease, distr_westernized,
                    labels = c("A","B"),
                    ncol =1,
                    nrow = 2)
figure


```

## figure for manuscript

```{r fig.height=5, fig.width=15}

(distr_age + scale_y_continuous("prevalence per partition", label=scales::percent) |
distr_westernized + scale_y_continuous("prevalence per partition", label=scales::percent, limits = c(0,1)) |
distr_disease + scale_y_continuous("prevalence per partition", label=scales::percent, limits = c(0,1))) +
  plot_annotation(tag_levels = "A")

ggsave("figure_distributions.pdf")

```



## table de contingence enterotypes * bifidotypes

```{r,message=F,warning=F}

## library for chi squared test ##

library("gplots")

## enterotypes are specific to bifid, so we'll call it bifidotype, comparing it to enterotypes

contingency <- function(file){
  
  
  bifidotype <- enterotype

  enterotype_whole_microbiota <- fread(file) %>%
    data.frame(.) %>%
    select(2:3)
  
  colnames(enterotype_whole_microbiota) <- c("sample_id","bifidotype")

  contingeny_table <- bifidotype %>%
    merge(., enterotype_whole_microbiota, by="sample_id") %>%
    select(- sample_id)

  names(contingeny_table) <- c('bifidotypes','enterotypes')

  contingeny_table <- contingeny_table %>%
    table(.)

  ## représentation graphique 

  test_chisq <- chisq.test(contingeny_table)

  return(test_chisq)
  
  }


## contingency between current bifidotypes(predict or whole - predict) against enterotypes
test_chisq_enterotypes <- contingency("enterotypes_curated_v3_prediction.csv")

## contingency between current bifidotypes(predict or whole - predict) against whole curated bifidotypes

test_chisq_whole <- contingency("enterotype_DMM")

```

## enterotypes * bifidotypes plot with pearson residuals

```{r,message=F,warning=F, fig.height=4, fig.width = 7}
library(corrplot)

corrplot(test_chisq_enterotypes$residuals, is.cor = FALSE)
#mtext("Bifidotypes x Enterotypes", at=5, line=-9, cex=1.5, side= 3)

```

## enterotypes * bifidotypes plot with pearson residuals

```{r,message=F,warning=F, fig.height=6, fig.width = 7}
library(corrplot)

corrplot(test_chisq_whole$residuals, is.cor = FALSE)
#mtext("Bifidotypes x Whole Bifidotypes", at=5, line=-7, cex=1.5, side= 3)

```



## test health_healthy

```{r,message=F,warning=F}

#library(readxl)



## we select the "control", which are "healthy" in the database :

healthy <- map %>%
  select(sample_id, disease) %>%
  filter(disease == 'Control')

print(paste0('le nombre de healthy dans curation est : ',nrow(healthy)))

healthy_healthy <- readxl::read_excel('metadata_healthy_healthy.xlsx') %>%
  data.frame(.) %>%
  dplyr::rename('sample_id' = Sample_ID) %>%
  merge(., healthy, by='sample_id')
  
print(paste0('le nombre de healthy-healthy dans curation est : ',nrow(healthy_healthy)))

```

## "validation interne"

```{r,message=F,warning=F}

load("DMM_files/best_fit_DMM_6clusters")

# we will check the predict funtion on the learning set, is it coherent with the enterotype file ?

best_fit <- best

count_learning_set <- count

## we need to select the common species between count and best fit object (estimate)
best_fit@fit$Estimate <- best_fit@fit$Estimate %>% as.data.frame %>% filter(row.names(.) %in% colnames(count_learning_set)) %>% as.matrix


predict_df <- predict(best_fit, count_learning_set, assign=TRUE) %>%
  as.data.frame

colnames(predict_df) <- 1:6

## create a column which will check which bifidotype corresponds to each subject (maximum of fit score per row)

predict_df$assigned_bifidotype = colnames(predict_df)[apply(predict_df, 1, which.max)]


predict_df

predict_df <- predict_df %>%
  select(assigned_bifidotype) %>%
  rownames_to_column("sample_id")

## check the bifido file, join the two
enterotype

test_join <- predict_df %>%
  inner_join(enterotype, by = "sample_id")

test_join

if (setequal(test_join$assigned_bifidotype,test_join$bifidotype)){
  
  print("OK - prediction on learning set corresponds to clustering results")
  
}else{
  
  print("warning ! - prediction on learning set DOES NOT correspond to clustering results")
  
}

```


```{r}


table(test_join$assigned_bifidotype,test_join$bifidotype) 

```



## bifidotype per study


```{r}

sampleMetadata %>% 
  select(study_name, sample_id) %>% 
  merge(bifidotypes, by="sample_id") %>% 
  select(study_name,bifidotype) %>% 
  group_by(study_name) %>%
  mutate(n=n()) %>%
  filter(n>500) %>%
  select(study_name,bifidotype) %>%
  ungroup() %>%
  table() %>%
  prop.table(1) %>%
  reshape2::melt() %>%
  mutate(study_name = ifelse(study_name=="LifeLinesDeep_2016","ZhernakovaA_2016",as.character(study_name))) %>%
  ggplot() + 
  geom_bar(aes(x=as.character(bifidotype), y=value, fill=study_name), stat="identity", position = "dodge") +
  theme_classic() +
  scale_y_continuous("% of individuals",labels=scales::percent) +
  xlab("Bididobacterium-based gut microbiome partitions") +
  scale_fill_brewer("", type="qual", palette = 3)
  

ggsave("partitions_bif_studies.pdf")

```


# sequencing depth vs partitions

```{r}

sampleMetadata %>% 
  select(number_reads, sample_id) %>% 
  merge(bifidotypes, by="sample_id") %>% 
  select(bifidotype, number_reads) %>%
  with(., pairwise.wilcox.test(.$number_reads, .$bifidotype, data=.)) %>%
  broom::tidy() %>%
  mutate(fdr=p.adjust(p.value,method="fdr")) %>%
  readr::write_tsv(file = "bif_number_reads_wilcox.tsv")
  



```