R/rate.map.R

In RRphylo: Phylogenetic Ridge Regression Methods for Comparative Studies

#' @title Mapping rate and direction of phenotypic change on 3D surfaces.
#' @description Given vectors of RW (or PC) scores, the function selects the
#'   RW(PC) axes linked to highest (and lowest) evolutionary rate values and
#'   reconstructs the morphology weighted on them. In this way, \code{rate.map}
#'   shows where and how the phenotype changed the most between any pair of
#'   taxa.
#' @usage rate.map(x, RR, PCscores, pcs, mshape, out.rem = TRUE,
#'   shape.diff=FALSE, show.names=TRUE)
#' @param x the species/nodes to be compared; it can be a single species, or a
#'   vector containing two species or a species and any of its parental nodes.
#' @param RR an object generated by using the \code{\link{RRphylo}} function
#' @param PCscores PC scores.
#' @param pcs RW (or PC) vectors (eigenvectors of the covariance matrix) of all
#'   the samples.
#' @param mshape the Consensus configuration.
#' @param out.rem logical: if \code{TRUE} mesh triangles with outlying area
#'   difference are removed.
#' @param shape.diff logical: if \code{TRUE}, the mesh area differences are
#'   displayed in an additional 3d plot.
#' @param show.names logical: if \code{TRUE}, the names of the species are
#'   displayed in the 3d plot.
#' @details After selecting the PC axes, \code{rate.map} automatically builds a
#'   3D mesh on the mean shape calculated from the Relative Warp Analysis (RWA)
#'   or Principal Component Analysis (PCA) (\cite{Schlager 2017}) by applying
#'   the function \code{\link[Rvcg]{vcgBallPivoting}} (\pkg{Rvcg}). Then, it
#'   compares the area differences between corresponding triangles of the 3D
#'   surfaces reconstructed for the species and surface of the mrca. Finally,
#'   \code{rate.map} returns a 3D plot showing such comparisons displayed on
#'   shape of the mrca used as the reference.The colour gradient goes from blue
#'   to red, where blue areas represent expansion of the mesh, while the red
#'   areas represent contractions of the mesh triangles. In the calculation of
#'   the differences of areas we supply the possibility to find and remove
#'   outliers from the vectors of areas calculated on the reference and target
#'   surfaces. We suggest considering this possibility if the mesh may contain
#'   degenerate facets. Additionally, \code{rate.map} allows to investigate the
#'   pure morphological comparison of shapes by excluding the evolutionary rate
#'   component by setting the argument \code{show.diff  = TRUE}. In this case, a
#'   second 3D plot will be displayed highlighting area differences in terms of
#'   expansion (green) and contraction (yellow).
#' @export
#' @seealso \href{../doc/RRphylo.html}{\code{RRphylo} vignette} ;
#'   \code{\link[Morpho]{relWarps}} ; \code{\link[Morpho]{procSym}}
#' @return The function returns a list including:
#'   \itemize{\item\strong{$selected} a list of PCs axes selected for higher
#'   evolutionary rates for each species. \item\strong{$surfaces} a list of
#'   reconstructed coloured surfaces of the given species and of the most recent
#'   common ancestor.}
#' @author Marina Melchionna, Antonio Profico, Silvia Castiglione, Gabriele
#'   Sansalone, Pasquale Raia
#' @references Schlager, S. (2017). \emph{Morpho and Rvcg-Shape Analysis in R:
#'   R-Packages for geometric morphometrics, shape analysis and surface
#'   manipulations.} In: Statistical shape and deformation analysis. Academic
#'   Press.
#'   Castiglione, S., Melchionna, M., Profico, A., Sansalone, G.,
#'   Modafferi, M., Mondanaro, A., Wroe, S., Piras, P., & Raia, P. (2021). Human
#'   face-off: a new method for mapping evolutionary rates on three-dimensional
#'   digital models. Palaeontology. doi:10.1111/pala.12582
#' @examples
#'   \dontrun{
#'   data(DataSimians)
#'   DataSimians$pca->pca
#'   DataSimians$tree->tree
#'   dato<-pca$PCscores
#'   cc<- 2/parallel::detectCores()
#'
#'   RRphylo(tree,dato,clus=cc)->RR
#'
#'   Rmap<-rate.map(x=c("Pan_troglodytes","Gorilla_gorilla"),RR=RR, PCscores=dato,
#'                  pcs=pca$PCs, mshape=pca$mshape, shape.diff = TRUE)
#'
#'   }

