Notebook source code: notebooks/23_real_world_applications__yeast_bhv_trees.ipynb
Run it yourself on binder

Yeast gene trees in BHV space

One of the biggest challenges in phylogenetics is that different data sets based on different genes can yield different phylogenetic trees. In an article published in Nature in 2003, Rokas et al. found that when data sets contain a small number of genes, it is exceedingly likely that the resulting trees will have conflicting topologies. One possible resolution to this is doing a majority-rule tree; another possibility is to adopt a Bayesian posterior over possible trees. Here, we use the geomstats implementation of BHV space (Billera et al. 2001) to look at how the trees differ based on gene sets considered, and how as the gene sets become larger and larger the resulting trees converge to both the Frechet mean of the single-gene trees and the all-gene tree.

We will be working in BHV(8), with 8 leaves. Each leaf represents one of the 8 yeast species from Rokas et al. 2003:

• S. cerevisiae (Scer), S. paradoxus (Spar), S. mikatae (Smik),

• S. kudriavzevii (Skud), S. bayanus (Sbay), S. castellii (Scas),

• S. kluyveri (Sklu), and C. albicans (Calb, outgroup).

The data we have provided in datasets/data/yeast/yeast_concatenated.fasta is concatenated alignment of all 106 considered genes (127026 basepairs). We then divide this concatenated gene situation into 106 equal-sized windows (of around 1198 basepairs each). Note: These windows do not exactly correspond to genes (gene boundaries do not agree with our fixed window size). We use ‘single-gene’ for the plot! (If you want true single-gene trees, feel free to get the alignment partitions from the authors of the paper, or get the individual sequences from the National Center for Biotechnology Information.) We build one Neighbor-Joining tree (using the biopython python package) for each window, giving us 106 trees that live in BHV(8). Each tree represents the inferred phylogenetic relationships for the 8 yeast species considering the alignment in that window. These trees will likely disagree, as different windows and different genes capture different evolutionary signals, which can lead to different inferred topologies. This is also what Rokas et al. quantified in their paper. We then increase our window size to yield 15 windows building more trees. We hopefully show that these are converging to a stable representation of what should be the underlying species tree. We confirm this by finding the Fréchet mean of the single-gene trees, and the variance around that mean. Then we find the Fréchet mean of the multi-gene trees and find the variance around that. We then plot everything using MultiDimensional Scaling to put it in 2d, and we indeed see convergence. Hoorah!

Data note: We have provided for you the processed data in data/yeast/yeast_concatenated.fasta. Some species of scientists appear to be allergic to Python, so in order to get this data, we use the phangorn R package. If you are desperate to have the R code to get it yourself, please see the code below the References section at the very bottom of this notebook.

0. Imports

from collections import Counter

import matplotlib.pyplot as plt
from sklearn.manifold import MDS

import geomstats.backend as gs
import geomstats.datasets.utils as data_utils
from geomstats.learning.frechet_mean import FrechetMean
from geomstats.metric_geometry.bhv_space import TreeSpace

gs.random.seed(66)

