diff --git a/.github/workflows/release-file-system-image.yml b/.github/workflows/release-file-system-image.yml
new file mode 100644
index 00000000..87bf570d
--- /dev/null
+++ b/.github/workflows/release-file-system-image.yml
@@ -0,0 +1,17 @@
+# Workflow derived from https://github.com/r-wasm/actions/tree/v1/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  release:
+    # Must republish release to update assets
+    types: [ published ]
+
+name: Build and deploy wasm R package image
+
+jobs:
+  release-file-system-image:
+    uses: r-wasm/actions/.github/workflows/release-file-system-image.yml@v1
+    permissions:
+      # For publishing artifact files to the release
+      contents: write
+      # To download GitHub Packages within action
+      repository-projects: read
diff --git a/.github/workflows/rhub.yaml b/.github/workflows/rhub.yaml
new file mode 100644
index 00000000..74ec7b05
--- /dev/null
+++ b/.github/workflows/rhub.yaml
@@ -0,0 +1,95 @@
+# R-hub's generic GitHub Actions workflow file. It's canonical location is at
+# https://github.com/r-hub/actions/blob/v1/workflows/rhub.yaml
+# You can update this file to a newer version using the rhub2 package:
+#
+# rhub::rhub_setup()
+#
+# It is unlikely that you need to modify this file manually.
+
+name: R-hub
+run-name: "${{ github.event.inputs.id }}: ${{ github.event.inputs.name || format('Manually run by {0}', github.triggering_actor) }}"
+
+on:
+  workflow_dispatch:
+    inputs:
+      config:
+        description: 'A comma separated list of R-hub platforms to use.'
+        type: string
+        default: 'linux,windows,macos'
+      name:
+        description: 'Run name. You can leave this empty now.'
+        type: string
+      id:
+        description: 'Unique ID. You can leave this empty now.'
+        type: string
+
+jobs:
+
+  setup:
+    runs-on: ubuntu-latest
+    outputs:
+      containers: ${{ steps.rhub-setup.outputs.containers }}
+      platforms: ${{ steps.rhub-setup.outputs.platforms }}
+
+    steps:
+    # NO NEED TO CHECKOUT HERE
+    - uses: r-hub/actions/setup@v1
+      with:
+        config: ${{ github.event.inputs.config }}
+      id: rhub-setup
+
+  linux-containers:
+    needs: setup
+    if: ${{ needs.setup.outputs.containers != '[]' }}
+    runs-on: ubuntu-latest
+    name: ${{ matrix.config.label }}
+    strategy:
+      fail-fast: false
+      matrix:
+        config: ${{ fromJson(needs.setup.outputs.containers) }}
+    container:
+      image: ${{ matrix.config.container }}
+
+    steps:
+      - uses: r-hub/actions/checkout@v1
+      - uses: r-hub/actions/platform-info@v1
+        with:
+          token: ${{ secrets.RHUB_TOKEN }}
+          job-config: ${{ matrix.config.job-config }}
+      - uses: r-hub/actions/setup-deps@v1
+        with:
+          token: ${{ secrets.RHUB_TOKEN }}
+          job-config: ${{ matrix.config.job-config }}
+      - uses: r-hub/actions/run-check@v1
+        with:
+          token: ${{ secrets.RHUB_TOKEN }}
+          job-config: ${{ matrix.config.job-config }}
+
+  other-platforms:
+    needs: setup
+    if: ${{ needs.setup.outputs.platforms != '[]' }}
+    runs-on: ${{ matrix.config.os }}
+    name: ${{ matrix.config.label }}
+    strategy:
+      fail-fast: false
+      matrix:
+        config: ${{ fromJson(needs.setup.outputs.platforms) }}
+
+    steps:
+      - uses: r-hub/actions/checkout@v1
+      - uses: r-hub/actions/setup-r@v1
+        with:
+          job-config: ${{ matrix.config.job-config }}
+          token: ${{ secrets.RHUB_TOKEN }}
+      - uses: r-hub/actions/platform-info@v1
+        with:
+          token: ${{ secrets.RHUB_TOKEN }}
+          job-config: ${{ matrix.config.job-config }}
+      - uses: r-hub/actions/setup-deps@v1
+        with:
+          job-config: ${{ matrix.config.job-config }}
+          token: ${{ secrets.RHUB_TOKEN }}
+      - uses: r-hub/actions/run-check@v1
+        with:
+          job-config: ${{ matrix.config.job-config }}
+          token: ${{ secrets.RHUB_TOKEN }}
diff --git a/NEWS b/NEWS
index 9416a4ee..a4aa1058 100644
--- a/NEWS
+++ b/NEWS
@@ -1,4 +1,4 @@
-    CHANGES in PHANGORN VERSION 3.0.0
+    CHANGES in PHANGORN VERSION 2.12.0
 
 NEW FEATURES
 
@@ -18,15 +18,21 @@ NEW FEATURES
 
       character vector for the node argument and not only integers.
 
-    o Ancestral state reconstructionwas in parts rewritten. Output is now an
+    o Ancestral state reconstruction was rewritten. anc_pml and anc_pars are now
 
-      object of class ancestral. E.g. ancestral states of constant sites are now
+      prefered over ancestral.pml and ańcestral.pars.
 
-      always this state.
+      Output is now an object of class ancestral, which contains tree and
+
+      reconstructions.
+
+      Joint reconstruction for ML is now available.
+
+      Ancestral states of constant sites are now always this state.
 
     o Improvements to several plot functions to get nicer plots out of the box.
 
-    o Nicer defaults for plot.pml, mainly for rooted trees.
+      Nicer defaults for plot.pml, mainly for rooted trees.
 
     o plotAnc got an argument scheme allowing to use different color schemes,
 
@@ -38,8 +44,7 @@ NEW FEATURES
 
     o new function composition_test comparing indicating possible problems
 
-      with base / state composition.
-
+      with base / state composition as in iqtree.
 
 OTHER CHANGES
 
