diff --git a/DESCRIPTION b/DESCRIPTION
index 454285a..a74ff78 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,12 @@
 Package: MicEco
 Title: Various functions for microbial community data
-Version: 0.9.4
+Version: 0.9.5
 Authors@R: person("Jakob", "Russel", email = "russel2620@gmail.com", role = c("aut", "cre"))
 Description: Collection of functions for microbiome analyses. E.g. fitting neutral models and standardized effect sizes of phylogenetic beta diversities, and much more.
 Depends: R (>= 3.2.5)
-Imports: stats, utils, phyloseq, foreach, doSNOW, picante, vegan, snow, bbmle, Hmisc, abind
+Imports: stats, utils, phyloseq, foreach, doSNOW, picante, vegan, snow, bbmle, Hmisc, abind, reshape2
 License: GPL (>= 3) | file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 6.1.1
+RoxygenNote: 7.0.2
 Remotes: bioc::release/phyloseq
diff --git a/NAMESPACE b/NAMESPACE
index 4af906a..c1afa76 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -8,6 +8,7 @@ export(community_rrna)
 export(neutral.fit)
 export(neutral.rand)
 export(proportionality)
+export(ps_refactor)
 export(rarefy_rrna)
 export(rarefy_rrna.matrix)
 export(rarefy_rrna.phyloseq)
@@ -30,6 +31,7 @@ import(picante)
 import(utils)
 import(vegan)
 importFrom(abind,abind)
+importFrom(reshape2,dcast)
 importFrom(stats,as.dist)
 importFrom(stats,cov)
 importFrom(stats,dnorm)
diff --git a/R/ps_refactor.R b/R/ps_refactor.R
new file mode 100644
index 0000000..849800f
--- /dev/null
+++ b/R/ps_refactor.R
@@ -0,0 +1,18 @@
+#' Relevel the Sample variable in a psmelted phyloseq object 
+#' 
+#' Relevel the Sample variable in a psmelted phyloseq object, 
+#' such that similar samples are plotted together with ggplot barcharts 
+#' 
+#' @param psmelted A phyloseq object melted into a data.frame with psmelt
+#' @param ... Arguments passed to hclust
+#' @importFrom reshape2 dcast
+#' @export
+ps_refactor <- function(psmelted, ...){
+    if(!is(psmelted, "data.frame")){
+        stop("Input should be a phyloseq object melted into a data.frame with psmelt")
+    }
+    mat <- reshape2::dcast(OTU ~ Sample, value.var = "Abundance", data = psmelted)   
+    hc <- hclust(dist(t(mat[, -1])), ...)
+    psmelted$Sample <- factor(psmelted$Sample, levels = hc$labels[hc$order])
+    return(psmelted)
+}    
\ No newline at end of file
diff --git a/README.md b/README.md
index adbb802..d578143 100644
--- a/README.md
+++ b/README.md
@@ -20,6 +20,10 @@ Calculate the unbiased effect size estimation (partial) omega-squared for adonis
 
 Rarefaction curve (theoretical and fast) from a phyloseq object. Output ready for plotting in ggplot2
 
+#### ps_refactor
+
+Relevel the Sample variable in a psmelted phyloseq object, such that similar samples are plotted together with ggplot barcharts.
+
 #### UniFrac.multi
 
 With unrooted phylogenies UniFrac sets the root randomly on the tree. 
