rer_converge.R 2.04 KB
Newer Older
1 2 3 4 5
suppressPackageStartupMessages({
  library(tidyverse)
  library(RERconverge)
})

6
master_tree <- read.tree(file="<<<dep master_tree>>>")
7
gene_trees <- readTrees("<<<dep gene_tree_set>>>", max.read=<<<max_read>>>, masterTree=master_tree, minSpecs=<<<min_specs>>>)
8
is_continuous = <<<template_bool continuous>>>
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

relative_rates <- getAllResiduals(
  gene_trees,
  transform = "<<<template_transform transform>>>",
  weighted = <<<template_bool weighted>>>,
  scale = <<<template_bool scale>>>
)

load_species_file <- function(path) {
    read_table(path, col_names = "species") %>% pull(species)
}

load_traits_file <- function(path, is_continuous) {
  if (is_continuous)
    read_table(path, col_names = c('species', 'value')) %>% deframe()
  else
    load_species_file(path)
}

28 29 30 31
phenotypes <- load_traits_file("<<<dep phenotypes>>>", is_continuous)

build_traits_paths <-
  function(phenotypes,
32 33
           gene_trees,
           is_continuous,
34
           clade)
35 36
  {
    if (is_continuous) {
37
      char2Paths(phenotypes, gene_trees)
38
    } else {
39
      foreground2Tree(phenotypes,
40 41 42 43 44 45 46 47
                      gene_trees,
                      clade = clade,
                      plot=F # FIXME : plot=T raises cryptic error
                      ) %>%
        tree2Paths(gene_trees)
    }
  }

48 49 50
clades = c("ancestral", "terminal", "all")

results = lapply(clades, FUN=function(clade){
51

52 53 54
  traits_paths <- build_traits_paths(
    phenotypes,
    gene_trees,
55 56
    is_continuous,
    clade
57
  )
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77

  res <- if (is_continuous) {
    correlateWithContinuousPhenotype(
      relative_rates,
      traits_paths,
    )
  } else {
    correlateWithBinaryPhenotype(
      relative_rates,
      traits_paths,
    )
  }

  table <- res %>% 
    mutate(gene_id=names(gene_trees$trees)) %>% 
    select(gene_id, everything())

  list(
    table=table,
    traits_paths=traits_paths
78
  )
79
})
80

81 82 83 84 85 86 87 88 89
names(results) = clades

results = results %>% 
  list_modify(
    relative_rates = relative_rates,
    master_tree = master_tree,
    gene_trees = gene_trees,
    phenotypes = phenotypes
  )
90

91
saveRDS(results, file="<<<dest>>>")