Heatmaps.Rmd

---
title: "Anlyse différentielle"
author: "ruben"
date: "10/03/2022"
output: html_document
---


## Librairies

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

library(devtools)
library(microbiome)
library(phyloseq)
library(ggpubr)
library(dplyr)
library(SpiecEasi)  
library(data.table)
library(microbiomeutilities)
library(viridis)
library(RColorBrewer)
library(fantaxtic)
library(tibble)
library("pheatmap")

source(file = "functions.R")

```


## c

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

load("curated_v3_otu_tax.rda")

otus <- data.frame(OTU)

taxo <- data.frame(as.matrix(TAX))


```

## Filtre données

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


## on prend les ID en commun

common.ids <- intersect(t(colnames(otus)), sampleMetadata$sample_id)


map <- sampleMetadata


otu <- otus[,common.ids]

tax <- taxo

map <- sampleMetadata %>%
  filter(sample_id %in% common.ids)



# si un patient a plusieurs échantillons, on prend celui qui aura le nombre de reads le plus élevé

map <- unique(setDT(map)[order(sample_id, -number_reads)], by = "sample_id")

nrow(map)

ncol(otu)


# tout est ok !
```
## Regrouper les maladies

```{r,message=F,warning=F}
# On va regrouper les facteurs dans la colonne "disease car il y a beaucoup de facteurs"

map$disease <- factor(map$disease)

# tout d'abord on simplifie les facteurs car il y a en a 23

# tous les facteurs contenant "adenoma" :
Adenoma <- grep("adenoma", levels(map$disease),fixed=TRUE,state.name, value = TRUE)

# tous les facteurs contenant "CRC" :
Colorectal <- c("metastases",grep("CRC", levels(map$disease),fixed=TRUE,state.name, value = TRUE))

# les maladies métaboliques (T2D, hypercholestérolémie etc etc)

Metabolic <- c("ACVD", "hypercholesterolemia", "hypertension", "hypertension;metastases", "IGT", "T2D")

# enfin le bipolar disorder et 

Bowel <- c("IBD")
  
Arthritis <- c("RA")

levels(map$disease) = list(
  "Control" = c("healthy"),
  "Adenoma" =  Adenoma,
  "Colorectal" = Colorectal,
  "Metabolic" = Metabolic,
  "Bowel" = Bowel,
  "Arthritis" = Arthritis,
  "BD" = c('BD'))


```


### Processing

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

# définition du nom des lignes pour map

mapp <- map %>% 
  tibble::column_to_rownames("sample_id") 

# création de l'object phyloseq

## matrices

otu_df <- otu

otu <- as.matrix(otu)


tax <- as.matrix(tax)

##

```

## Objet Phyloseq

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

otu <- otu_table(otu, taxa_are_rows = TRUE)
mapp <- sample_data(mapp)
taxo <- tax_table(tax)

phy <- phyloseq(otu, mapp, taxo)

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

physeq_norm <- prune_samples(sample_sums(physeq_norm) > 0, physeq_norm)
  
physeq_norm  <- transform_sample_counts(physeq_norm, function(x) x / sum(x))

## add the westernized column

physeq_norm@sam_data <- sample_data(transform_westernized(data.frame(physeq_norm@sam_data)))


```


## sous sélections des taxas et samples

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

heatmap_phylo <- function(category, x, order, label_title){
  
  # on remove les samples qui ont que des 0 partout
  
  phy_normalized <- merge_samples_mean(physeq_norm, category)

  #phy_normalized  <- transform_sample_counts(phy_normalized, function(x) log10(x))

   #function(x) log10(x/mean(x))
  
  phy_normalized_ten <- get_top_taxa(phy_normalized, 10, discard_other = TRUE)

  ##
  
  p <- plot_heatmap(phy_normalized_ten,
                    taxa.label = "Species", taxa.order = "Species", 
                    trans=NULL, low="beige", high="red", na.value="beige",sample.order = order, cluster_rows = TRUE)
  

  p$labels$fill <- label_title
  p$labels$x <- x
  
  p
  
  taxa <- plot_taxa_heatmap(phy_normalized_ten,
                            subset.top = 10,
                            transformation = "log10",
                            VariableA = "disease",
                            heatcolors = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(100),
                            )
  
  taxa
  
}

category_age <- c("newborn","child","schoolage","adult","senior")
category_disease <- c("Control","Adenoma","Colorectal","Metabolic", "Bowel", "Arthritis", "BD")
category_gender <- c("male","female")
category_westernized <- c("no","yes")


heatmap_phylo("age_category", "Age Category", category_age, expression(log[10]*'(Relative Mean Abundance)'))

#heatmap_phylo("disease", "Disease", category_disease, expression(log[10]*'(Relative Mean Abundance)'))

#heatmap_phylo("gender", "Gender", category_gender, expression(log[10]*'(Relative Mean Abundance)'))

#heatmap_phylo("non_westernized", "Non Westernized",  category_westernized, expression(log[10]*'(Relative Mean Abundance)'))

#heatmap_phylo("Enterotypes_id", "Enterotype Clusters", as.character(1:24), expression(log[10]*'(Relative Mean Abundance)'))

```


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

