2023-10-17.Rmd

---
title: "Untitled"
author: "kim soyeon"
date: "2023-10-17"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


https://www.microbiomestat.wiki/setting-up-microbiomestat-installation-and-data-preparation/kick-start-your-journey-the-microbiomestat-installation-guide
```{r }
library(phyloseq)
library(tidyverse)
library(devtools)
```


```{r}
devtools::install_github("cafferychen777/MicrobiomeStat")


# install.packages("MicrobiomeStat")
# 
# packages_to_install <- c(
#   "rlang",
#   "tibble",
#   "ggplot2",
#   "matrixStats",
#   "lmerTest",
#   "foreach",
#   "modeest",
#   "dplyr",
#   "pheatmap",
#   "tidyr",
#   "ggh4x",
#   "GUniFrac",
#   "stringr",
#   "rmarkdown",
#   "knitr",
#   "pander",
#   "tinytex",
#   "vegan",
#   "scales",
#    "ape",
#   "ggrepel",       # Used only in the linda.plot function
#   "parallel",      # Utilized only when setting parallel in linda
#   "ggprism",       # Active only when theme.choice is set to “prism”
#   "aplot",         # Operates exclusively in generate_beta_ordination_single
#   "philentropy",   # Used when dist.name is set to "JS"
#   "forcats",       # Comes into play only within barplot and areaplot
#   "yaml",          # Employed solely in mStat_import_qiime2_as_data_obj
#   "biomformat",    # Used in mStat_import_biom_as_data_obj
#   "Biostrings"     # Required for mStat_import_qiime2_as_data_obj and mStat_import_dada2_as_data_obj
# )
# 
# # Installing packages
# install.packages(packages_to_install)

```

## import phyloseq to mStat
```{r}
# Required package
library(microbiome)
library(MicrobiomeStat)
# Load the dataset
data(peerj32)
peerj32.phy <- peerj32$phyloseq

# Convert the Phyloseq object to a MicrobiomeStat data object
data.obj <- mStat_convert_phyloseq_to_data_obj(peerj32.phy)

data.obj$meta.dat$sex <- factor(as.factor(data.obj$meta.dat$sex), levels = c("male", "female"))

head(data.obj$meta.dat)
head(data.obj$feature.tab)
head(data.obj$feature.ann)




## alpha 
alpha <- generate_alpha_test_single(
  data.obj = data.obj,
  alpha.obj = NULL, 
  alpha.name = c("shannon", "observed_species"),
  depth = NULL,
  time.var = "time",
  t.level = "2",
  group.var = "group",
  adj.vars = "sex")
alpha

alpha.p <- generate_alpha_boxplot_single(
  data.obj = data.obj,
  alpha.obj = NULL,
  alpha.name = c("shannon"),
  depth = NULL,
  subject.var = "subject",
  time.var = NULL,
  t.level = NULL,
  group.var = "group",
  strata.var = NULL,
  base.size = 16,
  theme.choice = "bw",
  palette = NULL,
  pdf = TRUE,
  file.ann = NULL,
  pdf.wid = 11,
  pdf.hei = 8.5)

alpha.p

alpha.p2 <- generate_alpha_boxplot_single(
  data.obj = data.obj,
  alpha.obj = NULL,
  alpha.name = c("shannon"),
  depth = NULL,
  subject.var = "subject",
  time.var = "time",
  t.level = "2",
  group.var = "group",
  strata.var = "sex",
  base.size = 16,
  theme.choice = "bw",
  palette = NULL,
  pdf = TRUE,
  file.ann = NULL,
  pdf.wid = 11,
  pdf.hei = 8.5)
alpha.p2


## beta 
beta.df <- generate_beta_test_single(
  data.obj = data.obj,
  dist.obj = NULL,
  time.var = "time",
  t.level = "2",
  group.var = "group", 
  adj.vars = "sex",
  dist.name = c('BC', 'Jaccard') 
)
beta.df$p.tab
beta.df$aov.tab



library(aplot)
generate_beta_ordination_single(
  data.obj = data.obj,
  dist.obj = NULL,
  pc.obj = NULL,
  subject.var = "subject",
  time.var = "time",
  t.level = "2",
  group.var = "group",
  strata.var = "sex",
  adj.vars = "sex",
  dist.name = c("BC"),
  base.size = 20,
  theme.choice = "bw",
  custom.theme = NULL,
  palette = NULL,
  pdf = TRUE,
  file.ann = NULL,
  pdf.wid = 11,
  pdf.hei = 8.5
)

p <- generate_beta_ordination_single(
  data.obj = peerj32.obj,
  dist.obj = NULL,
  pc.obj = NULL,
  subject.var = "subject",
  time.var = "time",
  t.level = "2",
  group.var = "group",
  strata.var = "sex",
  adj.vars = "sex",
  dist.name = c("BC"),
  base.size = 20,
  theme.choice = "bw",
  custom.theme = NULL,
  palette = NULL,
  pdf = TRUE,
  file.ann = NULL,
  pdf.wid = 11,
  pdf.hei = 8.5
)

# Add a p-value annotation to the plot. The p-value is calculated from a PERMANOVA test.
p <- p + annotate("text", x = 0.3, y = 0.8, 
  label = paste("italic(p) == ", 
  format(pvalue, digits = 2)), 
  parse = TRUE, 
  size = 5)

p

```





