curated metagenomics data V3 enterotying

library(purrr)
library(dplyr)
library(magrittr)
library(DirichletMultinomial)

load("shannon_curated_v3.rda")

## group all bootstrap within one file
fit_genus_bootstrap_curated = vector("list",2)
load("fit_genus_list_1_curated_v3.rda")
fit_genus_bootstrap_curated[[1]] = fit_genus_list
load("fit_genus_list_2_curated_v3.rda")
fit_genus_bootstrap_curated[[2]] = fit_genus_list
load("fit_genus_list_3_curated_v3.rda")
fit_genus_bootstrap_curated[[3]] = fit_genus_list
load("fit_genus_list_4_curated_v3.rda")
fit_genus_bootstrap_curated[[4]] = fit_genus_list
load("fit_genus_list_5_curated_v3.rda")
fit_genus_bootstrap_curated[[5]] = fit_genus_list

extract BIC and Laplace coefficients from boostrap models

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


BIC_df = NULL

for(i in 1:5) {
  for(j in 1:5){
    
    BIC_df=cbind(BIC_df,BIC_list[[i]][[j]])
    
  }
}



lplc_df = NULL

for(i in 1:5) {
  for(j in 1:5){
    
    lplc_df=cbind(lplc_df,lplc_list[[i]][[j]])
    
  }
}

plot BIC and laplace results

nb_seed=5

cowplot::plot_grid(

BIC_df %>% as.data.frame %>% 
  mutate(`Dirichlet component`=k) %>%
  reshape2::melt(id.vars="Dirichlet component") %>% 
  mutate(seed=rep(sort(rep(1:nb_seed, length(k))),5)) %>%
  mutate(boot=sort(rep(paste0("boot_",1:5), length(k)*nb_seed))) %>%
  ggplot() + geom_line(aes(x=`Dirichlet component`,y=value,group=variable,color=boot)) +
  ggtitle("BIC optimun") + #xlab("Dirichlet component") +
  theme_bw() + scale_x_continuous(limits = c(1,101), breaks = scales::pretty_breaks(n = 20)),

lplc_df %>% as.data.frame %>% 
  mutate(`Dirichlet component`=k) %>%
  reshape2::melt(id.vars="Dirichlet component") %>% 
  mutate(seed=rep(sort(rep(1:nb_seed, length(k))),5)) %>%
  mutate(boot=sort(rep(paste0("boot_",1:5), length(k)*nb_seed))) %>%
  ggplot() + geom_line(aes(x=`Dirichlet component`,y=value,group=variable,color=boot)) +
  ggtitle("Laplace optimun") + #xlab("Dirichlet component") +
  theme_bw() + scale_x_continuous(limits = c(1,101), breaks = scales::pretty_breaks(n = 20)),



BIC_df %>% as.data.frame %>% 
  mutate(`Dirichlet component`=k) %>%
  reshape2::melt(id.vars="Dirichlet component") %>% 
  mutate(seed=rep(sort(rep(1:nb_seed, length(k))),5)) %>%
  mutate(boot=sort(rep(paste0("boot_",1:5), length(k)*nb_seed))) %>%
  group_by(boot,seed) %>%
  filter(value==min(value)) %>%
  ggplot() + geom_bar(aes(x=`Dirichlet component`)) + 
  ggtitle("BIC optimun") + 
  xlab("Dirichlet components")+
  theme_bw() + scale_x_continuous(limits = c(1,30), breaks = scales::pretty_breaks(n = 30)),


lplc_df %>% as.data.frame %>% 
  mutate(`Dirichlet component`=k) %>%
  reshape2::melt(id.vars="Dirichlet component") %>% 
  mutate(seed=rep(sort(rep(1:nb_seed, length(k))),5)) %>%
  mutate(boot=sort(rep(paste0("boot_",1:5), length(k)*nb_seed))) %>%
  group_by(boot,seed) %>%
  filter(value==min(value)) %>%
  ggplot() + geom_bar(aes(x=`Dirichlet component`)) + 
  ggtitle("Laplace optimun") + 
  xlim(1,30) +  xlab("Dirichlet components")+
  theme_bw() + scale_x_continuous(limits = c(1,30), breaks = scales::pretty_breaks(n = 30)),
labels = "AUTO"
)

Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing
scale.

ggsave(filename = "figures/figS3.pdf")

Saving 12 x 6 in image

BIC and Laplace have a global optimun among differents seeds and bootstraps. All bootstraps found at least one model optimal fit at k=24. Majority of best model based on Laplace were k=24. For BIC, k=24 was also found as an optimum. For downstream, analysis we then choose k=24.

boot=1
seed=1
fit_genus_list = fit_genus_bootstrap_curated[[boot]]
best_genus_lplc = 24 
best_genus_lplc_idx = which(sapply(fit_genus_list[[seed]], function(x) attr(x,"mixture")$Weight %>% length)==best_genus_lplc) # 19th items from list corresponded to k=24 here


