R/plot_orthotree.R

In neurogenomics/orthogene: Interspecies gene mapping

#' Create a phylogenetic tree of shared orthologs
#' 
#' Automatically creates a phylogenetic tree plot annotated with metadata
#' describing how many orthologous genes each species shares with the
#'  \code{reference_species} ("human" by default). 
#'  
#' @param tree A phylogenetic tree of class \link[ape]{phylo}. If no tree 
#' is provided (\code{NULL}) a 100-way multiz tree will be imported from  
#' \href{http://hgdownload.soe.ucsc.edu/goldenPath/hg38/multiz100way}{
#' UCSC Genome Browser}.
#' @param orth_report An ortholog report from one or more species 
#' generated by \link[orthogene]{report_orthologs}.
#' @param species Species to include in the final plot.
#'  If \code{NULL}, then all species from the given database (\code{method}) 
#'  will be included (via \link[orthogene]{map_species}), 
#'  so long as they also exist in the \code{tree}. 
#' @param clades A named list of clades each containing 
#' a character vector of species used to define the respective clade 
#' using \link[ggtree]{MRCA}. 
#' @param clades_rotate A list of clades to rotate (via \link[ape]{rotate}), 
#' each containing a character vector of species used to define the respective
#'  clade using \link[ggtree]{MRCA}. 
#' @param show_plot Whether to print the final tree plot.
#' @param save_paths Paths to save plot to.
#' @param height Saved plot height.
#' @param width Saved plot width. 
#' @param mc.cores Number of cores to parallelise different steps with.  
#' @inheritParams report_orthologs
#' @inheritParams convert_orthologs 
#' @inheritParams map_species 
#' @inheritParams ggtree_plot 
#' @inheritParams prepare_tree 
#' @inheritParams ggplot2::ggsave 
#' 
#' @source \href{https://yulab-smu.top/treedata-book/chapter8.html#phylopic}{
#' ggtree tutorial}
#' 
#' @returns A list containing:  
#' \itemize{
#' \item{plot : }{Annotated ggtree object.}
#' \item{tree : }{The pruned, standardised phylogenetic tree used in the plot.}
#' \item{orth_report : }{Ortholog reports for each \code{species} against
#'  the \code{reference_species}}.
#' \item{metadata : }{Metadata used in the plot,
#'  including silhouette PNG ids from \href{http://phylopic.org/}{phylopic}.}
#' \item{clades : }{Metadata used for highlighting clades.}
#' \item{method : }{\code{method} used.}
#' \item{reference_species : }{\code{reference_species} used.}
#' \item{save_paths : }{\code{save_paths} to plot.}
#' } 
#' @export
#' @importFrom methods show
#' @examples
#' if(require("ape")){
#' orthotree <- plot_orthotree(species = c("human","monkey","mouse"))
#' }
plot_orthotree <- function(tree = NULL,
                           orth_report = NULL,
                           species = NULL,
                           method = c("babelgene",
                                      "homologene",
                                      "gprofiler"),
                           tree_source = "timetree",
                           non121_strategy = "drop_both_species",
                           reference_species = "human", 
                           clades = list("Primates"=c("Homo sapiens",
                                                      "Macaca mulatta"),
                                         "Eutherians"=c("Homo sapiens",
                                                        "Mus musculus",
                                                        "Bos taurus"
                                         ), 
                                         "Mammals"=c("Homo sapiens",
                                                     "Mus musculus",
                                                     "Bos taurus",
                                                     "Ornithorhynchus anatinus",
                                                     "Monodelphis domestica"
                                                     ), 
                                         "Tetrapods"=c("Homo sapiens",
                                                       "Mus musculus",
                                                       "Gallus gallus",
                                                       "Anolis carolinensis",
                                                       "Xenopus tropicalis"),
                                         "Vertebrates"=c("Homo sapiens",
                                                       "Mus musculus",
                                                       "Gallus gallus",
                                                       "Anolis carolinensis",
                                                       "Xenopus tropicalis",
                                                       "Danio rerio"),
                                         "Invertebrates"=
                                             c("Drosophila melanogaster",
                                               "Caenorhabditis elegans")
                                         ),
                           clades_rotate = list(),
                           scaling_factor = NULL,
                           show_plot = TRUE,
                           save_paths = c(
                               tempfile(fileext = ".ggtree.pdf"),
                               tempfile(fileext = ".ggtree.png")
                           ),
                           width = 15, 
                           height = width,
                           mc.cores = 1,
                           verbose = TRUE){
    # devoptera::args2vars(plot_orthotree, reassign = TRUE)
    requireNamespace("ggtree")
    
    #### species ####
    if(is.null(species)){
        species <- map_species(method = method)$scientific_name 
    } 
    #### Standardise reference species ####
    reference_species_og <- reference_species
    reference_species <- map_species(species = reference_species,
                                     method = method,
                                     verbose = FALSE)
    #### Generate orth_report ####
    if(is.null(orth_report)){ 
        orth_report <- report_orthologs(
            target_species = species, 
            reference_species = reference_species,
            method_all_genes =  method,
            method_convert_orthologs = method,
            correct_intraspecies = TRUE,
            non121_strategy = non121_strategy,
            mc.cores = mc.cores, 
            verbose = verbose)$report 
        species <- orth_report$target_species
    }
    #### tree ####
    if(is.null(tree)){
        tr <- prepare_tree(species = species,
                           tree_source = tree_source,
                           method = method,
                           run_map_species = c(TRUE, TRUE),
                           show_plot = FALSE,
                           verbose = verbose)
    } else { tr <- tree} 
    #### Rotate clades ####
    tr <- rotate_clades(tr = tr,
                        clades = clades_rotate) 
    #### silhouettes ####
    ## The first human silhouette is a weird handprint. Skip that one.
    which <- rep(1,length(tr$tip.map))
    for(s in c("homo sapiens","caenorhabditis elegans",
               "rattus norvegicus","ornithorhynchus anatinus",
               "felis catus")){
        bool <- tolower(names(tr$tip.map))==s
        if(sum(bool)>0){
            which[bool] <- 2
        } 
    } 
    uids <- get_silhouettes(species = names(tr$tip.map),
                            which = which,
                            mc.cores = mc.cores,
                            verbose = verbose) 
    #### Make clades metadata ####
    clades <- prepare_clades(tree = tr,
                             clades = clades,
                             verbose = verbose)
    #### Merge metadata ####
    d <- plot_orthotree_metadata(orth_report = orth_report,
                                 uids = uids,
                                 tr = tr,
                                 verbose = verbose)
    #### Determine scaling factor ####
    if(is.null(scaling_factor)){
        scaling_factor <- 35*nrow(d)
    }
    #### plot ####
    p <- ggtree_plot(tr = tr,
                     d = d,
                     scaling_factor = scaling_factor,
                     clades = clades,
                     reference_species = reference_species_og, 
                     verbose = verbose)
    #### Show plot ####
    if(isTRUE(show_plot)) methods::show(p) 
    #### Save plot ####
    if(length(save_paths)>0){
        for(path in save_paths){
            dir.create(dirname(path),showWarnings = FALSE, recursive = TRUE)
            messager("Saving plot ==>",path,v=verbose)
            ggplot2::ggsave(path, p, 
                            width = width, 
                            height = height)
        } 
        # file.copy(save_paths,"~/Downloads/",overwrite = TRUE)
    } 
    return(list(plot=p,
                tree=tr,
                orth_report=orth_report,
                metadata=d,
                clades=clades,
                method=method,
                reference_species=reference_species,
                save_paths=save_paths))
}