```{r}
library(vegan)
library(glue)
library(ggtext)
library(aplot)
library(ggplotify)
library(ggpubr)
library(patchwork)


ps <- readRDS("./ps.rds")
ps.rel  <- phyloseq::transform_sample_counts(ps, function(x) x/sum(x))
meta <- sample_data(ps.rel)
meta$SampleID <- rownames(meta)
sample_data(ps.rel) <- meta

## 
phyloseq = ps.rel
index = "bray"
seed = 42
type = "body.site"
plot = "PCoA"
SampleID = "SampleID"
shap = NULL
type_col = col3 = c("#E31A1C", "#1F78B4", "#4D4D4D", "green")
## 
set.seed(seed)
x.dist <- phyloseq::distance(phyloseq, method = index)

myfunc <- function(v1) {
deparse(substitute(v1))
}
dist <- myfunc(x.dist)
a <- adonis2(as.formula (glue("{dist} ~ {type}")), data=data.frame(sample_data(phyloseq)),
             permutations=9999, method=index)
Perm.p <- a$`Pr(>F)`[1]

set.seed(seed)
ord <- ordinate(phyloseq, plot, index)


mat <- ord$vectors[,1:2] %>% as.data.frame()
mat[,SampleID] <- rownames(mat)
mat<- mat %>% dplyr::arrange_(SampleID)
meta <- sample_data(phyloseq) %>% data.frame() %>% dplyr::arrange_(SampleID)
pcoa_df <- merge(meta, mat, by = SampleID)
PC1 <- round(ord$values["Relative_eig"][1,]*100, 1)
PC2 <- round(ord$values["Relative_eig"][2,]*100, 1)

Title <- switch(index,bray="Bray-curtis",
                jaccard="Jaccard",
                unifrac="Unweighted UniFrac",
                wunifrac = "Weighted UniFrac")
##
main.plot <- pcoa_df %>% 
  ggplot2::ggplot(aes(x = Axis.1, y=Axis.2)) + 
  ggplot2::geom_vline(xintercept = 0, colour = "grey80") + 
  ggplot2::geom_hline(yintercept = 0, colour = "grey80") + 
  ggplot2::geom_point(aes_string(shape = shap, color=type),alpha = 0.7, size=2.5) +
  # ggplot2::geom_rug(aes_string(fill=type),, col=rgb(0,0,0,alpha=.2)) + 
  lims(x = c(-0.8, 0.8), y = c(-0.8, 0.8)) + 
  ggplot2::stat_ellipse(aes_string(color= type) ) + 
  ggplot2::theme_test() +
  ggplot2::scale_color_manual(values = type_col) +
  ggplot2::labs(
    # title = Title,
    y =paste0("PCoA2 (", PC2, "%)"),
    x =paste0("PCoA1 (", PC1, "%)")) +
  ggplot2::theme(plot.caption = element_text(hjust = 0)) +
  ggplot2::theme(plot.caption = element_markdown(), 
                 aspect.ratio=1,
                 legend.position = "bottom", 
                 legend.title = element_blank())+
  ggtext::geom_richtext(label.color = NA, size = 4, fill = NA, 
                         hjust = 0, vjust =0, x = -Inf, y = -Inf, 
                         label= paste0( "**PERMANOVA** *p*-value=", Perm.p ) )
main.plot

## 
pcoa1.boxplot <-
  ggplot2::ggplot(pcoa_df) +
  ggplot2::geom_boxplot(ggplot2::aes(x = !!sym(type),y = Axis.1, fill = !!sym(type)), color = "black", alpha = 0.5, show.legend = F) +
  ggplot2::scale_fill_manual(values = col3) +coord_flip() + 
  lims(y = c(-0.8, 0.8)) +  
  theme_void() 

pcoa1.boxplot

## 
pcoa2.boxplot <-
  ggplot2::ggplot(pcoa_df) +
  ggplot2::geom_boxplot(ggplot2::aes(x = !!sym(type),y = Axis.2, fill = !!sym(type)), color = "black", alpha = 0.5, show.legend = F) +
  ggplot2::scale_fill_manual(values = col3) +
  lims(y = c(-0.8, 0.8))  + 
   theme_void() 
pcoa2.boxplot

## patchwork
pcoa1.boxplot + plot_spacer() + main.plot + pcoa2.boxplot + 
  plot_layout(ncol = 2, nrow = 2, widths = c(4, 0.5), heights = c(0.5, 4))

ggsave("./patchwork.png", device = png)
```