Plotting results across H0, HA, H0bGC

Loading

library(ggplot2)
library(ggpubr)

Useful functions

#d<-read.table("Data/Test/youpi.tsv", h=T, sep="\t")
#d1<-read.table("Data/Test/tralala.txt", h=T, sep="\t")
readData <- function(h0ha, h0gBGC) {
  d<-read.table(h0ha, h=T, sep="\t")
  d1<-read.table(h0gBGC, h=T, sep="\t")
  d1$Oracle <- 2
  d<-rbind(d, d1)
  d$Oracle <- as.factor(d$Oracle)
  levels(d$Oracle)[levels(d$Oracle)=="0"] <- "H0"
  levels(d$Oracle)[levels(d$Oracle)=="1"] <- "HA"
  levels(d$Oracle)[levels(d$Oracle)=="2"] <- "H0bGC"
  return (d)
}

plotOrderedFactors <- function(d, method, maxorder, decreasing, name, plot_y) {
  # Get the three colors that I want
  colors_3 <- scale_color_brewer(palette="Dark2")$palette(3)
  ord <- order(d[which(colnames(d)==method)], decreasing=decreasing)
  d2 <- d[ord, ]
  if (plot_y) {
  plot(c(0.5, 2, 3.5),c (1, 50, maxorder), ylim=rev(c(1, maxorder)), type="n", xaxt='n', xlab="", ylab="")
  title(ylab = "Rank", cex.lab = 1.5, line = 1.7)
  }
  else {
  plot(c(0.5, 2, 3.5),c (1, 50, maxorder), ylim=rev(c(1, maxorder)), type="n", xaxt='n', xlab="", ylab="", yaxt='n')
  }
  title(xlab = name, cex.lab = 1.5, line = 0.5)
  colors_n <- sapply(as.numeric(d2$Oracle), function(x) {return (colors_3[x])})
  points(as.numeric(d2$Oracle) + rnorm(n=length(d2$Oracle), mean=0, sd=0.1), 1:length(ord), pch = as.numeric(d2$Oracle), cex=1, col= colors_n, lwd=2)
  return(d2)
} 

plotAllMethods <- function (maxorder, nrows, ncols, withPCAndOC=TRUE) {
  if (nrows > 1) {
    print("ERROR : more than 1 row")
    return()
  }
  if (withPCAndOC) {
  layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8), nrows, ncols, byrow = TRUE), widths=rep(3, 8), heights=rep(0.1, 8))
  }
  else {
      layout(matrix(c(1, 2, 3, 4, 5, 6), nrows, ncols, byrow = TRUE), widths=rep(3, 6), heights=rep(0.1, 6))
  }
  #par(oma=c(0,0,0,0))  # all sides have 3 lines of space  
  par(mar = c(2,3,2,0) + 0.1) ## default is c(5,4,4,2) + 0.1
  temp = plotOrderedFactors(d, "Identical_LG08", maxorder, TRUE, "Identical", plot_y=TRUE)
  temp = plotOrderedFactors(d, "Topological_LG08", maxorder, TRUE, "Topological", plot_y=FALSE)
  temp = plotOrderedFactors(d, "PCOC", maxorder, TRUE, "PCOC", plot_y=FALSE)
  if (withPCAndOC) {
    temp = plotOrderedFactors(d, "PC", maxorder, TRUE, "PC", plot_y=FALSE)
    temp = plotOrderedFactors(d, "OC", maxorder, TRUE, "OC", plot_y=FALSE)
  }
  temp = plotOrderedFactors(d, "Tdg09_1MinusLRT", maxorder, TRUE, "Tdg09", plot_y=FALSE)
  temp = plotOrderedFactors(d, "Mutinomial_1MinusLRT", maxorder, TRUE, "Multinomial", plot_y=FALSE)
  temp = plotOrderedFactors(d, "Diffsel_max", maxorder, TRUE, "Diffsel", plot_y=FALSE)
}

