In [ ]:

    

library(plyr,  warn.conflicts=F)
library(dplyr, warn.conflicts=F)
library(tidyr, warn.conflicts=F)
library(ggplot2)
library(caTools)
library(gtable)


    

In [ ]:

    

metadata = read.delim("chlamy_meta.tab")


    

In [ ]:

    

read.distmat =  function (filename) {
    dm = as.matrix(read.delim(filename, header=T, row.names=1))
    idxs = match(metadata$Run, row.names(dm))
    return(dm[idxs, idxs])
}


    

In [ ]:

    

full = read.distmat("kwip/full_wip.dist")
metadata = metadata[match(row.names(full), metadata$Run),]


    

In [ ]:

    

c = cmdscale(as.dist(full))


    

In [ ]:

    

# From https://gist.github.com/baptiste/6652815
gtable_arrange <- function(..., grobs=list(), as.table=TRUE,
                           top = NULL, bottom = NULL, 
                           left = NULL, right = NULL, draw=TRUE){
  require(gtable)
  # alias
  gtable_add_grobs <- gtable_add_grob
  
  dots <- list(...)
  params <- c("nrow", "ncol", "widths", "heights",
              "respect", "just", "z") # TODO currently ignored
  
  layout.call <- intersect(names(dots), params)
  params.layout <- dots[layout.call]
  
  if(is.null(names(dots)))
    not.grobnames <- FALSE else
      not.grobnames <- names(dots) %in% layout.call
  
  if(!length(grobs))
  grobs <- dots[! not.grobnames ]
  
  ## figure out the layout
  n <- length(grobs)
  nm <- n2mfrow(n)
  
  if(is.null(params.layout$nrow) & is.null(params.layout$ncol)) 
  {
    params.layout$nrow = nm[1]
    params.layout$ncol = nm[2]
  }
  if(is.null(params.layout$nrow))
    params.layout$nrow = ceiling(n/params.layout$ncol)
  if(is.null(params.layout$ncol))
    params.layout$ncol = ceiling(n/params.layout$nrow)
  
  if(is.null(params.layout$widths))
    params.layout$widths <- unit(rep(1, params.layout$ncol), "null")
  if(is.null(params.layout$heights))
    params.layout$heights <- unit(rep(1,params.layout$nrow), "null")

  positions <- expand.grid(row = seq_len(params.layout$nrow), 
                           col = seq_len(params.layout$ncol))
  if(as.table) # fill table by rows
    positions <- positions[order(positions$row),]
  
  positions <- positions[seq_along(grobs), ] # n might be < ncol*nrow
  
  ## build the gtable, similar steps to gtable_matrix
  
  gt <- gtable(name="table")
  gt <- gtable_add_cols(gt, params.layout$widths)
  gt <- gtable_add_rows(gt, params.layout$heights)
  gt <- gtable_add_grobs(gt, grobs, t = positions$row, 
                            l = positions$col)
  
  ## titles given as strings are converted to text grobs
  if (is.character(top)) 
    top <- textGrob(top)
  if (is.character(bottom)) 
    bottom <- textGrob(bottom)
  if (is.character(right)) 
    right <- textGrob(right, rot = -90)
  if (is.character(left)) 
    left <- textGrob(left, rot = 90)
  
  if(!is.null(top)){
    gt <- gtable_add_rows(gt, heights=grobHeight(top), 0)
    gt <- gtable_add_grobs(gt, top, t=1, l=1, r=ncol(gt))
  }
  if(!is.null(bottom)){
    gt <- gtable_add_rows(gt, heights=grobHeight(bottom), -1)
    gt <- gtable_add_grobs(gt, bottom, t=nrow(gt), l=1, r=ncol(gt))
  }
  if(!is.null(left)){
    gt <- gtable_add_cols(gt, widths=grobWidth(left), 0)
    gt <- gtable_add_grobs(gt, left, t=1, b=nrow(gt), l=1, r=1)
  }
  if(!is.null(right)){
    gt <- gtable_add_cols(gt, widths=grobWidth(right), -1)
    gt <- gtable_add_grobs(gt, right, t=1, b=nrow(gt), l=ncol(gt), r=ncol(gt))
  }
  
  if(draw){
   grid.newpage()
   grid.draw(gt)
  }
  invisible(gt)
}


    

In [ ]:

    

coverages = c("0.01x", "0.1x", "0.5x", "1x", "2x", "4x", "8x", "12x", "15x", "25x",
              "50x", "75x", "100x", "150x", "200x", "full")

plot_covs = c("0.1x","1x", "2x",
              "4x", "8x", "15x",
              "50x", "150x", "full")
plots = list()
pdf("all-pcoas.pdf")
for (coverage in coverages) {
    fname = paste0("kwip/", coverage, "_wip.dist")
    mat = read.distmat(fname)
    mds = cmdscale(mat, k=2, eig=T, x.ret=T)
    eigs = mds$eig
    pct.contrib = round(eigs / sum(eigs) * 100)
    
    # Invert axes to match the paper (Flowers et al.) figure.
    # The sample here is one of the two red ones in the top right corner.
    if (mds$points["SRR1734600", 1] < 0) {
        mds$points[,1] = mds$points[,1] * -1
    }
    if (mds$points["SRR1734600", 2] < 0) {
        mds$points[,2] = mds$points[,2] * -1
    }
    
    pts.df = as.data.frame(mds$points)
    pts.df$Group = metadata$origin
    
    cols = c("light blue", "blue", "dark green", "red" )
    p = ggplot(pts.df, aes(x=V1, y=V2, colour=Group)) + 
        geom_point(size=2) + 
        scale_color_manual(values = cols) +
        xlab(paste0("PC 1 (", pct.contrib[1], "%)")) +
        ylab(paste0("PC 2 (", pct.contrib[2], "%)")) +
        ggtitle(paste(coverage, "fold subset")) + 
        theme_classic() +
        theme(panel.border=element_rect(colour = "black", fill=NA),
              legend.position="bottom")
    print(p)
    if (coverage  %in% plot_covs) {
        p = ggplot(pts.df, aes(x=V1, y=V2, colour=Group)) + 
            geom_point(size=2) + 
            scale_color_manual(values = cols) +
            ggtitle(coverage) + 
            theme_classic() +
            theme(panel.border=element_rect(colour = "black", fill=NA),
                  legend.position="none",
                  axis.title.x=element_blank(),
                  axis.title.y=element_blank(),
                  axis.ticks=element_blank(),
                  axis.text.x = element_blank(),
                  axis.text.y = element_blank()
                 )
        plots = c(plots, list(p))
    }
}
dev.off()


    

In [ ]:

    

pdf("subset-pcoa-matrix.pdf", height = 5, width = 5)
gtable_arrange(ncol=3, grobs = lapply(plots, ggplotGrob), left="PC1", bottom="PC2")
dev.off()


    

In [ ]: