bak/pcoa_result_plot.R

In KidultXJT/Kitult: Kitult

pcoa_result_plot <- function(
  input,
  output,
  outputFormat = "png",
  main = NULL,
  group = NULL,
  group_colour = NULL,
  select = NULL,
  point_size  = 8,
  point_alpha = 1,
  point_text_size = 4,
  width  = 9,
  height = 8
){
  # library(ggplot2)
  # library(ape)
  # library(xlsx)
  # library(readxl)
  
  #input = "~/test/04.Beta/total/unweighted_unifrac_pcoa.txt"
  infile <- read.table(input, sep = "\n",check.names = F)
  sample = as.character(infile[c(-(1:6),-nrow(infile),-(nrow(infile)-1)),])
  sample = sapply(1:length(sample),function(i) strsplit(sample, "\t")[[i]])[1,]
  
  #######################
  ##      Eigv.m        ##
  ########################
  Eigv.m = matrix(as.numeric(strsplit(as.character(infile$V1[2]), "\t")[[1]]),
                  nrow = 1)
  colnames(Eigv.m) = c(sample)
  
  ########################
  ## Proportion Explain ##
  ########################
  Explain.m = matrix(as.numeric(strsplit(as.character(infile$V1[4]), "\t")[[1]]),
                     nrow = 1)
  colnames(Explain.m) = c(sample)
  
  
  xlab = paste("PC1(",as.numeric(sprintf("%.3f",Explain.m[1,1]))*100,"%)", sep = "")
  ylab = paste("PC2(",as.numeric(sprintf("%.3f",Explain.m[1,2]))*100,"%)", sep = "")
  
  
  ########################
  ##     Sites.df       ##
  ########################
  sites.df = matrix()
  for(i in 7:(6+length(sample))){
    sites.df = data.frame(sites.df, as.character(infile$V1[i]))
  }
  
  sites.df <- as.matrix(t(sites.df)[2:dim(sites.df)[2]])
  sites.df = as.list(sites.df)
  
  sites.df = data.frame(
    do.call('rbind',
            strsplit(as.character(sites.df),'\t',fixed=TRUE)
    ))
  
  sites.rowname = sites.df$X1#c(samples[[1]])
  sites.df      = sites.df[,2:(length(sample))]
  rownames(sites.df) = sites.rowname
  
  PC1 <- as.numeric(as.character(sites.df$X2))
  PC2 <- as.numeric(as.character(sites.df$X3))
  sites.df2 <- data.frame(PC1,PC2)
  
  if(!is.null(group)) {groups <- unlist(strsplit(group,","))} else groups = group
  sites.df2$Groups = groups
  sites.df2$Name  = c(sample)
  
  ########################
  ####     Plot       ####
  ########################
  ggpcoa = ggplot(sites.df2,aes(PC1,PC2))
  
  if(!is.null(main)) ggpcoa <- ggpcoa + labs(title = main)
  
  ggpcoa = ggpcoa +
    xlab(xlab)
  ggpcoa = ggpcoa +
    ylab(ylab)
  
  ggpcoa = ggpcoa +
    geom_hline(yintercept=0,linetype=2,color="grey",size=.5) + 
    geom_vline(xintercept=0,linetype=2,color="grey",size=.5)
  
  if(!is.null(group)) {ggpcoa <- ggpcoa +
    geom_point(data  = sites.df2,
               size  =  point_size,
               alpha =  point_alpha, 
               aes(color=Groups)) + 
    geom_text(data  = sites.df2,
              size  = point_text_size,
              aes(label = sample),
              color = "grey56")}else{
              # nudge_x = .02, 
              # nudge_y = .02
                
                ggpcoa <- ggpcoa +
                  geom_point(data  = sites.df2,
                             # colour = colour,
                             size  = point_size,
                             alpha = point_alpha) + 
                  geom_text(data  = sites.df2,
                            size  = point_text_size,
                            aes(label = sample),
                            color = "grey56")
                            # nudge_x = .02, 
                            # nudge_y = .02)
              }
  
  
  colours <- c("#F4BD6D","#87CCC5","olivedrab3","brown2","goldenrod1","deepskyblue3","maroon2","mediumpurple4",
               "cadetblue4","darkblue","hotpink3","indianred1","gray59",
               "darkorange1","slategray3","grey95","mediumorchid3","paleturquoise1",
               "orangered4","grey64","gray62","royalblue4","antiquewhite3")
  cols <- colours[1:dim(as.data.frame(table(sites.df2$Groups)))[1]]
  ggpcoa <- ggpcoa + 
    scale_colour_manual(values = cols)
  
  if(!is.null(group_colour)) {col <- unlist(strsplit(group_colour,",")); ggpcoa = ggpcoa + scale_colour_manual(values = col)}
  
  ggpcoa <- ggpcoa + theme_bw()
  
  ggpcoa <- ggpcoa + theme(
    legend.title = element_text(angle = 0,  
                                hjust =.95, 
                                vjust =.95,
                                size  = 12),
    legend.position   = "right",
    legend.background = element_rect(
      colour = "white", 
      fill = "white", 
      size = 1, 
      linetype = 1),#'dashed'
    legend.key.size   = unit(.5, "cm"),
    legend.key.height = unit(.5, "cm"),
    legend.key.width  = unit(.5, "cm"),
    legend.text       = element_text(angle = 0,  
                                     hjust =.95, 
                                     vjust =.95,
                                     size  = 10),
    legend.justification = "top",
    panel.background = element_rect(fill = "white",colour = "white",size = 1.5))
  # legend.justification = "top")
  
  ##########################################################################
  
  # ggpcoasp <- ggpcoa +
  #   geom_point(data   = sp_df,
  #              color  ="royalblue4", 
  #              shape  = 0,
  #              alpha  = .6,
  #              size   = 4)
  # 
  # ggpcoaspt <- ggpcoa +
  #   geom_text(data    = sp_df,
  #             aes(label = name),
  #             color   = sp_colour, 
  #             alpha   = .7,
  #             size    = sp_size)
  
  Format   = strsplit(outputFormat,",")[[1]]
  sapply(Format,function(x){
    if(x == "pdf") temp = paste0("cairo_pdf(file = paste0(output,\'.\',x), bg = \'transparent\', width = width, height = height)")
    else temp = paste0(x,"(file = paste0(output,\'.\',x), bg = \'transparent\', width = width*72, height = height*72, type = \'cairo\')")
    eval(parse(text = temp))
    # png(file = paste0(output_hc,".",x), bg = "transparent", width = width, height = height)
    plot(ggpcoa)
    dev.off()
    
  })
  return(ggpcoa)
  # if(!is.null(output)) ggsave(output,
  #                             ggpcoa,
  #                             width  = 9,
  #                             height = 8)
  
}