plotViolinPlots <- function (nrows, ncols, withPCAndOC=TRUE) {
    if (nrows > 1) {
    print("ERROR : more than 1 row")
    return()
  }
  par(mar = c(2,3,2,0) + 0.1) ## default is c(5,4,4,2) + 0.1
  p1 <- ggplot(d, aes(x=Oracle, y=PCOC, color=Oracle, fill=Oracle)) + 
    geom_violin(trim=TRUE) + 
    scale_color_brewer(palette="Dark2") + 
    scale_fill_brewer(palette="Dark2")+ 
    rremove("x.text")  + 
    theme(legend.position="none") +
    labs(y = "")
  p2 <- ggplot(d, aes(x=Oracle, y=PC, color=Oracle, fill=Oracle)) + 
    geom_violin(trim=TRUE) + 
    scale_color_brewer(palette="Dark2") + 
    scale_fill_brewer(palette="Dark2")+ rremove("x.text")  + theme(legend.position="none")+ rremove("y.text")  +
    labs(y = "")
  p3 <- ggplot(d, aes(x=Oracle, y=OC, color=Oracle, fill=Oracle)) + 
    geom_violin(trim=TRUE) + 
    scale_color_brewer(palette="Dark2") + 
    scale_fill_brewer(palette="Dark2")+ rremove("x.text")  + theme(legend.position="none")+ rremove("y.text")  +
    labs(y = "")
  p4 <- ggplot(d, aes(x=Oracle, y=Identical_LG08, color=Oracle, fill=Oracle)) + 
    geom_violin(trim=TRUE) + 
    scale_color_brewer(palette="Dark2") + 
    scale_fill_brewer(palette="Dark2")+ rremove("x.text")  + theme(legend.position="none")+ rremove("y.text")  +
    labs(y = "")
  p5 <- ggplot(d, aes(x=Oracle, y=Topological_LG08, color=Oracle, fill=Oracle)) + 
    geom_violin(trim=TRUE) + 
    scale_color_brewer(palette="Dark2") + 
    scale_fill_brewer(palette="Dark2")+ rremove("x.text")  + theme(legend.position="none")+ rremove("y.text")  +
    labs(y = "")
  p6 <- ggplot(d, aes(x=Oracle, y=Tdg09_1MinusLRT, color=Oracle, fill=Oracle)) + 
    geom_violin(trim=TRUE) + 
    scale_color_brewer(palette="Dark2") + 
    scale_fill_brewer(palette="Dark2")+ rremove("x.text")  + theme(legend.position="none")+ rremove("y.text")  +
    labs(y = "")
  p7 <- ggplot(d, aes(x=Oracle, y=Mutinomial_1MinusLRT, color=Oracle, fill=Oracle)) + 
    geom_violin(trim=TRUE) + 
    scale_color_brewer(palette="Dark2") + 
    scale_fill_brewer(palette="Dark2")+ rremove("x.text")  + theme(legend.position="none")+ rremove("y.text")  +
    labs(y = "")
  p8 <- ggplot(d, aes(x=Oracle, y=Diffsel_max, color=Oracle, fill=Oracle)) + 
    geom_violin(trim=TRUE) + 
    scale_color_brewer(palette="Dark2") + 
    scale_fill_brewer(palette="Dark2")+ rremove("x.text")  + theme(legend.position="none")+ rremove("y.text")  +
    labs(y = "")
  if(withPCAndOC) {
    ggarrange(p4, p5, p1, p2, p3, p6, p7, p8 , 
          labels = c("Identical", "Topological", "PCOC", "PC", "OC", "Tdg09", "Multinomial", "Diffsel"),
          ncol = ncols, nrow = nrows)
  }
  else {
    ggarrange(p4, p5, p1, p6, p7, p8 , 
          labels = c("Identical", "Topological", "PCOC", "Tdg09", "Multinomial", "Diffsel"),
          ncol = ncols, nrow = nrows)
  }
  }

By method, count the numbers of sites of each type from rank 1 to the last

