Data from: Phlag: Scalable detection of genomics regions with unexplained phylogenetic heterogeneity
Data files
Jul 13, 2026 version files 8.81 GB
-
analysis-avian.tar.gz
439.02 MB
-
analysis-mammals.tar.gz
427.86 MB
-
avian-simulations-tree_sequences.tar.gz
7.94 GB
-
README.md
10.42 KB
Abstract
Phylogenetic analyses of entire genomes (phylogenomics) have revealed abundant heterogeneity of evolutionary histories. While much has been done to model this heterogeneity and to infer species trees despite it, the current toolkit has a limitation. Most methods assume that gene trees across the genome differ but are all sampled from the same distribution, defined by models such as the multi-species coalescent (MSC), and parametrized consistently across the genome. Empirical data strongly suggest this assumption is often violated because the species tree, its parameters, or the process generating the gene trees can all change across the genome. Errors in the data can further compound this heterogeneity. To address this challenge, we define the problem of detecting what segments of the genome are inconsistent with a putative species tree, even after allowing discordance according to MSC. We model gene trees not as a set, but rather as a series (a realization of a stochastic process) along genomic positions. We propose a Hidden Markov Model (HMM) approach applied to quartet statistics measured from gene trees and tie the model to MSC using simulations. The combined use of these three ideas leads to a scalable method called Phlag. On simulated and real data, we show that Phlag can detect many cases of change in underlying evolutionary processes, including reduced recombination rates, population size changes, and admixture, all using the same algorithm.
Dataset DOI: 10.5061/dryad.g79cnp65h
This repository contains the Phlag benchmarking results, gene trees, species trees, and supporting data for the avian dataset (experiment E3) and mammalian dataset (experiment E4), together with the simulated ARGs (experiments E1 and E2).
For the individual dataset (excluding ARGs), please see:
- gene trees simulated using msprime and simulation experiments: [phlag-avian-simulations] GitHub and Zenodo,
- analysis conducted on the mammalian phylogeny: [phlag-mammalian-analysis] GitHub and Zenodo,
- experiments on the Stiller2024 avian phylogeny: [phlag-avian-analysis] GitHub and Zenodo.
Mammalian analysis (analysis-mammals.tar.gz)
A dataset consisting of 19,465 gene trees inferred from chromosome 3 of the mammalian genome alignment by Foley et al., covering 241 species. Gene trees were estimated by selecting 1Kbp subalignments with minimum missing data from each 10Kbp segment and running IQ-TREE under GTR+G4. Phlag was applied to 136 internal branches under key mammalian orders (Carnivora, Chiroptera, Primates, Artiodactyla, and Rodentia), retaining branches that define a quadripartition in at least 90% of the gene trees. Each branch was analyzed individually (single focal branch) using the prior-updated mode with topology-order emissions.
Data
alltrees.tree.gz: Compressed gene trees (19,465 trees) inferred from chromosome 3 of the mammalian alignment by Foley et al. We estimated gene trees using IQ-TREE under the GTR+G4 model from 1Kbp subalignments selected with minimum missing data from each 10Kbp segment.labelled_species_tree.nwk: Species tree in Newick format with 241 mammalian taxa and internal nodes labelledI0–I239with branch lengths.ref.topology: Reference species tree topology.qqs.txt: Precomputed quartet-quartet site (QQS) frequencies for all gene trees and internal branches, used as input to Phlag.pos: Genomic positions (on the human chromosome 3 coordinate) for each gene tree window.order.txt: File identifiers for each gene tree locus/window.taxon_map.txt: Mapping from internal node labels (I0–I239) to taxonomic family names via NCBI taxonomy.taxon_map_order.txt: Mapping from internal node labels to taxonomic order names.taxdump.tar.gz: NCBI taxonomy database dump used bymap_lca.py.map_lca.py: Python script (using ete3) for resolving internal node labels to their lowest common ancestor in NCBI taxonomy.prep.sh: Preprocessing script that extracts Hellinger distances from Phlag prediction files and produces summary files.echo_cmd.sh: Script containing the Phlag commands used to generate predictions under different hyperparameter settings.
Phlag predictions
-
all_pred-{PARAMS}-chr3/: Directories containing Phlag output for chromosome 3 under different hyperparameter combinations. Each directory contains:distances_chr3.txt: Hellinger distance between the null and alternative emission distributions for each internal branch.pred-I{NODE}-{PARAMS}.txt: Per-branch prediction file containing the Phlag command, the modified species tree, decoded state predictions, and the distance metric.
The naming convention encodes hyperparameters as follows:
eap{beta}: expected number of anomalies (--expected-num-anomalies), e.g.,eap50: beta = 50.ep{1-rho}: expected anomaly proportion (--expected-anomaly-proportion), e.g.,ep005: 1-rho = 0.05.penalty{lambda}/npenalty{N}: prior penalty strength (lambda), withnprefix indicating a negative value.
Available parameter combinations:
all_pred-eap50_ep005_penalty15-chr3all_pred-eap50_ep002_penalty15-chr3all_pred-eap40_ep005_penalty15-chr3all_pred-eap100_ep005_penalty15-chr3(the main paper analysis)all_pred-eap100_ep010_penalty15-chr3
Avian analysis (analysis-avian.tar.gz)
We applied Phlag on 39,849 gene trees sampled from five macrochromosomes (chr1–chr5) of the 363-taxon avian dataset by Stiller et al. Gene trees were estimated by selecting 1Kbp subalignments with minimum missing data from each 10Kbp segment and running IQ-TREE under GTR+G4. We focused on 20 branches near the base of Neoaves, selected from 33 basal branches (60–68 Mya) highlighted by Stiller et al. that had quartet support below 0.5. Each branch was analyzed individually (single focal branch) using the prior-updated mode with topology-order emissions.
Data
sorted_genetrees/: Compressed gene trees for each chromosome, sorted by genomic position.gene_trees-Stiller2024-chr{CHR}-sorted.nwk.gz: Gene trees from the 363-taxon avian dataset by Stiller et al., selecting 1Kbp subalignments with minimum missing data from each 10Kbp segment. Gene trees were estimated using IQ-TREE under the GTR+G4 model with approximate Bayesian support.- Covers 28 autosomes (chr1–chr28) and the Z chromosome (chrZ).
qqs/: Precomputed quartet-quartet site (QQS) frequencies, used as input to Phlag.gene_trees-Stiller2024-chr{CHR}-sorted-qqs.txt.gz: Per-chromosome QQS frequencies.qqs-chr1-5.txt.gz: Concatenated QQS for the five macrochromosomes (chr1–chr5).qqs-chr1-5-Z.txt.gz: Concatenated QQS for the five macrochromosomes plus the Z chromosome.
coordinates/: Coordinate files for each chromosome assembly.coordinates-Stiller2024-chr{CHR}.txt: Genomic coordinates for gene tree loci.
positions-gene_trees-Stiller2024-chr{CHR}-sorted.txt: Genomic positions (chromosome and base pair coordinate) for each gene tree window.- Species trees:
main.tre: Main avian species tree with CU branch lengths (363 taxa; Stiller et al.).main_alternative_mod.tre: Modified species tree.main-num_generations.tre: Species tree with branch lengths in number of generations.63K_dated.tre/2023-04-dated.tre: Dated species trees with branch lengths in time (million years).castlespro_stiller.rooted.tre: Species tree with branch lengths in substitution units, estimated by CASTLES-pro.
- Flagged regions (Z chromosome analysis):
flagged_regions-Stiller2024-chrZ-species_tree.nwk: Species tree used for flagged Z chromosome analysis.flagged_gene_trees-Stiller2024-chrZ-sorted.nwk: Gene trees from flagged regions of the Z chromosome.main_species_tree-Stiller2024-chrZ_flagged.nwk.support: Species tree with support values from flagged Z chromosome gene trees.
mapping_stiller_fig2a.tsv/mapping_stiller_main.tsv/mapping_stiller_merged.tsv: Mapping from internal node labels (N{ID}) to branch numbers used in Stiller et al. Fig. 2a and in the paper figures.list_chr.txt: List of all 34 chromosome/scaffold identifiers in the assembly.plot.R: R script (using ggplot2) for generating visualization plots from Phlag predictions.
Phlag predictions
Phlag was run on individual chromosomes and on concatenations of gene trees across the five macrochromosomes (chr1–chr5), with and without the Z chromosome. Each prediction directory contains:
distances_chr{CHR}.txt/distances_concat.txt: Hellinger distance between the null and alternative emission distributions for each branch.preds-chr{CHR}-{PARAMS}.txt/preds-concat-{PARAMS}.txt: Tab-separated prediction matrix with one row per gene tree and one column per branch. Values are binary (0 = null, 1 = alternative, nan = branch not present in the gene tree).
The naming convention encodes hyperparameters as follows:
eap{beta}: expected number of anomalies (--expected-num-anomalies), e.g.,eap40: beta = 40.ep{1-rho}: expected anomaly proportion (--expected-anomaly-proportion), e.g.,ep005: 1-rho = 0.05.penalty{lambda}/npenalty{lambda}: prior penalty strength (lambda), withnprefix indicating a negative value.noprior: segment mode (no MSC-based prior on emissions).
Available parameter combinations (per-chromosome predictions):
eap40_ep005_penalty15/(prior-updated, the main paper analysis)eap40_ep005_noprior/(segment mode)eap40_ep005_npenalty15/(negative penalty)
Concatenated analyses:
preds-concat-eap100_ep005_penalty15-chr_1_5.txtanddistances_concat-chr_1_5.txt: Five macrochromosomes concatenated, beta = 100 (the main paper analysis).preds-concat-eap100_ep005_penalty15-chr_1_5_Z.txtanddistances_concat-chr_1_5_Z.txt: Five macrochromosomes + Z chromosome, beta = 100.preds-merged-eap100_ep005_penalty15.txtanddistances_concat_merged.txt: Merged analysis.eap40_ep005_noprior-concat_1_5/andeap40_ep005_noprior-concat_1_5_Z/: Concatenated analyses without prior.
ARGs from simulation experiments (avian-simulations-tree_sequences.tar.gz)
Tree sequence objects were simulated using msprime. We used Stiller et al. 2024 avian phylogeny as our demographic model.
default-{SIZE}-{INDEX}.ts: Only Neoaves, size is 250Kb.favian-{SIZE}-{INDEX}.ts: Full Aves species tree, size is set to 500Kb.anomaly_case-{BRANCH}_{CHANGE}-{SIZE}-{INDEX}.ts: ARGs deviating from the background (default-*orfavian-*) parameterization (such as population size change, recombination rate change, and admixture).
{SIZE} correspond to independent block lengths that we used to simulate ARGs independently. Resulting sequences were concatenated based on the {INDEX} value.
You can use the below Python snippet to read the tree sequence objects.
import tskit
# Load the tree sequence object
ts = tskit.load("path/to/your/file.ts")
# Print summary information
print(ts)
Code/software
Phlag is available on GitHub, and the version used in the paper can be found on Zenodo.
