# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .r
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: R
#     language: R
#     name: ir
# ---

# + deletable=true editable=true
library(tidyverse)
library(ggplot2)
library(vegan)

# + [markdown] deletable=true editable=true
# ## Metadata

# + deletable=true editable=true
sra.md = read.delim("metadata/sra2name.tab")
#head(sra.md)

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

# + deletable=true editable=true
# 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)

# + [markdown] deletable=true editable=true
# ## kWIP

# + deletable=true editable=true
kwip = read.delim("kwip/greenhouse_wip.dist", row.names=1) %>%
       as.matrix()

# + deletable=true editable=true
kwip.sras = rownames(kwip)
kwip.names = metadata$name[match(kwip.sras, metadata$runid)]
rownames(kwip) = colnames(kwip) = kwip.names

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

# + [markdown] deletable=true editable=true
# ## Weighted UF

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

# + deletable=true editable=true
all(rownames(wuf) == rownames(kwip))

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

# + [markdown] deletable=true editable=true
# ## UUF

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

# + deletable=true editable=true
all(rownames(uuf) == rownames(kwip))

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

# + [markdown] deletable=true editable=true
# # Plotting 
#

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

# + deletable=true editable=true
all(rownames(kwip.mds) == rownames(wuf.mds))

# + deletable=true editable=true
all(rownames(kwip.mds) == rownames(uuf.mds))

# + deletable=true editable=true
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)

# + deletable=true editable=true
# 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)

# + deletable=true editable=true
# 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)

# + deletable=true editable=true
# 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)

# + [markdown] deletable=true editable=true
# # Distance Matrix Correlation
#
#

# + deletable=true editable=true
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)
    
}

# + deletable=true editable=true
distcor(kwip, uuf)

# + deletable=true editable=true
distcor(kwip, wuf)

# + deletable=true editable=true
distcor(uuf, wuf)

# + [markdown] deletable=true editable=true
# # Summary of results
#
# #### Axis pct contributions:
#
# |metric | PC1 | PC2 |
# |-------|-----|-----|
# |kWIP   |22.6 |15.8 |
# |WUF    |46.4 |11.5 |
# |UUF    |18.1 |14.9 | 
#
#
# #### Correlations between metrics:
#
# - kwip-> WUF: 0.88
# - kwip-> UUF: 0.90
# - WUF-> UUF: 0.73
# -



StarcoderdataJupyterScriptsDedupFiltered