phylogenetic analyses

Phylogenetic comparative analysis uses the concept of shared ancestry, as described by a phylogenetic tree, to account for the non-independence of related species in trait data, measuring phylogenetic signal, testing correlations between traits with independent contrasts, and estimating ancestral trait values. This analysis can be used to answer questions such as “Do closely related species in this tree have more similar values of this trait than distantly related species?”, “Are these two traits correlated after accounting for shared ancestry?”, or “What value of this trait is estimated for the ancestor at this node?”


Once a tree is built, it can be used to ask how traits evolved. Three questions come up whenever chemistry is compared across species: do close relatives have similar trait values (phylogenetic signal), are two traits related once shared ancestry is taken into account (independent contrasts), and what trait values did the ancestors probably have (ancestral trait reconstruction)? These questions matter because related species are not independent data points: two sister species may share a compound simply because both inherited it from the same ancestor. This chapter answers all three with the chemical bloom data and the angiosperm tree from the previous chapter. We start with a wrapper function, runPhylogeneticAnalyses(), that runs these phylogenetic comparative analyses in one call, then take each analysis in turn.

chemical_bloom_tree <- buildTree(
  scaffold_type = "newick",
  scaffold_in_path = "http://thebustalab.github.io/data/angiosperms.newick",
  members = unique(chemical_blooms$label)
)

runPhylogeneticAnalyses(
    traits = pivot_longer(chemical_blooms[,1:4], cols = c(3:4), names_to = "trait", values_to = "value"),
    column_w_names_of_tiplabels = "label",
    column_w_names_of_traits = "trait",
    column_w_values_for_traits = "value",
    tree = chemical_bloom_tree
)
## # A tibble: 310 × 17
##    parent  node branch.length label isTip     x     y branch
##     <int> <dbl>         <dbl> <chr> <lgl> <dbl> <dbl>  <dbl>
##  1     80     1         290.  Gink… TRUE   352.     1   207.
##  2     80     1         290.  Gink… TRUE   352.     1   207.
##  3     81     2         267.  Pice… TRUE   352.     2   219.
##  4     81     2         267.  Pice… TRUE   352.     2   219.
##  5     81     3         267.  Cupr… TRUE   352.     3   219.
##  6     81     3         267.  Cupr… TRUE   352.     3   219.
##  7     84     4         135.  Eryt… TRUE   352.     5   285.
##  8     84     4         135.  Eryt… TRUE   352.     5   285.
##  9     86     5          16.2 Iris… TRUE   352.     6   344.
## 10     86     5          16.2 Iris… TRUE   352.     6   344.
## # ℹ 300 more rows
## # ℹ 9 more variables: angle <dbl>, trait <chr>,
## #   value <dbl>, trait_type <chr>,
## #   phylogenetic_signal_k_value <dbl>,
## #   phylogenetic_signal_k_p_value <dbl>,
## #   phylogenetic_signal_lambda_value <dbl>,
## #   phylogenetic_signal_lambda_p_value <dbl>, pic <dbl>

For all the below, there are some structural requirements: (i) the tree needs to be a phylo object, (ii) the traits need to be a long-format data.frame, with one column of species names, one column naming the trait and one column holding its value (pivot_longer() makes this from a wide table, as in the examples below), and (iii) the species names must exactly match the tip labels of the tree (though they don’t have to be in the same order), for example:

chemical_bloom_tree <- buildTree(
  scaffold_type = "newick",
  scaffold_in_path = "http://thebustalab.github.io/data/angiosperms.newick",
  members = unique(chemical_blooms$label)
)
## IMPORTANT: Some species substitutions or removals were made as part of buildTree. Run build_tree_substitutions() to see them all.
## Pro tip: most tree read/write functions reset node numbers. Fortify your tree and save it as a csv file to preserve node numbering.

phylogeneticSignal

Phylogenetic signal is a measure of the degree to which related species share similar trait values. It is used to determine whether a trait has evolved in a manner that is consistent with the species’ evolutionary history. phylochemistry provides the phylogeneticSignal function, which can be used to calculate phylogenetic signal for a given set of traits and a phylogenetic tree. Here is an example:

phylogeneticSignal(
  traits = pivot_longer(chemical_blooms, cols = c(2:10), names_to = "compound", values_to = "value"),
  column_w_names_of_tiplabels = "label",
  column_w_names_of_traits = "compound",
  column_w_values_for_traits = "value",
  tree = chemical_bloom_tree
)
##             trait trait_type n_species number_of_levels
## 1         Alkanes continuous        78               NA
## 2    Sec_Alcohols continuous        78               NA
## 3          Others continuous        78               NA
## 4     Fatty_acids continuous        78               NA
## 5        Alcohols continuous        78               NA
## 6   Triterpenoids continuous        78               NA
## 7         Ketones continuous        78               NA
## 8 Other_compounds continuous        78               NA
## 9       Aldehydes continuous        78               NA
##   evolutionary_transitions_observed
## 1                                NA
## 2                                NA
## 3                                NA
## 4                                NA
## 5                                NA
## 6                                NA
## 7                                NA
## 8                                NA
## 9                                NA
##   median_evolutionary_transitions_in_randomization
## 1                                               NA
## 2                                               NA
## 3                                               NA
## 4                                               NA
## 5                                               NA
## 6                                               NA
## 7                                               NA
## 8                                               NA
## 9                                               NA
##   minimum_evolutionary_transitions_in_randomization
## 1                                                NA
## 2                                                NA
## 3                                                NA
## 4                                                NA
## 5                                                NA
## 6                                                NA
## 7                                                NA
## 8                                                NA
## 9                                                NA
##   evolutionary_transitions_in_randomization
## 1                                        NA
## 2                                        NA
## 3                                        NA
## 4                                        NA
## 5                                        NA
## 6                                        NA
## 7                                        NA
## 8                                        NA
## 9                                        NA
##   phylogenetic_signal_k_value phylogenetic_signal_k_p_value
## 1                 0.060905743                         0.186
## 2                 2.045611108                         0.001
## 3                 0.029719595                         0.711
## 4                 0.067575546                         0.318
## 5                 0.053729903                         0.353
## 6                 0.566806846                         0.001
## 7                 0.239488396                         0.030
## 8                 0.018420367                         0.868
## 9                 0.008596969                         0.928
##   phylogenetic_signal_lambda_value
## 1                           0.4522
## 2                           0.9999
## 3                           0.0001
## 4                           0.0001
## 5                           0.0001
## 6                           0.9999
## 7                           0.7874
## 8                           0.0001
## 9                           0.0001
##   phylogenetic_signal_lambda_p_value
## 1                             0.3748
## 2                             0.0000
## 3                             1.0000
## 4                             1.0000
## 5                             1.0000
## 6                             0.0000
## 7                             0.0450
## 8                             1.0000
## 9                             1.0000

independentContrasts

Phylogenetic independent contrasts are a method for analyzing the relationship between two or more traits while taking into account the evolutionary history of the species being studied. This method involves transforming the data into “independent contrasts” to remove the effects of shared ancestry, allowing for more accurate analysis of the relationship between traits. phylochemistry provides the independentContrasts function to calculate phylogenetic independent contrasts for a given set of traits and a phylogenetic tree. Here is an example of calculating independent contrasts for an example dataset.

contrasts <- independentContrasts(
  traits = pivot_longer(chemical_blooms, cols = c(2:10), names_to = "compound", values_to = "value"),
  column_w_names_of_tiplabels = "label",
  column_w_names_of_traits = "compound",
  column_w_values_for_traits = "value",
  tree = chemical_bloom_tree
)

# buildLinearModel(
#   data = contrasts,
#   formula = "Fatty_acids = Alkanes + 0"
# ) -> model
# 
# ggplot(model$data) +
#   geom_point(aes(x = input_x, y = input_y)) +
#   geom_line(aes(x = model_x, model_y))

ancestralTraits