def count_and_plot_topologies(trees, tree_str, n_common=5):
    """Count and plot the k most common different topologies.

    Parameters
    ----------
    trees : list of Tree
        Trees to look at.
    tree_str : str
        Just for printing/plotting convenience.
    n_common : int
        Max number of most common topologies to plot
    """
    all_topologies = [tree.topology.topology_key() for tree in trees]
    unique_topologies = set(all_topologies)

    counts = Counter(all_topologies)
    print(f"{len(unique_topologies)} unique topologies across {tree_str} trees")
    print("Top 5 (topology_id, count):")
    trees_to_plot = []
    for i, (topology_key, count) in enumerate(counts.most_common(n_common)):
        print(f"\ttopology #{i + 1}: {count}")
        trees_to_plot.append(single_gene_trees[all_topologies.index(topology_key)])

    _, axes = plt.subplots(1, len(trees_to_plot), figsize=(n_common, n_common // 2 + 1))
    axes = axes.flatten()
    for i, tree in enumerate(trees_to_plot):
        tree.plot(root_id=0, ax=axes[i])
    plt.suptitle("Top 5 single-gene yeast phylogenetic tree topologies in BHV(8)")
    plt.tight_layout()
    plt.show()

1. Load data

single_gene_trees = data_utils.load_yeast(n_windows=106)
print(f"Built {len(single_gene_trees)} 'single-gene' trees in BHV(8)")
print(f"Example: {single_gene_trees[0]}\n")
count_and_plot_topologies(single_gene_trees, "single-gene", n_common=5)

Built 106 'single-gene' trees in BHV(8)
Example: (({0, 1}|{2, 3, 4, 5, 6, 7}, {0, 1, 2}|{3, 4, 5, 6, 7}, {0, 1, 2, 3}|{4, 5, 6, 7}, {0, 1, 2, 3, 4, 5}|{6, 7}, {0, 1, 2, 3, 4}|{5, 6, 7});[0.02698943 0.01330342 0.01215568 0.02141068 0.05232679])

9 unique topologies across single-gene trees
Top 5 (topology_id, count):
        topology #1: 36
        topology #2: 30
        topology #3: 13
        topology #4: 12
        topology #5: 8

Great, let’s make some multi-gene trees!

We expect there will be more agreement among these.

multi_gene_trees = data_utils.load_yeast(n_windows=15)
print(f"Built {len(multi_gene_trees)} 'multi-gene' trees in BHV(8)")
print(f"Example: {multi_gene_trees[0]}\n")
count_and_plot_topologies(multi_gene_trees, "multi-gene", n_common=5)

Built 15 'multi-gene' trees in BHV(8)
Example: (({0, 1}|{2, 3, 4, 5, 6, 7}, {0, 1, 2, 5, 6, 7}|{3, 4}, {0, 1, 2}|{3, 4, 5, 6, 7}, {0, 1, 2, 3, 4, 5}|{6, 7}, {0, 1, 2, 3, 4}|{5, 6, 7});[0.01956877 0.00959495 0.01005255 0.01562352 0.06409424])

3 unique topologies across multi-gene trees
Top 5 (topology_id, count):
        topology #1: 7
        topology #2: 7
        topology #3: 1

2. Let’s calculate a mean!

First just of the single-gene boys.

bhv8 = TreeSpace(n_labels=8)
mean_estimator = FrechetMean(
    bhv8, sample_method="stochastic", max_iter=10000, epsilon=1e-4
)
mean_estimator = mean_estimator.fit(single_gene_trees)
single_gene_mean = mean_estimator.estimate_

print(
    f"Average distance from mean: {bhv8.metric.dist(single_gene_trees, single_gene_mean).mean():0.4f}"
)

Average distance from mean: 0.0244

Now the multi-gene boys.

mean_estimator = FrechetMean(
    bhv8, sample_method="stochastic", max_iter=10000, epsilon=1e-4
)
mean_estimator = mean_estimator.fit(multi_gene_trees)
multi_gene_mean = mean_estimator.estimate_

print(
    f"Average distance from mean: {bhv8.metric.dist(multi_gene_trees, multi_gene_mean).mean():0.4f}"
)

Average distance from mean: 0.0113

Are they the same?

_, axes = plt.subplots(1, 2, figsize=(5, 3))
single_gene_mean.plot(root_id=0, ax=axes[0])
axes[0].set_title("single-gene mean")
multi_gene_mean.plot(root_id=0, ax=axes[1])
axes[1].set_title("multi-gene mean")
plt.suptitle(
    f"Dist. between mean trees: {bhv8.metric.dist(single_gene_mean, multi_gene_mean):0.4f}"
)
plt.tight_layout()
plt.show()

if single_gene_mean.topology.topology_key() == multi_gene_mean.topology.topology_key():
    print("Single- and multi-gene means are same topology.")
else:
    print("Single- and multi-gene means are NOT same topology.")

Single- and multi-gene means are same topology.

3. Some visualisation with multi-dimensional scaling (MDS)

We embed all trees (106 single-gene trees + 15 multi-gene trees + Fréchet single-gene mean + Fréchet multi-gene mean) into two dimensions using metric MDS with the distance matrix.

Distance matrix

all_trees = single_gene_trees + multi_gene_trees + [single_gene_mean, multi_gene_mean]
n_single, n_multi, n_all = len(single_gene_trees), len(multi_gene_trees), len(all_trees)

# dist_matrix = bhv8.metric.dist(all_trees, all_trees)
dist_matrix = gs.zeros((n_all, n_all))
for i, tree in enumerate(all_trees):
    dist_matrix[i, :] = bhv8.metric.dist(tree, all_trees)

MDS plotting

coords = MDS(
    n_components=2, metric="precomputed", random_state=66, init="random", n_init=1
).fit_transform(dist_matrix)

single_coords = coords[:n_single]
multi_coords = coords[n_single : n_single + n_multi]
single_mean_coord = coords[-2]
multi_mean_coord = coords[-1]

fig, ax = plt.subplots(figsize=(7, 6))


sc = ax.scatter(
    single_coords[:, 0],
    single_coords[:, 1],
    color="lightblue",
    s=35,
    label="single-gene trees",
)
ax.scatter(
    multi_coords[:, 0],
    multi_coords[:, 1],
    marker="^",
    s=250,
    color="pink",
    edgecolors="k",
    zorder=5,
    label="multi-gene trees",
)
ax.scatter(
    *single_mean_coord,
    marker="*",
    s=350,
    color="blue",
    edgecolors="k",
    zorder=6,
    label="single-gene Fréchet mean",
)
ax.scatter(
    *multi_mean_coord,
    marker="*",
    s=350,
    color="magenta",
    edgecolors="k",
    zorder=7,
    label="multi-gene Fréchet mean",
)

ax.legend(fontsize=9, loc="lower left")
ax.set_title("BHV(8): yeast phylogenetic trees")
ax.set_xlabel("MDS1")
ax.set_ylabel("MDS2")
plt.tight_layout()
plt.show()

4. Convergence: dare we do something quantitative

We will divvy up the basepairs into an increasing number of genes. Then with that increasing information we will build NJ trees and measure the variance of these trees (as the average of squared distance from their Fréchet mean). We will hopefully show that the tree estimate stabilises (decreasing variance) as the number of genes included in the calculation grows.

window_counts = [106, 50, 40, 30, 20, 5, 1]
dists_from_single_mean = []
variances = []
for k in window_counts:
    trees = data_utils.load_yeast(n_windows=k)
    mean_tree = mean_estimator.fit(trees).estimate_

    dist_from_single_mean = bhv8.metric.dist(single_gene_mean, trees).mean()
    dists_from_single_mean.append(dist_from_single_mean)

    variance = bhv8.metric.squared_dist(mean_tree, trees).mean()
    variances.append(variance)

num_genes = [106 / k for k in window_counts]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
ax1.plot(num_genes, variances, c="blue")
ax2.plot(num_genes, dists_from_single_mean, c="orange")

ax1.set_xlabel("Number of 'genes' in tree")
ax2.set_xlabel("Number of 'genes' in tree")
ax1.set_ylabel("BHV variance")
ax2.set_ylabel("Dist to single-gene Fréchet mean")
ax1.set_yscale("log")
ax2.set_yscale("log")
ax1.grid(True)
ax2.grid(True)
plt.suptitle("More genes in tree computation converge toward species tree in BHV space")
plt.tight_layout()
plt.show()

References

Rokas, A., Williams, B., King, N. et al. Genome-scale approaches to resolving incongruence in molecular phylogenies. Nature 425, 798–804 (2003). https://doi.org/10.1038/nature02053

Billera, Louis J., Susan P. Holmes, and Karen Vogtmann. “Geometry of the space of phylogenetic trees.” Advances in Applied Mathematics 27.4 (2001): 733-767. https://doi.org/10.1006/aama.2001.0759

R code

Here is the R code. Put it in an appropriately named file in the notebooks directory and run it by Rscript save_yeast_fasta.R.

# Run once to export the Rokas et al. 2003 yeast gene alignments from the
# phangorn R package as a FASTA file.

# Install phangorn if needed:
# install.packages("phangorn", repos = "https://cloud.r-project.org")
library(phangorn)

data(yeast)

outdir <- "../datasets/data/yeast"
dir.create(outdir, showWarnings = FALSE)

# single concatenated phyDat
cat("Detected: single concatenated phyDat\n")
cat("Sequences:", length(yeast), "\n")
cat("Total sites:", attr(yeast, "nr"), "\n")
cat("Species:", paste(names(yeast), collapse = ", "), "\n\n")

out_path <- file.path(outdir, "yeast_concatenated.fasta")
write.phyDat(yeast, file = out_path, format = "fasta")
cat(sprintf("Wrote concatenated alignment to %s\n", out_path))