In [ ]:

    

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


    

In [ ]:

    

kwip_mds = function (filename, pdfname=NULL, title=NULL, dim.mult=c(1, 1)) {
    dm = as.matrix(read.delim(filename, header=T, row.names=1))
    metadata = read.delim("chlamy_meta.tab")
    metadata = metadata[match(row.names(dm), metadata$Run),]
    mds = cmdscale(dm, k=2, eig=T, x.ret=T)
    eigs = mds$eig
    pct.contrib = round(eigs / sum(eigs) * 100)
    #plot(cumsum(pct.contrib), type='l', ylim=c(0,100))
    
    pts = mds$points
    for (i in 1:ncol(pts)) {
        pts[,i] = pts[,i] * dim.mult[i]
    }
    pts.df = as.data.frame(pts)
    pts.df$Group = metadata$origin

    cols = c("Laboratory"="light blue", "Northeast"= "blue", "Southeast"="dark green", "West"="red" )
    p = ggplot(pts.df, aes(x=V1, y=V2, colour=Group)) + 
        geom_point(size=2) + 
        scale_color_manual(values = cols) +
        xlab(paste0("Dim 1 (", pct.contrib[1], "%)")) +
        ylab(paste0("Dim 2 (", pct.contrib[2], "%)")) +
        theme_classic() +
        theme(panel.border=element_rect(colour = "black", fill=NA),
              legend.position="bottom")
    if (!is.null(title)) {
        p = p + ggtitle(title)
    }
    print(p)

    if (!is.null(pdfname)) {
        pdf(pdfname, width=3.5, height=3.5)
        print(p)
        dev.off()
    }
}


    

In [ ]:

    

kwip_mds("kwip/flowers_wip.dist", title="All", dim.mult = c(-1,1), pdfname = "flowers_all.pdf")


    

In [ ]:

    

kwip_mds("kwip/flowers_wild_wip.dist", title="Wild", pdfname = "flowers_wild.pdf")


    

In [ ]:

    

kwip_mds("kwip/flowers_lab_wip.dist", title = "Lab", pdfname = "flowers_lab.pdf")


    

In [ ]: