In [ ]:

    

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


    

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 [ ]:

    

metadata = metadata[match(row.names(dm), metadata$Run),]


    

In [ ]:

    

full = read.distmat("kwip/full_wip.dist")


    

In [ ]:

    

coverages = c("0.5x", "1x", "2x", "4x", "8x", "12x", "15x", "25x", "50x", "100x", "full")
matricies = list()
for (coverage in coverages) {
    fname = paste0("kwip/", coverage, "_wip.dist")
    mat = read.distmat(fname)
    matricies = c(matricies, list(mat))
    mds = cmdscale(mat, k=2, eig=T, x.ret=T)
    eigs = mds$eig
    pct.contrib = round(eigs / sum(eigs) * 100)
    
    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("PC1") +
        # ylab("PC2") +
        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)
}


    

In [ ]:

    

# The negative here is to match the Flowers et al. paper
md.pts = mds$points


    

In [ ]:

    

pts.df = as.data.frame(md.pts)
pts.df$Group = metadata$origin


    

In [ ]:

    

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("PC1") +
    ylab("PC2") +
    # xlab(paste0("PC 1 (", pct.contrib[1], "%)")) +
    # ylab(paste0("PC 2 (", pct.contrib[2], "%)")) +
    theme_classic() +
    theme(panel.border=element_rect(colour = "black", fill=NA),
          legend.position="bottom")
print(p)


    

In [ ]: