In [ ]:
library(plyr, warn.conflicts=F)
library(dplyr, warn.conflicts=F)
library(tidyr, warn.conflicts=F)
library(ggplot2)
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))
}
In [ ]:
dm = as.matrix(read.delim("kwip/field-expt_wip.dist", header=T, row.names=1))
metadata = read.delim("nw-data/rice_metagenome_field_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("field-expt.pdf")
print(p)
dev.off()
print(p)
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)
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)