Build Your Own Tree

This page walks through the full pipeline: from downloading sequences off GenBank to a finished tree in R. All four scripts are in the GitHub repo and run on any laptop with R 4.0+.

What You Need

R 4.0+ from r-project.org. Four packages: ape and phangorn from CRAN; msa and Biostrings from Bioconductor. Run this once, then comment it out.

install.packages(c("ape", "phangorn"))
if (!requireNamespace("BiocManager")) install.packages("BiocManager")
BiocManager::install(c("msa", "Biostrings"))

For Bayesian analysis you also need MrBayes (runs in the terminal, outside R). Mac: brew install mrbayes. Other systems: binaries at mrbayes.sourceforge.net. Confirm by parsing mb in a terminal.

The four scripts (00_prepare_alignment.R, 01_tree_NJ.R, 02_tree_ML.R, 03_tree_Bayesian.R) are in the R/ folder of the repo. Run them from that folder as your working directory. Output folders (out/rdata/, out/trees/) are created automatically.

Getting Sequences from GenBank

All sequences come from NCBI Nucleotide. The gene is COI, the standard DNA barcode for animal species, ~650 bp long. Search with: [species name][ORGN] AND COI[GENE] AND 600:700[SLEN]. Download 2–3 sequences per species in FASTA format. Combine them: cat *.fasta > combined.fasta. Edit headers to clean short names like >Red_Drum_1.

Accession numbers used in this project:

Script 00: Alignment and Setup

00_prepare_alignment.R is the shared setup all three tree scripts depend on. Run it once (~30 seconds). It reads sequences, aligns them, and writes the NEXUS file for MrBayes.

# Read unaligned sequences directly into Biostrings format
dna_strings <- Biostrings::readDNAStringSet(
  "FASTA/gulf_fish_COI_combined.fasta"
)

# Align with MUSCLE
msa_result <- msa(dna_strings, method = "Muscle")
cat("Alignment length:", ncol(as.matrix(msa_result)), "bp\n")

# Convert to ape (for NJ distances) and phangorn (for ML)
aln_dnabin <- msaConvert(msa_result, type = "ape::DNAbin")
aln_phydat <- msaConvert(msa_result, type = "phangorn::phyDat")

# Save both objects plus the color helper for downstream scripts
save(aln_dnabin, aln_phydat, family_colors, tip_color,
     file = "out/rdata/alignment_objects.RData")

Output: FASTA/gulf_fish_COI_aligned.fasta, MrBayes/gulf_fish_COI_aligned.nex (for Script 03), and out/rdata/alignment_objects.RData, the main handoff to all three tree scripts, containing the aligned matrix in both DNAbin and phyDat formats plus the shared color palette.

Script 01: Neighbor-Joining

01_tree_NJ.R builds two NJ trees (rooted on Bull Shark and Flathead Mullet) with 1,000 bootstrap replicates each. NJ is the fastest method and serves as both a sanity check and a starting topology for ML.

# Compute K80 pairwise distances
dist_k80 <- dist.dna(aln_dnabin, model = "K80")

# Build unrooted NJ tree, then root on Bull Shark
tree_nj_raw   <- nj(dist_k80)
outgrp_shark  <- grep("Bull_Shark", tree_nj_raw$tip.label, value = TRUE)
tree_nj_shark <- root(tree_nj_raw, outgroup = outgrp_shark,
                       resolve.root = TRUE)

# Bootstrap: resample columns 1000 times, rebuild NJ each time
boot_nj_shark <- boot.phylo(
  tree_nj_shark, aln_dnabin,
  function(x) root(
    nj(dist.dna(x, model = "K80")),
    outgroup = grep("Bull_Shark", rownames(x), value = TRUE),
    resolve.root = TRUE
  ),
  B = 1000, quiet = TRUE
)

When reading the output: check that same-species sequences cluster together (validates labels), that all four Sciaenidae group together (tests family monophyly), and that the outgroup has a long branch while closely related ingroup species have short ones. Bootstrap at the Sciaenidae node should be near 100%; below 70% at internal nodes just means branching order is unresolved, not an error.

Script 02: Maximum Likelihood

02_tree_ML.R runs model selection, fits a GTR+Gamma+I tree starting from the NJ topology, and runs 1,000 bootstrap replicates. Plan for 10–20 minutes on a laptop.

modelTest() from phangorn ranks models by AIC. Regardless of the winner, the script uses GTR+Gamma+I: all six substitution rates free, rate variation across sites modeled with gamma, invariant sites estimated.

# Start ML search from the NJ topology (faster convergence)
fit_start <- pml(tree_nj_start, data = aln_phydat)
fit_ml    <- optim.pml(
  fit_start,
  model         = "GTR",
  optGamma      = TRUE,        # estimate Gamma shape parameter
  optInv        = TRUE,        # estimate proportion invariant
  rearrangement = "stochastic",# thorough topology search
  control       = pml.control(trace = 0)
)

# 1000 bootstrap replicates (~10-20 min)
bs_ml      <- bootstrap.pml(fit_ml, bs = 1000, optNni = TRUE,
                              control = pml.control(trace = 0))
tree_ml_bs <- plotBS(fit_ml$tree, bs_ml, type = "n")

rearrangement = "stochastic" adds random perturbations to avoid local optima. Starting from the NJ tree speeds up convergence. Bootstrap scale: ≥95% = strong; 70–94% = moderate; 50–69% = weak; <50% = not usually reported. Weak internal support within a recently diversified family like Sciaenidae is normal for a single 650 bp gene.

Script 03: Bayesian Inference

03_tree_Bayesian.R runs in two steps. First run: appends the MrBayes block to the NEXUS file and exits with instructions. Second run (after MrBayes finishes): reads the consensus tree and produces the figure.

What MrBayes is doing

MrBayes uses MCMC to explore tree space, randomly proposing and accepting/rejecting trees based on posterior probability. After a burn-in period, the surviving sampled trees form a consensus. This run: 500,000 generations, sampled every 500th (1,000 total), first 250 discarded as burn-in, leaving 750 trees for the consensus. Key settings: nst=6 (GTR), rates=gamma, ngen=500000. Increase ngen to 2,000,000+ for larger datasets or poor convergence.

begin mrbayes;
  set autoclose=yes;
  outgroup Bull_Shark_1;
  lset nst=6 rates=gamma;
  mcmc ngen=500000 samplefreq=500 printfreq=1000;
  sumt burnin=250;
  sump burnin=250;
end;

Running MrBayes

After the first R run, switch to a terminal and run:

cd R/MrBayes
mb gulf_fish_COI_aligned.nex

Takes ~5–10 minutes for this dataset. Output goes to R/MrBayes/, including .con.tre (consensus tree) and .p (parameter log). Then run 03_tree_Bayesian.R a second time to produce the figure.

Checking MCMC convergence

Load the .p file into Tracer (beast.community/tracer). A good run shows a "fuzzy caterpillar" trace with no drift, and ESS > 200 for all parameters. Low ESS or a drifting trace means the chain hasn't converged; re-run with larger ngen.

PP = 0.95 means 95% of sampled trees contained that clade, a direct probability statement. PP values run numerically higher than ML bootstrap and aren't directly comparable. Use them as complementary evidence.

Extending This Analysis

Adding your own species

Download 2–3 COI sequences, clean the headers (e.g., >Spotted_Seatrout_1), and append them to FASTA/gulf_fish_COI_combined.fasta. Re-run from Script 00. If the new species is a new family, add it to family_colors in Script 00 first.

Using different genes

COI is the default, but 16S rRNA (more conserved, better for deeper splits) and cytochrome b (large GenBank database) both work well. Multi-gene analyses are more powerful but require sequences from the same individuals; don't mix accessions across genes from different specimens.

Scaling up with TAMU-HPRC

For 50+ species or Bayesian runs with ngen=5,000,000+, a laptop isn't practical. TAMU-HPRC is free for students and faculty and can finish the same MrBayes run in hours. Get started at hprc.tamu.edu; create an account, complete the orientation, and submit via SLURM.