HumannANCOMBC.qmd

---
title: "Tripobiome Humman ANCON-BC"
format: html
editor: visual
---

```{r}
library(tidyverse)
library(vegan)
library(ggstatsplot)
library(BiodiversityR)
library(mixOmics)
library(ggsci)
library(ANCOMBC)
```

```{r message=FALSE, warning=FALSE, error=FALSE}

pathabundance_files <- dir("/storage2/tripobiome/original/", recursive = T, pattern = "pathabundance",full.names = T )
pathcoverage_files <- dir("/storage2/tripobiome/original/", recursive = T, pattern = "pathcoverage",full.names = T)



metadata <- read_csv("/storage2/tripobiome/Metadata.csv")
metadata <- metadata %>% dplyr::select(Sample,tto_previo_no_completo,tto_parasiticida_correcto,GRUPO,sexo,bristol) %>% 
  mutate(GRUPO2 = ifelse(GRUPO == "Control no infectado","Control","Chagas")) %>% 
  replace_na(list(tto_parasiticida_correcto = "No")) %>% 
  mutate(tto_parasiticida_correcto = ifelse(GRUPO == "Control no infectado","Control",tto_parasiticida_correcto))
```

Transform to tables

```{r message=FALSE, warning=FALSE, error=FALSE, results='hide'}


#load("/storage2/tripobiome/HumannANCOMSession.rds")

pathabundance_table <- tibble()

for(i in 1:length(pathabundance_files))
{
  name = basename(pathabundance_files[i])
  tmp_table <- vroom::vroom(pathabundance_files[i])
  tmp_table <- tmp_table %>% dplyr::rename(Pathway = 1, Abundance = 2) %>% mutate(Sample =name)

  pathabundance_table <- bind_rows(pathabundance_table,tmp_table)

}
pathabundance_table <- pathabundance_table %>% mutate(Sample = str_replace(Sample,"_pathabundance.tsv",""))



pathcoverage_table <- tibble()

for(i in 1:length(pathcoverage_files))
{
  name = basename(pathcoverage_files[i])
  tmp_table <- vroom::vroom(pathcoverage_files[i])
  tmp_table <- tmp_table %>% dplyr::rename(Pathway = 1, Coverage = 2) %>% mutate(Sample =name)

  pathcoverage_table <- bind_rows(pathcoverage_table,tmp_table)

}
pathcoverage_table <- pathcoverage_table %>% mutate(Sample = str_replace(Sample,"_pathcoverage.tsv",""))

path_table <- dplyr::full_join(pathabundance_table,pathcoverage_table)

path_table <- path_table %>% separate(Sample,c("Sample","Position","Lane","Read","number"), sep="_") %>% group_by(Pathway,Sample) %>% summarise(Abundance = sum(Abundance), Coverage = max(Coverage))




```

```{r}

path_table %>% 
  group_by(Sample) %>% 
  summarise(Total = sum(Abundance)) %>% 
  ggplot(aes(x = fct_reorder(Sample,Total), y = Total)) + 
  geom_col() + 
  coord_flip() + 
  labs(title = "Total Reads Pathway")

descartes <- c("S3","S27","S2","S7","S75")


```

# Pathway

```{r}
path_table_summ <- path_table %>% 
  separate(Pathway, c("Pathway","Taxonomy"), sep = "\\|") %>% 
  group_by(Pathway,Sample) %>% 
  summarise(Abundance = sum(Abundance), Coverage = sum(Coverage))

path_table_summ <- path_table_summ %>% 
  filter(!(Sample %in% descartes)) %>% 
  group_by(Sample) %>% 
  mutate(TotalAbundance = sum(Abundance)) %>% 
  ungroup()

min_abundance <- round(min(path_table_summ$TotalAbundance))

path_table_summ_matrix <- path_table_summ %>% 
  filter(Coverage > 0.5) %>% 
  group_by(Sample,Pathway) %>% 
  summarise(Abundance = round(sum(Abundance))) %>% 
  pivot_wider(names_from = "Pathway", values_from = "Abundance", values_fill = 0)

```

## pre-processing data

```{r}
library(mia)
library(ANCOMBC)

path_mat <- path_table_summ %>% 
  mutate(Sample = paste0("S",Sample)) %>% 
  group_by(Pathway,Sample) %>% 
  summarise(Abundance = sum(Abundance)) %>% 
  ungroup() %>% 
  filter(!(Sample %in% descartes)) %>% 
  group_by(Pathway) %>% 
  mutate(vari = var(Abundance), N = n()) %>% 
  filter(N > 5) %>% 
  ungroup() %>% 
  dplyr::select(Sample,Pathway,Abundance) %>% 
  pivot_wider(names_from = Sample, values_from = Abundance, values_fill = 0) %>% 
  column_to_rownames("Pathway") %>% as.matrix()


assays = SimpleList(counts = path_mat)

smd <- metadata %>% 
  filter(Sample %in% colnames(path_mat)) %>% 
  arrange(match(Sample,colnames(path_mat))) %>% 
  column_to_rownames("Sample")

smd <- DataFrame(smd)

tse <- TreeSummarizedExperiment(assays = assays,
                                colData = smd)
```