rate.map<-function (x, RR, PCscores, pcs, mshape,
                    out.rem = TRUE, shape.diff=FALSE, show.names=TRUE) {

  if (!requireNamespace("inflection", quietly = TRUE)) {
    stop("Package \"inflection\" needed for this function to work. Please install it.",
         call. = FALSE)
  }

  if (!requireNamespace("rgl", quietly = TRUE)) {
    stop("Package \"rgl\" needed for this function to work. Please install it.",
         call. = FALSE)
  }

  if (!requireNamespace("Rvcg", quietly = TRUE)) {
    stop("Package \"Rvcg\" needed for this function to work. Please install it.",
         call. = FALSE)
  }

  if (!requireNamespace("Morpho", quietly = TRUE)) {
    stop("Package \"Morpho\" needed for this function to work. Please install it.",
         call. = FALSE)
  }

  if (!requireNamespace("ddpcr", quietly = TRUE)) {
    stop("Package \"ddpcr\" needed for this function to work. Please install it.",
         call. = FALSE)
  }
  rates<-RR$multiple.rates
  tree<-RR$tree
  aces<-RR$aces
  phen<-rbind(aces,PCscores)

  if(all(x%in%tree$tip.label)) {
    if(length(x)==1){
      mrca<-getMommy(tree,which(tree$tip.label==x))[1]
      rates_sum<-rates[match(x,rownames(rates)),,drop=FALSE]
    } else {
      mrca<-getMRCA(tree,x)
      path1<-c(x[1],getMommy(tree,which(tree$tip.label==x[1])))
      path2<-c(x[2],getMommy(tree,which(tree$tip.label==x[2])))
      rates1<-apply(rates[match(path1[1:(which(path1==mrca)-1)],rownames(rates)),,drop=FALSE],2,sum)
      rates2<-apply(rates[match(path2[1:(which(path2==mrca)-1)],rownames(rates)),,drop=FALSE],2,sum)
      rates_sum<-rbind(rates1,rates2)
    }
  } else {
    mrca<-x[-which(x%in%tree$tip.label)]
    x<-x[which(x%in%tree$tip.label)]
    path1<-c(x,getMommy(tree,which(tree$tip.label==x)))
    if(!mrca%in%path1) stop("the node is not along the species path")
    rates1<-apply(rates[match(path1[1:(which(path1==mrca)-1)],rownames(rates)),,drop=FALSE],2,sum)
    rates_sum<-t(as.matrix(rates1))
  }

  sele<-list()
  sur<-list()
  rates.list<-list()


  for(i in 1:nrow(rates_sum)){
    cutter <- inflection::ede(seq(1:dim(rates_sum)[2]),rates_sum[i,order(rates_sum[i,],decreasing = TRUE)], 0)

    if (cutter[1] == 1 & cutter[2] == dim(rates_sum)[2]) {
      cutter2 <- inflection::ede(seq(1:(dim(rates_sum)[2]-2)),
                                 rates_sum[i,order(rates_sum[i,],decreasing = TRUE)][-c(1,dim(rates_sum)[2])], 0)
      cutter[1:2] <- cutter2[1:2] + 1
    } else if (cutter[1] == 1) {
      cutter2 <- inflection::ede(seq(1:(dim(rates_sum)[2]-1)),
                                 rates_sum[i,order(rates_sum[i,],decreasing = TRUE)][-1], 0)
      cutter[1] <- cutter2[1] + 1
    } else if (cutter[2] == dim(rates_sum)[2]) {
      cutter2 <- inflection::ede(seq(1:(dim(rates_sum)[2]-1)),
                                 rates_sum[i,order(rates_sum[i,],decreasing = TRUE)][-dim(rates_sum)[2]], 0)
      cutter[2] <- cutter2[2]
    }

    up=seq(1:cutter[1])
    dw=seq(cutter[2],dim(rates_sum)[2])

    while ((length(up)+length(dw)) > 0.5 * dim(rates_sum)[2]){
      if (length(up)>=length(dw)) up[-length(up)]->up else dw[-1]->dw
    }

    rates.up<-rates_sum[i,order(rates_sum[i,],decreasing = TRUE)][up]
    rates.dw<-rates_sum[i,order(rates_sum[i,],decreasing = TRUE)][dw]

    sele[[i]]<-match(c(names(rates.up),names(rates.dw)),colnames(rates_sum))
    names(sele[[i]])<-c(names(rates.up),names(rates.dw))

    sur[[i]] <- dosur(scores = phen[match(as.character(x[i]),rownames(phen)),], pcs = pcs,
                      sel = sele[[i]],
                      mshape = mshape, radius = 0)
    rates.list[[i]]<-c(rates.up,rates.dw)
  }

  names(sele)<-x
  names(sur)<-x
  names(rates.list)<-x

  sur.mshape <- Rvcg::vcgBallPivoting(mshape, radius = 0)
  sur.ref <- dosur(scores = aces[match(mrca,rownames(aces)),], pcs = pcs, sel = NULL,
                   mshape = mshape, radius = 0)

  meshes<-list()
  diffs<-list()
  colmap_tot <- colorRampPalette(c("darkred","red","orange","white","lightblue","blue","darkblue"))

  if (length(x)==1){

    ddpcr::quiet(meshes[[1]]<-localmeshdiff(sur[[1]],sur.ref, ploton = 2, n.int = 10000,
                                     paltot=c("darkred","red","orange","white","lightblue","blue","darkblue"),
                                     from = NULL,to= NULL,out.rem = out.rem,fact = 1.5,visual = 1,scale01 = TRUE, plot = TRUE)$mesh)
    title(main = x[1])
    if (show.names) rgl::mtext3d(text = paste(x[which(x%in%tree$tip.label)]), font = 2, edge = "x", line = 3)

  } else {

    rgl::mfrow3d(1,2,sharedMouse = TRUE)
    ddpcr::quiet(meshes[[1]]<-localmeshdiff(sur[[1]],sur.ref, ploton = 2, n.int = 10000,
                                     paltot=c("darkred","red","orange","white","lightblue","blue","darkblue"),
                                     from = NULL,to= NULL,out.rem = out.rem,fact = 1.5,visual = 1,scale01 = FALSE, densityplot = TRUE, plot = FALSE)$mesh)
    rgl::shade3d(meshes[[1]], specular = "black")

    title(main = x[1])

    if (show.names) rgl::mtext3d(text = paste(x[1]), font = 2,edge = "x", line = 3)
    rgl::next3d()
    ddpcr::quiet(meshes[[2]]<-localmeshdiff(sur[[2]],sur.ref, ploton = 2, n.int = 10000,
                                     paltot=c("darkred","red","orange","white","lightblue","blue","darkblue"),
                                     from = NULL,to= NULL,out.rem = out.rem,fact = 1.5,visual = 1,scale01 = FALSE, densityplot = TRUE, plot = FALSE)$mesh)
    title(main = x[2])
    rgl::shade3d(meshes[[2]], specular = "black")
    if (show.names) rgl::mtext3d(text = paste(x[2]), font = 2, edge = "x", line = 3)


  }

  names(meshes)<-x

  if (shape.diff==TRUE){

    surX1<-dosur(scores=phen[match(x[1],rownames(phen)),],pcs=pcs,mshape=mshape)
    mapD.surX1 <- areadiff(surX1, sur.mshape, out.rem = out.rem,fact = 1.5)


    if (all(x%in%tree$tip.label)&length(x)>1) {
      surX2<-dosur(scores=phen[match(x[2],rownames(phen)),],pcs=pcs,mshape=mshape)
      mapD.surX2 <- areadiff(surX2, sur.mshape, out.rem = out.rem,fact = 1.5)

      msurD <- mapD.surX1[[3]] - mapD.surX2[[3]]

      rgl::open3d()
      ddpcr::quiet(shape_diff<-localmeshdiff(surX1,surX2,ploton = 1,out.rem = out.rem,vec = msurD,scale01 = TRUE,
                                      paltot = rev(c("#004251","#248e69","#5fff99","white","gold1","darkgoldenrod3","#663c14")),
                                      plot=FALSE,densityplot = TRUE)$mesh)
      title(main = paste(x,collapse = " vs. "))

      rgl::shade3d(shape_diff, specular = "black")
    } else {
      rgl::open3d()
      ddpcr::quiet(shape_diff<-localmeshdiff(surX1,sur.mshape,ploton = 1,out.rem = out.rem,scale01 = TRUE,
                                      paltot = rev(c("#004251","#248e69","#5fff99","white","gold1","darkgoldenrod3","#663c14")),
                                      plot=FALSE,densityplot = TRUE)$mesh)
      title(main = paste(x[1],"vs. consensus shape"))
      rgl::shade3d(shape_diff, specular = "black")
    }
  }

  return(list(selected=rates.list,surfaces=c(mrca=list(sur.ref),meshes)))

}