diff --git a/R/branchLengthMetrics.R b/R/branchLengthMetrics.R
index 553c9f9..c8959a1 100644
--- a/R/branchLengthMetrics.R
+++ b/R/branchLengthMetrics.R
@@ -29,7 +29,14 @@ branchLengthMetrics <- function(tree) {
   # Dists with long branches are normally distributed
   # with clusters there is a peak below the mean which looks like a cutoff
   distances <- cophenetic.phylo(tree)[upper.tri(cophenetic.phylo(tree))]
-  ls$clustsShapiro <- shapiro.test(distances)$statistic
+
+  # shapiro-wilk test depends on sample size, and isn't implemented for s > 5000; >~70 tips
+  shap_dists <- distances
+  if (length(distances) > 400)
+  {
+    shap_dists <- sample(distances, size = 400)
+  }
+  ls$clustsShapiro <- as.numeric(shapiro.test(shap_dists)$statistic)
 
   # Larger for positive skew -> close tips; 1 for symmetric -> equally distant tips
   quantiles <- quantile(distances, probs = c(0.25,0.5,0.75))
diff --git a/R/scoreTree.R b/R/scoreTree.R
new file mode 100644
index 0000000..755926d
--- /dev/null
+++ b/R/scoreTree.R
@@ -0,0 +1,52 @@
+#' Scores an observed tree against a known simulation. Scores are the quantile the observed tree metric is in
+#'
+#' \code{scoreTree} returns the scoring of an observed tree
+#' Each row is a sequential pruning, each column a metric from allMetrics.
+#'
+#' @param obs_tree a phylogenetic tree to score (as a \code{phylo} object)
+#' @param pruned a set of trees from prunePermute to score against (as a \code{list} of \code{phylo} objects)
+#' @return An \code{array} of scores.
+#' @author John Lees (\email{lees.john6@@gmail.com})
+#' @export
+#' @examples
+#' sim_tree <- simulate()
+#' sim_tree2 <- simulate()
+#' pruned <- prunePermute(sim_tree$tree)
+#' obs_tree <- sim_tree2$tree
+#' scoreTree(obs_tree, pruned)
+
+scoreTree <- function(obs_tree, pruned_trees)
+{
+  # metrics calculated for simulated and observed passed trees
+  # TODO will the sim_metrics be pre-calculated and passed? Hopefully, but
+  # assume not for now
+  sim_metrics <- lapply(pruned$seqPrunes, lapply, allMetrics)
+  real_metrics <- lapply(timeprune(obs_tree)$trees, allMetrics)
+
+  # Loop over each metric within each pruning
+  scores <- array(dim=c(length(real_metrics),length(real_metrics[[1]])))
+  for (prunedTree in 1:length(real_metrics))
+  {
+    for (metric in 1:length(real_metrics[[1]]))
+    {
+      # This gives the quantile the observed score is in
+      scoreRange <- rep(NA, length=length(sim_metrics[[1]]))
+      for (permutation in 1:length(sim_metrics[[1]]))
+      {
+        scoreRange[permutation] <- sim_metrics[[prunedTree]][[permutation]][[metric]]
+      }
+      scoreRange <- sort(scoreRange)
+      scores[prunedTree,metric] <-
+        sum(scoreRange < real_metrics[[prunedTree]][[metric]])/length(scoreRange)
+    }
+  }
+
+  return(scores)
+}
+
+# Wrapper for branchLengthMetrics and imbalanceMetrics
+allMetrics <- function(tree)
+{
+  all_metrics <- c(imbalanceMetrics(tree), branchLengthMetrics(tree))
+  return(all_metrics)
+}