## DA

```{r warning=FALSE, message=FALSE}
out_grupo2 <- ancombc2(data = tse, fix_formula = "GRUPO2")
out_tto <- ancombc2(data = tse, fix_formula = "tto_parasiticida_correcto",)
```

## Results

```{r}

out_grupo2$res %>% dplyr::rename(Pathway = 1) %>% 
  filter(diff_GRUPO2Control ==TRUE) %>% 
    ggplot(aes(x = reorder(Pathway,lfc_GRUPO2Control), y = lfc_GRUPO2Control, fill = lfc_GRUPO2Control)) + 
    geom_bar(stat = "identity", width = 0.7, color = "black", 
             position = position_dodge(width = 0.4)) +
    geom_errorbar(aes(ymin = lfc_GRUPO2Control - se_GRUPO2Control, ymax = lfc_GRUPO2Control + se_GRUPO2Control), 
                  width = 0.2, position = position_dodge(0.05), color = "black") +
    scale_fill_gsea(name= "Log fold Change") +
  labs(x = NULL, y = "Log fold change", 
         title = "Log fold changes as one unit increase of Chagas vs Control") + 
  coord_flip() +
  theme_bw() + 
    theme(plot.title = element_text(hjust = 0.5),
          panel.grid.minor.y = element_blank(),
          axis.text.x = element_text(angle = 60, hjust = 1))


out_tto$res %>% dplyr::rename(Pathway = 1) %>% 
  filter(diff_tto_parasiticida_correctoSi ==TRUE) %>% 
    ggplot(aes(x = reorder(Pathway,lfc_tto_parasiticida_correctoSi), y = lfc_tto_parasiticida_correctoSi, fill = lfc_tto_parasiticida_correctoSi)) + 
    geom_bar(stat = "identity", width = 0.7, color = "black", 
             position = position_dodge(width = 0.4)) +
    geom_errorbar(aes(ymin = lfc_tto_parasiticida_correctoSi - se_tto_parasiticida_correctoSi, ymax = lfc_tto_parasiticida_correctoSi + se_tto_parasiticida_correctoSi), 
                  width = 0.2, position = position_dodge(0.05), color = "black") +
  scale_fill_gsea(name= "Log fold Change") +
  labs(x = NULL, y = "Log fold change", 
         title = "Log fold changes as one unit increase of Tratamiento Parasiticida Correcto") + 
  coord_flip() +
  theme_bw() + 
    theme(plot.title = element_text(hjust = 0.5),
          panel.grid.minor.y = element_blank(),
          axis.text.x = element_text(angle = 60, hjust = 1))

```

## Alpha y Beta basados en la normalización de ANCOMBC

```{r}
bc_counts <- out_grupo2$bias_correct_log_table

bc_counts[is.na(bc_counts)] = 0

bc_counts = (10^bc_counts)
```

```{r}
diversities <- estimateR(round(bc_counts) %>% t(), parallel(40))%>% 
   t() %>% 
  as.data.frame() %>% 
  rownames_to_column("Sample") %>% 
  dplyr::inner_join(metadata)

ggbetweenstats(diversities, x = GRUPO2, y = S.obs, type = "robust")
ggbetweenstats(diversities, x = GRUPO2, y = S.chao1, type = "robust")
ggbetweenstats(diversities, x = GRUPO2, y = se.chao1, type = "robust")
ggbetweenstats(diversities, x = GRUPO2, y = S.ACE, type = "robust")
ggbetweenstats(diversities, x = GRUPO2, y = se.ACE, type = "robust")

ggbetweenstats(diversities, x = tto_parasiticida_correcto, y = S.obs, type = "robust")
ggbetweenstats(diversities, x = tto_parasiticida_correcto, y = S.chao1, type = "robust")
ggbetweenstats(diversities, x = tto_parasiticida_correcto, y = se.chao1, type = "robust")
ggbetweenstats(diversities, x = tto_parasiticida_correcto, y = S.ACE, type = "robust")
ggbetweenstats(diversities, x = tto_parasiticida_correcto, y = se.ACE, type = "robust")
```