diff --git a/man/comdist.par.Rd b/man/comdist.par.Rd
index 6351473..2a591a0 100644
--- a/man/comdist.par.Rd
+++ b/man/comdist.par.Rd
@@ -4,8 +4,7 @@
 \alias{comdist.par}
 \title{Inter-community mean pairwise distance}
 \usage{
-comdist.par(comm, dis, abundance.weighted = FALSE, cores = 1,
-  progress = TRUE)
+comdist.par(comm, dis, abundance.weighted = FALSE, cores = 1, progress = TRUE)
 }
 \arguments{
 \item{comm}{Community data matrix with samples as rows}
diff --git a/man/comdistnt.par.Rd b/man/comdistnt.par.Rd
index 5d9843f..dd6946c 100644
--- a/man/comdistnt.par.Rd
+++ b/man/comdistnt.par.Rd
@@ -4,8 +4,14 @@
 \alias{comdistnt.par}
 \title{Inter-community mean nearest taxon distance}
 \usage{
-comdistnt.par(comm, dis, abundance.weighted = FALSE,
-  exclude.conspecifics = FALSE, cores = 1, progress = TRUE)
+comdistnt.par(
+  comm,
+  dis,
+  abundance.weighted = FALSE,
+  exclude.conspecifics = FALSE,
+  cores = 1,
+  progress = TRUE
+)
 }
 \arguments{
 \item{comm}{Community data matrix with samples as rows}
diff --git a/man/neutral.rand.Rd b/man/neutral.rand.Rd
index a9aeadb..f9fe8a2 100644
--- a/man/neutral.rand.Rd
+++ b/man/neutral.rand.Rd
@@ -4,8 +4,15 @@
 \alias{neutral.rand}
 \title{Fit Sloan et al. (2006) Neutral Model several times}
 \usage{
-neutral.rand(data, n = NULL, s = NULL, rRNA = NULL, rn = NULL,
-  cores = 1, naming = NULL)
+neutral.rand(
+  data,
+  n = NULL,
+  s = NULL,
+  rRNA = NULL,
+  rn = NULL,
+  cores = 1,
+  naming = NULL
+)
 }
 \arguments{
 \item{data}{A phyloseq object}
diff --git a/man/ps_refactor.Rd b/man/ps_refactor.Rd
new file mode 100644
index 0000000..d1913e8
--- /dev/null
+++ b/man/ps_refactor.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ps_refactor.R
+\name{ps_refactor}
+\alias{ps_refactor}
+\title{Relevel the Sample variable in a psmelted phyloseq object}
+\usage{
+ps_refactor(psmelted, ...)
+}
+\arguments{
+\item{psmelted}{A phyloseq object melted into a data.frame with psmelt}
+
+\item{...}{Arguments passed to hclust}
+}
+\description{
+Relevel the Sample variable in a psmelted phyloseq object, 
+such that similar samples are plotted together with ggplot barcharts
+}
diff --git a/man/rarefy_rrna.Rd b/man/rarefy_rrna.Rd
index 1cd84c6..d0c5fc9 100644
--- a/man/rarefy_rrna.Rd
+++ b/man/rarefy_rrna.Rd
@@ -6,8 +6,7 @@
 \alias{rarefy_rrna.phyloseq}
 \title{Rarefy and normalize based on 16S rRNA copy numbers}
 \usage{
-rarefy_rrna(x, reads, copy.database = "v13.5", seed = NULL,
-  trim = FALSE)
+rarefy_rrna(x, reads, copy.database = "v13.5", seed = NULL, trim = FALSE)
 
 rarefy_rrna.matrix(x, reads, copy.database, seed = NULL, trim)
 
diff --git a/man/rcurve.Rd b/man/rcurve.Rd
index 81ce9ee..3dd49ab 100644
--- a/man/rcurve.Rd
+++ b/man/rcurve.Rd
@@ -5,8 +5,7 @@
 \title{Rarefaction curve on phyloseq object
 Theoretical richness with vegan's rarefy function}
 \usage{
-rcurve(physeq, subsamp = 10^c(1:5), trim = TRUE,
-  add_sample_data = TRUE)
+rcurve(physeq, subsamp = 10^c(1:5), trim = TRUE, add_sample_data = TRUE)
 }
 \arguments{
 \item{physeq}{phyloseq object}
diff --git a/man/ses.UniFrac.Rd b/man/ses.UniFrac.Rd
index 82e028c..fef28da 100644
--- a/man/ses.UniFrac.Rd
+++ b/man/ses.UniFrac.Rd
@@ -4,10 +4,20 @@
 \alias{ses.UniFrac}
 \title{Standardized effect size of unifrac}
 \usage{
-ses.UniFrac(physeq, method = "taxa.labels", fixedmar = "both",
-  shuffle = "both", strata = NULL, mtype = "count", burnin = 0,
-  thin = 1, weighted = TRUE, normalized = TRUE, runs = 99,
-  cores = 1)
+ses.UniFrac(
+  physeq,
+  method = "taxa.labels",
+  fixedmar = "both",
+  shuffle = "both",
+  strata = NULL,
+  mtype = "count",
+  burnin = 0,
+  thin = 1,
+  weighted = TRUE,
+  normalized = TRUE,
+  runs = 99,
+  cores = 1
+)
 }
 \arguments{
 \item{physeq}{phyloseq-class, containing at minimum a phylogenetic tree and otu table}
diff --git a/man/ses.comdist.Rd b/man/ses.comdist.Rd
index 00d9ce2..49a2b18 100644
--- a/man/ses.comdist.Rd
+++ b/man/ses.comdist.Rd
@@ -4,10 +4,16 @@
 \alias{ses.comdist}
 \title{Standardized effect size of inter-community MPD (betaMPD, betaNRI)}
 \usage{
-ses.comdist(samp, dis, null.model = c("taxa.labels", "richness",
-  "frequency", "sample.pool", "phylogeny.pool", "independentswap",
-  "trialswap"), abundance.weighted = FALSE, runs = 999,
-  iterations = 1000, cores = 1)
+ses.comdist(
+  samp,
+  dis,
+  null.model = c("taxa.labels", "richness", "frequency", "sample.pool",
+    "phylogeny.pool", "independentswap", "trialswap"),
+  abundance.weighted = FALSE,
+  runs = 999,
+  iterations = 1000,
+  cores = 1
+)
 }
 \arguments{
 \item{samp}{Community data matrix with samples as rows}
diff --git a/man/ses.comdist2.Rd b/man/ses.comdist2.Rd
index 0d31b1b..5194d9b 100644
--- a/man/ses.comdist2.Rd
+++ b/man/ses.comdist2.Rd
@@ -4,9 +4,20 @@
 \alias{ses.comdist2}
 \title{Standardized effect size of inter-community MPD (betaMPD, betaNRI)}
 \usage{
-ses.comdist2(samp, dis, method = "quasiswap", fixedmar = "both",
-  shuffle = "both", strata = NULL, mtype = "count", burnin = 0,
-  thin = 1, abundance.weighted = FALSE, runs = 999, cores = 1)
+ses.comdist2(
+  samp,
+  dis,
+  method = "quasiswap",
+  fixedmar = "both",
+  shuffle = "both",
+  strata = NULL,
+  mtype = "count",
+  burnin = 0,
+  thin = 1,
+  abundance.weighted = FALSE,
+  runs = 999,
+  cores = 1
+)
 }
 \arguments{
 \item{samp}{Community data matrix with samples as rows}
diff --git a/man/ses.comdistnt.Rd b/man/ses.comdistnt.Rd
index bd0f505..325e45a 100644
--- a/man/ses.comdistnt.Rd
+++ b/man/ses.comdistnt.Rd
@@ -4,11 +4,17 @@
 \alias{ses.comdistnt}
 \title{Standardized effect size of inter-community MNTD (betaMNTD, betaNTI)}
 \usage{
-ses.comdistnt(samp, dis, null.model = c("taxa.labels", "richness",
-  "frequency", "sample.pool", "phylogeny.pool", "independentswap",
-  "trialswap"), abundance.weighted = FALSE,
-  exclude.conspecifics = FALSE, runs = 999, iterations = 1000,
-  cores = 1)
+ses.comdistnt(
+  samp,
+  dis,
+  null.model = c("taxa.labels", "richness", "frequency", "sample.pool",
+    "phylogeny.pool", "independentswap", "trialswap"),
+  abundance.weighted = FALSE,
+  exclude.conspecifics = FALSE,
+  runs = 999,
+  iterations = 1000,
+  cores = 1
+)
 }
 \arguments{
 \item{samp}{Community data matrix with samples as rows}
diff --git a/man/ses.comdistnt2.Rd b/man/ses.comdistnt2.Rd
index 366f902..4ce72a3 100644
--- a/man/ses.comdistnt2.Rd
+++ b/man/ses.comdistnt2.Rd
@@ -4,10 +4,21 @@
 \alias{ses.comdistnt2}
 \title{Standardized effect size of inter-community MNTD (betaMNTD, betaNTI)}
 \usage{
-ses.comdistnt2(samp, dis, method = "quasiswap", fixedmar = "both",
-  shuffle = "both", strata = NULL, mtype = "count", burnin = 0,
-  thin = 1, abundance.weighted = FALSE, exclude.conspecifics = FALSE,
-  runs = 999, cores = 1)
+ses.comdistnt2(
+  samp,
+  dis,
+  method = "quasiswap",
+  fixedmar = "both",
+  shuffle = "both",
+  strata = NULL,
+  mtype = "count",
+  burnin = 0,
+  thin = 1,
+  abundance.weighted = FALSE,
+  exclude.conspecifics = FALSE,
+  runs = 999,
+  cores = 1
+)
 }
 \arguments{
 \item{samp}{Community data matrix with samples as rows}
diff --git a/man/ses.mntd.par.Rd b/man/ses.mntd.par.Rd
index 202c1b8..5f1293c 100644
--- a/man/ses.mntd.par.Rd
+++ b/man/ses.mntd.par.Rd
@@ -4,10 +4,16 @@
 \alias{ses.mntd.par}
 \title{Standardized effect size of MNTD}
 \usage{
-ses.mntd.par(samp, dis, null.model = c("taxa.labels", "richness",
-  "frequency", "sample.pool", "phylogeny.pool", "independentswap",
-  "trialswap"), abundance.weighted = FALSE, runs = 999,
-  iterations = 1000, cores = 1)
+ses.mntd.par(
+  samp,
+  dis,
+  null.model = c("taxa.labels", "richness", "frequency", "sample.pool",
+    "phylogeny.pool", "independentswap", "trialswap"),
+  abundance.weighted = FALSE,
+  runs = 999,
+  iterations = 1000,
+  cores = 1
+)
 }
 \arguments{
 \item{samp}{Community data matrix with samples as rows}
diff --git a/man/ses.mpd.par.Rd b/man/ses.mpd.par.Rd
index 1181358..4110eb7 100644
--- a/man/ses.mpd.par.Rd
+++ b/man/ses.mpd.par.Rd
@@ -4,10 +4,16 @@
 \alias{ses.mpd.par}
 \title{Standardized effect size of MPD}
 \usage{
-ses.mpd.par(samp, dis, null.model = c("taxa.labels", "richness",
-  "frequency", "sample.pool", "phylogeny.pool", "independentswap",
-  "trialswap"), abundance.weighted = FALSE, runs = 999,
-  iterations = 1000, cores = 1)
+ses.mpd.par(
+  samp,
+  dis,
+  null.model = c("taxa.labels", "richness", "frequency", "sample.pool",
+    "phylogeny.pool", "independentswap", "trialswap"),
+  abundance.weighted = FALSE,
+  runs = 999,
+  iterations = 1000,
+  cores = 1
+)
 }
 \arguments{
 \item{samp}{Community data matrix with samples as rows}