enterotypes_curated =

  fit_genus_list[[seed]][[best_genus_lplc_idx]] %>%
  mixture(assign=TRUE) %>% as.data.frame %>% set_colnames(c("Enterotypes_id"))

###################################################
### 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))))

       m1        m2        m3        m4        m5        m6        m7        m8        m9       m10 
0.5150830 0.6794278 0.7418365 1.1821217 1.1740381 0.5320536 0.4839814 0.5281410 0.3535401 0.5849210 
      m11       m12       m13       m14       m15       m16       m17       m18       m19       m20 
0.6168881 1.0230422 1.1717559 0.7744601 0.5246573 1.5678386 0.7076409 0.4040203 1.6276143 1.0186103 
      m21       m22       m23       m24 
1.4157758 1.0864123 0.9955640 1.2710302 

sum(meandiff)

[1] 20.98045

###################################################
### 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_curated_v3 = 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:24)[m_order]) )  -> genus_alpha_weight_curated

save(genus_alpha_weight_curated, file="genus_alpha_weight_curated.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"))

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:24)[m_order]) ) %>% #pull(variable)
  #group_by(taxa) %>%
  #mutate(total = sum(value)) %>%
  #mutate(value=value/total) %>%
  ggplot() + 
  geom_tile(aes(x=variable,y=taxa,fill=log10(value))) + 
  scale_fill_viridis_c("DMM model alpha") +
  ylab("taxa") +
  xlab("DMM components")

NA
NA

correlation between DMM alpha_weight and shannon alpha diversity

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:24)[m_order]) ) %>%
  
  merge(
    
    shannon_curated_v3 %>%
  merge(enterotypes_curated, by="row.names") %>%
  mutate(Enterotypes_id=paste0("m",Enterotypes_id)) %>%
  group_by(Enterotypes_id) %>%
  summarise(m=median(shannon)),
  
  by.x="variable", by.y="Enterotypes_id"
    
    
    
  ) %>%
  group_by(taxa) %>%
  do(broom::tidy(cor.test(.$value, .$m, method="spearman"))) %>%
  ungroup() %>%
  mutate(p.value = p.adjust(p.value, method = "fdr") %>% round(4)) %>%
  arrange(desc(abs(estimate))) %>%
  select(taxa,estimate,p.value)-> cor_alpha_table
  
write.csv2(cor_alpha_table, file="cor_alpha_table.csv")



# df_curated_v3 %>%
#   select(Mean,diff) %>%
#   merge(cor_alpha_table,.,by.x="taxa", by.y="row.names") %>%
#   ggplot(aes(x=abs(estimate),y=diff)) + geom_point()
#   #ggrepel::geom_label_repel(aes(label=taxa), max.overlaps = 20)

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") + 
  xlab("DMM partitions") + ylab("theta extracted from DMM modelling") +
  scale_x_continuous(breaks = 1:24) +
  cowplot::theme_cowplot() -> figS4A
  
figS4A  

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,8,11,30),] %>%
  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)) + 
  geom_rect() +      
  coord_polar(theta="y") + 
  facet_wrap(~variable) + 
  xlim(c(2, 4)) + theme_void() +
  scale_fill_brewer("",type="qual")

ggsave("curated_v3_microbiota_types.pdf")

Saving 15 x 8 in image

classify data with the choosen model

load("curated_v3_genus_count.rda")


predict(fit_genus_list[[1]][[best_genus_lplc_idx]], t(curated_v3_genus_count), 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]][[best_genus_lplc_idx]],"group")), "train", "predicted")) -> enterotypes_prediction

table(enterotypes_prediction$value > 0.50) %>% prop.table()

     FALSE       TRUE 
0.03160195 0.96839805 

table(enterotypes_prediction$value > 0.80) %>% prop.table()

    FALSE      TRUE 
0.2524794 0.7475206 

enterotypes_prediction %>%
  ggplot() + geom_histogram(aes(x=value), bins = 50) + 
  #facet_wrap(~set) + 
  scale_x_continuous(limits = c(0,1.05), labels = scales::percent) +
  xlab("classification confidence") +
  ylab("number of samples") +
  cowplot::theme_cowplot() -> figS4B

figS4B

NA
NA

cowplot::plot_grid(figS4A,figS4B, labels = "AUTO", rel_widths = c(1.7,1))

Avis : Removed 2 rows containing missing values (geom_bar).

ggsave("figures/figS4.pdf", h=5, w=10)

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

combine CMD and AGP DMM weights

permanova en AGP/curated dmm partitions

anova(vegan::betadisper(microbiota_types_jsd, microbiota_annotation$dataset ))


vegan::adonis(microbiota_types_jsd~microbiota_annotation$dataset)