leffj / mctoolsr

Microbial community analysis tools in R
http://leffj.github.io/mctoolsr/
20 stars 8 forks source link

n_perm on calc_pairwise_permanovas() #27

Closed padpadpadpad closed 7 years ago

padpadpadpad commented 7 years ago

Hi

Firstly, thank you for making the package calc_pairwise_permanovas() is awesome! I just had a problem where only 999 permutations was giving me non-reproducible results for which pairs were significantly different from each other.

Was wondering if it would be worth allowing the user to specify more permutations if they wanted. Would take nearly hardly any changes. I have already done it but would be nice to just update the package here.

calc_pairwise_permanovas <- function (dm, metadata_map, compare_header, n_perm) 
{
  comp_var = metadata_map[, compare_header]
  comp_pairs = combn(levels(comp_var), 2)
  pval = c()
  R2 = c()
  for (i in 1:ncol(comp_pairs)) {
    pair = comp_pairs[, i]
    dm_w_map = list(dm_loaded = dm, map_loaded = metadata_map)
    dm_w_map$map_loaded$in_pair = comp_var %in% pair
    dm_w_map_filt = mctoolsr::filter_dm(dm_w_map, filter_cat = "in_pair", 
                              keep_vals = TRUE)
    m = vegan::adonis(dm_w_map_filt$dm_loaded ~ dm_w_map_filt$map_loaded[, compare_header], permutations = n_perm)
    pval = c(pval, m$aov.tab$`Pr(>F)`[1])
    R2 = c(R2, m$aov.tab$R2[1])
  }
  results = data.frame(t(comp_pairs), R2, pval)
  results$pvalBon = pval * length(pval)
  results$pvalFDR = round(pval * (length(pval)/rank(pval, ties.method = "average")), 
                          3)
  results
}

Let me know and I can do a pull request if thats easiest.

Cheers Dan

padpadpadpad commented 7 years ago

Merged.