R/net_dis_pcoa.R

In Guan06/mina: Microbial community dIversity and Network Analysis

##############################################################################

#' Visulization of spectra network distance as PCoA.
#'
#' @import ggplot2
#' @import stringr
#' @param x The folder with all egv files generated by net_dis_indi().
#' @return p The plotted figure.
#' @examples
#' \dontrun{
#' data(maize)
#' maize <- norm_tab(maize)
#' maize <- fit_tabs(maize)
#' maize <- get_rep(maize, top = 5)
#' maize <- bs_pm(maize, group = "Compartment", individual = TRUE, out_dir =
#' "./individual_bs_pm/")
#' maize <- net_dis_indi("./individual_bs_pm/", method = "spectra", egv = TRUE,
#' dir = "./egv_folder/")
#' p <- net_dis_pcoa("./egv_folder/")
#' }
#' @rdname net_dis_pcoa
#' @export

setMethod("net_dis_pcoa", signature("character"),
          function(x){

              bs_files <- list.files(x, pattern = "^spectra_bs_",
                                   full.names = TRUE)
              
              ## for bootstrap results and plot
              bs_vectors <- c()
              for (bs in bs_files) {
                  this <- readRDS(bs)
                  bs_vectors <- rbind(bs_vectors, this)
              }
              dis_bs <- as.matrix(parDist(bs_vectors, method = "euclidean"))

              rc_n <- paste0(rownames(dis_bs), "_", seq(1 : nrow(dis_bs)))
              rownames(dis_bs) <- colnames(dis_bs) <- rc_n
              ## get the group information
              idx <- seq(1,length(colnames(dis_bs)) * 2, by = 2)
              design <- data.frame(Net = colnames(dis_bs),
                  Group = unlist(strsplit(colnames(dis_bs), "_bs"))[idx])

              dis_bs_dmr <- cmdscale(dis_bs, k = 4, eig = T)
              eig <- dis_bs_dmr$eig
              eig[eig < 0] <- 0
              eig1 <- 100 * eig[1] / sum(eig)
              eig2 <- 100 * eig[2] / sum(eig)
              
              points <- as.data.frame(dis_bs_dmr$points[, 1:2])
              colnames(points) <- c("x", "y")
              points$Net <- rownames(points)
              points <- merge(points, design)
              
              p <- ggplot(points, aes(x, y, color = Group)) +
                  geom_point(alpha = 0.8) +
                  theme(panel.background = element_blank(),
                    panel.grid = element_blank(),
                    axis.line.x = element_line(color = "black"),
                    axis.line.y = element_line(color = "black"),
                    axis.ticks = element_line(color = "black"),
                    axis.text = element_text(colour = "black", size = 10),
                    text = element_text(family="sans"),
                    legend.position = "top") +
                  coord_fixed(ratio = 1) +
                  labs (x = paste0("PCo1 (", format(eig1, digits = 4), "%)"),
                        y = paste0("PCo2 (", format(eig2, digits = 4), "%)"))
          }
)