library("scales")

category <- "age_category"

heatmap_phylo_prevalence <- function(category, x, order, physeq_norm){
  
  map <- data.frame(physeq_norm@sam_data) %>%
    rownames_to_column(.,"sample_id") %>%
    select(category, sample_id)
  

  tax <- data.frame(physeq_norm@tax_table) %>%
    rownames_to_column(.,"tax") %>%
    select(tax)
    
  #rownames(prevalence_group_mean) <- gsub(".", "|", rownames(prevalence_group_mean))

  prevalence_group_mean <- data.frame(t(physeq_norm@otu_table)) %>%
    rownames_to_column(.,"sample_id") %>%
    merge(.,map, by="sample_id") %>%
    select(-sample_id) %>%
    group_by_at(category) %>%
    summarise_all(., funs(mean(.>0))) %>%
    column_to_rownames(var = category) %>%
    t(.) %>%
    data.frame()
  
  rownames(prevalence_group_mean) <- gsub("\\.", "\\|", rownames(prevalence_group_mean))

  
  prevalence_group_mean <- prevalence_group_mean %>%
    rownames_to_column(.,"tax") %>%
    merge(., tax, by="tax") %>%
    column_to_rownames(.,var='tax')

  
  ## on va regrouper

  
  rownames(map) <- NULL
  prevalence_group_mean <- otu_table(data.frame(prevalence_group_mean), taxa_are_rows = TRUE)
  mapp <- data.frame(category = order, fill = order) %>%
    column_to_rownames(.,var='category')
  

      ## change taxa species names Bifidobacterium -> B.
  taxa_renamed <- data.frame(physeq_norm@tax_table)
  taxa_renamed$Species <- gsub("Bifidobacterium", "B.",taxa_renamed$Species)
  taxa_renamed <- as.matrix(taxa_renamed)
  taxa_renamed <- tax_table(taxa_renamed)

  phy <- phyloseq(prevalence_group_mean, taxa_renamed, sample_data(mapp))
  

  phy <- get_top_taxa(phy, 10, discard_other = TRUE)

  p <- plot_heatmap(phy,
                    taxa.label = "Species", taxa.order = "Species", 
                    low="beige", high="red", na.value="beige",sample.order = order,
                    #trans = identity_trans(),
                    trans = trans_new(name = "test", transform = function(x) x,
                              inverse = function(x) x,
                              domain = c(0, 0.9)),
                    title = x)

  #p$labels$fill <- expression('Prevalence')
  p$labels$fill <- ""
  p$labels$x <- ""
  p$labels$y <- ""
  #p$data$Species <- gsub("Bifidobacterium", "B.",p$data$Species)

  
  p

}



heatmap_phylo_prevalence("age_category", "age category", category_age, physeq_norm = subset_samples(physeq_norm, !is.na(age_category)))


heat_age <- heatmap_phylo_prevalence("age_category", "age category", category_age, physeq_norm = subset_samples(physeq_norm, !is.na(age_category)))
heat_age

heat_disease <- heatmap_phylo_prevalence("disease", "disease", category_disease, physeq_norm = subset_samples(physeq_norm, !is.na(disease)))
heat_disease

heatmap_phylo_prevalence("gender", "Gender", category_gender, physeq_norm = subset_samples(physeq_norm, !is.na(gender)))

heat_western <- heatmap_phylo_prevalence("westernized", "westernized",  category_westernized, physeq_norm = subset_samples(physeq_norm, !is.na(non_westernized)))
heat_western

#heatmap_phylo_prevalence("Enterotypes_id", "Enterotype Clusters", as.character(1:24), physeq_norm = physeq_norm)

#physeq_norm = subset_samples(physeq_norm, !is.na(Enterotypes_id))

```

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

figure <- ggarrange(heat_western, heat_age, heat_disease,
                    labels = c("D","E","F"),
                    ncol =1,
                    nrow = 3,
                    common.legend = TRUE,
                    legend = "right")

annotate_figure(figure, top = text_grob("Bifidobacterium species prevalence", size = 18))
#figure

```