plotCumSum <- function (d, method, maxorder, decreasing, name, plot_y, legend=FALSE) {
  # Get the three colors that I want
  colors_3 <- scale_color_brewer(palette="Dark2")$palette(3)
  lwidth <- 3
  ord <- order(d[which(colnames(d)==method)], decreasing=decreasing)
  d2 <- d[ord, ]
  count_h0=ave(d2$Oracle=="H0", "H0", FUN=cumsum)[1:maxorder]
  count_ha=ave(d2$Oracle=="HA", "HA", FUN=cumsum)[1:maxorder]
  count_h0bGC=ave(d2$Oracle=="H0bGC", "H0bGC", FUN=cumsum)[1:maxorder]
  data <- data.frame(H0=count_h0, HA=count_ha, H0bGC=count_h0bGC)
  maximum_y = max(data)
  if (plot_y) {plot(data$H0, t="l", col=colors_3[1], lwd=lwidth, xlim=c(0, maxorder), ylim=c(0, maximum_y), xlab="", ylab="", main=name)
  }
  else {
    plot(data$H0, t="l", col=colors_3[1], lwd=lwidth, xlim=c(0, maxorder), ylim=c(0, maximum_y), xlab="", ylab="", main=name, yaxt='n')
  }
  lines(data$HA, col=colors_3[2], lwd=lwidth)
  lines(data$H0bGC, col=colors_3[3], lwd=lwidth)
  if (legend) {
    legend("topleft", col=colors_3, lwd=lwidth, legend=c("H0", "HA", "H0bGC"))
  }
  return(data)
}

plotCumSums <- function (nrows, ncols, maxorder, withPCAndOC = TRUE) {
  if (nrows > 1) {
    print("ERROR : more than 1 row")
    return()
  }
  if (withPCAndOC) {
    layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8), nrows, ncols, byrow = TRUE), widths=rep(3, 8), heights=rep(0.1, 8))
  }
  else {
    layout(matrix(c(1, 2, 3, 4, 5, 6), nrows, ncols, byrow = TRUE), widths=rep(3, 6), heights=rep(0.1, 6))
  }
  par(mar = c(2,3,2,0) + 0.1) ## default is c(5,4,4,2) + 0.1
  temp = plotCumSum(d, "Identical_LG08", maxorder, TRUE, "Identical", plot_y=TRUE)
  temp = plotCumSum(d, "Topological_LG08", maxorder, TRUE, "Topological", plot_y=FALSE)
  temp = plotCumSum(d, "PCOC", maxorder, TRUE, "PCOC", plot_y=FALSE)
  if (withPCAndOC) {
    temp = plotCumSum(d, "PC", maxorder, TRUE, "PC", plot_y=FALSE)
    temp = plotCumSum(d, "OC", maxorder, TRUE, "OC", plot_y=FALSE)
  }
  temp = plotCumSum(d, "Tdg09_1MinusLRT", maxorder, TRUE, "Tdg09", plot_y=FALSE)
  temp = plotCumSum(d, "Mutinomial_1MinusLRT", maxorder, TRUE, "Multinomial", plot_y=FALSE)
  temp = plotCumSum(d, "Diffsel_max", maxorder, TRUE, "Diffsel", plot_y=FALSE)
}

First, we analyze the results on the Amaranth phylogeny

#d <- readData("Data/Test/youpi.tsv", "Data/Test/tralala.txt")
d <- readData("Data/h0ha_amaranth.tsv", "Data/h0_gbgc_amaranth.tsv")

All sites

plotAllMethods(2000, 1, 6, withPCAndOC = FALSE)

Zooming in

plotAllMethods(150, 1, 8)

plotAllMethods(150, 1, 6, withPCAndOC = FALSE)

Violin plots

plotViolinPlots(1, 8)

plotViolinPlots(1, 6, withPCAndOC = FALSE)

Cumulative plots

tmp <- plotCumSums (1, 8, 150 ) 

tmp <- plotCumSums (1, 6, 150, withPCAndOC = FALSE ) 

Second, we analyze the results on the Cyp phylogeny

#d <- readData("Data/Test/youpi.tsv", "Data/Test/tralala.txt")
d <- readData("Data/h0ha_cyp.tsv", "Data/h0_gbgc_cyp.tsv")

All sites

plotAllMethods(2000, 1, 8)

plotAllMethods(2000, 1, 6, withPCAndOC = FALSE)

Zooming in

plotAllMethods(150, 1, 8)

plotAllMethods(150, 1, 6, withPCAndOC = FALSE)

Violin plots

plotViolinPlots(1, 8)

plotViolinPlots(1, 6, withPCAndOC = FALSE)

Cumulative plots

tmp <- plotCumSums (1, 8, 150 ) 

tmp <- plotCumSums (1, 6, 150, withPCAndOC = FALSE )