In [ ]:

    

library(tidyverse)
library(ggplot2)
library(vegan)


    

In [ ]:

    

sra.md = read.delim("metadata/sra2name.tab")
#head(sra.md)


    

In [ ]:

    

edw.md = read.delim("edwards-data/greenhouse_metadata.txt")
#head(edw.md)


    

In [ ]:

    

# Join and match metadata
metadata = sra.md %>%
    mutate(name=gsub('_', '.', name)) %>%
    inner_join(edw.md, by=c("name"="SampleID")) %>%
    arrange(runid) %>%
    mutate(Compartment=factor(Compartment, levels=c("Bulk Soil", "Rhizosphere", "Rhizoplane", "Endosphere")))
head(metadata)


    

In [ ]:

    

kwip = read.delim("kwip/greenhouse_wip.dist", row.names=1) %>%
       as.matrix()


    

In [ ]:

    

kwip.sras = rownames(kwip)
kwip.names = metadata$name[match(kwip.sras, metadata$runid)]
rownames(kwip) = colnames(kwip) = kwip.names


    

In [ ]:

    

kwip.mds = cmdscale(as.dist(kwip), eig=T)
kwip.var = (kwip.mds$eig / sum(kwip.mds$eig))[1:2]
round(kwip.var * 100, 2)


    

In [ ]:

    

wuf = read.delim("edwards-data//weighted.gh.unifrac", row.names=1, sep=" ") %>% as.matrix()


    

In [ ]:

    

all(rownames(wuf) == rownames(kwip))


    

In [ ]:

    

wuf.mds = cmdscale(as.dist(wuf),eig=T)
wuf.var = (wuf.mds$eig / sum(wuf.mds$eig))[1:2]
round(wuf.var * 100, 2)


    

In [ ]:

    

uuf = read.delim("edwards-data/unweighted.gh.unifrac", row.names=1, sep=" ") %>% as.matrix()


    

In [ ]:

    

all(rownames(uuf) == rownames(kwip))


    

In [ ]:

    

uuf.mds = cmdscale(as.dist(uuf),eig=T)
uuf.var = (uuf.mds$eig / sum(uuf.mds$eig))[1:2]
round(uuf.var * 100, 2)


    

In [ ]:

    

edw.colours = c("#E41A1C", "#984EA3", "#4DAF4A", "#377EB8")


    

In [ ]:

    

all(rownames(kwip.mds) == rownames(wuf.mds))


    

In [ ]:

    

all(rownames(kwip.mds) == rownames(uuf.mds))


    

In [ ]:

    

plot.dat = data.frame(
        PC1.kwip = kwip.mds$points[,1],
        PC2.kwip = kwip.mds$points[,2],
        PC1.wuf = wuf.mds$points[,1],
        PC2.wuf = wuf.mds$points[,2],
        PC1.uuf = uuf.mds$points[,1],
        PC2.uuf = uuf.mds$points[,2],
        name = rownames(kwip.mds$points)
)

plot.dat = left_join(plot.dat, metadata)


    

In [ ]:

    

# kWIP
p = ggplot(plot.dat, aes(x=PC1.kwip, y=PC2.kwip, colour=Compartment, shape=Site)) +
    geom_point(alpha=0.75, size=2) +
    scale_color_manual(values = edw.colours) +
    xlab("PC1") + 
    ylab("PC2") +
    theme_bw() +
    ggtitle("kWIP") +
    theme(panel.grid=element_blank())

svg("gh_kwip.svg", width=4, height = 3)
print(p)
dev.off()
print(p)


    

In [ ]:

    

# WUF
p = ggplot(plot.dat, aes(x=PC1.wuf, y=PC2.wuf, colour=Compartment, shape=Site)) +
    geom_point(alpha=0.75, size=2) +
    scale_color_manual(values = edw.colours) +
    xlab("PC1") + 
    ylab("PC2") +
    theme_bw() +
    ggtitle("Weighted UniFrac") +
    theme(panel.grid=element_blank())

svg("gh_wuf.svg", width=4, height = 3)
print(p)
dev.off()
print(p)


    

In [ ]:

    

# UUF
p = ggplot(plot.dat, aes(x=PC1.uuf, y=PC2.uuf, colour=Compartment, shape=Site)) +
    geom_point(alpha=0.75, size=2) +
    scale_color_manual(values = edw.colours) +
    xlab("PC1") + 
    ylab("PC2") +
    theme_bw() +
    ggtitle("UniFrac") +
    theme(panel.grid=element_blank())

svg("gh_uuf.svg", width=4, height = 3)
print(p)
dev.off()
print(p)


    

In [ ]:

    

distcor = function(a, b, method="spearman") {
    a = as.matrix(a)
    a = a[lower.tri(a)]
    b = as.matrix(b)
    b = b[lower.tri(b)]
    cor(a, b, method=method)
    
}


    

In [ ]:

    

distcor(kwip, uuf)


    

In [ ]:

    

distcor(kwip, wuf)


    

In [ ]:

    

distcor(uuf, wuf)


    

In [ ]: