Evolutionary Phylogeny

This walkthrough follows one alignment from input checks through tree inference and interpretation. The commands are templates: verify each program version and analysis setting locally, then keep the resulting files and logs with your project.

Learning objectives and prerequisites

By the end of this lesson you will be able to:

• Produce a reproducible multiple sequence alignment (MSA) for homologous sequences using MAFFT.
• Build a phylogenetic tree from that alignment using MEGA-CC (command-line MEGA) or alternative CLI tools (IQ-TREE, FastTree).
• Record and interpret key analysis choices: alignment file, evolutionary model or method, rooting strategy, and support values.
• Critically evaluate what the resulting tree does and does not demonstrate about evolutionary relationships.

Prerequisites:

• Familiarity with biological sequence files (FASTA), basic UNIX command line, and the concept of homologous characters and multiple sequence alignment.
• Access to a Unix-like environment (Linux, macOS, or WSL) with the MAFFT aligner and one selected tree-inference program installed from their official instructions. The tutorial uses file paths as examples; adapt them to your workspace.

Why evolutionary phylogeny matters (biological and computational context)

Phylogenetic analysis reconstructs relationships among genes, proteins, or organisms by comparing homologous characters (nucleotides or amino acids). Trees summarize hypotheses of shared ancestry: topology (branching pattern) suggests relative relationships, branch lengths often represent estimated genetic change, and support values indicate statistical confidence in clades under the chosen model and data. Phylogenies are used in systematics, comparative genomics, functional inference, and epidemiology, but they are hypotheses, not definitive histories; results depend on the input alignment, model choice, taxon sampling, and methodological assumptions.

To infer phylogenies robustly, perform two main preprocessing steps: (1) generate a reliable MSA with MAFFT that aligns homologous positions across sequences, and (2) choose a phylogenetic method (distance, maximum likelihood, Bayesian, or approximate methods) and explicitly state the evolutionary model and the support strategy (bootstraps, posterior probabilities). This tutorial emphasizes reproducible command-line workflows and careful recording of all settings.

Reproducible workflow (sequential)

Overview: work in a project directory containing raw sequences, intermediate files, configuration files, and results. These examples use paths relative to the project directory. Supply homologous, unaligned sequences in a single FASTA file; records should represent the same locus or a defensible set of orthologous characters, not merely sequences with similar names. Create the results directory before analysis. Validate program flags against installed versions.

1) Prepare a project folder and raw input

# Create a project directory (safe operation)
mkdir -p ~/phylogeny-project/results
cd ~/phylogeny-project

# Place your own unaligned FASTA file here; no example biological data are supplied.
test -s sequences.fasta && echo "Input file present" || { echo "STOP: add a non-empty sequences.fasta" >&2; exit 1; }

Use one of the equivalent structural preflight routes below on sequences.fasta. Both require at least three records, unique first-token IDs, and non-empty sequences. They do not determine whether records are homologous, the right orthologs, or the right biological sequence type. The input is your unaligned FASTA file; no R package is required for the primary check.

R · base R (primary route)

read_fasta <- function(path) {
  lines <- readLines(path, warn = FALSE)
  records <- list()
  for (line in trimws(lines)) {
    if (!nzchar(line)) next
    if (startsWith(line, ">")) {
      header <- trimws(sub("^>", "", line))
      tokens <- strsplit(header, "\\s+")[[1]]
      id <- if (length(tokens)) tokens[[1]] else ""
      if (!nzchar(id)) stop("STOP: empty FASTA identifier")
      records[[length(records) + 1L]] <- list(id = id, seq = "")
    } else {
      if (!length(records)) stop("STOP: sequence text before first FASTA header")
      records[[length(records)]]$seq <- paste0(records[[length(records)]]$seq,
                                                gsub("\\s+", "", line, perl = TRUE))
    }
  }
  ids <- vapply(records, `[[`, character(1), "id")
  seqs <- vapply(records, `[[`, character(1), "seq")
  if (length(records) < 3L || anyDuplicated(ids) || any(!nzchar(seqs)) ||
      any(!nzchar(ids))) {
    stop("STOP: require >=3 records, unique IDs, and non-empty sequences")
  }
  list(ids = ids, seqs = seqs)
}
fasta <- read_fasta("sequences.fasta")
cat("FASTA structure OK:", length(fasta$ids),
    "records; confirm sequence type and homology biologically\n")

Save as check_fasta.R and run Rscript check_fasta.R from the project directory. The expected success line reports the record count; invalid input stops with a diagnostic. This is a structural check, not a scientific validation.

Python · optional alternative

from pathlib import Path

records = []
for line in Path("sequences.fasta").read_text().splitlines():
    line = line.strip()
    if not line:
        continue
    if line.startswith(">"):
        records.append([line[1:].split()[0], ""])
        if not records[-1][0]:
            raise SystemExit("STOP: empty FASTA identifier")
    elif not records:
        raise SystemExit("STOP: sequence text before first FASTA header")
    else:
        records[-1][1] += "".join(line.split())
ids = [name for name, _ in records]
if len(records) < 3 or len(set(ids)) != len(ids) or any(not seq for _, seq in records):
    raise SystemExit("STOP: require >=3 records, unique IDs, and non-empty sequences")
print(f"FASTA structure OK: {len(records)} records; confirm sequence type and homology biologically")

Save as check_fasta.py and run python3 check_fasta.py from the project directory. Dependencies: Python 3 standard library only. With the same input it checks the same three conditions as the R route.

Files you should have at this stage:

• sequences.fasta: raw unaligned nucleotide or protein sequences in FASTA format.
• A writable results/ directory for outputs.

2) Install tools

Install MAFFT and one tree-inference program using their official instructions. MEGA-CC distribution, licensing, and command-line support depend on the MEGA release and platform; IQ-TREE and FastTree are alternatives, not additional required steps. Verify each executable is available and record its version before analysis.

Installation is intentionally not given as a single Conda command: package names, channels, and availability differ among tools and platforms. Install each program from its official instructions. Check MEGA documentation for the supported MEGA-CC package and license for your environment.

3) Multiple sequence alignment with MAFFT

MAFFT aligns multiple nucleotide or amino-acid sequences. Choose the input type that matches the biological question and inspect the result. For coding sequences, protein-guided alignment can be appropriate when codon structure matters, but this lesson does not perform back-translation or codon-aware inference. --auto selects an MAFFT strategy based on input; it does not determine homology or make a poor input biologically valid.

# Native MAFFT CLI: unaligned FASTA input to aligned FASTA output
mafft --auto sequences.fasta > results/aligned_sequences.fasta

This writes the alignment to results/aligned_sequences.fasta. The result must be reviewed for homology, alignment quality, expected sequence type, and problematic regions before inference.

4) Inspect the alignment

# Quick text-based inspection; view header and first alignment block
head -n 200 results/aligned_sequences.fasta

# Re-run the base-R structural check shown below after MAFFT completes
Rscript check_alignment.R

Run one of these equivalent checks on results/aligned_sequences.fasta. They require at least three records with equal, non-zero sequence lengths. They do not assess alignment quality or biological homology. The R option again uses base R only.

R · base R (primary route)

read_aligned_fasta <- function(path) {
  lines <- readLines(path, warn = FALSE)
  seqs <- character()
  for (line in trimws(lines)) {
    if (!nzchar(line)) next
    if (startsWith(line, ">")) {
      seqs <- c(seqs, "")
    } else {
      if (!length(seqs)) stop("STOP: sequence text before first FASTA header")
      seqs[length(seqs)] <- paste0(seqs[length(seqs)], gsub("\\s+", "", line, perl = TRUE))
    }
  }
  lengths <- nchar(seqs)
  if (length(seqs) < 3L || any(lengths == 0L) || length(unique(lengths)) != 1L)
    stop("STOP: inspect alignment; records must share non-zero length")
  lengths
}
aligned_lengths <- read_aligned_fasta("results/aligned_sequences.fasta")
cat("Alignment structure OK:", length(aligned_lengths), "records x",
    aligned_lengths[[1]], "columns\n")

Save as check_alignment.R and run Rscript check_alignment.R from the project directory after MAFFT completes. Success reports record count and aligned columns; a malformed or unequal-length file stops.

Python · optional alternative

from pathlib import Path

seqs = []
for line in Path("results/aligned_sequences.fasta").read_text().splitlines():
    line = line.strip()
    if line.startswith(">"):
        seqs.append("")
    elif line:
        if not seqs:
            raise SystemExit("STOP: sequence text before first FASTA header")
        seqs[-1] += "".join(line.split())
lengths = {len(seq) for seq in seqs}
if len(seqs) < 3 or 0 in lengths or len(lengths) != 1:
    raise SystemExit("STOP: inspect the alignment; records must share non-zero length")
print(f"Alignment structure OK: {len(seqs)} records x {next(iter(lengths))} columns")

Save as check_alignment.py and run python3 check_alignment.py from the project directory. Dependencies: Python 3 standard library only. It checks the same count and length properties as the R version, not alignment quality.

Inspect these features:

• Do sequence names match expectations and are unique?
• Are there long stretches of gaps or sequences with excessive missing data?
• Has the alignment preserved expected conserved motifs (if known)?

5) Build a tree with MEGA-CC (command-line MEGA)

MEGA-CC analyses use a MEGA analysis-options (.mao) file. Create and validate the settings with the matching MEGA version and consult its current command-line documentation before running: the example below is a template, not a version-independent contract.

# Example file locations:
# results/aligned_sequences.fasta  (input alignment)
# ML_Tree_Settings.mao             (MEGA options file you prepared)

# Run MEGA-CC (adjust path to megacc binary if not on PATH)
megacc -a ML_Tree_Settings.mao -d results/aligned_sequences.fasta -o results/megacc_output

Expected MEGA outputs (look in results/megacc_output/):

• *.nwk: Newick-format tree files (may be rooted/unrooted depending on settings).
• Log files and result summaries that include substitution model, likelihood, and support values (if bootstraps were requested).

Record the exact .mao file, the version reported by the installed MEGA-CC executable using its documented version option, and the full command line in your project notes. Do not assume megacc --version is accepted across releases; check local help/documentation.

6) Alternative: build with IQ-TREE (explicit model selection & bootstraps)

IQ-TREE can select a substitution model with ModelFinder and infer a maximum-likelihood tree. The example requests 1,000 ultrafast bootstrap replicates; this is a resampling-based support procedure, not a guarantee that the model or alignment is adequate.

# Run model selection plus tree inference, with 1,000 ultrafast bootstrap replicates.
# -pre sets the output prefix; -wbt asks IQ-TREE to write bootstrap trees.
iqtree -s results/aligned_sequences.fasta -m MFP -bb 1000 -wbt -pre results/iqtree_run

# To show the version for provenance
iqtree --version

IQ-TREE will produce:

• results/iqtree_run.treefile: ML tree with branch lengths.
• results/iqtree_run.ufboot: bootstrap trees, because this command includes -wbt.
• results/iqtree_run.log and results/iqtree_run.iqtree: model selection output, log-likelihoods, and parameter estimates.

7) FastTree for very large datasets (approximate ML)

# For nucleotide data, request the GTR model explicitly with -gtr and set nucleotide mode with -nt
FastTree -gtr -nt results/aligned_sequences.fasta > results/fasttree.nwk

FastTree writes a Newick tree to standard output, redirected here to a file. For nucleotide data this command uses GTR+CAT. FastTree reports local support values on a 0-to-1 scale; consult its documentation before comparing them with bootstrap percentages from another program.

Output inspection and interpretation

After producing trees, inspect both the tree file contents and the associated logs to interpret meaningfully:

• Open the Newick file in a tree viewer (FigTree, iTOL, Dendroscope) or use text tools to display structure. Example: cat results/iqtree_run.treefile.
• Check support values: bootstrap percentages (IQ-TREE) or ultrafast bootstrap (UFboot) values near nodes. Support values are method-specific and should not be compared as though they were the same statistic. Interpret them with the analysis report and model assumptions; no universal cutoff establishes a clade as true.
• Branch lengths generally represent estimated substitutions per site under the fitted model. A long branch does not by itself imply elapsed time, and a short branch does not alone establish that a relationship is unresolved.
• Model choice: in IQ-TREE logs, ModelFinder reports the best-fitting substitution model and AIC/BIC scores. Record this model as part of your analysis provenance.
• Check the program documentation and tree file for root placement. A rooted tree gives an orientation for interpreting lineage order, but its root is an inference/choice, not a direct observation of an ancestor. Use a justified outgroup where appropriate or clearly label midpoint rooting; report the method. Verify the tree viewer has not changed it silently.

What a plausible result means:

• A supported clade means the sampled alignment and chosen inference procedure repeatedly recover that grouping under its assumptions. It is evidence about this dataset and analysis, not proof of organismal history or a direct measure of evolutionary distance.
• Branch-length patterns can suggest relative rates of evolution but do not alone prove absolute timing or direction of change without calibration.

What the tree does NOT prove:

• A single tree does not definitively establish true historical events; trees are model-dependent hypotheses.
• Low support nodes indicate unresolved relationships, not necessarily absence of relationship.
• Convergence, horizontal gene transfer, recombination, or alignment error can mislead inference; additional evidence is often required to support strong evolutionary claims.

Troubleshooting (symptoms, likely causes, fixes)

• Symptom: Alignment shows sequences with many terminal gaps or extremely short aligned length.
  Likely cause: Mixed sequence types (proteins and nucleotides) or truncated sequences.
  Fix: Verify sequence type and convert consistently (translate nucleotides to protein if aligning coding sequences by amino acids), remove obviously incomplete sequences, and re-run MAFFT.
• Symptom: MEGA-CC fails with a parsing error for the .mao file.
  Likely cause: The .mao file is corrupted or produced by an incompatible MEGA version.
  Fix: Recreate the .mao file using the same MEGA GUI version as your MEGA-CC installation, or use IQ-TREE if its methods and assumptions fit the analysis; inspect its log and report rather than relying on a command alone.
• Symptom: One or more branches have low support under the selected analysis.
  Likely cause: The data may contain little signal for that split, or alignment, sampling, or model assumptions may be inadequate; these are distinct possibilities.
  Fix: Inspect the affected alignment region and confirm taxon/locus orthology, then assess model fit and whether sampling addresses the biological question. Add taxa or loci only when justified and analyze concatenated data with an appropriate partition/model strategy. Higher support after changing data is not, by itself, validation of the relationship.
• Symptom: Unexpected topologies contrasting well-established relationships.
  Likely cause: Long-branch attraction, recombination, paralogy (mixing paralogous gene copies), or poor alignment.
  Fix: Validate sequence identity, orthology, alignment, and taxon sampling first; investigate recombination or gene-tree/species-tree discordance where relevant. Consider partitioning only when biologically justified partitions exist (for example, distinct genes or codon positions), and assess model adequacy; a more complex model alone does not guarantee the expected topology.
• Symptom: Tool installation issues or missing binaries.
  Likely cause: Channel mismatch or version conflicts in Conda/Mamba environment.
  Fix: Check the installation instructions for each tool and the channels available on your platform; resolve channel priority using the current Conda/Mamba documentation rather than assuming one command covers MAFFT and whichever single tree-inference tool you selected. MEGA-CC is an alternative, not an additional required installation.

Limitations, version and reference-data caveats, and responsible use

Limitations:

• Results depend on sequence sampling, alignment quality, model choice, and method. No single method guarantees the "true" tree.
• Bootstraps and posterior probabilities measure data consistency under chosen assumptions, not absolute truth.
• FastTree uses an approximate maximum-likelihood approach; its model and support summaries differ from other ML workflows, so state the method rather than treating it as interchangeable with a full alternative.

Version and reference-data caveats:

• Software behavior and available models change across versions. Record program versions using each installed tool's documented version command (for example, IQ-TREE documents iqtree --version); do not assume version flags are identical across tools.
• Substitution model sets differ between programs; the same model name may have different parameterizations. Report the exact model string and where it was selected.
• Reference data (taxon names, outgroup selection) can bias results. Document why each sequence was included and any sequences omitted.

Responsible use:

• Avoid overinterpreting trees. When making biological claims (e.g., taxonomic changes, epidemiological inferences), combine phylogenetic evidence with additional data and appropriate statistical tests.
• Follow data-sharing policies and privacy considerations for sequence data (e.g., human-derived viruses) and cite sequence sources.

Save an interpretable analysis record

1. Record provenance: save copies of aligned_sequences.fasta, the MEGA .mao file (if used), IQ-TREE logs, and the exact commands in command.txt.
2. Visualize the primary tree in a GUI (FigTree or iTOL) and annotate clades of interest.
3. If a branch has low support, diagnose signal, orthology, alignment, sampling, and model adequacy before considering more loci or justified partitioning; do not optimize support alone.
4. Prepare a figure caption that lists alignment method, program versions, substitution model, bootstrap method and replicates, and rooting choice.
5. Deposit alignment and tree files in an archive (e.g., Dryad, Zenodo) if publishing results to enable reproducibility.

References

• MAFFT: Multiple sequence alignment program, a general reference for MSA alternatives.
• MEGA-CC command-line guide: official invocation and .mao options-file instructions.
• IQ-TREE documentation: ModelFinder, ultrafast bootstraps, and command-line examples.
• FastTree documentation: approximate ML for very large alignments.
• MAFFT official documentation: installation and command-line usage.