```{r}



bc_dist <- vegdist(bc_counts %>% t(), method = "bray")


data_plot_bc <- metaMDS(bc_dist)

data_plot_bc$points %>% 
  as_tibble(rownames ="Sample") %>% 
  dplyr::inner_join(metadata) %>% 
  ggplot(aes(x = MDS1, y = MDS2, label = Sample , fill = GRUPO2)) + geom_label()

data_plot_bc$points %>% 
  as_tibble(rownames ="Sample") %>% 
  dplyr::inner_join(metadata) %>% 
  ggplot(aes(x = MDS1, y = MDS2, label = Sample , fill = tto_parasiticida_correcto)) + geom_label()

D_meta_bc <- bc_dist %>% as.matrix() %>% as_tibble(rownames= "Sample") %>% dplyr::inner_join(metadata)

D_meta_bc_dist = D_meta_bc %>% 
  dplyr::select(Sample:S9) %>% 
  as.data.frame() %>% 
  column_to_rownames("Sample") %>% 
  as.dist()

D_meta_bc_data <- D_meta_bc %>% 
  dplyr::select(Sample,GRUPO2,tto_parasiticida_correcto,GRUPO,sexo,bristol) %>% 
  mutate(GRUPO2 = as.factor(GRUPO2),
         tto_parasiticida_correcto = as.factor(tto_parasiticida_correcto)) %>% column_to_rownames("Sample")

adonis2(D_meta_bc_dist ~ GRUPO2 + tto_parasiticida_correcto, data = D_meta_bc_data)

anosim(D_meta_bc_dist, D_meta_bc_data$GRUPO2) %>% 
  plot(ylab = "Rank Similarities", xlab = "Grupo2")
anosim(D_meta_bc_dist, D_meta_bc_data$tto_parasiticida_correcto)%>% 
  plot(ylab = "Rank Similarities", xlab = "Tratamiento Correcto")


plsda(D_meta_bc_dist,D_meta_bc_data$GRUPO2) %>% plotIndiv(ellipse = TRUE)
plsda(D_meta_bc_dist,D_meta_bc_data$tto_parasiticida_correcto) %>% plotIndiv(ellipse = TRUE)

  
```

# Correlations

# Load data

## OTUS Ruminococcaceae, Faecalibacterium, Prevotella y Bacteroides, Parabacteroides spp, Enterococcus hirae, Prevotella

```{r}
source("/storage/megalodon/R/import_kraken.R")
source("/storage/megalodon/R/kraken_heat_tree.R")

files = dir("/storage2/tripobiome/taxonomia/",pattern = "mic.kreport",full.names = T)

OTUs <- import_kraken(files,threads = 40)
OTUs <- OTUs %>% mutate(Sample = gsub(".mic.kreport","",Sample))

OTU_rare <- OTUs %>%
   filter(grepl("Bacteria",Fullname)) %>%
   filter(grepl("S",TaxRank)) %>%
   dplyr::select(Name2,Reads,Sample) %>% 
   distinct() %>% 
   pivot_wider(names_from = Name2, values_from = Reads, values_fill = 0) %>%
   column_to_rownames("Sample") 

OTU_rare <- rrarefy(OTU_rare, sample = min(rowSums(OTU_rare)))

OTU_rare <- OTU_rare %>% 
  as.data.frame() %>% 
  rownames_to_column("Sample") %>%
  pivot_longer(names_to = "OTU", values_to = "Abundance", -Sample) %>% 
  group_by(OTU) %>% 
  mutate(Var = var(Abundance), Mean = mean(Abundance)) %>% 
  filter(Var > 0.25*Mean, Mean > 10) %>%
  ungroup() %>% 
  dplyr::select(Sample,OTU,Abundance) %>% 
  pivot_wider(names_from = OTU, values_from = Abundance, values_fill = 0) %>%
  column_to_rownames("Sample")
```


# Pathway vs Taxonomy
```{r}
commons <- path_mat %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("Sample") %>% 
  dplyr::select(Sample) %>% 
  semi_join(OTU_rare %>%
              as.data.frame() %>% 
              rownames_to_column("Sample") %>% 
              dplyr::select(Sample))

path_mat_comm <- path_mat %>%
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("Sample") %>% 
  filter(Sample %in% commons$Sample) %>% 
  arrange(Sample) %>% 
  column_to_rownames("Sample") %>% 
  as.matrix()

OTU_matrix <- OTU_rare %>%
  as.data.frame() %>% 
  rownames_to_column("Sample") %>% 
  filter(Sample %in% commons$Sample) %>% 
  arrange(Sample) %>% 
  column_to_rownames("Sample") %>% 
  as.matrix()

OTU_vs_Path<- cor(path_mat_comm,OTU_matrix, method = "spearman")


library(igraph)
library(ggraph)
library(tidygraph)
```


```{r fig.width= 20, fig.height=20}
OTU_vs_Path %>% 
  as.data.frame() %>% 
  rownames_to_column("Pathway") %>% 
  pivot_longer(names_to = "OTU", values_to = "Correlation",-Pathway) %>% 
  filter(Correlation > 0.45 ) %>% 
  rename(from = Pathway, to = OTU) %>% 
  as_tbl_graph(nodes = nodes) %>%
  ggraph(layout = "gem") + 
  geom_node_point()+
  geom_node_text(aes(label = name),repel = TRUE) +
  geom_edge_link() +
  theme_void()
```