library(vegan)
set.seed(1)


load('clean_data/Figure5/Female.RData')

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


#sample name data
X = data.frame('P'=rep(0,nrow(names)),
               'G'=rep(0,nrow(names)),
               'W'=rep(0,nrow(names)),
               'F'=rep(0,nrow(names)),
               'Diet'=rep('',nrow(names)))

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$F = (sapply(strsplit(names$Sample_name, 'F'), function(x) x[2]))
X$Age = names$Weeks
X$names = names$Sample_name
X$matchmouse = paste0(X$G,X$F)

#remove those that do not have three observations for three timepoints
ind = which(table(X$matchmouse) != 3)
remove = names(table(X$matchmouse)[ind])
X = X[!X$matchmouse %in% remove, ]

#match the order of the subjects
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))]


#exploratory pca plot
u_svd = svd(u_dm)
plot(u_svd$v[,1] ~ u_svd$v[,2], col = X$Age)
w_svd = svd(w_dm)
plot(w_svd$v[,1] ~ w_svd$v[,2], col = X$Age)

#adonis
F_res_w = adonis(w_dm ~ Age, X, strata = X$matchmouse, permutations = 1e4)
F_res_u = adonis(u_dm ~ Age, X, strata = X$matchmouse, permutations = 1e4)

#total number of permutations : 3!^6 = 46656

#only compare 0 vs 12
ind = which(X$Age %in% c('week0', 'week12'))
tmpX = X[ind, ]
tmp_w_dm = w_dm[ind, ind]
tmp_u_dm = u_dm[ind, ind]
F_res_w_0_12 = adonis(tmp_w_dm ~ Age, tmpX, strata = tmpX$matchmouse, permutations = 1e4)
F_res_u_0_12 = adonis(tmp_u_dm ~ Age, tmpX, strata = tmpX$matchmouse, permutations = 1e4)

#only compare 12 vs 24
ind = which(X$Age %in% c('week12', 'week24'))
tmpX = X[ind, ]
tmp_w_dm = w_dm[ind, ind]
tmp_u_dm = u_dm[ind, ind]
F_res_w_12_24 = adonis(tmp_w_dm ~ Age, tmpX, strata = tmpX$matchmouse, permutations = 1e5)
F_res_u_12_24 = adonis(tmp_u_dm ~ Age, tmpX, strata = tmpX$matchmouse, permutations = 1e5)

result = list(F_res_w = F_res_w,
     F_res_u = F_res_u,
     F_res_w_0_12 = F_res_w_0_12,
     F_res_u_0_12 = F_res_u_0_12,
     F_res_w_12_24 = F_res_w_12_24,
     F_res_u_12_24 = F_res_u_12_24)

save(result, file = 'Female_adonis_results.Rdata')