library(reshape2)
library(ggplot2)
library(vegan)
set.seed(1)


#### Males ####
load('clean_data/Figure3/Male.RData')

#renaming data..
names = Diet_Male_sample_name;
names$Sample_name = as.character(names$Sample_name)
u_dm = unweighted_unifrac_dm
u_pc1 = unweighted_unifrac_pc1[,-(1:2)]
u_pc2 = unweighted_unifrac_pc2
w_dm = weighted_unifrac_dm
w_pc1 = weighted_unifrac_pc1[,-(1:2)]
rm(Diet_Male_sample_name,
   unweighted_unifrac_pc1,
   unweighted_unifrac_pc2,
   unweighted_unifrac_dm,
   weighted_unifrac_pc1,
   weighted_unifrac_pc2)


#sample name data
X = data.frame('P'=rep(0,113),
               'G'=rep(0,113),
               'W'=rep(0,113),
               'M'=rep(0,113),
               'Diet'=rep('',113))

X$P = substr(sapply(strsplit(names$Sample_name, 'P'), function(x) x[2]), 1, 1)
X$G = substr(sapply(strsplit(names$Sample_name, 'G'), function(x) x[2]), 1, 1)
X$W = substr(sapply(strsplit(names$Sample_name, 'W'), function(x) x[2]), 1, 1)
X$M = (sapply(strsplit(names$Sample_name, 'M'), function(x) x[2]))
X$Diet = names$Diet
X$names = names$Sample_name
X$matchmouse = paste0(X$G,X$M)
ind = which(table(X$matchmouse) != 2)
X = X[-ind, ]


#rearrange the samples order
u_pc1 = u_pc1[, match(X$names, colnames(u_pc1))]
w_pc1 = w_pc1[, match(X$names, colnames(w_pc1))]
u_dm = u_dm[match(X$names, rownames(u_dm)), match(X$names, colnames(u_dm))]
w_dm = w_dm[match(X$names, rownames(w_dm)), match(X$names, colnames(w_dm))]

# recreate pc plots
# plot(as.numeric(u_pc1[1,-(1:2)]), as.numeric(u_pc1[2,-(1:2)]),
#      col = as.factor(X$Diet))
# plot(as.numeric(w_pc1[1,-(1:2)]) ~ as.numeric(w_pc1[2,-(1:2)]),
#      col = as.factor(X$Diet))

u_svd = svd(u_dm)
plot(u_svd$v[,1] ~ u_svd$v[,2], col = X$Diet)

w_svd = svd(w_dm)
plot(w_svd$v[,1] ~ w_svd$v[,2], col = X$Diet)

# heat(u_dm, X)
# heat(w_dm, X)

#### adonis ####
res_w_m = adonis(w_dm ~ Diet, X, strata = X$matchmouse, permutations = 1e4)
res_u_m = adonis(u_dm ~ Diet, X, strata = X$matchmouse, permutations = 1e4)


adonis_result = list(res_w_f = res_w,
                     res_w_m = res_w_m,
                     res_u_f = res_u,
                     res_u_m = res_u_m)

#save(adonis_result, file = 'adonis_result.Rdata')