

In [ ]:

    
library(tidyverse,  warn.conflicts=F)
library(ggplot2)

Setup

sra2name maps SRRs (from the distance matrices) to sample names (that the metadata use)



In [ ]:

    
sra2name = read.delim("sra2name.tab")

domds runs PCoA (metric mds) over a distance matrix, returning an nicer data structure



In [ ]:

    
domds = function (dm, k=2) {
    mds = cmdscale(as.dist(dm), k, eig=T, x.ret=T)
    eigs = mds$eig
    pct.contrib = round(eigs / sum(eigs) * 100, 1)
    return(list(points=as.data.frame(mds$points, col.names=c("PC1", "PC2")),
                pct.contrib=pct.contrib, mds.x=mds$x, eigs=eigs))
}

Greenhouse data

kWIP



In [ ]:

    
dm = as.matrix(read.delim("kwip/greenhouse_wip.dist", header=T, row.names=1))
metadata = read.delim("nw-data/rice_metagenome_greenhouse_metadata.txt")

# Add the run id to metadata so we can match w/ dist matrices
metadata$runid = sra2name[match(metadata$SampleID, sra2name$name), 1]



In [ ]:

    
mds = domds(dm)
plot_df = mds$points
plot_df$Site = metadata$Site
plot_df$Compartment = metadata$Compartment
str(plot_df)



In [ ]:

    
p = ggplot(plot_df, aes(x=V1, y=V2, colour=Site, shape=Compartment)) + 
    geom_point(size=2) + 
    scale_color_brewer(palette="Paired") +
    xlab(paste0("PC 1 (", mds$pct.contrib[1], "%)")) +
    ylab(paste0("PC 2 (", mds$pct.contrib[2], "%)")) +
    theme_classic() +
    theme(panel.border=element_rect(colour = "black", fill=NA),
          legend.position='bottom')

pdf("greenhouse.pdf")
print(p)
dev.off()

print(p)

Timecourse

including the samples from the greenhouse



In [ ]:

    
dm = as.matrix(read.delim("kwip/timecourse_ghdavis_wip.dist", header=T, row.names=1))
metadata = read.delim("nw-data/rice_metagenome_TC_and_DavisGH_metadata.txt")

# Add the run id to metadata so we can match w/ dist matrices
metadata$runid = sra2name[match(metadata$SampleID, sra2name$name), 1]



In [ ]:

    
mds = domds(dm)
plot_df = mds$points
plot_df$Compartment = metadata$Compartment
plot_df$Days = as.factor(as.character(metadata$Days))
str(plot_df)



In [ ]:

    
p = ggplot(plot_df, aes(x=V1, y=V2, colour=Days, shape=Compartment)) + 
    geom_point(size=2) + 
    scale_color_brewer(palette="Paired") +
    xlab(paste0("PC 1 (", mds$pct.contrib[1], "%)")) +
    ylab(paste0("PC 2 (", mds$pct.contrib[2], "%)")) +
    theme_classic() +
    theme(panel.border=element_rect(colour = "black", fill=NA),
          legend.position='bottom')

pdf("timecourse.pdf")
print(p)
dev.off()
print(p)



In [ ]:



In [ ]: