hierarchical enterotyping

load genus abundance

genus_path = system.file("data-raw/qiime/generated-files-20190512/taxa/genus.qza", package = "gutzilla")

genus = qiime2R::read_qza(genus_path)$data %>% as.data.frame %>% tibble::rownames_to_column("taxa")  %>% as_tibble()

load metadata (cleaned) and add united nation subregions

country_info = readr::read_tsv(system.file(package="gutzilla", "data-raw/countryInfo.txt"), comment = "#", col_names = FALSE, na = "") %>%
  select(X5,X9) %>%
  dplyr::rename(country="X5", continent="X9" )

Avis : One or more parsing issues, call `problems()` on your data frame for details, e.g.:
  dat <- vroom(...)
  problems(dat)Rows: 252 Columns: 14── Column specification ──────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (12): X1, X2, X3, X4, X5, X6, X9, X10, X11, X12, X13, X14
dbl  (2): X7, X8
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

data("UNSD_countries")

metadata = readr::read_csv2(system.file("data-raw/Metadata_10317_20191022-112414_curatedv4_VSv1.csv", package = "gutzilla"))

ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.
Rows: 24251 Columns: 278── Column specification ──────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ";"
chr  (251): SAMPLE_NAME, SAMPLE_TYPE, SEX, AGE_CAT, BMI_CAT, COUNTRY, ALCOHOL_TYPES_BEERCIDER, ALCOHOL_TYPES_R...
dbl    (7): AGE_YEARS, HEIGHT_CM, WEIGHT_KG, BIRTH_YEAR, HOST_TAXID, TAXON_ID, VIOSCREEN_NBITEMS
num    (3): BMI, VIOSCREEN_KCAL, VIOSCREEN_TIME
lgl   (16): HOST, ALCOHOL_TYPES, ALLERGIC_TO, FERMENTED_CONSUMED, FERMENTED_PRODUCE_COMMERCIAL, FERMENTED_PROD...
dttm   (1): COLLECTION_TIMESTAMP
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#to do: check issue parsing IBD variables

colnames(metadata) = stringr::str_to_lower(colnames(metadata))



metadata %>%
  select(country_of_birth) %>%
  merge(UNSD_countries %>% 
          select(`Country or Area`, `Sub-region Name`, `Region Name`) %>% 
          mutate(`Country or Area` =  gsub("United Kingdom of Great Britain and Northern Ireland", "United Kingdom", `Country or Area`)) %>% 
          mutate(`Country or Area` =  gsub("United States of America", "USA", `Country or Area`)) %>%
          mutate(`Sub-region Name` =  ifelse(`Sub-region Name` %in% c("Central Asia", "Southern Asia"), "Central and Southern Asia",`Sub-region Name`)) %>%
          mutate(`Sub-region Name` =  ifelse(`Sub-region Name` %in% c("Melanesia", "Micronesia","Polynesia"), "Oceania (others)",`Sub-region Name`)), 
        by.x="country_of_birth", by.y="Country or Area") %>% 
  group_by(`Region Name`,`Sub-region Name`) %>% summarise(n=n()) %>%
  ungroup() %>%
  ggplot() + geom_bar(aes(y=n, x=`Sub-region Name`), stat="identity") + scale_y_log10() +
  facet_wrap(~`Region Name`, scales = "free_x", nr=1) + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  xlab("Region of Birth\n(United Nations-m49)") +
  ylab("Number of stool sample")

`summarise()` has grouped output by 'Region Name'. You can override using the `.groups` argument.

  
ggsave("Region_of_birth.pdf")

Saving 7.29 x 4.5 in image

metadata <- metadata %>%
  merge(UNSD_countries %>% 
          select(`Country or Area`, `Sub-region Name`, `Region Name`) %>% 
          mutate(`Country or Area` =  gsub("United Kingdom of Great Britain and Northern Ireland", "United Kingdom", `Country or Area`)) %>% 
          mutate(`Country or Area` =  gsub("United States of America", "USA", `Country or Area`)) %>%
          mutate(`Sub-region Name` =  ifelse(`Sub-region Name` %in% c("Central Asia", "Southern Asia"), "Central and Southern Asia",`Sub-region Name`)) %>%
          mutate(`Sub-region Name` =  ifelse(`Sub-region Name` %in% c("Melanesia", "Micronesia","Polynesia"), "Oceania (others)",`Sub-region Name`)), 
        by.x="country_of_birth", by.y="Country or Area")
  

