Running UPGMA on Real Data in Python

Code along video

Learning objectives

In this code along, we will complete main.py so that it reads a distance matrix from a file given on the command line, runs UPGMA, and writes the resulting tree to file in Newick format. Then we will run our code on eight real biological datasets.

Code along summary

Setup

We will continue working in the upgma directory. You can move the code for our sample tree from the previous lessons into its own function (or delete it), since main() will now be devoted to running UPGMA on real data.

Command line arguments

Our program will take two command line arguments: the location of the input data, and the directory where we should write the Newick tree. Recall that sys.argv[0] is always the name of the program, so the two arguments are sys.argv[1] and sys.argv[2], and if sys.argv does not have exactly three elements, then we raise an error.

Once we have the name of the input file, we call read_matrix_from_file() from io_util.py. Its type hints tell us exactly what it does: it takes the name of a file and returns a tuple containing a list of strings (the species names) and a distance matrix. Then we pass the matrix and the species names into upgma(); remember that upgma() takes the matrix first.

Finally, we need a name for the output file. If our input file is Data/Cetaceans/cetaceans.mtx, then we would like to write our tree to Output/cetaceans.nwk. The Path class from Python’s built-in pathlib module makes this easy: Path(input_file).stem strips away both the directories at the start of the file name and its extension, leaving just cetaceans, to which we add .nwk. We then call write_newick_to_file(), which takes the tree, the output directory, and the output file name. Here is our complete main.py.

import sys
from pathlib import Path

from functions import upgma
from io_util import read_matrix_from_file, write_newick_to_file


def main() -> None:
    if len(sys.argv) != 3:
        raise ValueError("Error: incorrect number of command line arguments. Usage: main.py <input_file> <output_dir>")

    input_file = sys.argv[1]
    output_dir = sys.argv[2]

    species_names, matrix = read_matrix_from_file(input_file)

    tree = upgma(matrix, species_names)

    output_filename = Path(input_file).stem + ".nwk"
    write_newick_to_file(tree, output_dir, output_filename)
    print(f"Tree written to {output_dir}/{output_filename}")


if __name__ == "__main__":
    main()

Running UPGMA on real data

Each dataset lives in its own subfolder of Data, so the path to a matrix includes that subfolder. (In the video, the first attempt fails with a No such file or directory error for exactly this reason.) Navigate into python/src/upgma and run the following command to build a tree for our cetaceans dataset and write it to the Output directory. If you are on Windows, use python in place of python3.

python3 main.py Data/Cetaceans/cetaceans.mtx Output

If all goes well, you will see a message that the tree was written to Output/cetaceans.nwk. Open this file, and you will see a Newick string with lots of parentheses and branch lengths and exactly one semicolon at the end. Now let’s do the same for the rest of our datasets.

python3 main.py Data/Cytochrome-C/cytochrome_c.mtx Output
python3 main.py Data/Giant-Panda-Bears/giant_panda_bears.mtx Output
python3 main.py Data/Great-Apes/great_apes.mtx Output
python3 main.py Data/HBA1/hemoglobin.mtx Output
python3 main.py Data/HIV-Subtypes/hiv_subtypes.mtx Output
python3 main.py Data/Mitochondrial-Haplogroups/mitochondrial_haplogroups.mtx Output
python3 main.py Data/UK-SARS-CoV-2/SARS_spike.mtx Output

Datasets

Each subfolder of Data contains a note.txt file explaining where its data came from. Most datasets use Levenshtein (edit) distances between protein or nucleotide sequences, which count the minimum number of insertions, deletions, and substitutions needed to transform one sequence into another; the datasets built from complete mitochondrial genomes use Jukes-Cantor corrected distances.

• Cetaceans. Complete mitochondrial genomes for 12 mammals, including whales, dolphins, the hippopotamus, and several hoofed mammals. This dataset concerns one of the most surprising results in molecular phylogenetics: whales and dolphins evolved from within the even-toed ungulates.
• Cytochrome-C. Cytochrome c protein sequences from 16 diverse eukaryotes. Cytochrome c evolves slowly under strong functional constraint, which makes it a canonical target for UPGMA. (We do not visualize this tree in the next lesson.)
• Giant-Panda-Bears. Complete mitochondrial genomes for 10 species in and around the bear family. Where does the giant panda belong?
• Great-Apes. Complete mitochondrial genomes for great apes and other primates, along with a mouse and a zebra finch as outgroups.
• HBA1. Hemoglobin subunit alpha protein sequences from 108 animals.
• HIV-Subtypes. Gag polyprotein sequences from reference strains of HIV-1, HIV-2, and simian immunodeficiency viruses (SIVs) from other primates.
• Mitochondrial-Haplogroups. Sequences from a hypervariable region of human mitochondrial DNA representing major haplogroups from around the world.
• UK-SARS-CoV-2. Spike protein sequences from 261 SARS-CoV-2 genomes sampled in the UK every two weeks between November 2020 and March 2024. Each sample is labeled with its date.

Looking ahead

Our Output directory now contains a Newick file for each dataset. But a Newick file is not much to look at. In the next lesson, we will use R to turn these files into pictures of trees, and then we will take a tour of what they tell us about evolution.