ConstrainedNJ()

.github/workflows/R-CMD-check.yml

@@ -84,10 +84,10 @@ jobs:
           remotes::install_cran("rcmdcheck")
         shell: Rscript {0}
 
-      - name: Install vdiffr for checks (macOS)
+      - name: Install covr, vdiffr for checks (macOS)
         if: runner.os == 'macOS'
         run: |
-          remotes::install_cran("vdiffr", type = 'source')
+          remotes::install_cran(c("covr", "vdiffr"))
           library('vdiffr')
         shell: Rscript {0}
 

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: TreeSearch
 Title: Phylogenetic Tree Search Using Custom Optimality Criteria
-Version: 0.4.3.9008
+Version: 0.4.3.9009
 Authors@R: c(person("Martin R.", 'Smith',
                   email="martin.smith@durham.ac.uk", 
                   role=c("aut", "cre", "cph"), 

NAMESPACE

@@ -19,6 +19,7 @@ export(ClusteringConcordance)
 export(ConcordantInfo)
 export(ConcordantInformation)
 export(Consistency)
+export(ConstrainedNJ)
 export(DoNothing)
 export(EasyTrees)
 export(EasyTreesy)

NEWS.md

@@ -1,10 +1,14 @@
-# TreeSearch 0.4.3.9008 (development)
+# TreeSearch 0.4.3.9009 (development)
+
+## Changes in functionality
 
  * v......9007 introduces a pre-ratchet TBR search and adjusts default settings.
  * v0.4.3.9005 fixes some bugs that may affect tree scoring under IW & PP.
 
+## New functions
  - `EasyTrees()` 'shiny' graphical user interface for tree search
- - `PlotCharacter()` reconstructs character distributions on trees.
+ - `PlotCharacter()` reconstructs character distributions on trees
+ - `ConstrainedNJ()` constructs starting trees that respect a constraint
  - `SiteConcordance()` calculates exact site concordance
  - `ConcordantInformation()` evaluates signal:noise of dataset implied by a
    given tree.
@@ -14,7 +18,9 @@
  - `Consistency()` calculates consistency and retention 'indices'
  - `MinimumLength()` and `MaximumLength()` calculate range of possible lengths
    of characters in a dataset on any tree
- - `TreeLength()` supports lists of trees.
+
+## Improvements
+ - `TreeLength()` supports lists of trees
  - Set handling of 'gap' token (-) when creating Morphy object with `gap = `
  - Label nodes with split frequencies using `JackLabels(plot = FALSE)`
  - Support for topological constraints during tree search

R/MaximizeParsimony.R

@@ -30,7 +30,8 @@
 #' @param tree (optional) A bifurcating tree of class \code{\link{phylo}},
 #' containing only the tips listed in `dataset`, from which the search
 #' should begin.
-#' If unspecified, a neighbour-joining tree will be generated from `dataset`.
+#' If unspecified, a neighbour-joining tree will be generated from `dataset`,
+#' respecting any supplied `constraint`.
 #' Edge lengths are not supported and will be deleted.
 #' @param ratchIter Numeric specifying number of iterations of the 
 #' parsimony ratchet (Nixon 1999) to conduct.
@@ -61,9 +62,8 @@
 #' Setting to larger values will include trees suboptimal by up to `tolerance`
 #' in search results, which may improve the accuracy of the consensus tree
 #' (at the expense of resolution) (Smith 2019).
-#' @param constraint Either `NULL` or an object of class `phyDat`. Trees that
-#' are not perfectly compatible with each character in `constraint` will not
-#' be considered during search.
+#' @param constraint An object of class `phyDat`; returned trees will be
+#' perfectly compatible with each character in `constraint`.
 #' See [vignette](https://ms609.github.io/TreeSearch/articles/inapplicable.html)
 #' for further examples.
 #' @param verbosity Integer specifying level of messaging; higher values give
@@ -157,7 +157,7 @@
 #' \insertRef{Smith2019}{TreeSearch}
 #' @encoding UTF-8
 #' @export
-MaximizeParsimony <- function (dataset, tree = NJTree(dataset),
+MaximizeParsimony <- function (dataset, tree,
                                ratchIter = 6L,
                                tbrIter = 2,
                                startIter = 2L, finalIter = 1L,
@@ -166,7 +166,7 @@ MaximizeParsimony <- function (dataset, tree = NJTree(dataset),
                                quickHits = 1 / 3,
                                concavity = Inf,
                                tolerance = sqrt(.Machine$double.eps),
-                               constraint = NULL,
+                               constraint,
                                verbosity = 2L, session = NULL) {
   # Define functions
   .Message <- if (is.null(session)) function (level, ...) {
@@ -463,8 +463,23 @@ MaximizeParsimony <- function (dataset, tree = NJTree(dataset),
     class = 'multiPhylo')
   }
 
+
+  # Define constants
+  epsilon <- tolerance #sqrt(.Machine$double.eps)
+  profile <- .UseProfile(concavity)
+  iw <- is.finite(concavity)
+  constrained <- !missing(constraint)
+  startTime <- Sys.time()
+  stopTime <- startTime + as.difftime(maxTime, units = 'mins')
+
   # Initialize tree
-  if (inherits(tree, 'multiPhylo')) {
+  if (missing(tree)) {
+    if (constrained) {
+      tree <- ConstrainedNJ(dataset, constraint)
+    } else {
+      tree <- NJTree(dataset)
+    }
+  } else if (inherits(tree, 'multiPhylo')) {
     .Message(1L, "Starting search from `tree[[1]]`.")
     tree <- tree[[1]]
   }
@@ -487,13 +502,6 @@ MaximizeParsimony <- function (dataset, tree = NJTree(dataset),
     outgroup <- NA
   }
 
-  # Define constants
-  epsilon <- tolerance #sqrt(.Machine$double.eps)
-  profile <- .UseProfile(concavity)
-  iw <- is.finite(concavity)
-  constrained <- !is.null(constraint)
-  startTime <- Sys.time()
-  stopTime <- startTime + as.difftime(maxTime, units = 'mins')
 
   # Initialize constraints
   if (constrained) {
@@ -859,7 +867,7 @@ Resample <- function (dataset, tree = NJTree(dataset), method = 'jack',
                       ratchIter = 1L, tbrIter = 8L, finalIter = 3L,
                       maxHits = 12L, concavity = Inf,
                       tolerance = sqrt(.Machine$double.eps),
-                      constraint = NULL,
+                      constraint,
                       verbosity = 2L, session = NULL,
                       ...) {
   if (!inherits(dataset, 'phyDat')) {
@@ -897,6 +905,33 @@ Resample <- function (dataset, tree = NJTree(dataset), method = 'jack',
                     verbosity = verbosity, session = session, ...) 
 }
 
+#' Constrained neighbour-joining tree
+#' 
+#' Constructs a neighbour-joining tree such that the tree is consistent with a
+#' constraint.
+#' 
+#' @param weight Numeric specifying degree to upweight characters in
+#' `constraint`.
+#' 
+#' @return `ConstrainedNJ()` returns a tree of class `phylo`.
+#' @importFrom TreeTools NJTree
+#' @inheritParams MaximizeParsimony
+#' @examples
+#' dataset <- TreeTools::MatrixToPhyDat(matrix(
+#'   c(0, 1, 1, 1, 0, 1,
+#'     0, 1, 1, 0, 0, 1), ncol = 2,
+#'   dimnames = list(letters[1:6], NULL)))
+#' constraint <- TreeTools::MatrixToPhyDat(
+#'   c(a = 0, b = 0, c = 0, d = 0, e = 1, f = 1))
+#' plot(ConstrainedNJ(dataset, constraint))
+#' @template MRS
+#' @export
+ConstrainedNJ <- function (dataset, constraint, weight = 12345) {
+  conData <- c(constraint, dataset)
+  attr(conData, 'weight')[seq_len(attr(constraint, 'nr'))] <- weight
+  NJTree(conData)
+}
+
 #' Launch tree search graphical user interface
 #' 
 #' @rdname MaximizeParsimony

tests/testthat/test-MaximizeParsimony.R

@@ -35,6 +35,17 @@ test_that("Constraints work", {
 
 })
 
+test_that("Constrained NJ trees work", {
+  dataset <- MatrixToPhyDat(matrix(
+   c(0, 1, 1, 1, 0, 1,
+     0, 1, 1, 0, 0, 1), ncol = 2,
+   dimnames = list(letters[1:6], NULL)))
+  constraint <- MatrixToPhyDat(c(a = 0, b = 0, c = 0, d = 0, e = 1, f = 1))
+  expect_equal(ape::read.tree(text = "(a, (d, ((c, b), (e, f))));"),
+               ConstrainedNJ(dataset, constraint))
+  expect_equal(NJTree(dataset), ConstrainedNJ(dataset, dataset))
+})
+
 test_that("Inconsistent constraints fail", {
   constraint <- MatrixToPhyDat(matrix(
     c(0, 1, 1, 1, 0, 0,