table(metadata$`Sub-region Name`)

      Australia and New Zealand       Central and Southern Asia                    Eastern Asia 
                            632                             356                             140 
                 Eastern Europe Latin America and the Caribbean                 Northern Africa 
                            328                             286                              39 
               Northern America                 Northern Europe                Oceania (others) 
                          14244                            5734                              11 
             South-eastern Asia                 Southern Europe              Sub-Saharan Africa 
                            227                             267                             163 
                   Western Asia                  Western Europe 
                            105                             778 

cowplot::plot_grid(

metadata %>%
  mutate(age_cat = age_cat %>% forcats::fct_relevel(c("baby","child","teen"), after=0 )) %>%
  filter(!is.na(age_cat)) %>%
  ggplot() + geom_bar(aes(x=age_cat)) + theme_bw(),

metadata %>%
  filter(!is.na(age_cat), !is.na(sex)) %>%
  ggplot() + geom_bar(aes(x=sex)) + theme_bw(), rel_widths = c(2,1))

NA
NA

import outliers to filter them out and select the top 30 genus based on read mass.

outliers = readLines(con="outliers_samples.txt")

n=30 # n samples per group

metadata %>% 
  filter(!(sample_name %in% outliers)) %>%
  #filter(!(`#SampleID` %in% outliers)) %>%
  merge(country_info, by.x="country_of_birth", by.y="country") %>%
  select(sample_name,age_cat,sex, `Sub-region Name`,country_of_birth) %>% 
  #select(`#SampleID`,age_cat,sex, continent,country_of_birth) %>% 
  filter(sex %in% c("male","female")) %>%
  filter(!is.na(age_cat) & age_cat != "Not provided") %>%
  group_by(age_cat,sex,`Sub-region Name`) %>% 
  sample_n(if(n() < n) n() else n) %>%
  #sample_n(size=30, replace=TRUE) %>%
  #unique() %>% 
  pull(sample_name) -> samples_id_select

Error in h(simpleError(msg, call)) : 
  erreur d'�valuation de l'argument 'y' lors de la s�lection d'une m�thode pour la fonction 'merge' : object 'country_info' not found

load DMM model fit from boostrap

 load("fit_genus_bootstrap.rda")

Avis : cannot open compressed file 'fit_genus_bootstrap.rda', probable reason 'No such file or directory'Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection

extract BIC and Laplace coefficients from boostrap models

BIC_list = map(fit_genus_bootstrap, function(z) map(z, function(y) sapply(y, function(x){attr(x,"goodnessOfFit")[["BIC"]]}))) 
lplc_list = map(fit_genus_bootstrap, function(z) map(z, function(y) sapply(y, function(x){attr(x,"goodnessOfFit")[["Laplace"]]})))
k = map(fit_genus_bootstrap, function(z) map(z, function(y) sapply(y, function(x){attr(x,"mixture")$Weight %>% length})))[[1]][[1]]