diff --git a/R/Coalescent.R b/R/Coalescent.R
index 21925ad1..d572c43c 100644
--- a/R/Coalescent.R
+++ b/R/Coalescent.R
@@ -64,7 +64,7 @@ comp2 <- function(x, y) {
 #' @param sTree A species tree which fixes the topology.
 #' @return The function returns an object of class \code{phylo}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com} Emmanuel Paradies
-#' @seealso \code{\link{speciesTree}}
+#' @seealso \code{\link[ape]{speciesTree}}
 #' @references Liu, L., Yu, L. and Pearl, D. K. (2010) Maximum tree: a
 #' consistent estimator of the species tree. \emph{Journal of Mathematical
 #' Biology}, \bold{60}, 95--106.
diff --git a/R/Densi.R b/R/Densi.R
index 8708fe78..ec48e18a 100644
--- a/R/Densi.R
+++ b/R/Densi.R
@@ -97,8 +97,8 @@ add_tiplabels <- function(xy, tip.label, direction, adj, font, srt = 0, cex = 1,
 #' @param \dots further arguments to be passed to plot.
 #' @returns \code{densiTree} returns silently x.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{plot.phylo}}, \code{\link{plot.networx}},
-#' \code{\link{jitter}}, \code{\link{rtt}}
+#' @seealso \code{\link[ape]{plot.phylo}}, \code{\link{plot.networx}},
+#' \code{\link{jitter}}, \code{\link[ape]{rtt}}
 #' @references densiTree is inspired from the great
 #' \href{https://www.cs.auckland.ac.nz/~remco/DensiTree/}{DensiTree} program of
 #' Remco Bouckaert.
diff --git a/R/addConfidences.R b/R/addConfidences.R
index 37b33788..6f6d1c08 100644
--- a/R/addConfidences.R
+++ b/R/addConfidences.R
@@ -15,7 +15,7 @@
 #' @return The object \code{x} with added bootstrap / MCMC support values.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{as.splits}}, \code{\link{as.networx}},
-#' \code{\link{RF.dist}}, \code{\link{plot.phylo}}
+#' \code{\link{RF.dist}}, \code{\link[ape]{plot.phylo}}
 #' @references Schliep, K., Potts, A. J., Morrison, D. A. and Grimm, G. W.
 #' (2017), Intertwining phylogenetic trees and networks.
 #' \emph{Methods Ecol Evol}.\bold{8}, 1212--1220. doi:10.1111/2041-210X.12760
diff --git a/R/ancestral.R b/R/ancestral.R
index ce2a45d9..96ec2b48 100644
--- a/R/ancestral.R
+++ b/R/ancestral.R
@@ -29,7 +29,7 @@
 #' @param data an object of class phyDat
 #' @param type method used to assign characters to internal nodes, see details.
 #' @param cost A cost matrix for the transitions between two states.
-## @param return return a \code{phyDat} object or matrix of probabilities.
+#' @param return return a \code{phyDat} object or matrix of probabilities.
 ##  @param x an object of class ancestral.
 #' @param \dots Further arguments passed to or from other methods.
 #' @return An object of class ancestral. This is a list containing the tree with
@@ -42,7 +42,7 @@
 ## returned.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{pml}}, \code{\link{parsimony}}, \code{\link[ape]{ace}},
-#' \code{\link{plotAnc}}, \code{\link{ltg2amb}}, \code{\link{latag2n}},
+#' \code{\link{plotAnc}}, \code{\link{ltg2amb}}, \code{\link[ape]{latag2n}},
 #' \code{\link{gap_as_state}}, \code{\link[ape]{root}},
 #' \code{\link[ape]{makeNodeLabel}}
 #' @references Felsenstein, J. (2004). \emph{Inferring Phylogenies}. Sinauer
@@ -61,13 +61,8 @@
 #' # generate node labels to ensure plotting will work
 #' tree <- makeNodeLabel(tree)
 #' fit <- pml(tree, Laurasiatherian)
-#' anc.ml <- ancestral.pml(fit, type = "ml")
-#' anc.p <- ancestral.pars(tree, Laurasiatherian)
-## \dontrun{
-## require(seqLogo)
-## seqLogo( t(subset(anc.ml, 48, 1:20)[[1]]), ic.scale=FALSE)
-## seqLogo( t(subset(anc.p, 48, 1:20)[[1]]), ic.scale=FALSE)
-## }
+#' anc.ml <- anc_pml(fit)
+#' anc.p <- anc_pars(tree, Laurasiatherian)
 #' # plot ancestral sequences at the root
 #' plotSeqLogo( anc.ml, 48, 1, 20)
 #' plotSeqLogo( anc.p, 48, 1, 20)
@@ -78,9 +73,10 @@
 #'
 #' @rdname ancestral.pml
 #' @export
-ancestral.pml <- function(object, type = "marginal", ...) {
+ancestral.pml <- function(object, type = "marginal", return = "prob", ...) {
   call <- match.call()
   pt <- match.arg(type, c("marginal", "ml", "bayes")) # "joint",
+  rt <- match.arg(return, c("prob", "phyDat", "ancestral"))
   tree <- object$tree
   INV <- object$INV
   inv <- object$inv
@@ -107,6 +103,9 @@ ancestral.pml <- function(object, type = "marginal", ...) {
   tree$node.label <- node_label
   tmp <- length(data)
 
+  joint <- TRUE
+  if(length(w) > 1 || object$inv > 0) joint <- FALSE
+
   eig <- object$eig
 
   bf <- object$bf
@@ -165,17 +164,34 @@ ancestral.pml <- function(object, type = "marginal", ...) {
       result2[[k]][ind] <- data[[1]][ind]
     }
   }
+#  browser()
   attrib$names <- node_label
   attributes(result2) <- attrib
+  if(rt == "phyDat") return(rbind(data, result2))
+  if(rt == "prob") {
+    tmp <- c(unclass(new2old.phyDat(data)), result)
+    attrib$names <- c(tree$tip.label, node_label)
+    attributes(tmp) <- attrib
+    return(tmp)
+  }
   attributes(result) <- attrib
-  result <- list2df_ancestral(result, result2)
-  result2 <- compress.phyDat(result2)
+  if(joint) result2 <- joint_pml(object)
+#  result <- list2df_ancestral(result, result2)
+#  result2 <- compress.phyDat(result2)
   erg <- list(tree=tree, data=data, prob=result, state=result2)
   class(erg) <- "ancestral"
   erg
 }
 
 
+#' @rdname ancestral.pml
+#' @export
+anc_pml <- function(object, type = "marginal", ...) {
+  ancestral.pml(object, type=type, return="ancestral")
+}
+
+
+
 #' Export and convenience functions for  ancestral reconstructions
 #'
 #' \code{write.ancestral} allows to export ancestral reconstructions. It writes
@@ -192,21 +208,22 @@ ancestral.pml <- function(object, type = "marginal", ...) {
 #' @examples
 #' data(Laurasiatherian)
 #' fit <- pml_bb(Laurasiatherian[,1:100], "JC", rearrangement = "none")
-#' anc_ml <- ancestral.pml(fit)
+#' anc_ml <- anc_pml(fit)
 #' write.ancestral(anc_ml)
 #' # Can be also results from iqtree
-#' align <- read.phyDat("ancestral_align.fasta")
+#' align <- read.phyDat("ancestral_align.fasta", format="fasta")
 #' tree <- read.tree("ancestral_tree.nwk")
 #' df <- read.table("ancestral.state", header=TRUE)
-#' anc_ml_disc <- ancestral(tree, align, df)
+#' anc_ml_disc <- as.ancestral(tree, align, df)
 #' plotAnc(anc_ml_disc, 20)
 #' unlink(c("ancestral_align.fasta", "ancestral_tree.nwk", "ancestral.state"))
 #' @rdname write.ancestral
 #' @export
 write.ancestral <- function(x, file="ancestral"){
   stopifnot(inherits(x, "ancestral"))
-  write.phyDat(x$data, file=paste0(file, "_align.fasta"))
-  write.table(x$prob, file=paste0(file, ".state"), quote=FALSE, row.names=FALSE,
+  write.phyDat(x$data, file=paste0(file, "_align.fasta"), format="fasta")
+  tab <- as.data.frame(x)
+  write.table(tab, file=paste0(file, ".state"), quote=FALSE, row.names=FALSE,
               sep="\t")
   write.tree(x$tree, file=paste0(file, "_tree.nwk"))
   invisible(x)
@@ -220,13 +237,14 @@ write.ancestral <- function(x, file="ancestral"){
 #' @importFrom utils head
 #' @rdname write.ancestral
 #' @export
-ancestral <- function(tree, align, prob){
+as.ancestral <- function(tree, align, prob){
   stopifnot(inherits(tree, "phylo"))
   stopifnot(inherits(align, "phyDat"))
   stopifnot(inherits(prob, "data.frame"))
   if(is.null(tree$node.label))stop("tree needs node.label")
   state <- extract_states(prob, attr(align, "type"),
                           levels=attr(align, "levels"))
+  prob <- extract_prob(prob, align, tree$node.label)
   erg <- list(tree=tree, data=align[tree$tip.label], prob=prob,
               state=state[tree$node.label])
   class(erg) <- "ancestral"
@@ -242,6 +260,21 @@ extract_states <- function(x, type, levels=NULL){
   phyDat(y, type=type, levels=levels)
 }
 
+extract_prob <- function(x, align, node_label){
+  index <- attr(align, "index")
+  idx <- which(!duplicated(index))
+  att <- attributes(align)
+  res <- vector("list", length(node_label))
+  for(i in seq_along(node_label)){
+     tmp <- x[x$"Node"==node_label[i], -c(1:3)] |> as.matrix()
+     dimnames(tmp) <- NULL
+     res[[i]] <- as.matrix(tmp)[idx,]
+  }
+  att$names <- node_label
+  attributes(res) <- att
+  res
+}
+
 #' @rdname write.ancestral
 #' @export
 print.ancestral <- function(x, ...){
@@ -249,8 +282,8 @@ print.ancestral <- function(x, ...){
   print(x$tree)
   cat("\n")
   print(x$data)
-  cat("\n")
-  print(head(x$prob))
+#  cat("\n")
+#  print(head(x$prob))
 }
 
 #' @rdname ancestral.pml
@@ -311,6 +344,46 @@ list2df_ancestral <- function(x, y=NULL, ...) {
   res
 }
 
+#' @rdname ancestral.pml
+#' @export
+as.data.frame.ancestral <- function(x, ...) {
+  stopifnot(inherits(x, "ancestral"))
+  y <- rbind(x$data, x$state)
+  contrast <- attr(x$data, "contrast")
+  Y <- unlist(as.data.frame(y))
+  n_node <- Nnode(x$tree)
+  n_tip <- Ntip(x$tree)
+  l <- length(y)
+  nr <- attr(x$data, "nr")
+  nc <- attr(x$data, "nc")
+  index <- attr(x$data, "index")
+  nr <- length(index)
+  nam <- names(y)
+  X <- matrix(0, l*length(index), nc)
+  j <- 0
+  for(i in seq_len(n_tip)){
+    X[(j+1):(j+nr), ] <- contrast[x$data[[i]][index], ]
+    j <- j + nr
+  }
+  for(i in seq_len(n_node)){
+    X[(j+1):(j+nr), ] <- x$prob[[i]][index, ]
+    j <- j + nr
+  }
+#  if(!is.null(y) & inherits(y, "phyDat")){
+#    y <- y[names(x)]
+#    Y <- unlist(as.data.frame(y))
+  res <- data.frame(Node=rep(nam, each=nr), Site=rep(seq_len(nr), l), Y, X)
+  colnames(res) <- c("Node", "Site", "State", paste0("p_", attr(x$data, "levels")))
+#  } else{
+#    res <- data.frame(Node=rep(nam, each=nr), Site=rep(seq_len(nr), l), X)
+#    colnames(res) <- c("Node", "Site", paste0("p_", attr(x, "levels")))
+#  }
+  rownames(res) <- NULL
+  res
+}
+
+
+
 
 fitchCoding2ambiguous <- function(x, type = "DNA") {
   y <- c(1L, 2L, 4L, 8L, 8L, 3L, 5L, 9L, 6L, 10L, 12L, 7L, 11L, 13L,
@@ -319,46 +392,59 @@ fitchCoding2ambiguous <- function(x, type = "DNA") {
 }
 
 
+
 #' @rdname ancestral.pml
 #' @export
 ancestral.pars <- function(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"),
-                           cost = NULL, ...) {
-  #, return = "prob"
+                           cost = NULL, return = "prob") {
   call <- match.call()
   type <- match.arg(type)
-  tree$node.label <- makeAncNodeLabel(tree, ...)
-  contrast <- attr(data, "contrast")
-  data <- data[tree$tip.label,]
+  tree$nodel.label <- makeAncNodeLabel(tree)
   if (type == "ACCTRAN" || type=="POSTORDER") {
-    #, return = return
-    res <- ptree(tree, data, acctran=(type == "ACCTRAN"))
-#    attr(res, "call") <- call
+    res <- ptree(tree, data, return = return, acctran=(type == "ACCTRAN"))#, tips=TRUE)
+    attr(res, "call") <- call
   }
   if (type == "MPR") {
-    res <- mpr(tree, data, cost = cost)
-    #, return = return
-#    attr(res, "call") <- call
+    res <- mpr(tree, data, cost = cost, return = return) #, tips=TRUE)
+    attr(res, "call") <- call
   }
-  result <- res[[1]]
-  result2 <- res[[2]]
+  res
+}
+
+
+
+#' @rdname ancestral.pml
+#' @export
+anc_pars <- function(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"),
+                           cost = NULL, ...) {
+  #
+  call <- match.call()
+  type <- match.arg(type)
+  tree$node.label <- makeAncNodeLabel(tree, ...)
+  contrast <- attr(data, "contrast")
+  data <- data[tree$tip.label,]
+
+  prob <- ancestral.pars(tree, data, type, cost)
+  joint <- joint_sankoff(tree, data, cost)
   ind <- identical_sites(data)
   if(length(ind)>0){
     for(k in seq_len(Nnode(tree))){
-      result[[k]][ind,] <- contrast[data[[1]][ind],]
-      result2[[k]][ind] <- data[[1]][ind]
+      prob[[k]][ind,] <- contrast[data[[1]][ind],]
+      joint[[k]][ind] <- data[[1]][ind]
     }
   }
 #  attrib$names <- node_label
 #  attributes(result2) <- attrib
 #  attributes(result) <- attrib
-  result <- list2df_ancestral(result, result2)
-  result2 <- compress.phyDat(result2) # needed?
-  erg <- list(tree=tree, data=data, prob=result, state=result2)
+#  result <- list2df_ancestral(result, result2)
+#  result2 <- compress.phyDat(result2) # needed?
+  erg <- list(tree=tree, data=data, prob=prob, state=joint, call=call)
   class(erg) <- "ancestral"
   erg
 }
 
 
+
 #' @rdname ancestral.pml
 #' @export
 pace <- ancestral.pars
@@ -386,11 +472,12 @@ mpr.help <- function(tree, data, cost = NULL) {
   res
 }
 
-# , return = "prob"
-mpr <- function(tree, data, cost = NULL, ...) {
+
+mpr <- function(tree, data, cost = NULL, return="prob", tips=FALSE, ...) {
   data <- subset(data, tree$tip.label)
   att <- attributes(data)
-  att$names <- tree$node.label
+  if(tips) att$names <- c(tree$tip.label, tree$node.label)
+  else att$names <- tree$node.label
   type <- att$type
   nr <- att$nr
   nc <- att$nc
@@ -410,31 +497,31 @@ mpr <- function(tree, data, cost = NULL, ...) {
   }
 #  for (i in 1:ntips) res[[i]] <- contrast[data[[i]], , drop = FALSE]
   for (i in (ntips + 1):m) res[[i]][] <- as.numeric(res[[i]] < RM)
-  res <- res[(ntips + 1):m]
-#  if (return == "prob") {
-    #        for(i in 1:ntips) res[[i]] <- contrast[data[[i]],,drop=FALSE]
-    res_prob <- lapply(res, fun)
-    attributes(res_prob) <- att
-    class(res_prob) <- c("ancestral", "phyDat")
-#  }
+  if(tips)for(i in seq_len(ntips)) res[[i]] <- contrast[data[[i]],,drop=FALSE]
+  else  res <- res[(ntips + 1):m]
+  res_prob <- lapply(res, fun)
+  attributes(res_prob) <- att
+  if(return=="prob") return(res_prob)
+#    class(res_prob) <- c("ancestral", "phyDat")
   #    else res[1:ntips] <- data[1:ntips]
-  fun2 <- function(x) {
-    x <- p2dna(x)
-    fitchCoding2ambiguous(x)
-  }
+  #fun2 <- function(x) {
+  #  x <- p2dna(x)
+  #  fitchCoding2ambiguous(x)
+  #}
 #  if (return != "prob") {
-  if (type == "DNA") {
-    res_state <- lapply(res, fun2)
-    attributes(res_state) <- att
-  }
-  else {
-    attributes(res) <- att
-    res_state <- highest_state(res)
-    attributes(res_state) <- att
-  }
+#  if (type == "DNA") {
+#    res_state <- lapply(res, fun2)
+#    attributes(res_state) <- att
+#  }
+#  else {
+  attributes(res) <- att
+  res_state <- highest_state(res)
+  attributes(res_state) <- att
+  if(tips)res_state[seq_len(ntips)] <- data
+#  }
 #    res[1:ntips] <- data
 #  }
-  list(res_prob, res_state)
+  res_state # list(res_prob, res_state)
 }
 
 
@@ -483,12 +570,13 @@ acctran <- function(tree, data) {
 }
 
 #, return = "prob"
-ptree <- function(tree, data, acctran=TRUE, ...) {
+ptree <- function(tree, data, acctran=TRUE, return = "prob", tips=FALSE, ...) {
   tree <- reorder(tree, "postorder")
   data <- subset(data, tree$tip.label)
   edge <- tree$edge
   att <- attributes(data)
-  att$names <- tree$node.label
+  att$names <- c(tree$tip.label, tree$node.label)
+  #else att$names <- tree$node.label
   nr <- att$nr
   type <- att$type
   m <- max(edge)
@@ -500,33 +588,44 @@ ptree <- function(tree, data, acctran=TRUE, ...) {
   tmp <- tmp[tmp[,2]>Ntip(tree),]
   if(length(tmp)>0 && acctran==TRUE)f$acctran_traverse(tmp)
   res <- res_state <- vector("list", nNode)
-#  res <- vector("list", m)
-  att$names <- tree$node.label #makeAncNodeLabel(tree, ...)
+  res <- vector("list", m)
+  att$names <- c(tree$tip.label, tree$node.label) #makeAncNodeLabel(tree, ...)
 #  else {
     fun <- function(X) {
       rs <- rowSums(X)
       X / rs
     }
     contrast <- att$contrast
-#    for(i in seq_len(nTip)) res[[i]] <- contrast[data[[i]], , drop=FALSE]
+    for(i in seq_len(nTip)) res[[i]] <- contrast[data[[i]], , drop=FALSE]
 #    for(i in (nTip+1):m) res[[i]] <- f$getAnc(i)[1:nr, , drop=FALSE]
-    for(i in seq_len(nNode)) res[[i]] <- f$getAnc(i+nTip)[seq_len(nr), , drop=FALSE]
+    for(i in seq_len(nNode)) res[[i+nTip]] <- f$getAnc(i+nTip)[seq_len(nr), , drop=FALSE]
     res <- lapply(res, fun)
     attributes(res) <- att
-    class(res) <- c("ancestral", "phyDat")
+#    class(res) <- c("ancestral", "phyDat")
 #  }
-  if(type=="DNA"){
-    indx <- c(1, 2, 6, 3, 7, 9, 12, 4, 8, 10, 13, 11, 14, 15, 16)
-    #    res[1:nTip] <- data[1:nTip]
-    for(i in seq_len(nNode)) #  (nTip+1):m)
-      res_state[[i]] <- indx[f$getAncAmb(i+nTip)[1:nr]]
-    attributes(res_state) <- att
-    #    return(res)
-  } else{
-    res_state <- highest_state(res)
-    class(res_state) <- "phyDat"
+  if(!tips) res <- res[tree$node.label]
+
+  if(return=="prob"){
+#    if(tips)
+    return(res)
+#    else return(res[tree$node.label])
   }
-  list(res, res_state)
+
+#  if(type=="DNA"){
+#    indx <- c(1, 2, 6, 3, 7, 9, 12, 4, 8, 10, 13, 11, 14, 15, 16)
+#    res_state[seq_len(nTip)] <- data
+#    for(i in (nTip+1):m)
+#      res_state[[i]] <- indx[f$getAncAmb(i+nTip)[1:nr]]
+#    attributes(res_state) <- att
+#    #    return(res)
+ # } else{
+    res_state <- highest_state(res)
+    attributes(res_state) <-  attributes(res)
+#    class(res_state) <- "phyDat"
+#  }
+#  if(tips) return(res_state)
+#  res_state[tree$node.label]
+    res_state
 }
 
 
diff --git a/R/baseFreq.R b/R/baseFreq.R
index 8e8629d7..fd32dac5 100644
--- a/R/baseFreq.R
+++ b/R/baseFreq.R
@@ -17,7 +17,7 @@
 #'
 #' @return  \code{baseFreq} returns a named vector and \code{glance} a one row
 #' \code{data.frame}.
-#' @seealso \code{\link{phyDat}, \link{base.freq}, \link{glance}}
+#' @seealso \code{\link{phyDat}, \link[ape]{base.freq}, \link{glance}}
 #' @author Klaus Schliep
 #' @examples
 #'
diff --git a/R/bootstrap.R b/R/bootstrap.R
index b0d6ac67..6c2105d1 100644
--- a/R/bootstrap.R
+++ b/R/bootstrap.R
@@ -29,8 +29,8 @@
 #' if not supplied the tree with labels supplied in the \code{node.label} slot.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{optim.pml}}, \code{\link{pml}},
-#' \code{\link{plot.phylo}}, \code{\link{maxCladeCred}}
-#' \code{\link{nodelabels}},\code{\link{consensusNet}} and
+#' \code{\link[ape]{plot.phylo}}, \code{\link{maxCladeCred}}
+#' \code{\link[ape]{nodelabels}},\code{\link{consensusNet}} and
 #' \code{\link{SOWH.test}} for parametric bootstrap
 #' @references Felsenstein J. (1985) Confidence limits on phylogenies. An
 #' approach using the bootstrap. \emph{Evolution} \bold{39}, 783--791
diff --git a/R/cladePar.R b/R/cladePar.R
index fc28b477..103652f8 100644
--- a/R/cladePar.R
+++ b/R/cladePar.R
@@ -14,7 +14,7 @@
 #' @param \dots Further arguments passed to or from other methods.
 #' @return A list containing the information about the edges and tips.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{plot.phylo}}
+#' @seealso \code{\link[ape]{plot.phylo}}
 #' @keywords plot
 #' @examples
 #'
diff --git a/R/codon.R b/R/codon.R
index 64769a3e..e2fa28e4 100644
--- a/R/codon.R
+++ b/R/codon.R
@@ -45,7 +45,7 @@
 #' @return The functions return an object of class \code{phyDat}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @references \url{https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=cgencodes}
-#' @seealso \code{\link{trans}}, \code{\link{phyDat}} and the chapter 4 in the
+#' @seealso \code{\link[ape]{trans}}, \code{\link{phyDat}} and the chapter 4 in the
 #' \code{vignette("phangorn-specials", package="phangorn")}
 #' @keywords cluster
 #' @examples
diff --git a/R/draw_CI.R b/R/draw_CI.R
index 06b89042..b33f9437 100644
--- a/R/draw_CI.R
+++ b/R/draw_CI.R
@@ -58,8 +58,9 @@ edge_length_matrix <- function(tree, trees, rooted=TRUE){
 ##'               add_edge_length(bs)
 ##' plot(tree_compat)
 ##' add_boxplot(tree_compat, bs, boxwex=.7)
-##' @seealso \code{\link{node.depth.edgelength}}, \code{\link{consensus}},
-##' \code{\link{maxCladeCred}}, \code{\link{add_boxplot}}
+##' @seealso \code{\link[ape]{node.depth}},
+##' \code{\link[ape]{consensus}}, \code{\link{maxCladeCred}},
+##' \code{\link{add_boxplot}}
 ##' @keywords aplot
 ##' @export
 add_edge_length <- function(tree, trees, fun=\(x)median(na.omit(x)),
@@ -109,7 +110,7 @@ add_edge_length <- function(tree, trees, fun=\(x)median(na.omit(x)),
 ##' add_ci(tree, trees)
 ##' plot(tree, direction="downwards")
 ##' add_boxplot(tree, trees, boxwex=.7)
-##' @seealso \code{\link{plot.phylo}}, \code{\link{plotBS}},
+##' @seealso \code{\link[ape]{plot.phylo}}, \code{\link{plotBS}},
 ##' \code{\link{add_edge_length}}, \code{\link{maxCladeCred}}
 ##' @keywords aplot
 ##' @rdname add_ci
diff --git a/R/gap_as_state.R b/R/gap_as_state.R
index c6917a11..b7d1555e 100644
--- a/R/gap_as_state.R
+++ b/R/gap_as_state.R
@@ -10,8 +10,9 @@
 #' @param ambiguous a character which codes for the ambiguous state
 #' @return The functions return an object of class \code{phyDat}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{phyDat}}, \code{\link{ltg2amb}}, \code{\link{latag2n}},
-#' \code{\link{ancestral.pml}}, \code{\link{gap_as_state}}
+#' @seealso \code{\link{phyDat}}, \code{\link{ltg2amb}},
+#' \code{\link[ape]{latag2n}}, \code{\link{ancestral.pml}},
+#' \code{\link{gap_as_state}}
 #' @keywords cluster
 #' @examples
 #' data(Laurasiatherian)
diff --git a/R/image_phyDat.R b/R/image_phyDat.R
index 8a686477..33ea8156 100644
--- a/R/image_phyDat.R
+++ b/R/image_phyDat.R
@@ -7,7 +7,7 @@
 #' @param x	 an object containing sequences, an object of class \code{phyDat}.
 #' @param ... further arguments passed to or from other methods.
 #' @returns Nothing. The function is called for plotting.
-#' @seealso \code{\link{image.DNAbin}}, \code{\link{image.AAbin}}
+#' @seealso \code{\link[ape]{image.DNAbin}}, \code{\link[ape]{image.AAbin}}
 #' @examples
 #' data("chloroplast")
 #' image(chloroplast[, 1:50], scheme="Clustal", show.aa = TRUE)
diff --git a/R/joint_ASR.R b/R/joint_ASR.R
index e279b46d..df332b1d 100644
--- a/R/joint_ASR.R
+++ b/R/joint_ASR.R
@@ -70,3 +70,70 @@ joint_pml <- function(x){
   res
 }
 
+
+# Joint reconstruction for parsimony
+joint_sankoff <- function(tree, data, cost=NULL){
+  stopifnot(inherits(data, "phyDat"))
+  stopifnot(inherits(tree, "phylo"))
+  tree <- reorder(tree, "postorder")
+  edge <- tree$edge
+  ntip <- Ntip(tree)
+  if (is.null(cost)) {
+    levels <- attr(data, "levels")
+    l <- length(levels)
+    cost <- matrix(1, l, l)
+    cost <- cost - diag(l)
+  }
+  P <- cost
+
+  desc <- Descendants(tree, type="children")
+  nr <- attr(data, "nr")
+  nc <- attr(data, "nc")
+  levels <- attr(data, "levels")
+  allLevels <- attr(data, "allLevels")
+
+  l <- nrow(tree$edge)
+  C <- array(NA, c(nr, nc, max(tree$edge)))
+  tmp <- sankoff(tree, data, cost, "data")
+  pscore <- tmp$pscore
+  L <- tmp$dat
+  nn <- unique(edge[,1])
+  pa <- 0
+  root <- edge[l, 1]
+  lnr <- seq_len(nr)
+
+  for(i in seq_len(length(nn))){
+    pa <- nn[i]
+    if(pa==root) break()
+    for(j in 1:nc){
+      pp <- L[[pa]] + rep(P[j,], each=nr)
+      pos <- max.col(-pp)
+      C[,j,pa] <- pos
+    }
+  }
+  pp <- L[[pa]]
+  pos <- max.col(-pp)
+  C[,,pa] <- pos
+  tree <- reorder(tree)
+  if(is.null(tree$node.label)) tree <- makeNodeLabel(tree)
+  res <- vector("list", length(tree$node.label))
+  names(res) <- tree$node.label
+  res[[1]] <- pos
+  att <- attributes(data)
+  att$names <- tree$node.label
+  labels <- c(tree$tip.label, tree$node.label)
+  edge <- tree$edge
+  nrw <- seq_len(nr)
+  for(i in seq_along(edge[,1])){
+    ch_i <- edge[i,2]
+    pa_i <- edge[i,1]
+    if(ch_i > ntip){
+      pos <-res[[labels[pa_i]]]
+      res[[labels[ch_i]]] <- C[cbind(nrw,pos,ch_i)]
+    }
+  }
+  ind <- match(levels, allLevels)
+  for(i in length(res))  res[[i]] <- ind[res[[i]]]
+  attributes(res) <- att
+  res
+}
diff --git a/R/ltg2amb.R b/R/ltg2amb.R
index 09c5aa1a..7b984f72 100644
--- a/R/ltg2amb.R
+++ b/R/ltg2amb.R
@@ -7,7 +7,7 @@
 #' @param amb character of the ambiguous state t replace the gaps.
 #' @param gap gap parameter to replace.
 #' @returns returns an object of class \code{phyDat}.
-#' @seealso \code{\link{latag2n}}, \code{\link{ancestral.pml}},
+#' @seealso \code{\link[ape]{latag2n}}, \code{\link{ancestral.pml}},
 #' \code{\link{gap_as_state}}
 #' @examples
 #' x <- phyDat(matrix(c("-", "A", "G", "-", "T", "C"), 2, 3))
diff --git a/R/maxCladeCred.R b/R/maxCladeCred.R
index 42a3a0c4..ab17d088 100644
--- a/R/maxCladeCred.R
+++ b/R/maxCladeCred.R
@@ -24,10 +24,10 @@
 #' @return a tree (an object of class \code{phylo}) with the highest clade
 #' credibility or a numeric vector of clade credibilities for each tree.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{consensus}}, \code{\link{consensusNet}},
-#' \code{\link{prop.part}}, \code{\link{bootstrap.pml}}, \code{\link{plotBS}},
-#' \code{\link{transferBootstrap}}, \code{\link{add_edge_length}},
-#' \code{\link{add_boxplot}}
+#' @seealso \code{\link[ape]{consensus}}, \code{\link{consensusNet}},
+#' \code{\link[ape]{prop.part}}, \code{\link{bootstrap.pml}},
+#' \code{\link{plotBS}}, \code{\link{transferBootstrap}},
+#' \code{\link{add_edge_length}}, \code{\link{add_boxplot}}
 #' @keywords cluster
 #' @importFrom fastmatch fmatch
 #' @examples
diff --git a/R/phyDat.R b/R/phyDat.R
index 64431b65..935cca99 100644
--- a/R/phyDat.R
+++ b/R/phyDat.R
@@ -177,6 +177,12 @@ rbind.phyDat <- function(...){
   for(i in seq_along(x)){
     res[(mcs[i]+1):mcs[i+1], ] <- as.character(x[[i]])
   }
+  if(types[1]=="USER"){
+    contrast <- attr(x[[1]], "contrast")
+    dimnames(contrast) <- list(attr(x[[1]], "allLevels"),
+                               attr(x[[1]], "levels"))
+    return(phyDat(res, type="USER", contrast=contrast))
+  }
   phyDat(res, type=types[1])
 }
 
diff --git a/R/phyDat_conversion.R b/R/phyDat_conversion.R
index de664302..b46792c7 100644
--- a/R/phyDat_conversion.R
+++ b/R/phyDat_conversion.R
@@ -15,6 +15,7 @@
 #' @aliases
 #' as.phyDat.character as.phyDat.data.frame as.phyDat.matrix
 #' as.MultipleAlignment as.MultipleAlignment.phyDat
+#' as.StringSet as.StringSet.phyDat
 #' acgt2ry phyDat2MultipleAlignment
 #' @param data An object containing sequences.
 #' @param x An object containing sequences.
@@ -28,9 +29,9 @@
 #' @param ... further arguments passed to or from other methods.
 #' @return The functions return an object of class \code{phyDat}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso [DNAbin()], [as.DNAbin()],
+#' @seealso \code{\link[ape]{DNAbin}}, \code{\link[ape]{as.DNAbin}},
 #' \code{\link{baseFreq}}, \code{\link{glance.phyDat}},
-#' \code{\link{read.dna}}, \code{\link{read.nexus.data}}
+#' \code{\link[ape]{read.dna}}, \code{\link[ape]{read.nexus.data}}
 #' and the chapter 1 in the \code{vignette("phangorn-specials",
 #' package="phangorn")} and the example of \code{\link{pmlMix}} for the use of
 #' \code{allSitePattern}
@@ -158,6 +159,31 @@ as.phyDat.DNAStringSet <- function(x, ...){
 }
 
 
+# @rdname phyDat
+#' @export
+as.StringSet <- function (x, ...){
+  if (inherits(x, "StringSet")) return(x)
+  UseMethod("as.StringSet")
+}
+
+
+
+#' @rdname as.phyDat
+#' @export
+as.StringSet.phyDat <- function(x, ...){
+  if (requireNamespace("Biostrings")){
+    z <- as.character(x)
+    type <- attr(x, "type")
+    seq <- switch(type,
+                  DNA = tolower(apply(z, 1, paste, collapse="")),
+                  AA = toupper(apply(z, 1, paste, collapse="")))
+    if(type=="DNA") return(Biostrings::DNAStringSet(seq))
+    if(type=="AA") return(Biostrings::AAStringSet(seq))
+  }
+  return(NULL)
+}
+
+
 # @rdname phyDat
 #' @export
 as.MultipleAlignment <- function (x, ...){
diff --git a/R/phylo.R b/R/phylo.R
index 9f3fa581..ae41ff78 100644
--- a/R/phylo.R
+++ b/R/phylo.R
@@ -1247,7 +1247,8 @@ pml.fit <- function(tree, data, bf = rep(1 / length(levels), length(levels)),
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{pml_bb}}, \code{\link{bootstrap.pml}},
 #' \code{\link{modelTest}}, \code{\link{pmlPart}}, \code{\link{pmlMix}},
-#' \code{\link{plot.phylo}}, \code{\link{SH.test}}, \code{\link{ancestral.pml}}
+#' \code{\link[ape]{plot.phylo}}, \code{\link{SH.test}},
+#' \code{\link{ancestral.pml}}
 #' @references Felsenstein, J. (1981) Evolutionary trees from DNA sequences: a
 #' maximum likelihood approach. \emph{Journal of Molecular Evolution},
 #' \bold{17}, 368--376.
@@ -1340,7 +1341,7 @@ pml <- function(tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
                 rate = 1, model = NULL, site.rate = "gamma", ASC = FALSE, ...) {
   call <- match.call()
   extras <- match.call(expand.dots = FALSE)$...
-  pmla <- c("wMix", "llMix", "dnds", "tstv")
+  pmla <- c("wMix", "llMix", "dnds", "tstv", "w", "g")
   existing <- match(pmla, names(extras))
   wMix <- ifelse(is.na(existing[1]), 0,
     eval(extras[[existing[1]]], parent.frame()))
@@ -1348,14 +1349,18 @@ pml <- function(tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
 #    eval(extras[[existing[2]]], parent.frame()))
   llMix <- 0
   if(!is.na(existing[2])) llMix <- eval(extras[[existing[2]]], parent.frame())
-
-  # allow
   dnds <- tstv <- 1
   dnds <- ifelse(is.na(existing[3]), 1,
     eval(extras[[existing[3]]], parent.frame()))
   tstv <- ifelse(is.na(existing[4]), 1,
     eval(extras[[existing[4]]], parent.frame()))
-
+  w <- g <- NULL
+  if(!is.na(existing[5]) & !is.na(existing[6])) {
+    w <- eval(extras[[existing[5]]], parent.frame())
+    g <- eval(extras[[existing[6]]], parent.frame())
+    stopifnot(length(g) == length(w))
+    k <- length(w)
+  }
   if (!inherits(tree, "phylo")) stop("tree must be of class phylo")
   if (is.null(tree$edge.length)) stop("tree must have edge weights")
   if (any(duplicated(tree$tip.label))) stop("tree must have unique labels!")
@@ -1419,22 +1424,23 @@ pml <- function(tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
     Q <- rep(1, length(levels) * (length(levels) - 1) / 2)
   m <- 1
   eig <- edQt(bf = bf, Q = Q)
-  if(site.rate=="free_rate"){
-    w <- rep(1/k, k)
-    g <- rep(1, k)
-  }
-  else{
-    rw <- rates_n_weights(shape, k, site.rate)
-    g <- rw[, 1]
-    w <- rw[, 2]
-  }
-  if (inv > 0){
-    w <- (1 - inv) * w
-    g <- g / (1 - inv)
-  }
+  if(is.null(g) | is.null(w)){
+    if(site.rate=="free_rate"){
+      w <- rep(1/k, k)
+      g <- rep(1, k)
+    }
+    else{
+      rw <- rates_n_weights(shape, k, site.rate)
+      g <- rw[, 1]
+      w <- rw[, 2]
+    }
+    if (inv > 0){
+      w <- (1 - inv) * w
+      g <- g / (1 - inv)
+    }
 #  if (wMix > 0) w <- (1 - wMix) * w
-  g <- g * rate
-
+    g <- g * rate
+  }
   inv0 <- inv
   kmax <- k
   if (any(g < .gEps)) {
@@ -1460,7 +1466,7 @@ pml <- function(tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
   tmp <- pml.fit(tree, data, bf, shape = shape, k = k, Q = Q,
     levels = attr(data, "levels"), inv = inv, rate = rate, g = g, w = w,
     eig = eig, INV = INV, ll.0 = ll.0, llMix = llMix, wMix = wMix, site = TRUE,
-    ASC=ASC)
+    ASC=ASC, site.rate = site.rate)
   df <- ifelse(is.ultrametric(tree), tree$Nnode, length(tree$edge.length))
 
   df <- switch(type,
diff --git a/R/plotAnc.R b/R/plotAnc.R
index 9ea7fca4..c0edf625 100644
--- a/R/plotAnc.R
+++ b/R/plotAnc.R
@@ -65,7 +65,7 @@ getAncDF <- function(x){
 #' example(NJ)
 #' # generate node labels to ensure plotting will work
 #' tree <- makeNodeLabel(tree)
-#' anc.p <- ancestral.pars(tree, Laurasiatherian)
+#' anc.p <- anc_pars(tree, Laurasiatherian)
 #' # plot the third character
 #' plotAnc(anc.p, 3, pos="topright")
 #' plotSeqLogo(anc.p, node="Node10", 1, 25)
@@ -73,7 +73,7 @@ getAncDF <- function(x){
 #' data(chloroplast)
 #' tree <- pratchet(chloroplast,  maxit=10, trace=0)
 #' tree <- makeNodeLabel(tree)
-#' anc.ch <- ancestral.pars(tree, chloroplast)
+#' anc.ch <- anc_pars(tree, chloroplast)
 #' image(as.phyDat(anc.ch)[, 1:25])
 #' plotAnc(anc.ch, 21, scheme="Ape_AA", pos="topleft")
 #' plotAnc(anc.ch, 21, scheme="Clustal", pos="topleft")
@@ -88,7 +88,7 @@ plotAnc <- function(x, i = 1, type="phylogram", ..., col = NULL,
   type <- match.arg(type, c("phylogram", "cladogram", "fan", "unrooted",
                             "radial", "tidy"))
   phylo_clado <- type %in% c("phylogram", "cladogram")
-  df <- getAncDF(x)
+  df <- as.data.frame(x)
   data <- x$data
   tree <- x$tree
   subset <- df[,"Site"] == i
@@ -158,32 +158,36 @@ plotAnc <- function(x, i = 1, type="phylogram", ..., col = NULL,
 my_ggseqlogo <-function (data, facet = "wrap", scales = "free_x", ncol = NULL,
           nrow = NULL, start=NULL, end=NULL, ...)
 {
-  x <- geom_logo(data = data, ...)
-  x[[2]] <- scale_x_continuous(limits = c(start-0.5, end+.5) ,
+  x <- ggseqlogo::geom_logo(data = data, ...)
+  x[[2]] <- ggplot2::scale_x_continuous(limits = c(start-0.5, end+.5) ,
                                breaks=pretty(seq(start, end)))
-  p <- ggplot() + x + theme_logo()
+  p <- ggplot2::ggplot() + x + ggseqlogo::theme_logo()
   if (!"list" %in% class(data)) return(p)
   facet <- match.arg(facet, c("grid", "wrap"))
   if (facet == "grid") {
-    p <- p + facet_grid(~seq_group, scales = scales)
+    p <- p + ggplot2::facet_grid(~seq_group, scales = scales)
   }
   else if (facet == "wrap") {
-    p <- p + facet_wrap(~seq_group, scales = scales, nrow = nrow, ncol = ncol)
+    p <- p + ggplot2::facet_wrap(~seq_group, scales = scales, nrow = nrow,
+                                 ncol = ncol)
   }
   return(p)
 }
 
 
 #' @rdname plot.ancestral
-#' @importFrom ggplot2 scale_x_continuous ggplot facet_grid facet_wrap
-#' @importFrom ggseqlogo geom_logo theme_logo
+## @importFrom ggplot2 scale_x_continuous ggplot facet_grid facet_wrap
+## @importFrom ggseqlogo geom_logo theme_logo
 #' @returns \code{plotSeqLogo} returns a ggplot object.
 #' @export
-plotSeqLogo <- function(x, node=getRoot(x$tree), start=1, end=10, scheme="Ape_NT", ...){
+plotSeqLogo <- function(x, node=getRoot(x$tree), start=1, end=10,
+                        scheme="Ape_NT", ...){
   stopifnot(inherits(x, "ancestral"))
+  chk <- requireNamespace("ggseqlogo", quietly = TRUE)
+  if (!chk) stop("package ggseqlogo needs to be not installed!\n")
   type <- attr(x$data, "type")
   tree <- x$tree
-  df <- getAncDF(x)
+  df <- as.data.frame(x)
   nodes <- c(tree$tip.label, tree$node.label)
   if(is.numeric(node)) node <- nodes[node]
   subset <- df[,"Node"] == node
diff --git a/R/plotBS.R b/R/plotBS.R
index 5588f297..1d2ce7d7 100644
--- a/R/plotBS.R
+++ b/R/plotBS.R
@@ -75,10 +75,10 @@ support <- function(tree, trees, method="FBP", tol=2e-8, scale=TRUE){
 #' \code{trees} is optional and if not supplied the labels supplied
 #' in the \code{node.label} slot will be used.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso  \code{\link{plot.phylo}}, \code{\link{add_ci}},
-#' \code{\link{nodelabels}},
-#' \code{\link{prop.clades}}, \code{\link{maxCladeCred}},
-#' \code{\link{transferBootstrap}}, \code{\link{consensus}},
+#' @seealso  \code{\link[ape]{plot.phylo}}, \code{\link{add_ci}},
+#' \code{\link[ape]{nodelabels}},
+#' \code{\link[ape]{prop.clades}}, \code{\link{maxCladeCred}},
+#' \code{\link{transferBootstrap}}, \code{\link[ape]{consensus}},
 #' \code{\link{consensusNet}}
 #' @references Felsenstein J. (1985) Confidence limits on phylogenies. An
 #' approach using the bootstrap. \emph{Evolution} \bold{39}, 783--791
diff --git a/R/plot_networx.R b/R/plot_networx.R
index 70de35f9..a48d80a1 100644
--- a/R/plot_networx.R
+++ b/R/plot_networx.R
@@ -325,7 +325,8 @@ plot.networx <- function(x, type = "equal angle", use.edge.length = TRUE,
     plot2D(coord, x, show.tip.label = show.tip.label,
            show.edge.label = show.edge.label, edge.label = edge.label,
            show.node.label = show.node.label, node.label = node.label,
-           show.nodes = show.nodes, tip.color = tip.color, edge.color = edge.color,
+           #show.nodes = show.nodes,
+           tip.color = tip.color, edge.color = edge.color,
            edge.width = edge.width, edge.lty = edge.lty, font = font, cex = cex,
            cex.node.label = cex.node.label, cex.edge.label = cex.edge.label,
            col.node.label = col.node.label, col.edge.label = col.edge.label,
@@ -408,24 +409,27 @@ plot2D <- function(coords, net, show.tip.label = TRUE, show.edge.label = FALSE,
   xx <- coords[, 1]
   yy <- coords[, 2]
   nTips <- length(label)
+  offset <- max(nchar(label)) * 0.018 * cex * diff(range(xx))
   if(is.null(xlim)){
     xlim <- range(xx)
-    offset <- max(nchar(label)) * 0.018 * cex * diff(xlim)
+#    offset <- max(nchar(label)) * 0.018 * cex * diff(xlim)
     if (show.tip.label) xlim <- c(xlim[1] - offset, xlim[2] + offset)
   }
   if(is.null(ylim)){
     ylim <- range(yy)
     if (show.tip.label){
       if(direction=="axial"){
-        offset <- max(nchar(label)) * 0.018 * cex * diff(ylim)
+#        offset <- max(nchar(label)) * 0.018 * cex * diff(ylim)
         ylim <- c(ylim[1] - offset, ylim[2] + offset)
       } else ylim <- c(ylim[1] - 0.03 * cex * diff(ylim),
                      ylim[2] + 0.03 * cex * diff(ylim))
     }
   }
   if (!add) {
-    plot.new()
-    plot.window(xlim, ylim, asp = 1)
+#    plot.new()
+#    plot.window(xlim, ylim, asp = 1, ...)
+    plot.default(0, type = "n", xlim = xlim, ylim = ylim, xlab = "",
+                 ylab = "", axes = FALSE, asp = 1, ...)
   }
   cladogram.plot(edge, xx, yy, edge.color, edge.width, edge.lty)
   if (show.tip.label) {
diff --git a/R/plot_pml.R b/R/plot_pml.R
index 553eb864..66507283 100644
--- a/R/plot_pml.R
+++ b/R/plot_pml.R
@@ -14,8 +14,8 @@
 #' @return \code{plot.pml} returns invisibly a list with arguments dexcribing the plot.
 #' For further details see the \code{plot.phylo}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{plot.phylo}}, \code{\link{axisPhylo}},
-#' \code{\link{add.scale.bar}}
+#' @seealso \code{\link[ape]{plot.phylo}}, \code{\link[ape]{axisPhylo}},
+#' \code{\link[ape]{add.scale.bar}}
 #' @keywords plot
 #' @examples
 #' fdir <- system.file("extdata/trees", package = "phangorn")
diff --git a/R/pml_bb.R b/R/pml_bb.R
index 9e2d1a7a..aeb28734 100644
--- a/R/pml_bb.R
+++ b/R/pml_bb.R
@@ -36,8 +36,8 @@
 #' @param \dots Further arguments passed to or from other methods.
 #' @return \code{pml_bb} returns an object of class pml.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{optim.pml}}, \code{\link{modelTest}}, \code{\link{rtt}},
-#' \code{\link{pml.control}}
+#' @seealso \code{\link{optim.pml}}, \code{\link{modelTest}},
+#' \code{\link[ape]{rtt}}, \code{\link{pml.control}}
 #' @keywords cluster
 #' @examples
 #'
diff --git a/R/read.nexus.partitions.R b/R/read.nexus.partitions.R
index b23ada8d..f68e1442 100644
--- a/R/read.nexus.partitions.R
+++ b/R/read.nexus.partitions.R
@@ -48,7 +48,7 @@ read.nexus.charset <- function(file){
 #' @return a list where each element is a 'phyDat' object or an object of class
 #' 'multiphyDat'.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{read.nexus.data}}, \code{\link{read.phyDat}}
+#' @seealso \code{\link[ape]{read.nexus.data}}, \code{\link{read.phyDat}}
 #' @keywords cluster
 #' @examples
 #' tree <- rtree(10)
diff --git a/R/read.nexus.splits.R b/R/read.nexus.splits.R
index 9f0e9f57..528abb4b 100644
--- a/R/read.nexus.splits.R
+++ b/R/read.nexus.splits.R
@@ -30,7 +30,7 @@
 #' assumes that different co-variables are tab delimited and the bipartition
 #' are separated with white-space. Comments in square brackets are ignored.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{prop.part}}, \code{\link{lento}},
+#' @seealso \code{\link[ape]{prop.part}}, \code{\link{lento}},
 #' \code{\link{as.splits}}, \code{\link{as.networx}}
 #' @keywords cluster
 #' @examples
diff --git a/R/splits.R b/R/splits.R
index 88d7abc9..c1fefc63 100644
--- a/R/splits.R
+++ b/R/splits.R
@@ -26,7 +26,7 @@
 #' two splits are incompatible.
 #' @note The internal representation is likely to change.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{prop.part}}, \code{\link{lento}},
+#' @seealso \code{\link[ape]{prop.part}}, \code{\link{lento}},
 #' \code{\link{as.networx}}, \code{\link{distanceHadamard}},
 #' \code{\link{read.nexus.splits}}
 #' @keywords cluster
diff --git a/R/transferBootstrap.R b/R/transferBootstrap.R
index e01533bb..9675ddbf 100644
--- a/R/transferBootstrap.R
+++ b/R/transferBootstrap.R
@@ -13,7 +13,7 @@
 #' the node labels.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso  \code{\link{plotBS}}, \code{\link{maxCladeCred}},
-#' \code{\link{drawSupportOnEdges}}
+#' \code{\link[ape]{drawSupportOnEdges}}
 #' @references Lemoine, F., Entfellner, J. B. D., Wilkinson, E., Correia, D.,
 #' Felipe, M. D., De Oliveira, T., & Gascuel, O. (2018). Renewing Felsenstein’s
 #' phylogenetic bootstrap in the era of big data. \emph{Nature},
diff --git a/R/treeManipulation.R b/R/treeManipulation.R
index fd0e4837..a98541f2 100644
--- a/R/treeManipulation.R
+++ b/R/treeManipulation.R
@@ -430,7 +430,7 @@ addOneTree <- function(tree, subtree, i, node) {
 #' @param edge.length optional numeric vector with edge length
 #' @return an object of class phylo
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{bind.tree}}
+#' @seealso \code{\link[ape]{bind.tree}}
 #' @keywords cluster
 #' @examples
 #' tree <- rcoal(10)
diff --git a/R/upgma.R b/R/upgma.R
index 7ac5275d..5607696c 100644
--- a/R/upgma.R
+++ b/R/upgma.R
@@ -20,8 +20,8 @@
 #' @return A phylogenetic tree of class \code{phylo}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{hclust}}, \code{\link{dist.hamming}}, \code{\link{NJ}},
-#' \code{\link{as.phylo}}, \code{\link{fastme}}, \code{\link{nnls.tree}},
-#' \code{\link{rtt}}
+#' \code{\link[ape]{as.phylo}}, \code{\link[ape]{fastme}},
+#' \code{\link{nnls.tree}}, \code{\link[ape]{rtt}}
 #' @references Sneath, P. H., & Sokal, R. R. (1973). \emph{Numerical taxonomy.
 #' The principles and practice of numerical classification.}
 #'
diff --git a/inst/tinytest/test_ancestral.R b/inst/tinytest/test_ancestral.R
index bd01eb41..b0bc2fb1 100644
--- a/inst/tinytest/test_ancestral.R
+++ b/inst/tinytest/test_ancestral.R
@@ -15,9 +15,9 @@ fit <- pml(tree, dna)
 # dna tests differs from other data types as it may returns ambiguous data
 # test ancestral generics
 #    test.ml1 <- ancestral.pml(fit, type = "ml")
-test_ml <- ancestral.pml(fit, type = "ml")
-test_mpr <- ancestral.pars(tree, dna, "MPR")
-test_acctran <- ancestral.pars(tree, dna, "ACCTRAN")
+#test_ml <- ancestral.pml(fit, type = "ml")
+#test_mpr <- ancestral.pars(tree, dna, "MPR")
+#test_acctran <- ancestral.pars(tree, dna, "ACCTRAN")
 
 #expect_equal(as.character(test_ml), as.character(test_acctran))
 #expect_equal(as.character(test_ml), as.character(test_mpr))
@@ -29,17 +29,18 @@ test_acctran <- ancestral.pars(tree, dna, "ACCTRAN")
 
 
 data(Laurasiatherian)
-fit <- pml_bb(Laurasiatherian[,1:100], "JC", rearrangement = "none")
-anc_ml <- ancestral.pml(fit)
+fit <- pml_bb(Laurasiatherian[,1:100], "JC", rearrangement = "none",
+              control=pml.control(trace=0))
+anc_ml <- anc_pml(fit)
 write.ancestral(anc_ml)
-align <- read.phyDat("ancestral_align.fasta")
+align <- read.phyDat("ancestral_align.fasta", format = "fasta")
 tree <- read.tree("ancestral_tree.nwk")
 df <- read.table("ancestral.state", header=TRUE)
-anc_ml_disc <- ancestral(tree, align, df)
-expect_equal(anc_ml[[1]], anc_ml_disc[[1]])
-expect_equal(anc_ml[[2]], anc_ml_disc[[2]])
-expect_equal(anc_ml[[3]], anc_ml_disc[[3]])
-expect_equal(anc_ml[[4]], anc_ml_disc[[4]])
+anc_ml_disc <- as.ancestral(tree, align, df)
+expect_equal(anc_ml$tree, anc_ml_disc$tree)
+expect_equal(anc_ml$data[tree$tip.label], anc_ml_disc$data)
+#expect_equal(anc_ml[[3]], anc_ml_disc[[3]])
+expect_equal(anc_ml$state, anc_ml_disc$state)
 unlink(c("ancestral_align.fasta", "ancestral_tree.nwk", "ancestral.state"))
 
 
diff --git a/man/add.tips.Rd b/man/add.tips.Rd
index 7689364f..bdd80ab3 100644
--- a/man/add.tips.Rd
+++ b/man/add.tips.Rd
@@ -31,7 +31,7 @@ tree1 <- add.tips(tree, c("A", "B", "C"), c(1,2,15))
 plot(tree1)
 }
 \seealso{
-\code{\link{bind.tree}}
+\code{\link[ape]{bind.tree}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/addConfidences.Rd b/man/addConfidences.Rd
index 06ba9911..29c31797 100644
--- a/man/addConfidences.Rd
+++ b/man/addConfidences.Rd
@@ -65,7 +65,7 @@ Schliep, K., Potts, A. J., Morrison, D. A. and Grimm, G. W.
 }
 \seealso{
 \code{\link{as.splits}}, \code{\link{as.networx}},
-\code{\link{RF.dist}}, \code{\link{plot.phylo}}
+\code{\link{RF.dist}}, \code{\link[ape]{plot.phylo}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/add_ci.Rd b/man/add_ci.Rd
index 9a3c8535..83cf11a5 100644
--- a/man/add_ci.Rd
+++ b/man/add_ci.Rd
@@ -52,7 +52,7 @@ plot(tree, direction="downwards")
 add_boxplot(tree, trees, boxwex=.7)
 }
 \seealso{
-\code{\link{plot.phylo}}, \code{\link{plotBS}},
+\code{\link[ape]{plot.phylo}}, \code{\link{plotBS}},
 \code{\link{add_edge_length}}, \code{\link{maxCladeCred}}
 }
 \author{
diff --git a/man/add_edge_length.Rd b/man/add_edge_length.Rd
index 424cba2a..cfc75e67 100644
--- a/man/add_edge_length.Rd
+++ b/man/add_edge_length.Rd
@@ -38,8 +38,9 @@ plot(tree_compat)
 add_boxplot(tree_compat, bs, boxwex=.7)
 }
 \seealso{
-\code{\link{node.depth.edgelength}}, \code{\link{consensus}},
-\code{\link{maxCladeCred}}, \code{\link{add_boxplot}}
+\code{\link[ape]{node.depth}},
+\code{\link[ape]{consensus}}, \code{\link{maxCladeCred}},
+\code{\link{add_boxplot}}
 }
 \author{
 Klaus Schliep
diff --git a/man/ancestral.pml.Rd b/man/ancestral.pml.Rd
index b84f3a60..033b2afc 100644
--- a/man/ancestral.pml.Rd
+++ b/man/ancestral.pml.Rd
@@ -2,25 +2,38 @@
 % Please edit documentation in R/ancestral.R
 \name{ancestral.pml}
 \alias{ancestral.pml}
+\alias{anc_pml}
 \alias{as.phyDat.ancestral}
+\alias{as.data.frame.ancestral}
 \alias{ancestral.pars}
+\alias{anc_pars}
 \alias{pace}
 \title{Ancestral character reconstruction.}
 \usage{
-ancestral.pml(object, type = "marginal", ...)
+ancestral.pml(object, type = "marginal", return = "prob", ...)
+
+anc_pml(object, type = "marginal", ...)
 
 \method{as.phyDat}{ancestral}(x, ...)
 
+\method{as.data.frame}{ancestral}(x, ...)
+
 ancestral.pars(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"),
-  cost = NULL, ...)
+  cost = NULL, return = "prob")
 
-pace(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"), cost = NULL, ...)
+anc_pars(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"), cost = NULL,
+  ...)
+
+pace(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"), cost = NULL,
+  return = "prob")
 }
 \arguments{
 \item{object}{an object of class pml}
 
 \item{type}{method used to assign characters to internal nodes, see details.}
 
+\item{return}{return a \code{phyDat} object or matrix of probabilities.}
+
 \item{\dots}{Further arguments passed to or from other methods.}
 
 \item{x}{an object of class ancestral}
@@ -69,8 +82,8 @@ example(NJ)
 # generate node labels to ensure plotting will work
 tree <- makeNodeLabel(tree)
 fit <- pml(tree, Laurasiatherian)
-anc.ml <- ancestral.pml(fit, type = "ml")
-anc.p <- ancestral.pars(tree, Laurasiatherian)
+anc.ml <- anc_pml(fit)
+anc.p <- anc_pars(tree, Laurasiatherian)
 # plot ancestral sequences at the root
 plotSeqLogo( anc.ml, 48, 1, 20)
 plotSeqLogo( anc.p, 48, 1, 20)
@@ -92,7 +105,7 @@ Press, Oxford.
 }
 \seealso{
 \code{\link{pml}}, \code{\link{parsimony}}, \code{\link[ape]{ace}},
-\code{\link{plotAnc}}, \code{\link{ltg2amb}}, \code{\link{latag2n}},
+\code{\link{plotAnc}}, \code{\link{ltg2amb}}, \code{\link[ape]{latag2n}},
 \code{\link{gap_as_state}}, \code{\link[ape]{root}},
 \code{\link[ape]{makeNodeLabel}}
 }
diff --git a/man/as.phyDat.Rd b/man/as.phyDat.Rd
index afba3199..6bc1412d 100644
--- a/man/as.phyDat.Rd
+++ b/man/as.phyDat.Rd
@@ -7,6 +7,8 @@
 \alias{as.phyDat.matrix}
 \alias{as.MultipleAlignment}
 \alias{as.MultipleAlignment.phyDat}
+\alias{as.StringSet}
+\alias{as.StringSet.phyDat}
 \alias{acgt2ry}
 \alias{phyDat2MultipleAlignment}
 \alias{as.phyDat}
@@ -45,6 +47,8 @@ phyDat2alignment(x)
 
 \method{as.phyDat}{DNAStringSet}(x, ...)
 
+\method{as.StringSet}{phyDat}(x, ...)
+
 \method{as.MultipleAlignment}{phyDat}(x, ...)
 
 \method{as.character}{phyDat}(x, allLevels = TRUE, ...)
@@ -111,9 +115,9 @@ all.equal(Laurasiatherian, as.phyDat(LauraDNAbin))
 
 }
 \seealso{
-[DNAbin()], [as.DNAbin()],
+\code{\link[ape]{DNAbin}}, \code{\link[ape]{as.DNAbin}},
 \code{\link{baseFreq}}, \code{\link{glance.phyDat}},
-\code{\link{read.dna}}, \code{\link{read.nexus.data}}
+\code{\link[ape]{read.dna}}, \code{\link[ape]{read.nexus.data}}
 and the chapter 1 in the \code{vignette("phangorn-specials",
 package="phangorn")} and the example of \code{\link{pmlMix}} for the use of
 \code{allSitePattern}
diff --git a/man/as.splits.Rd b/man/as.splits.Rd
index 4f82b6b5..bccb480e 100644
--- a/man/as.splits.Rd
+++ b/man/as.splits.Rd
@@ -105,7 +105,7 @@ plot(as.networx(spl))
 
 }
 \seealso{
-\code{\link{prop.part}}, \code{\link{lento}},
+\code{\link[ape]{prop.part}}, \code{\link{lento}},
 \code{\link{as.networx}}, \code{\link{distanceHadamard}},
 \code{\link{read.nexus.splits}}
 }
diff --git a/man/baseFreq.Rd b/man/baseFreq.Rd
index 5ba3339e..822796a1 100644
--- a/man/baseFreq.Rd
+++ b/man/baseFreq.Rd
@@ -51,7 +51,7 @@ glance(chloroplast)
 composition_test(Laurasiatherian)[1:10,]
 }
 \seealso{
-\code{\link{phyDat}, \link{base.freq}, \link{glance}}
+\code{\link{phyDat}, \link[ape]{base.freq}, \link{glance}}
 }
 \author{
 Klaus Schliep
diff --git a/man/bootstrap.pml.Rd b/man/bootstrap.pml.Rd
index 1616bc4e..71613afc 100644
--- a/man/bootstrap.pml.Rd
+++ b/man/bootstrap.pml.Rd
@@ -100,8 +100,8 @@ trees. \emph{Molecular Biology and Evolution} \bold{3}, 403--417
 }
 \seealso{
 \code{\link{optim.pml}}, \code{\link{pml}},
-\code{\link{plot.phylo}}, \code{\link{maxCladeCred}}
-\code{\link{nodelabels}},\code{\link{consensusNet}} and
+\code{\link[ape]{plot.phylo}}, \code{\link{maxCladeCred}}
+\code{\link[ape]{nodelabels}},\code{\link{consensusNet}} and
 \code{\link{SOWH.test}} for parametric bootstrap
 }
 \author{
diff --git a/man/cladePar.Rd b/man/cladePar.Rd
index 0fc2d612..fbfacca8 100644
--- a/man/cladePar.Rd
+++ b/man/cladePar.Rd
@@ -42,7 +42,7 @@ cladePar(tree, 18, "blue", "blue", x=x, plot=TRUE)
 
 }
 \seealso{
-\code{\link{plot.phylo}}
+\code{\link[ape]{plot.phylo}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/coalSpeciesTree.Rd b/man/coalSpeciesTree.Rd
index 878be9ed..0af0610b 100644
--- a/man/coalSpeciesTree.Rd
+++ b/man/coalSpeciesTree.Rd
@@ -40,7 +40,7 @@ consistent estimator of the species tree. \emph{Journal of Mathematical
 Biology}, \bold{60}, 95--106.
 }
 \seealso{
-\code{\link{speciesTree}}
+\code{\link[ape]{speciesTree}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com} Emmanuel Paradies
diff --git a/man/densiTree.Rd b/man/densiTree.Rd
index a6d88b1a..7b4e63ad 100644
--- a/man/densiTree.Rd
+++ b/man/densiTree.Rd
@@ -117,8 +117,8 @@ Remco R. Bouckaert (2010) DensiTree: making sense of sets of phylogenetic
 trees \emph{Bioinformatics}, \bold{26 (10)}, 1372-1373.
 }
 \seealso{
-\code{\link{plot.phylo}}, \code{\link{plot.networx}},
-\code{\link{jitter}}, \code{\link{rtt}}
+\code{\link[ape]{plot.phylo}}, \code{\link{plot.networx}},
+\code{\link{jitter}}, \code{\link[ape]{rtt}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/dna2codon.Rd b/man/dna2codon.Rd
index f9077f4c..88e6c958 100644
--- a/man/dna2codon.Rd
+++ b/man/dna2codon.Rd
@@ -75,7 +75,7 @@ dna2codon(Laurasiatherian)
 \url{https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=cgencodes}
 }
 \seealso{
-\code{\link{trans}}, \code{\link{phyDat}} and the chapter 4 in the
+\code{\link[ape]{trans}}, \code{\link{phyDat}} and the chapter 4 in the
 \code{vignette("phangorn-specials", package="phangorn")}
 }
 \author{
diff --git a/man/gap_as_state.Rd b/man/gap_as_state.Rd
index 0cb61e83..76148104 100644
--- a/man/gap_as_state.Rd
+++ b/man/gap_as_state.Rd
@@ -36,8 +36,9 @@ rownames(contr) <- attr(tmp, "allLevels")
 contr
 }
 \seealso{
-\code{\link{phyDat}}, \code{\link{ltg2amb}}, \code{\link{latag2n}},
-\code{\link{ancestral.pml}}, \code{\link{gap_as_state}}
+\code{\link{phyDat}}, \code{\link{ltg2amb}},
+\code{\link[ape]{latag2n}}, \code{\link{ancestral.pml}},
+\code{\link{gap_as_state}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/image.phyDat.Rd b/man/image.phyDat.Rd
index 637c9e1b..3e7b1cd6 100644
--- a/man/image.phyDat.Rd
+++ b/man/image.phyDat.Rd
@@ -26,5 +26,5 @@ data("chloroplast")
 image(chloroplast[, 1:50], scheme="Clustal", show.aa = TRUE)
 }
 \seealso{
-\code{\link{image.DNAbin}}, \code{\link{image.AAbin}}
+\code{\link[ape]{image.DNAbin}}, \code{\link[ape]{image.AAbin}}
 }
diff --git a/man/ltg2amb.Rd b/man/ltg2amb.Rd
index 4e7b9655..63e9dc1e 100644
--- a/man/ltg2amb.Rd
+++ b/man/ltg2amb.Rd
@@ -28,7 +28,7 @@ image(x)
 image(y)
 }
 \seealso{
-\code{\link{latag2n}}, \code{\link{ancestral.pml}},
+\code{\link[ape]{latag2n}}, \code{\link{ancestral.pml}},
 \code{\link{gap_as_state}}
 }
 \keyword{cluster}
diff --git a/man/maxCladeCred.Rd b/man/maxCladeCred.Rd
index 1de7ead8..82b07ce9 100644
--- a/man/maxCladeCred.Rd
+++ b/man/maxCladeCred.Rd
@@ -81,10 +81,10 @@ maxCladeCred(c(tree, bs[[1]]), tree=FALSE, part = pp)
 
 }
 \seealso{
-\code{\link{consensus}}, \code{\link{consensusNet}},
-\code{\link{prop.part}}, \code{\link{bootstrap.pml}}, \code{\link{plotBS}},
-\code{\link{transferBootstrap}}, \code{\link{add_edge_length}},
-\code{\link{add_boxplot}}
+\code{\link[ape]{consensus}}, \code{\link{consensusNet}},
+\code{\link[ape]{prop.part}}, \code{\link{bootstrap.pml}},
+\code{\link{plotBS}}, \code{\link{transferBootstrap}},
+\code{\link{add_edge_length}}, \code{\link{add_boxplot}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/plot.ancestral.Rd b/man/plot.ancestral.Rd
index e6a5bcf2..9b699bf5 100644
--- a/man/plot.ancestral.Rd
+++ b/man/plot.ancestral.Rd
@@ -57,7 +57,7 @@ For further details see vignette("Ancestral").
 example(NJ)
 # generate node labels to ensure plotting will work
 tree <- makeNodeLabel(tree)
-anc.p <- ancestral.pars(tree, Laurasiatherian)
+anc.p <- anc_pars(tree, Laurasiatherian)
 # plot the third character
 plotAnc(anc.p, 3, pos="topright")
 plotSeqLogo(anc.p, node="Node10", 1, 25)
@@ -65,7 +65,7 @@ plotSeqLogo(anc.p, node="Node10", 1, 25)
 data(chloroplast)
 tree <- pratchet(chloroplast,  maxit=10, trace=0)
 tree <- makeNodeLabel(tree)
-anc.ch <- ancestral.pars(tree, chloroplast)
+anc.ch <- anc_pars(tree, chloroplast)
 image(as.phyDat(anc.ch)[, 1:25])
 plotAnc(anc.ch, 21, scheme="Ape_AA", pos="topleft")
 plotAnc(anc.ch, 21, scheme="Clustal", pos="topleft")
diff --git a/man/plot.pml.Rd b/man/plot.pml.Rd
index 35a212e5..ed0a8f57 100644
--- a/man/plot.pml.Rd
+++ b/man/plot.pml.Rd
@@ -48,8 +48,8 @@ plot(fit_unrooted, cex=.5)
 
 }
 \seealso{
-\code{\link{plot.phylo}}, \code{\link{axisPhylo}},
-\code{\link{add.scale.bar}}
+\code{\link[ape]{plot.phylo}}, \code{\link[ape]{axisPhylo}},
+\code{\link[ape]{add.scale.bar}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/plotBS.Rd b/man/plotBS.Rd
index 04d50ea2..4f80e832 100644
--- a/man/plotBS.Rd
+++ b/man/plotBS.Rd
@@ -91,10 +91,10 @@ Penny D. and Hendy M.D. (1986) Estimating the reliability of evolutionary
 trees. \emph{Molecular Biology and Evolution} \bold{3}, 403--417
 }
 \seealso{
-\code{\link{plot.phylo}}, \code{\link{add_ci}},
-\code{\link{nodelabels}},
-\code{\link{prop.clades}}, \code{\link{maxCladeCred}},
-\code{\link{transferBootstrap}}, \code{\link{consensus}},
+\code{\link[ape]{plot.phylo}}, \code{\link{add_ci}},
+\code{\link[ape]{nodelabels}},
+\code{\link[ape]{prop.clades}}, \code{\link{maxCladeCred}},
+\code{\link{transferBootstrap}}, \code{\link[ape]{consensus}},
 \code{\link{consensusNet}}
 }
 \author{
diff --git a/man/pml.Rd b/man/pml.Rd
index 9d79876b..451b3319 100644
--- a/man/pml.Rd
+++ b/man/pml.Rd
@@ -232,7 +232,8 @@ discrete morphological character data. \emph{Systematic Biology} \bold{50},
 \seealso{
 \code{\link{pml_bb}}, \code{\link{bootstrap.pml}},
 \code{\link{modelTest}}, \code{\link{pmlPart}}, \code{\link{pmlMix}},
-\code{\link{plot.phylo}}, \code{\link{SH.test}}, \code{\link{ancestral.pml}}
+\code{\link[ape]{plot.phylo}}, \code{\link{SH.test}},
+\code{\link{ancestral.pml}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/pml_bb.Rd b/man/pml_bb.Rd
index 333386f2..fb79d66e 100644
--- a/man/pml_bb.Rd
+++ b/man/pml_bb.Rd
@@ -69,8 +69,8 @@ fit_HKY_R2 <- pml_bb(woodmouse, model="HKY+R(2)")
 }
 }
 \seealso{
-\code{\link{optim.pml}}, \code{\link{modelTest}}, \code{\link{rtt}},
-\code{\link{pml.control}}
+\code{\link{optim.pml}}, \code{\link{modelTest}},
+\code{\link[ape]{rtt}}, \code{\link{pml.control}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/read.nexus.partitions.Rd b/man/read.nexus.partitions.Rd
index 468746c5..2000a44c 100644
--- a/man/read.nexus.partitions.Rd
+++ b/man/read.nexus.partitions.Rd
@@ -42,7 +42,7 @@ tmp
 unlink(zz)
 }
 \seealso{
-\code{\link{read.nexus.data}}, \code{\link{read.phyDat}}
+\code{\link[ape]{read.nexus.data}}, \code{\link{read.phyDat}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/read.nexus.splits.Rd b/man/read.nexus.splits.Rd
index 5b9d9660..87c75b78 100644
--- a/man/read.nexus.splits.Rd
+++ b/man/read.nexus.splits.Rd
@@ -75,7 +75,7 @@ write.splits(spl, print.labels = FALSE)
 
 }
 \seealso{
-\code{\link{prop.part}}, \code{\link{lento}},
+\code{\link[ape]{prop.part}}, \code{\link{lento}},
 \code{\link{as.splits}}, \code{\link{as.networx}}
 }
 \author{
diff --git a/man/transferBootstrap.Rd b/man/transferBootstrap.Rd
index fb714273..71989c32 100644
--- a/man/transferBootstrap.Rd
+++ b/man/transferBootstrap.Rd
@@ -45,7 +45,7 @@ phylogenetic bootstrap in the era of big data. \emph{Nature},
 }
 \seealso{
 \code{\link{plotBS}}, \code{\link{maxCladeCred}},
-\code{\link{drawSupportOnEdges}}
+\code{\link[ape]{drawSupportOnEdges}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/upgma.Rd b/man/upgma.Rd
index ddf6be05..94d94f7c 100644
--- a/man/upgma.Rd
+++ b/man/upgma.Rd
@@ -56,8 +56,8 @@ samples under the assumption of a molecular clock using serial-sample UPGMA.
 }
 \seealso{
 \code{\link{hclust}}, \code{\link{dist.hamming}}, \code{\link{NJ}},
-\code{\link{as.phylo}}, \code{\link{fastme}}, \code{\link{nnls.tree}},
-\code{\link{rtt}}
+\code{\link[ape]{as.phylo}}, \code{\link[ape]{fastme}},
+\code{\link{nnls.tree}}, \code{\link[ape]{rtt}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/write.ancestral.Rd b/man/write.ancestral.Rd
index 923ee026..2a0daa7c 100644
--- a/man/write.ancestral.Rd
+++ b/man/write.ancestral.Rd
@@ -2,13 +2,13 @@
 % Please edit documentation in R/ancestral.R
 \name{write.ancestral}
 \alias{write.ancestral}
-\alias{ancestral}
+\alias{as.ancestral}
 \alias{print.ancestral}
 \title{Export and convenience functions for  ancestral reconstructions}
 \usage{
 write.ancestral(x, file = "ancestral")
 
-ancestral(tree, align, prob)
+as.ancestral(tree, align, prob)
 
 \method{print}{ancestral}(x, ...)
 }
@@ -41,13 +41,13 @@ plotting capabilities of R.
 \examples{
 data(Laurasiatherian)
 fit <- pml_bb(Laurasiatherian[,1:100], "JC", rearrangement = "none")
-anc_ml <- ancestral.pml(fit)
+anc_ml <- anc_pml(fit)
 write.ancestral(anc_ml)
 # Can be also results from iqtree
-align <- read.phyDat("ancestral_align.fasta")
+align <- read.phyDat("ancestral_align.fasta", format="fasta")
 tree <- read.tree("ancestral_tree.nwk")
 df <- read.table("ancestral.state", header=TRUE)
-anc_ml_disc <- ancestral(tree, align, df)
+anc_ml_disc <- as.ancestral(tree, align, df)
 plotAnc(anc_ml_disc, 20)
 unlink(c("ancestral_align.fasta", "ancestral_tree.nwk", "ancestral.state"))
 }
diff --git a/src/dist.c b/src/dist.c
index e6d7c770..08a061ba 100644
--- a/src/dist.c
+++ b/src/dist.c
@@ -1,7 +1,7 @@
 /*
  * dist.c
  *
- * (c) 2008-2021 Klaus Schliep (klaus.schliep@gmail.com)
+ * (c) 2008-2024 Klaus Schliep (klaus.schliep@gmail.com)
  *
  *
  * This code may be distributed under the GNU GPL
diff --git a/src/fitch64.cpp b/src/fitch64.cpp
index 2b20bf2b..cf71a805 100644
--- a/src/fitch64.cpp
+++ b/src/fitch64.cpp
@@ -846,6 +846,7 @@ double pscore(Fitch* obj, const IntegerMatrix & orig){
 
 NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
   int i,j;
+  size_t ij;
   int states = obj->nStates;
   int nBits = obj->nBits;
   std::vector< std::vector<uint64_t> > vector = obj->X;
@@ -867,6 +868,7 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
   uint64_t ones = ~0ull;
   uint64_t tmp = 0ull;
   for(int k=0; k<nl; k+=2){
+    ij = (size_t) anc[k] - 1L;
     child1 = vector[desc[k] - 1L].data();
     child2 = vector[desc[k+1] - 1L].data();
     parent = vector[anc[k] - 1L].data();
@@ -884,7 +886,7 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
       parent += states;
       tmp = ~orvand & ones;
       for(int l=0; l<BIT_SIZE; ++l){
-        if((tmp >> l) & 1ull) pars[anc[k] - 1L]+= obj->weight[i*BIT_SIZE + l];
+        if((tmp >> l) & 1ull) pars[ij]+= obj->weight[i*BIT_SIZE + l];
       }
     }
     for (int i = obj->wBits; i < nBits; ++i){
@@ -897,12 +899,13 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
       child2 += states;
       parent += states;
       tmp = ~orvand & ones;
-      pars[anc[k] - 1L] += popcnt64(tmp);
+      pars[ij] += popcnt64(tmp);
     }
   }
   if(unrooted){
     child1 = vector[desc[nl] - 1].data();
     parent = vector[anc[nl] - 1].data();
+    ij = (size_t) anc[nl] - 1L;
     for (i = 0; i < obj->wBits; ++i){
       // OR the ANDs of all states
       uint64_t orvand = 0ull;
@@ -915,7 +918,7 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
       parent += states;
       uint64_t tmp = ~orvand & ones;
       for(int l=0; l<BIT_SIZE; ++l){
-        if((tmp >> l) & 1ull) pars[anc[nl] - 1L]+= obj->weight[i*BIT_SIZE + l];
+        if((tmp >> l) & 1ull) pars[ij]+= obj->weight[i*BIT_SIZE + l];
       }
     }
     for (int i = obj->wBits; i < nBits; ++i){
@@ -929,7 +932,7 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
       child1 += states;
       parent += states;
       uint64_t tmp = ~orvand & ones;
-      pars [anc[nl] - 1L]+= popcnt64(tmp);
+      pars [ij]+= popcnt64(tmp);
     }
   }
   return(pars);
diff --git a/src/ml.c b/src/ml.c
index e7678e04..7ed105d7 100644
--- a/src/ml.c
+++ b/src/ml.c
@@ -1,7 +1,7 @@
 /*
  * ml.c
  *
- * (c) 2008-2021  Klaus Schliep (klaus.schliep@gmail.com)
+ * (c) 2008-2024  Klaus Schliep (klaus.schliep@gmail.com)
  *
  *
  * This code may be distributed under the GNU GPL
diff --git a/src/phangorn_utils.cpp b/src/phangorn_utils.cpp
index de7e80d0..d8f3dde0 100644
--- a/src/phangorn_utils.cpp
+++ b/src/phangorn_utils.cpp
@@ -88,14 +88,14 @@ int countCycle_cpp(IntegerMatrix M){
 // [[Rcpp::export]]
 IntegerVector countCycle2_cpp(IntegerMatrix M){
   int tmp;
-  int l = M.nrow();
-  int m = M.ncol();
+  size_t l = M.nrow();
+  size_t m = M.ncol();
   IntegerVector res(l);
-  for (int i=0; i<l; i++) {
+  for (size_t i=0; i<l; i++) {
     tmp = 0L;
     if(M[i] != M[i + (m -1) * l])tmp=1L;
-    for (int j=1; j<m; j++) {
-      if(M[i + (j-1L)* l] != M[i + j * l])tmp++;
+    for (size_t j=1; j<m; j++) {
+      if(M[i + (j-1L)* l] != M[i + j * l])tmp++;    ///
     }
     res[i]=tmp;
   }
@@ -395,12 +395,12 @@ Fast cophenetic distance
 */
 
 
-void copheneticHelpCpp(std::vector<int> left, std::vector<int> right, int h, NumericVector nh, int nTips, NumericVector dm){
+void copheneticHelpCpp(std::vector<int> left, std::vector<int> right, size_t h, NumericVector nh, int nTips, NumericVector dm){
     int ind;
     for(std::size_t i=0; i<left.size(); i++){
         for(std::size_t j=0; j<right.size(); j++){
             ind = give_index3(left[i], right[j], nTips);
-            dm[ind] = 2.0*nh[h] - nh[left[i]-1L] - nh[right[j]-1L];
+            dm[ind] = 2.0*nh[h] - nh[left[i]-1] - nh[right[j]-1];
         }
     }
 }
@@ -416,7 +416,7 @@ NumericVector cophenetic_cpp(IntegerMatrix edge, NumericVector edge_length,
     std::vector< std::vector<int> > bip = bipCPP(edge, nTips);
     NumericVector dm( nTips * (nTips-1) /2 );
     int l=0, left, right;
-    for(int h=nNode; h<(nNode + nTips); h++){
+    for(size_t h=nNode; h<(nNode + nTips); h++){
         // changed from NumericVector to IntegerVector
         IntegerVector tmp_ch = ch[h];
         l = tmp_ch.size();
diff --git a/src/sankoff.c b/src/sankoff.c
index 0c1edbe4..214e8e37 100644
--- a/src/sankoff.c
+++ b/src/sankoff.c
@@ -1,7 +1,7 @@
 /*
  * dist.c
  *
- * (c) 2008-2023  Klaus Schliep (klaus.schliep@gmail.com)
+ * (c) 2008-2024  Klaus Schliep (klaus.schliep@gmail.com)
  *
  *
  * This code may be distributed under the GNU GPL
@@ -31,7 +31,7 @@ SEXP C_rowMin(SEXP sdat, SEXP sn, SEXP sk){
     return(result);
 }
 
-
+/*
 void rowMin2(double *dat, int n,  int k, double *res){
     int i, h;
     double x;
@@ -41,7 +41,7 @@ void rowMin2(double *dat, int n,  int k, double *res){
         res[i] = x;
     }
 }
-
+*/
 
 double get_ps(double *dat, int n,  int k, double *weight){
   int i, h;
diff --git a/tests/testthat/test_plot_ancestral.R b/tests/testthat/test_plot_ancestral.R
index 3b79e5b2..d384ca89 100644
--- a/tests/testthat/test_plot_ancestral.R
+++ b/tests/testthat/test_plot_ancestral.R
@@ -7,7 +7,7 @@ dat <- matrix(c("a", "a",
               dimnames = list(c("t1", "t2", "t3", "t4"), NULL))
 dna <- phyDat(dat)
 fit <- pml(tree, dna)
-test_ml <- ancestral.pml(fit, type = "ml")
+test_ml <- anc_pml(fit)
 
 
 test_that("Pie_plots works", {
diff --git a/vignettes/Ancestral.Rmd b/vignettes/Ancestral.Rmd
index caafacb8..ccaa0f72 100644
--- a/vignettes/Ancestral.Rmd
+++ b/vignettes/Ancestral.Rmd
@@ -38,37 +38,30 @@ parsimony(tree, primates)
 ```
 
 For parsimony analysis of the edge length represent the observed number of changes. Reconstructing ancestral states therefore defines also the edge lengths of a tree. However there can exist several equally parsimonious reconstructions or states can be ambiguous and therefore edge length can differ.
+In _phangorn_ all the ancestral reconstructions for parsimony used to be based on the fitch algorithm (below version 3.0) and needed bifurcating trees. However trees can get pruned afterwards using the function `multi2di` from _ape_. Recently we replaced the acctran routine with a method based on the sankoff algorithm adopting the algorithm for joint reconstruction [@Pupko2000] and breaking ties at random. This has the additional benefit that to infer phylogenies with multifurcations. 
+
 "MPR" reconstructs the ancestral states for each (internal) node as if the tree would be rooted in that node. However the nodes are not independent of each other. If one chooses one state for a specific node, this can restrict the choice of neighboring nodes (figures 2 and 3).
-The function acctran (accelerated transformation) assigns edge length and internal nodes to the tree [@Swofford1987].
+There is also an option "POSTORDER" which is only a one pass algotrithm, which we use for teaching purposes. 
+The function acctran (accelerated transformation) assigns edge length and internal nodes to the tree [@Swofford1987]. 
 ```{r ancestral_reconstruction}
-anc.acctran <- ancestral.pars(tree, primates, "ACCTRAN")
-anc.mpr <- ancestral.pars(tree, primates, "MPR")
+anc.pars <- anc_pars(tree, primates)
 ```
 
-All the ancestral reconstructions for parsimony are based on the fitch algorithm and so far only bifurcating trees are allowed. However trees can get pruned afterwards using the function `multi2di` from _ape_.
-
-The `plotSeqLogo` function is a wrapper around the from the _ggseqlogo_ function the the _ggseqlogo_ package [@ggseqlogo]  and provides a simple way to show proportions of a nucleotides of ancestral states (see figure 1). 
+The `plotSeqLogo` function is a wrapper around the from the _ggseqlogo_ function in the _ggseqlogo_ package [@ggseqlogo]  and provides a simple way to show proportions of a nucleotides of ancestral states (see figure 1). 
 
 ```{r seqLogo, fig.cap="Fig 1. Ancestral reconstruction for a node.", eval=TRUE}
-#library(seqLogo)
-#seqLogo( t(subset(anc.mpr, getRoot(tree), 1:20)[[1]]), ic.scale=FALSE)
-plotSeqLogo(anc.mpr, node=getRoot(tree), 1, 20)
+plotSeqLogo(anc.pars, node=getRoot(tree), 1, 20)
 ```
 ```{r MPR, fig.cap="Fig 2. Ancestral reconstruction using MPR."}
-plotAnc(anc.mpr, 17)
+plotAnc(anc.pars, 17)
 title("MPR")
 ```
-```{r ACCTRAN, fig.cap="Fig 3. Ancestral reconstruction using ACCTRAN."}
-plotAnc(anc.acctran, 17)
-title("ACCTRAN")
-```
-
 
 # Likelihood reconstructions
 
 _phangorn_ also offers the possibility to estimate ancestral states using ML.
 The advantages of ML over parsimony is that the reconstruction accounts for different edge lengths.
-So far only a marginal construction is implemented (see [@Yang2006][@Koshi1996]) and no joint reconstruction [@Pupko2000].
+Currently a marginal construction is implemented (see [@Yang2006][@Koshi1996]), but the joint reconstruction [@Pupko2000] only for models without rate variation (e.g. gamma models) or invariant sites.
 ```{r fit_ML}
 fit <- pml(tree, primates)
 fit <- optim.pml(fit, model="F81", control = pml.control(trace=0))
@@ -81,11 +74,9 @@ and the highest posterior probability ("bayes") criterion:
 \[
 P(x_r=A) =  \frac{\pi_A L(x_r=A)}{\sum_{k \in \{A,C,G,T\}}\pi_k L(x_r=k)},
 \]
-where $L(x_r)$ is the joint probability of states at the tips and the state at the root $x_r$ and $\pi_i$ are the estimated base frequencies of state $i$.
-Both methods agree if all states (base frequencies) have equal probabilities.
+where $L(x_r)$ is the joint probability of states at the tips and the state at the root $x_r$ and $\pi_i$ are the estimated base frequencies of state $i$. Both methods agree if all states (base frequencies) have equal probabilities.
 ```{r ML_reconstruction}
-anc.ml <- ancestral.pml(fit, "ml")
-anc.bayes <- ancestral.pml(fit, "bayes")
+anc.ml <- anc_pml(fit)
 ```
 The differences of the two approaches for a specific site (17) are represented in the following figures.
 ```{r plotML, fig.cap="Fig 4. Ancestral reconstruction the using the maximum likelihood."}
@@ -93,8 +84,8 @@ plotAnc(anc.ml, 17)
 title("ML")
 ```
 ```{r plotB, fig.cap="Fig 5. Ancestral reconstruction using (empirical) Bayes."}
-plotAnc(anc.bayes, 17)
-title("Bayes")
+#plotAnc(anc.bayes, 17)
+#title("Bayes")
 ```
 
 # Fitting for discrete comparative data
@@ -142,7 +133,7 @@ all.equal(fit_SYM$logLik, fit_ace$loglik+log(1/3))
 
 
 ```{r SYM_reconstruction}
-anc_SYM <- ancestral.pml(fit_SYM, "ml")
+anc_SYM <- anc_pml(fit_SYM)
 plotAnc(anc_SYM)
 ```