MI_MSP_exploration.Rmd

---
title: "MI MSP exploration"
output: html_notebook
---



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

library(dplyr)

```



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

load("data-raw/MilieuInterieur/MilieuInterieur_df_long.rda")


```

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

df_long %>% head()

```

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

microbiome_resources = "/lustre/workgroups/microbiome_resources/reference/IGC/annotation/"

msp_tax = readr::read_tsv(paste0(microbiome_resources,"1661_msps.gtdb_r95_taxonomy.tsv"))

msp_genes = readr::read_tsv(paste0(microbiome_resources,"msp.tsv"))

msp_genes = 
  msp_genes %>% 
  tidyr::separate(msp_name_module_name, into = c("msp","id","module"), extra = "merge", sep = "_") %>%
  mutate(msp = paste0(msp,"_",id)) %>%
  select(-id)


igc_eggnog = readr::read_tsv("/lustre/workgroups/microbiome_resources/reference/IGC/annotation_raw/eggNOG_v5.0/IGC.eggNOG_v5.0.tsv.gz",
                             skip = 2, col_names = TRUE)


eggnog_diff_healthy = readr::read_csv2("data-raw/association_eggnogs_metadata_grp_healthy.csv") %>% filter(p.value < 0.1) %>%
  filter(prevalence_unhealthy > prevalence_healthy)

head(eggnog_diff_healthy)

genes_richness = readr::read_tsv("data-raw/MilieuInterieur/MilieuInterieur_genes_richness.tsv")


```

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

msp_tax %>%
  filter(gtdb_classification %in% grep("bifidum", msp_tax$gtdb_classification, value = TRUE)) %>%
  merge(msp_genes, by.x="msp_name", by.y ="msp") %>% 
  merge(igc_eggnog, by.x="gene_name", by.y="query") %>%
  merge(eggnog_diff_healthy, by.x="eggNOG_OGs", by.y="eggnog") -> msp_healthy_associated_genes
  
# msp_tax %>%
#   filter(gtdb_classification %in% grep("Bifidobacterium", msp_tax$gtdb_classification, value = TRUE))


msp_healthy_associated_genes %>%
  filter(prevalence_unhealthy>prevalence_healthy) %>%
  filter(module != "core")


```

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

head(msp_healthy_associated_genes)

```

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

#fake_sample = data.frame(X1="fake", count=0, sample = df_long %>% pull(sample) %>% unique) # to keep all samples

df_healthy_associated_genes_cast <- df_long %>% 
  filter(X1 %in% msp_healthy_associated_genes$gene_name) %>%
  rbind(fake_sample) %>%
  reshape2::dcast(X1~sample, value.var = "count", fill=0) %>%
  filter(X1 != "fake") -> df_healthy_associated_genes_cast


head(df_long_filtered)

df_healthy_associated_genes_cast <- df_long_filtered %>%
  reshape2::dcast(X1~sample, value.var = "count", fill=0) %>%
  filter(X1 != "fake")


```

``````{r,message=F,warning=F}
df_long %>% 
  filter(X1 %in% msp_healthy_associated_genes$gene_name) %>%
  merge(genes_richness, by="sample") -> bifidum_genes_selected


bifidum_genes_selected %>%
  mutate(count = count/read_mapped_counts) %>%
  select(sample,X1,count) %>%
  rbind(fake_sample) %>%
  reshape2::dcast(X1~sample, value.var = "count", fill=0) %>%
  filter(X1 != "fake") -> df_healthy_associated_genes_cast_prop


head(df_healthy_associated_genes_cast)

write.csv2(df_healthy_associated_genes_cast, "df_healthy_msp.csv", quote = F, col.names = TRUE, row.names = TRUE)

```

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


df_healthy_associated_genes_cast_prop %>%
  tibble::column_to_rownames("X1") %>%
  as.matrix %>%
  heatmap()


```



```{r}

readr::read_tsv("data-raw/MilieuInterieur/MilieuInterieur_df_motus_long.tsv") %>%
  group_by(sample) %>%
  mutate(depth = sum(count)) %>%
  mutate(count=count/depth) %>%
  reshape2::dcast(`#consensus_taxonomy`~sample, value.var="count", fill=0) -> motus_prop_cast
  
motus_bifidum_name = "Bifidobacterium bifidum [ref_mOTU_v25_03116]"

motus_prop_cast %>%
  filter(`#consensus_taxonomy` == motus_bifidum_name) %>%
  tibble::column_to_rownames("#consensus_taxonomy") %>% 
  t %>%
  merge(df_healthy_associated_genes_cast_prop %>%
  tibble::column_to_rownames("X1") %>% t, by="row.names") %>%
  filter(`Bifidobacterium bifidum [ref_mOTU_v25_03116]` > 10^-4) %>%
  mutate_if(is.numeric, function(x) {ifelse(x>0,1,0)} ) %>%
  select(-2) %>%
  dplyr::rename(sample=Row.names) %>%
  reshape2::melt(id.vars="sample") %>%
  merge(genes_richness, by="sample") %>%
  nest(data = -variable) %>% 
  mutate(
    test = map(data, ~ cor.test(.x$value, .x$genes_richness_1M), method="spearman"), # S3 list-col
    tidied = map(test, tidy)
  ) %>% 
  unnest(tidied)
  
motus_prop_cast %>%
  filter(`#consensus_taxonomy` == motus_bifidum_name) %>%
  tibble::column_to_rownames("#consensus_taxonomy") %>% 
  t %>%
  merge(df_healthy_associated_genes_cast_prop %>%
  tibble::column_to_rownames("X1") %>% t, by="row.names") %>%
  filter(`Bifidobacterium bifidum [ref_mOTU_v25_03116]` > 10^-4) %>%
  mutate_if(is.numeric, function(x) {ifelse(x>0,1,0)} ) %>%
  select(-2) %>%
  dplyr::rename(sample=Row.names) %>%
  tibble::column_to_rownames("sample") %>%
  select(
    any_of(msp_healthy_associated_genes %>%
  filter(prevalence_unhealthy>prevalence_healthy) %>%
  filter(module != "core") %>% pull(gene_name) )) %>%
  as.matrix %>%
  heatmap(col = cm.colors(256), 
          scale = "none", 
          distfun = function(x) dist(x,method = 'binary'), 
          hclustfun = function(x) hclust(x,method = 'ward.D'))


motus_prop_cast %>%
  tibble::column_to_rownames("#consensus_taxonomy") %>%
  t %>%
  cor(.,df_healthy_associated_genes_cast_prop %>%
  tibble::column_to_rownames("X1") %>%
  as.matrix %>% t, method="spearman") -> motus_bifidum_cor 


#motus_bifidum_cor 


heatmap(motus_bifidum_cor, col = cm.colors(256))

```

```{r}

motus_bifidum_cor %>%
  reshape2::melt() %>%
  filter(value < -0.2)


```