Error in map(fit_genus_bootstrap, function(z) map(z, function(y) sapply(y,  : 
  subscript out of bounds

plot BIC and laplace results

BIC and Laplace have a global optimun among differents seeds and bootstraps. All bootstraps found at least one model optimal fit at k=19. Majority of best model based on Laplace were k=19. For BIC, k=19 arrived in second postition. For downstream, analysis we then choose k=19.

For downstream analysis, we choose models among bootstrap where BIC and Laplace showed k=19 as optimum. Among models, we choose the worst one to avoid model overfitting.

###################################################
### code chunk number 11: posterior-mean-diff
###################################################
p0 <- fitted(fit_genus_list[[seed]][[1]], scale=TRUE)     # scale by theta
pbest <- fitted(fit_genus_list[[seed]][[best_genus_lplc_idx]], scale=TRUE)
colnames(pbest) <- paste("m", 1:best_genus_lplc, sep="")
(meandiff <- colSums(abs(pbest - as.vector(p0))))
sum(meandiff)


###################################################
### code chunk number 12: table-1
###################################################
diff <- rowSums(abs(pbest - as.vector(p0)))
o <- order(diff, decreasing=TRUE)
cdiff <- cumsum(diff[o]) / sum(diff)
df <- head(cbind(Mean=p0[o], pbest[o,], diff=diff[o], cdiff), 30)
df %>% as.data.frame()

df_gutzilla = df %>% as.data.frame()

m_order = df %>% as.data.frame() %>% select(-Mean, -diff,-cdiff)  %>% BiotypeR::dist.JSD() %>% hclust(method="complete") %>% unclass %>% .$order
#m_order = df %>% as.data.frame() %>% select(-Mean, -diff,-cdiff)  %>% t %>% dist() %>% hclust(method="complete") %>% unclass %>% .$order 


df %>% as.data.frame() %>%
  tibble::rownames_to_column("taxa") %>%
  mutate(taxa=factor(taxa, levels=row.names(df))) %>% 
  reshape2::melt(id.vars="taxa") %>%
  filter(!variable %in% c("Mean","diff","cdiff")) %>%
  mutate(variable = variable %>% as.character) %>% 
  mutate(variable = variable %>% factor(levels=paste0("m",1:19)[m_order]) )  -> genus_alpha_weight_gutzilla

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

df %>% as.data.frame() %>%
  tibble::rownames_to_column("taxa") %>%
  mutate(taxa=factor(taxa, levels=row.names(df))) %>% 
  reshape2::melt(id.vars="taxa") %>%
  filter(!variable %in% c("Mean","diff","cdiff")) %>%
  #group_by(taxa) %>%
  #mutate(total = sum(value)) %>%
  #mutate(value=value/total) %>%
  ggplot() + 
  geom_tile(aes(x=variable,y=taxa,fill=value)) + 
  scale_fill_viridis_c("DMM model alpha") +
  ylab("taxa") +
  xlab("DMM components")
  

look at theta per component

#fitted(fit_genus_list[[1]][[best_genus_lplc]])

mixturewt(fit_genus_list[[seed]][[1]])

mixturewt(fit_genus_list[[seed]][[best_genus_lplc_idx]])

mixturewt(fit_genus_list[[seed]][[best_genus_lplc_idx]]) %>% 
  tibble::rowid_to_column("k") %>% 
  ggplot() + 
  geom_bar(aes(x=k,y=theta),stat="identity") + 
  geom_hline(yintercept = mixturewt(fit_genus_list[[seed]][[1]])$theta, col="red") +
  cowplot::theme_cowplot() +
  ylab("theta") +
  xlab("DMM partitions") -> figS8E

figS8E

metadata %>%
  merge(enterotypes, by.x="sample_name", by.y="row.names", all.y="TRUE") %>%
  group_by(`Sub-region Name`) %>%
  mutate(total=n()) %>%
  group_by(`Sub-region Name`,Enterotypes_id,total) %>%
  summarise(prop=n()) %>%
  mutate(prop=prop/total) %>%
  ggplot() + geom_tile(aes(x=Enterotypes_id %>% as.character(),y=`Sub-region Name`,fill=prop)) + scale_fill_viridis_c()


metadata %>%
  merge(enterotypes, by.x="sample_name", by.y="row.names", all.y="TRUE") %>%
  filter(age_cat != "baby") %>%
  group_by(age_cat) %>%
  mutate(total=n()) %>%
  group_by(age_cat,Enterotypes_id,total) %>%
  summarise(prop=n()) %>%
  mutate(prop=prop/total) %>%
  ungroup() %>%
  group_by(Enterotypes_id) %>%
  mutate(prop2=sum(prop)) %>%
  mutate(prop=prop/prop2) %>%
  ggplot() + geom_tile(aes(x=Enterotypes_id %>% as.character(),y=age_cat,fill=prop)) + scale_fill_viridis_c()


metadata %>%
  merge(enterotypes, by.x="sample_name", by.y="row.names", all.y="TRUE") %>%
  filter(age_cat != "baby") %>%
  filter(!is.na(bmi_cat), bmi_cat != "Not provided") %>%
  group_by(bmi_cat) %>%
  mutate(total=n()) %>%
  group_by(bmi_cat,Enterotypes_id,total) %>%
  summarise(prop=n()) %>%
  mutate(prop=prop/total) %>%
  ungroup() %>%
  group_by(bmi_cat) %>%
  mutate(prop2=sum(prop)) %>%
  mutate(prop=prop/prop2) %>%
  ggplot() + geom_tile(aes(x=Enterotypes_id %>% as.character(),y=bmi_cat,fill=prop)) + scale_fill_viridis_c()


metadata %>%
  merge(enterotypes, by.x="sample_name", by.y="row.names", all.y="TRUE") %>%
  filter(age_cat != "baby") %>%
  filter(!is.na(gluten), gluten != "Not provided") %>%
  group_by(gluten) %>%
  mutate(total=n()) %>%
  group_by(gluten,Enterotypes_id,total) %>%
  summarise(prop=n()) %>%
  mutate(prop=prop/total) %>%
  ungroup() %>%
  group_by(gluten) %>%
  mutate(prop2=sum(prop)) %>%
  mutate(prop=prop/prop2) %>%
  ggplot() + geom_tile(aes(x=Enterotypes_id %>% as.character(),y=gluten,fill=prop)) + scale_fill_viridis_c()

metadata %>%
  merge(enterotypes, by.x="sample_name", by.y="row.names", all.y="TRUE") %>%
  filter(age_cat != "baby") %>%
  filter(!is.na(types_of_plants), types_of_plants != "Not provided") %>%
  group_by(types_of_plants) %>%
  mutate(total=n()) %>%
  group_by(types_of_plants,Enterotypes_id,total) %>%
  summarise(prop=n()) %>%
  mutate(prop=prop/total) %>%
  ungroup() %>%
  group_by(types_of_plants) %>%
  mutate(prop2=sum(prop)) %>%
  mutate(prop=prop/prop2) %>%
  ggplot() + geom_tile(aes(x=Enterotypes_id %>% as.character(),y=types_of_plants,fill=prop)) + scale_fill_viridis_c()

metadata %>%
  merge(enterotypes, by.x="sample_name", by.y="row.names", all.y="TRUE") %>%
  filter(age_cat != "baby") %>%
  filter(!is.na(diet_type), diet_type != "Not provided") %>%
  group_by(diet_type) %>%
  mutate(total=n()) %>%
  group_by(diet_type,Enterotypes_id,total) %>%
  summarise(prop=n()) %>%
  mutate(prop=prop/total) %>%
  ungroup() %>%
  group_by(diet_type) %>%
  mutate(prop2=sum(prop)) %>%
  mutate(prop=prop/prop2) %>%
  ggplot() + geom_tile(aes(x=Enterotypes_id %>% as.character(),y=diet_type,fill=prop)) + scale_fill_viridis_c()

to do: compute enterotype proportion for all categories and do a PCA analysis with minimal spanning tree

metadata %>%
  merge(enterotypes, by.x="sample_name", by.y="row.names", all.y="TRUE") %>%
  select(sample_name,age_cat,sex,`Sub-region Name`,bmi_cat, Enterotypes_id) %>%
  reshape2::melt(id.vars=c("sample_name","Enterotypes_id")) %>%
  filter(!is.na(value), value != "Not provided") %>%
  select(-variable, -`sample_name`) %>%
  group_by(Enterotypes_id,value) %>%
  filter(!is.na(value)) %>%
  #mutate(total=n()) %>%
  #group_by(Enterotypes_id,value,total) %>%
  #summarise(prop=n()) %>%
  reshape2::dcast(Enterotypes_id~value) %>%
  mutate_at(-1,~./(sum(.))) -> et_repartition

et_repartition %>%
  tibble::column_to_rownames("Enterotypes_id") %>%
  apply(1,function(x) x/sum(x)) %>%
  t %>%
  as.data.frame %>%
  ade4::dudi.pca(.,scan=F, nf=3) -> et_metadata_pca


ade4::scatter(et_metadata_pca)

ade4::s.label(et_metadata_pca$li, yax=2)

df %>% as.data.frame() %>%
  select(-diff,-cdiff,-Mean) %>%
  mutate_at(1:19,~log10(.)) %>%
  t %>%
  as.data.frame %>%
  ade4::dudi.pca(scannf = F,nf=3) %>% .$li %>% ade4::s.label()




df %>% as.data.frame() %>%
  select(-diff,-cdiff,-Mean) %>%
  tibble::rownames_to_column("taxa") %>%
  #mutate(taxa2=taxa) %>%
  #mutate(taxa2=ifelse())
    as_tibble %>%
  .[c(1:3,7,11),] %>%
  reshape2::melt(id.vars="taxa") %>%
  group_by(variable) %>%
  mutate(ymax=cumsum(value)) %>%
  mutate(ymin=c(0, head(ymax, n=-1))) %>%
  #ggplot() + geom_bar(aes(y=value,x=variable,fill=taxa),stat="identity") 
  ggplot(aes(ymax=ymax, ymin=ymin, xmax=4, xmin=3, fill=taxa %>% strsplit(x=.,split=";") %>% sapply(function(x)x[6])  )) + 
  geom_rect() +      
  coord_polar(theta="y") + 
  facet_wrap(~variable) + 
  xlim(c(2, 4)) + theme_void() +
  scale_fill_brewer("",type="qual")

ggsave("agp_microbiota_types.pdf")

reclassify all dataset

# fit_genus_list[[1]][[best_genus_lplc]] %>%
#   mixture(assign=TRUE) %>% as.data.frame %>% set_colnames(c("Enterotypes_id"))

genus_select=rownames(fitted(fit_genus_bootstrap[[1]][[1]][[best_genus_lplc]]))
genus %>%
      select(-outliers) %>%
      tibble::column_to_rownames("taxa") %>% .[genus_select,]-> genus_df

genus_df %>% colnames() %>% data.frame(sample_name=.) -> enterotypes_clusters_df


 k=1
for(i in 1:5) {
  
  for(j in 1:5){
    

   
    
    predict(fit_genus_bootstrap[[i]][[j]][[best_genus_lplc]], t(genus_df), assign = TRUE) %>%
      reshape2::melt() %>%
      group_by(Var1) %>%
      filter(value==max(value)) %>%
      dplyr::rename(`#SampleID`="Var1", Enterotypes_id="Var2") %>%
      select(`#SampleID`,Enterotypes_id) -> enterotypes_prediction_max
    
    colnames(enterotypes_prediction_max) = c("sample_name",paste0("Enterotypes_id_",k))
    
    enterotypes_clusters_df=
    enterotypes_clusters_df %>%
      merge(enterotypes_prediction_max, by="sample_name")
    
    cat(k,"\n")
    k=k+1
   
    
  }
  
  
}

 enterotypes_clusters_dist = matrix(nr=25,nc=25)  

 for(i in 1:25) {
     for (j in 1:25) {
     
     enterotypes_clusters_dist[i,j] = mclust::adjustedRandIndex(enterotypes_clusters_df[-1][,i],enterotypes_clusters_df[-1][,j])
     
     
   }
 }
 
 

classify data with the choosen model

predict(fit_genus_list[[1]][[14]], t(genus_df), assign = TRUE) %>%
  reshape2::melt() %>%
  group_by(Var1) %>%
  filter(value==max(value)) %>%
  dplyr::rename(`sample_name`="Var1", Enterotypes_id="Var2") %>%
  mutate(set=ifelse(sample_name %in% row.names(attr(fit_genus_list[[1]][[14]],"group")), "train", "predicted")) -> enterotypes_prediction


enterotypes_prediction %>%
  ggplot() + geom_histogram(aes(x=value, fill=value>0.80), bins = 50) + 
  #facet_wrap(~set) + 
  scale_x_continuous(limits = c(0,1.05), labels = scales::percent)

enterotypes_prediction %>%
  ggplot() + 
  geom_histogram(aes(x=value), bins = 50) + 
  #facet_wrap(~set) + 
  scale_x_continuous(limits = c(0,1.05), labels = scales::percent) +
  cowplot::theme_cowplot() -> figS8F

export enterotypes with outliers

enterotypes_prediction %>%
  as.data.frame %>%
  rbind(data.frame(sample_name = outliers,Enterotypes_id=NA,value=NA,set="outliers", stringsAsFactors = FALSE)) %>%
  write.csv2(file="enterotypes_prediction_outliers.csv")

Enterotypes vs Age and Richness

Age

metadata %>%
  merge(enterotypes, by.x="sample_name", by.y="row.names", all.y="TRUE") %>%
  select(Enterotypes_id, age_years) %>% ggplot() + geom_boxplot(aes(x=Enterotypes_id%>% as.character, y=age_years))
  group_by(Enterotypes_id) %>%
  summarise(age_years=mode(age_years))
  

richness

shannon_path = system.file("data-raw/qiime/generated-files-20190512/alpha/shannon.qza", package = "gutzilla")
shannon = qiime2R::read_qza(shannon_path)$data %>% as.data.frame

merge(enterotypes_prediction,shannon,by.x="sample_name", by.y="row.names")  -> enterotypes_alpha_div
  
enterotypes_alpha_div %>% 
  ggplot()   + 
  geom_boxplot(aes(x=paste0("M",Enterotypes_id %>% as.character),y=shannon, fill=set)) +
  xlab("DMM partitions") + ylab("alpha-diversity (shannon)") +
  cowplot::theme_cowplot() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) -> figS8G
figS8G

figure supp 8

(figS8ABCD/(figS8E+figS8F+figS8G + plot_layout(widths = c(1,1,1.5)))) + plot_layout(heights = c(2,1)) + plot_annotation(tag_levels = 'a') -> figS8

figS8
ggsave("figS8.pdf")

check fibres

fibre_path = system.file("data-raw/vioscreen/FICHIER_JULIEN_FIBER_HOSTSUBJECTID_2.xlsx", package = "agp")

readxl::read_xlsx(fibre_path) %>%
  select(-FIBER...6) %>%
  dplyr::rename(FIBER="FIBER...4") -> fibre_df
  

colnames(fibre_df) = colnames(fibre_df) %>% stringr::str_to_lower()

n=1275 with fiber and enterotypes data

fibre_df %>%
  merge(enterotypes_alpha_div, by="sample_name") %>% 
  ggplot() + 
  geom_boxplot(aes(y=shannon,x=paste0("Q",quartile_fiber),fill=paste0("Q",quartile_fiber))) +
  xlab("fiber\nquartiles") + ylab("alpha-diversity (shannon)") +
  scale_fill_brewer("fiber\nquartiles",type="div",palette = 2)


fibre_df %>%
  merge(enterotypes_alpha_div, by="sample_name") %>% 
  ggplot() + 
  geom_boxplot(aes(y=shannon,x=paste0("m",Enterotypes_id %>% as.character),fill=paste0("Q",quartile_fiber))) +
  xlab("Microbiota DMM based clusters") + ylab("alpha-diversity (shannon)") +
  scale_fill_brewer("fiber\nquartiles",type="div",palette = 2)


fibre_df %>%
  merge(enterotypes_alpha_div, by="sample_name") %>% 
  group_by(Enterotypes_id) %>%
  do(with(.,kruskal.test(.$shannon,.$quartile_fiber, method="spearman") %>% broom::tidy() ))

fibre_df %>%
  merge(enterotypes_alpha_div, by="sample_name") %>% 
  group_by(Enterotypes_id) %>%
  summarise(n=n())

fibre_df %>%
  merge(enterotypes_alpha_div, by="sample_name") %>% 
  filter(Enterotypes_id %in% c(2,5,6,7,13)) %>%
  ggplot() + 
  geom_boxplot(aes(y=shannon,x=paste0("m",Enterotypes_id %>% as.character),fill=paste0("Q",quartile_fiber))) +
  xlab("Microbiota DMM based clusters") + ylab("alpha-diversity (shannon)") +
  scale_fill_brewer("fiber\nquartiles",type="div",palette = 2)

check type of plant

metadata %>%
  select(sample_name,types_of_plants) %>%
   mutate(types_of_plants=types_of_plants%>% forcats::fct_rev() %>% forcats::fct_relevel("More than 30", after=Inf) %>% forcats::fct_relevel("21 to 30", after=3)) %>%
  na.omit() %>%
  merge(enterotypes_alpha_div, by="sample_name") %>% 
  ggplot() + 
  geom_boxplot(aes(y=shannon,x=types_of_plants,fill=types_of_plants)) +
  xlab("nb type of plants per week") + ylab("alpha-diversity (shannon)") +
  scale_fill_brewer("Type of plants",type="div",palette = 2) +
  xlab("Type of plants") + ylab("alpha-diversity (shannon)") + ggtitle("all pop")

metadata %>%
  select(sample_name,types_of_plants) %>%
  merge(fibre_df, by="sample_name") %>%
  mutate(types_of_plants=types_of_plants%>% forcats::fct_rev() %>% forcats::fct_relevel("More than 30", after=Inf) %>% forcats::fct_relevel("21 to 30", after=3)) %>%
  na.omit() %>%
  merge(enterotypes_alpha_div, by="sample_name") %>% 
  ggplot() + 
  geom_boxplot(aes(y=shannon,x=types_of_plants,fill=types_of_plants)) +
  xlab("nb type of plants per week") + ylab("alpha-diversity (shannon)") +
  scale_fill_brewer("Type of plants",type="div",palette = 2) +
  xlab("Type of plants ") + ylab("alpha-diversity (shannon)") + ggtitle("with vioscreened only")



metadata %>%
  select(sample_name,types_of_plants) %>%
   mutate(types_of_plants=types_of_plants%>% forcats::fct_rev() %>% forcats::fct_relevel("More than 30", after=Inf) %>% forcats::fct_relevel("21 to 30", after=3)) %>%
  na.omit() %>%
  merge(enterotypes_alpha_div, by="sample_name") %>% 
  ggplot() + 
  geom_boxplot(aes(y=shannon,x=paste0("m",Enterotypes_id %>% as.character),fill=types_of_plants)) +
  xlab("Microbiota DMM based clusters") + ylab("alpha-diversity (shannon)") +
  scale_fill_brewer("Type of\nplant",type="div",palette = 2) + ggtitle("all pop")



metadata %>%
  select(sample_name,types_of_plants) %>%
  merge(fibre_df, by="sample_name") %>%
   mutate(types_of_plants=types_of_plants%>% forcats::fct_rev() %>% forcats::fct_relevel("More than 30", after=Inf) %>% forcats::fct_relevel("21 to 30", after=3)) %>%
  na.omit() %>%
  merge(enterotypes_alpha_div, by="sample_name") %>% 
  ggplot() + 
  geom_boxplot(aes(y=shannon,x=paste0("m",Enterotypes_id %>% as.character),fill=types_of_plants)) +
  xlab("Microbiota DMM based clusters") + ylab("alpha-diversity (shannon)") +
  scale_fill_brewer("Type of\nplant",type="div",palette = 2) + ggtitle("with vioscreened only")

metadata %>%
  select(sample_name,types_of_plants) %>%
  merge(fibre_df, by="sample_name") %>%
   mutate(types_of_plants=types_of_plants%>% forcats::fct_rev() %>% forcats::fct_relevel("More than 30", after=Inf) %>% forcats::fct_relevel("21 to 30", after=3)) %>%
  na.omit() %>%
  merge(enterotypes_alpha_div, by="sample_name") %>%
  filter(Enterotypes_id %in% c(1,2,5,7,13,15)) %>%
  ggplot() + geom_point(aes(x=fiber,y=shannon,col=types_of_plants)) + 
  scale_color_brewer("Type of\nplant",type="div",palette = 2) + 
  facet_wrap(~paste0("m",Enterotypes_id), scales = "free_y") + theme_dark() +
  ggtitle("example of enterotypes")

rescale enterotypes using world demographic data

we will use UNSD subregion to rescale enterotypes proportion observed

data("UNSD_pop_subregion")

UNSD_pop_subregion %>%
          mutate(Subregion = ifelse(Subregion %in% c("Melanesia","Micronesia","Polynesia"), "Oceania (others)",Subregion)) %>%
          mutate(Subregion =  ifelse(Subregion %in% c("Central Asia", "Southern Asia"), "Central and Southern Asia",Subregion)) %>%
          mutate(Subregion =  ifelse(Subregion %in% c("South-Eastern Asia"), "South-eastern Asia",Subregion)) %>%
  ggplot() + geom_bar(aes(y=n_pop_2020, x=Subregion), stat="identity") +
  theme(axis.text.x = element_text(angle = -90, vjust=0.5))



metadata %>%
  merge(enterotypes_prediction, by.x="sample_name", all.y="TRUE") %>%
  select(Enterotypes_id,`Sub-region Name`) %>%
  na.omit %>%
  group_by(`Sub-region Name`) %>%
  summarise(total=n()) %>%
  ggplot() + geom_bar(aes(y=total, x=`Sub-region Name`), stat="identity") +
  theme(axis.text.x = element_text(angle = -90, vjust=0.5))

metadata %>%
  merge(enterotypes_prediction, by.x="sample_name",  all.y="TRUE") %>%
  select(`Sub-region Name`) %>%
  na.omit %>%
  group_by(`Sub-region Name`) %>%
  summarise(total=n()) %>%
  merge(UNSD_pop_subregion %>%
          mutate(Subregion = ifelse(Subregion %in% c("Melanesia","Micronesia","Polynesia"), "Oceania (others)",Subregion)) %>%
          mutate(Subregion =  ifelse(Subregion %in% c("Central Asia", "Southern Asia"), "Central and Southern Asia",Subregion)) %>%
          mutate(Subregion =  ifelse(Subregion %in% c("South-Eastern Asia"), "South-eastern Asia",Subregion)) %>%
                   group_by(Subregion) %>%
                   summarise(n_pop_2020 = sum(n_pop_2020)),
  by.x="Sub-region Name", by.y="Subregion", all.x=TRUE) %>%
  mutate_if(is.numeric,function(x) x/sum(x) ) %>%
  dplyr::rename(n_pop_agp_birth="total") -> region_n_pop_agp_unsd

  region_n_pop_agp_unsd %>%
  reshape2::melt(id.vars="Sub-region Name") %>%
  ggplot() + 
  geom_bar(aes(x=`Sub-region Name`, y=value, fill=variable), stat="identity", position = "dodge") +
  theme(axis.text.x = element_text(angle = -90, vjust=0.5)) 

  
  
  

metadata %>%
  merge(enterotypes_prediction, by.x="sample_name",  all.y="TRUE") %>%
  select(`Sub-region Name`,Enterotypes_id,set) %>%
  ggplot() + geom_bar(aes(x=Enterotypes_id,fill=set)) 

metadata %>%
  merge(enterotypes_prediction, by.x="sample_name",  all.y="TRUE") %>%
  select(`Sub-region Name`,Enterotypes_id,set) %>%
  na.omit %>%
  group_by(`Sub-region Name`,Enterotypes_id) %>%
  summarise(n_per_et=n()) %>%
  ungroup() %>%
  group_by(`Sub-region Name`) %>%
  mutate(n_per_region=sum(n_per_et)) %>%
  ungroup() %>%
  #mutate(n=n %>% (function(x) x / sum(x)) ) %>%
  merge(UNSD_pop_subregion %>%
          filter(!(Subregion == "Latin America and the Caribbean" & Type =="Region")) %>%
          mutate(Subregion = ifelse(Subregion %in% "Oceania (Excluding Australia and New Zealand)", "Oceania (others)",Subregion)) %>%
          #mutate(Subregion =  ifelse(Subregion %in% c("Central Asia", "Southern Asia"), "Central and Southern Asia",Subregion)) %>%
          mutate(Subregion =  ifelse(Subregion %in% c("South-Eastern Asia"), "South-eastern Asia",Subregion)) ,
                   
                   #summarise(n_pop_2020 = sum(n_pop_2020)),
  by.x="Sub-region Name", by.y="Subregion", all.x=TRUE) %>%
  mutate(n_estimated_per_pop = (n_per_et * n_pop_2020)/n_per_region) %>%
  select(-n_per_region,-Type,-n_pop_2020) %>%
  mutate(p_per_et=n_per_et/sum(n_per_et), p_estimated_per_pop=n_estimated_per_pop/sum(n_estimated_per_pop))  -> enterotypes_by_UNSD_regions

write.csv2(enterotypes_by_UNSD_regions, file="enterotypes_by_UNSD_regions.csv")
  
  
enterotypes_by_UNSD_regions %>%
  group_by(Enterotypes_id) %>%
  summarise_if(is.numeric,sum) %>%
  select(Enterotypes_id,p_per_et,p_estimated_per_pop) %>%
  reshape2::melt(id.vars="Enterotypes_id") %>%
  ggplot() + geom_bar(aes(x=Enterotypes_id,y=value,fill=variable), stat="identity", position = "dodge")