library("stringr")
library("dplyr")
library("reshape2")
library("vegan")
dir.wrk <- str_replace(getwd(), "/scripts", "")
dir.data <- file.path(dir.wrk, "data/data_permutation_withReplacement")
dir.output <- file.path(dir.wrk, "data/data_processed")
TransformLog <- function(dat) {
dml <- vegan::decostand(dat[, 5:10], "log")
dml <- cbind(dat[, -c(5:10)], dml)
return(dml)
}
getPermanova <- function(dat){
# GET DISTANCE MATRIX ---
dist_dml <- vegan::vegdist(x=as.matrix(dat[,5:10]), method="euclidean", binary=FALSE, diag=TRUE, upper=TRUE, na.rm=FALSE)
# PERMANOVA TEST ---
set.seed(12345)
y_permanova <- vegan::adonis(dist_dml ~ IncomeGroup+Ethnicity+EducationLevel+GeoRegion,
data=dat, permutations=999, method="euclidean", parallel=8)
return(y_permanova)
}
nperm <- 1000
for (i in 1:nperm) {
# DEFINE FILE ---
file.dat <- file.path(dir.data, paste("abundance_permmute_", i, ".tsv", sep = ""))
# LOAD DATA ---
dat <- read.delim(file.dat, header = TRUE, stringsAsFactors = FALSE)
# LOG TRANSFORM ---
dml <- TransformLog(dat)
# PERFORM PERMANOVA TEST ---
pnova <- getPermanova(dml)
dm <- as.data.frame(pnova$aov.tab)
# WRITE OUTPUT ---
file.output <- file.path(dir.data, paste("stats_permanova_permmute_", i, ".tsv",
sep = ""))
write.table(dm, file.output, sep = "\t", row.names = TRUE, col.names = NA, quote = FALSE)
}