Ancestral trait reconstruction is a method to infer the characteristics (or “traits”) of ancestral organisms based on the traits of their modern descendants. By examining the traits of present-day species and using phylogenetic trees, we can estimate or “reconstruct” the traits of common ancestors. This method can be applied to various types of traits, including continuously varying and discrete traits. Ancestral trait reconstruction helps us gain insights into the evolutionary processes and the historical transitions that led to current biodiversity. phylochemistry provides the function ancestralTraits to perform these operations. Note that ancestralTraits is different from buildTree’s “ancestral_states”. “ancestral_states” estimates ancestral sequence states at phylogeny nodes, while ancestralTraits will estimate the traits of an ancestor, given the traits of extant species that are present on the leaves of a phylogeny. Here is an example. Please note that ancestralTraits accepts data in a long-style data frame.

anc_traits_tree <- ancestralTraits(
  traits = pivot_longer(chemical_blooms, cols = -1),
  column_w_names_of_tiplabels = "label",
  column_w_names_of_traits = "name",
  column_w_values_for_traits = "value",
  tree = chemical_bloom_tree
)
head(anc_traits_tree)
## # A tibble: 6 × 11
##   parent  node branch.length label  isTip     x     y branch
##    <int> <dbl>         <dbl> <chr>  <lgl> <dbl> <dbl>  <dbl>
## 1     80     1          290. Ginkg… TRUE   352.     1   207.
## 2     80     1          290. Ginkg… TRUE   352.     1   207.
## 3     80     1          290. Ginkg… TRUE   352.     1   207.
## 4     80     1          290. Ginkg… TRUE   352.     1   207.
## 5     80     1          290. Ginkg… TRUE   352.     1   207.
## 6     80     1          290. Ginkg… TRUE   352.     1   207.
## # ℹ 3 more variables: angle <dbl>, trait <chr>, value <dbl>

In addition to providing ancestral state estimations, there is also a function for plotting those estimations on a phylogeny: geom_ancestral_pie. Here is an example. pie_size is the size of the pie chart that will be plotted at each node. geom_ancestral_pie relies on having columns in its input called trait and value, such as those output by ancestralTraits. Note that an object passed to ggtree() with duplicate node names requires the distinct function to remove the duplicates, otherwise geom_ancestral_pie will get confused about where to place the pies.

ggtree(
  distinct(anc_traits_tree, node, .keep_all = TRUE)
) +
  geom_ancestral_pie(
    data = filter(anc_traits_tree, isTip == FALSE),
    pie_size = 0.1, pie_alpha = 1
  ) +
  geom_tiplab(offset = 20, align = TRUE) +
  scale_x_continuous(limits = c(0,650)) +
  theme_void()

chapter summary

Download this summary as a one-page PDF

The four analysis functions in this chapter take the same inputs: a tree as a phylo object, and a long table of traits with one column of species names that match the tree’s tip labels exactly. Choosing among the analyses is a matter of choosing the question. Phylogenetic signal asks whether related species share similar trait values. Independent contrasts remove the effects of shared ancestry before two traits are compared. Ancestral trait reconstruction estimates the trait values of ancestors at the tree’s nodes. It is not the same as buildTree()’s ancestral_states in the phylogenies chapter, which computes ancestral sequences.

Function Question it answers Settings students vary
runPhylogeneticAnalyses() the comparative analyses in one call traits, column_w_names_of_tiplabels, column_w_names_of_traits, column_w_values_for_traits, tree
phylogeneticSignal() do close relatives have similar trait values? the same five arguments
independentContrasts() are two traits related once shared ancestry is removed? the same five arguments
ancestralTraits() what trait values did the ancestors have? the same five arguments
geom_ancestral_pie() pie charts of the ancestral estimates at each node of a ggtree() plot data = filter(anc_traits_tree, isTip == FALSE), pie_size, pie_alpha

The long trait table is made with pivot_longer(), and a tree passed to ggtree() for ancestral pies first needs duplicate nodes removed with distinct(anc_traits_tree, node, .keep_all = TRUE).