Unraveling Genomic Introgression: A Hidden Markov Model Framework for Comparative Genomics

Abstract

This article provides a comprehensive exploration of Hidden Markov Models (HMMs) and their transformative role in detecting and analyzing introgression in comparative genomic studies. Aimed at researchers, scientists, and drug development professionals, it covers foundational concepts from the basic principles of HMMs and the evolutionary significance of adaptive introgression to advanced methodological implementations in tools like Ancestry_HMM-S and PhyloNet-HMM. The scope extends to troubleshooting common analytical challenges, optimizing model parameters, and validating findings through empirical case studies and performance benchmarks against other methods. By synthesizing theoretical knowledge with practical application, this review serves as a critical resource for leveraging HMMs to uncover adaptive genetic variation with implications for evolution, disease resistance, and biomedical discovery.

The Core Concepts: Understanding Introgression and Hidden Markov Models

Defining Adaptive Introgression and Its Evolutionary Impact

Adaptive introgression (AI) is the natural process by which beneficial genetic material is transferred from a donor species to a recipient species through hybridization and repeated backcrossing, followed by selection on the introgressed alleles [1]. This evolutionary mechanism allows recipient species to rapidly acquire adaptive traits that enhance their fitness and survival in changing environments. Historically, introgression was viewed primarily as a homogenizing force that could hinder adaptation by introducing maladaptive alleles or leading to genetic swamping [1]. However, genomic studies have established that introgression can serve as a critical adaptive driving force, often enabling faster adaptation than reliance solely on de novo mutations [1] [2].

The study of adaptive introgression has been revolutionized by the genomic revolution, which has provided tools for detecting and analyzing introgressed regions across diverse taxonomic groups, from bacteria to mammals [1] [3]. This application note details the methodologies for investigating adaptive introgression, with a specific focus on the application of Hidden Markov Models (HMMs) within a comparative genomic framework.

Computational Detection: Hidden Markov Models

Hidden Markov Models (HMMs) provide a powerful statistical framework for identifying introgressed genomic regions by modeling the underlying genealogical states that give rise to observed genetic variation. These models are particularly effective because they can distinguish signals of introgression from other evolutionary processes, such as Incomplete Lineage Sorting (ILS), which can create similar patterns of topological incongruence in gene trees [4].

Core Principle of PhyloNet-HMM

The PhyloNet-HMM framework integrates phylogenetic networks with HMMs to scan aligned genomes and infer the local genealogical history at each site [4]. The model operates on the principle that genomic regions of introgressive descent will exhibit genealogies that are discordant with the predominant species tree, clustering the introgressed individual with the donor species rather than its closest relative.

The following diagram illustrates the core analytical workflow of a PhyloNet-HMM analysis for detecting introgressed genomic regions.

Key HMM-Based Tools and Protocols

Protocol 1: Detecting Introgression with PhyloNet-HMM

This protocol is designed for analyzing whole-genome alignment data from at least three species or populations to detect specific introgressed regions [4].

• Input Preparation: Generate a multiple whole-genome sequence alignment from the recipient species, donor species, and an outgroup species.
• Model Specification: Define the set of possible parental species trees (phylogenetic networks) that include the hypothesized introgression event.
• Model Training: Use an Expectation-Maximization (EM) algorithm to train the HMM on the genomic data, estimating transition and emission probabilities.
• Path Decoding: Apply the Viterbi algorithm to compute the most likely sequence of hidden states (genealogies) across the genome.
• Annotation: Identify contiguous genomic regions where the inferred genealogy is discordant with the species tree and indicates a donor-recipient relationship.

Protocol 2: Local Ancestry Inference with an HMM for Pileup Data

This protocol is suited for inferring local ancestry directly from sequencing read pileup data, which is useful for non-model organisms or low-coverage sequencing projects [5].

• Data Input: Use aligned sequencing reads (BAM files) from admixed individuals and genotype likelihoods or allele counts from reference populations.
• HMM Definition:
  • Hidden States: The possible ancestry assignments (e.g., Population A, Population B) for a genomic segment.
  • Observations: The aligned sequencing reads supporting different alleles at each site.
  • Transition Probabilities: Governed by the recombination rate and the time since admixture.
• Parameter Estimation: Simultaneously estimate local ancestry and the time since admixture by modeling the decay of ancestry tract lengths.
• Ancestry Painting: Output the probability of each ancestry type along each chromosome of the admixed individual.

Table 1: Key HMM-Based Software for Introgression Analysis

Software	Methodology	Key Features	Applicability	Reference
PhyloNet-HMM	Phylogenetic Network + HMM	Teases apart introgression from ILS; works with genome alignments	Diverged species, ancient introgression	[4]
HMM for Pileup Data	HMM on read counts	Works directly on NGS reads; estimates admixture time; any ploidy	Non-model organisms, pooled sequencing	[5]
Enhanced HMMs	Improved state/emission models	Increased accuracy for archaic introgression inference	Complex demographic histories	[6]

Performance Evaluation of Detection Methods

The performance of adaptive introgression detection methods varies significantly based on evolutionary parameters. A 2025 benchmark study evaluated three modern methods—VolcanoFinder, Genomatnn, and MaLAdapt—alongside the summary statistic Q95(w, y) [7].

Table 2: Impact of Evolutionary Parameters on Method Performance (Based on [7])

Evolutionary Parameter	Impact on Detection Power & False Positive Rate	Recommendation
Divergence Time	Power increases with longer divergence times between species.	Methods perform best on deeply diverged lineages.
Migration Time	Power is higher for more recent migration/introgression events.
Selection Coefficient	Stronger selection (higher s) leads to easier detection of AI.	Methods are robust to detect AI under strong selection.
Recombination Hotspots	Can break down AI signals, reducing power.	Use methods accounting for linked selection.
Adjacent Windows	Including them in training data is critical to correctly pinpoint the selected locus and avoid misclassification.	Ensure analytical framework models hitchhiking effects.

The study concluded that methods based on the Q95 statistic are currently the most efficient for exploratory studies of adaptive introgression, though HMM-based approaches remain powerful for full-genome annotation [7].

Empirical Evidence and Applications

Case Studies Across the Tree of Life

Adaptive introgression has been documented as a key evolutionary mechanism across a wide spectrum of life.

• Mammals (Mice): The Vkorc1 gene, which confers resistance to rodent poison, was adaptively introgressed from an Algerian mouse population into the European house mouse [4]. PhyloNet-HMM analysis of mouse chromosome 7 estimated that approximately 9% of sites (covering 13 Mbp and over 300 genes) were of introgressive origin, far exceeding the previously known single locus [4].
• Plants (Galápagos Tomato): Introgression of carotenoid biosynthesis loci from the endemic Solanum cheesmaniae (orange fruit) into the invasive Solanum pimpinellifolium (red fruit) has led to convergent evolution of orange fruits on the islands, a trait likely favored by selection [8].
• Trees (Cottonwoods): A 31-year common garden experiment with Populus trees demonstrated that introgression from the low-elevation P. fremontii into the high-elevation P. angustifolia enhanced survival and biomass accumulation in a warmer, drier environment. Specific genetic markers (e.g., RFLP-1286) were significantly associated with increased survival, providing direct evidence for climate adaptation via introgression [9].
• Bacteria: Although asexual, bacteria also experience introgression defined as gene flow in core genomes. A systematic analysis of 50 bacterial lineages found an average of 2% of core genes were introgressed, with rates up to 14% in Escherichia–Shigella, shaping their evolution without completely blurring species borders [10].

Protocol for Validating Adaptive Introgression

Protocol 3: From Genomic Scan to Functional Validation

This integrated protocol outlines the steps from computational prediction to functional validation of an adaptively introgressed allele [2] [8] [9].

• Genome-Wide Scan: Use an HMM-based or Q95-based method to scan genomes and identify candidate introgressed regions.
• Phenotypic Association: Correlate the introgressed haplotype with a specific adaptive trait (e.g., fruit color, disease resistance, climate tolerance). In a crop context, this leverages farmer knowledge about trait importance [2].
• Population Genetic Tests: Perform tests for selection (e.g., reduced diversity, extended haplotype homozygosity) on the candidate region to distinguish neutral from adaptive introgression.
• Gene Function Analysis: Annotate the region to identify candidate genes and use functional assays (e.g., gene expression, CRISPR knockout) to confirm the gene's role in the adaptive phenotype.
• Fitness Assessment: Measure the fitness consequence of the introgressed allele versus the native allele through common garden or field experiments, as demonstrated in the Populus study [9].

The following diagram maps this multi-stage validation workflow.

The Scientist's Toolkit

Table 3: Essential Research Reagents and Solutions for Introgression Studies

Reagent / Resource	Function and Application in AI Research
Reference Genomes	High-quality genomes for the recipient, donor, and outgroup species are essential for alignment and phylogenetic inference.
Whole-Genome Alignment	Serves as the primary input for methods like PhyloNet-HMM to identify genealogical discordance.
Annotated GenBank Files	Provide gene annotations to quickly identify functional elements within candidate introgressed regions.
RFLP or SNP Markers	Used for genotyping and tracking the presence of specific introgressed haplotypes in experimental populations (e.g., common gardens) [9].
Pileup Data from NGS	Raw read alignment data (BAM/CRAM files) used for local ancestry inference without prior genotype calling [5].
Common Garden Collections	Living collections of plants or organisms from different populations/species grown in a controlled environment to assess the fitness and phenotypic effects of introgressed alleles [9].

Adaptive introgression is a fundamental evolutionary mechanism that facilitates rapid adaptation across the tree of life. Hidden Markov Models provide a robust computational framework for detecting these events by systematically analyzing genomic landscapes and distinguishing introgression from confounding processes. The integration of advanced HMMs with population genetic theory and empirical validation through controlled experiments and functional assays creates a powerful paradigm for understanding how gene flow contributes to adaptation. As methodological improvements continue, particularly in handling complex demographics and diverse data types, HMMs will remain indispensable for deciphering the functional and evolutionary impact of adaptive introgression.

Fundamental Principles of Hidden Markov Models (HMMs)

Hidden Markov Models (HMMs) are powerful statistical frameworks for modeling sequential data where the underlying system is assumed to be a Markov process with unobservable (hidden) states. In computational biology, HMMs have become indispensable for analyzing genomic sequences, protein families, and evolutionary processes [11]. The model operates on the principle that a sequence of observations is generated by a sequence of hidden states, which themselves form a first-order Markov chain—meaning each state depends only on the previous state [11]. This mathematical formalism provides a computationally straightforward yet powerful approach for solving time-series problems in bioinformatics.

For comparative genomic introgression research, HMMs offer particular value by enabling researchers to distinguish true introgression signatures from spurious signals that arise due to other evolutionary processes like incomplete lineage sorting (ILS) [4]. The PhyloNet-HMM framework exemplifies this application by combining phylogenetic networks with HMMs to simultaneously capture reticulate evolutionary history and dependencies within genomes [4]. This integration allows systematic scanning of genomes for regions of introgressive descent while accounting for point mutations, recombination, and ancestral polymorphism.

Key Components and Mathematical Framework

A standard HMM consists of two interconnected components: a hidden state sequence following a first-order Markov chain, and an observable sequence where each observation depends on the current hidden state [11]. The following table summarizes the core components and their mathematical representations:

Table 1: Fundamental Components of Hidden Markov Models

Component	Mathematical Representation	Description	Role in Genomic Analysis
Hidden States (Z)	( Z = {z1, z2, ..., z_N} )	Unobservable states forming Markov chain	Evolutionary histories (e.g., species tree vs. introgression) [4]
Observations (X)	( X = {x1, x2, ..., x_N} )	Measurable data dependent on hidden states	Genomic sequence data or alignment patterns [4] [11]
Transition Probability (A)	( a{ij} = P(z{t+1} = j | z_t = i) )	Probability of moving from state i to j	Evolutionary change between different genealogies [4]
Emission Probability (B)	( bj(k) = P(xt = k | z_t = j) )	Probability of observation k given state j	Likelihood of genomic data under specific evolutionary history [4]
Initial State Distribution (π)	( πi = P(z1 = i) )	Probability of starting in state i	Prior probability of evolutionary scenarios at sequence start [4]

The joint probability of an observation sequence and hidden state sequence is given by:

[ P(X,Z) = π{z1} \prod{t=1}^{T-1} a{zt z{t+1}} \prod{t=1}^{T} b{zt}(xt) ]

This factorization enables efficient algorithms for computational inference, making HMMs practical for genome-scale analyses [11].

HMM Topologies and Variants for Genomic Analysis

Different HMM topologies have been developed to address specific modeling needs in computational biology. The topology refers to the permitted transitions between states in the underlying Markov chain [11]. The major HMM variants used in genomic research include:

Table 2: HMM Architectures and Their Genomic Applications

HMM Type	Topology	Key Features	Genomic Applications
Standard HMM	Fully connected or left-right	First-order Markov chain, geometric state duration	Basic gene finding, profile searches [11]
Generalized HMM (GHMM)	Left-right with explicit duration modeling	Non-geometric state duration distributions	Gene structure prediction (e.g., GENSCAN) [11]
Pair HMM (PHMM)	Three-state (match, insert, delete)	Generates pairs of aligned sequences	Pairwise sequence alignment, homology detection [11]
Generalized Pair HMM (GPHMM)	Hybrid of GHMM and PHMM	Models sequences of different lengths	Cross-species gene finding, synteny analysis [11]
Profile HMM	Linear left-right (match, insert, delete)	Models conserved protein domains	Protein family classification, remote homology detection [11]
PhyloNet-HMM	Network-structured	Incorporates phylogenetic networks	Comparative genomics, introgression detection [4]

For introgression research, the PhyloNet-HMM represents a significant advancement over standard architectures. It models the evolutionary history of genomic regions as evolving within the branches of parental species trees, enabling detection of introgressed regions by identifying topological incongruence in local genealogies [4]. This framework can distinguish introgression from incomplete lineage sorting, a major confounding factor in comparative genomics.

Application Notes: HMMs for Introgression Detection

Experimental Protocol: PhyloNet-HMM for Introgression Scanning

Purpose: To detect genomic regions of introgressive descent in aligned genomes while accounting for ILS and dependencies across loci.

Input Requirements:

• Multiple sequence alignments from the studied species
• A set of candidate parental species trees representing potential evolutionary histories
• Genomic coordinates and annotation data

Methodology:

• Data Preparation and Alignment

  • Obtain whole-genome sequences from target species and outgroups
  • Perform multiple sequence alignment using appropriate tools (e.g., MAFFT, MUSCLE)
  • Verify alignment quality and handle missing data appropriately

• Model Configuration

  • Define hidden states corresponding to different evolutionary histories (e.g., species tree vs. introgression scenarios)
  • Specify emission probabilities based on evolutionary models (e.g., GTR, HKY)
  • Set transition probabilities between different genealogical states

• Parameter Estimation

  • Use forward-backward algorithm to compute posterior probabilities
  • Apply Baum-Welch algorithm for parameter optimization
  • Implement multivariate optimization heuristics for convergence [4]

• Introgression Detection

  • Scan genomic positions using sliding window approach
  • Compute ( P(\text{Parental Tree}i | \text{Alignment Site}j) ) for each site j and parental tree i [4]
  • Identify regions with high probability of introgression ancestry

• Validation and Interpretation

  • Compare detected regions with known functional elements
  • Perform statistical significance testing
  • Integrate with additional evidence (e.g., gene annotations, functional data)

Expected Output:

• Probability of introgression for each genomic site
• Delineation of introgressed genomic regions
• Estimates of introgression timing and donor species

Workflow Visualization

Figure 1: PhyloNet-HMM workflow for introgression detection

Successful implementation of HMM-based introgression analysis requires both computational tools and biological data resources. The following table details essential components for establishing this research pipeline:

Table 3: Research Reagent Solutions for HMM-Based Introgression Studies

Category	Specific Resource	Function/Purpose	Implementation Notes
Software Libraries	PhyloNet [4]	Phylogenetic network analysis	Open-source platform for PhyloNet-HMM implementation
	HMMER [11]	Profile HMM construction	Protein family modeling, homology detection
	BioHMM [11]	General biological HMM applications	Gene finding, genomic segmentation
Statistical Frameworks	Forward-Backward Algorithm [11]	Posterior probability calculation	Efficient computation using dynamic programming
	Baum-Welch Algorithm [11]	Parameter estimation	Expectation-Maximization for HMM training
	Viterbi Algorithm [11]	Most probable path finding	Optimal state sequence identification
Data Resources	Multiple Genome Alignments [4]	Input observational data	UCSC Genome Browser, ENSEMBL comparisons
	Species Phylogenies [4]	Evolutionary framework constraints	Time-calibrated trees from literature
	Functional Annotations [4]	Biological interpretation	Gene ontology, pathway databases
Computational Infrastructure	High-Performance Computing Cluster	Genome-scale analysis	Parallel processing for large datasets
	Adequate Storage Solutions	Data management	Secure backup for genomic sequences and results

Advanced Protocol: Multi-Species Introgression Analysis

Purpose: To identify and characterize introgression events across multiple eukaryotic species using genome-scale data.

Specialized Reagents:

• Whole-genome sequences from minimum of three closely related species
• Reference genome annotation files (GFF/GTF format)
• Pre-computed species trees with divergence time estimates

Detailed Procedure:

• Hypothesis Formulation

  • Define candidate introgression scenarios based on biological knowledge
  • Specify phylogenetic networks representing possible evolutionary histories
  • Establish statistical thresholds for significance

• Model Training

  • Initialize parameters using empirical Bayes approach
  • Implement training with expectation-conditional maximization
  • Assess convergence using likelihood trajectory monitoring

• Genome Scanning

  • Execute PhyloNet-HMM across all chromosomes
  • Compute posterior probabilities for each introgression scenario
  • Apply multiple testing correction (e.g., Benjamini-Hochberg FDR control)

• Results Integration

  • Annotate detected regions with gene information
  • Test for functional enrichment (e.g., GO term analysis)
  • Correlate with phenotypic data when available

Troubleshooting Guide:

• Poor convergence: Adjust initialization parameters, increase iteration limits
• Computational limitations: Implement chromosome-wise analysis, then combine results
• Ambiguous signals: Incorporate additional genomic features or validation experiments

Interpretation Framework and Validation Strategies

The output from HMM-based introgression analyses requires careful biological interpretation. Key considerations include:

Statistical Confidence Metrics:

• Posterior probabilities for introgression at each site
• Genome-wide false discovery rates
• Confidence intervals for parameter estimates

Biological Validation Approaches:

• Independent verification using alternative methods (e.g., D-statistics)
• Experimental validation through functional assays
• Comparison with previously published introgression events

Evolutionary Interpretation:

• Estimation of introgression timing and directionality
• Assessment of adaptive significance through selection tests
• Correlation with biogeographic and ecological data

The application of PhyloNet-HMM to mouse chromosome 7 data demonstrates the power of this approach, successfully detecting a previously reported adaptive introgression event involving the rodent poison resistance gene Vkorc1, in addition to novel introgressed regions covering approximately 13 Mbp and over 300 genes [4]. This illustrates how HMM frameworks can reveal previously unrecognized evolutionary events at genome-wide scale.

Figure 2: HMM data flow from input to biological interpretation

Why HMMs are Uniquely Suited for Genomic Sequence Analysis

Hidden Markov Models (HMMs) are powerful statistical tools that have become fundamental to modern genomic sequence analysis. As a doubly-embedded stochastic process, HMMs describe systems with two layers: an invisible process of hidden states and a visible process of observable symbols [12]. In genomic terms, the hidden states often represent functional elements (such as exons, introns, or protein domains), while the observed symbols correspond to nucleotide or amino acid sequences [12] [13]. This dual structure makes HMMs exceptionally well-suited for genomics, where biological function (hidden state) must be inferred from molecular sequences (observed data).

The mathematical foundation of HMMs makes them ideal for capturing the complex dependencies in biological sequences. An HMM is completely specified by three probability measures: transition probabilities between states, emission probabilities of symbols from states, and initial state probabilities [12] [13]. This probabilistic framework allows researchers to decode the most likely functional architecture underlying observed sequence data, enabling predictions about gene structure, protein domains, and evolutionary relationships that would be difficult to ascertain by direct examination alone.

Core Strengths of HMMs for Sequence Analysis

Capabilities and Advantages

HMMs possess several intrinsic properties that make them uniquely appropriate for genomic applications:

• Ability to model domain correlations: HMMs effectively capture correlations between adjacent sequence elements, which is crucial for identifying features like codons in protein-coding genes or base-paired regions in structural RNAs [12]
• Handling of variable-length features: The Markov property enables natural modeling of genomic elements with variable lengths, such as exons and introns, without fixed-length constraints [12]
• Integration of evolutionary information: Profile HMMs can incorporate position-specific conservation information from multiple sequence alignments, making them sensitive detectors of remote homologs [14] [15]
• Probabilistic framework: Unlike rigid pattern-matching approaches, HMMs provide probability scores that quantify uncertainty in predictions, allowing researchers to assess confidence levels [12] [15]

Comparison to Alternative Methods

Table 1: Comparison of genomic sequence analysis methods

Method Type	Strengths	Limitations	Best Suited Applications
HMM-based	Probabilistic; handles variable lengths; incorporates evolutionary information	Computationally intensive; complex parameter estimation	Gene finding, remote homology detection, structural RNA identification
BLAST/Sequence similarity	Fast; simple interpretation	Limited sensitivity for divergent sequences; binary outcomes	Quick searches against databases; identifying close homologs
Regular expressions	Computationally efficient; simple implementation	Rigid patterns; poor handling of degeneracy	Simple motif finding; restriction enzyme sites
Machine learning (CNN/RNN)	Can capture complex patterns; minimal feature engineering	Requires large training datasets; black box nature	Pattern recognition in high-throughput data

Key HMM Architectures in Genomics

Profile Hidden Markov Models

Profile HMMs represent one of the most successful applications of HMMs in bioinformatics. These models are built from multiple sequence alignments of protein families and capture position-specific conservation patterns [14]. A profile HMM typically includes match states (representing conserved columns), insertion states (accommodating variable regions), and deletion states (handling gaps) [15]. The emission probabilities at each match state reflect the amino acid or nucleotide preferences at that position, while transition probabilities model the likelihood of insertions and deletions.

The power of profile HMMs lies in their ability to detect remote homologous relationships that evade simpler similarity searches. By accumulating evidence from weakly conserved positions across entire domains, profile HMMs can identify family members with sequence identities below 20-30% [14]. This sensitivity has made them indispensable tools for automated genome annotation and protein family classification in projects like Pfam and InterPro.

Pair-HMMs for Sequence Alignment

Pair-HMMs provide a principled probabilistic framework for pairwise sequence alignment by explicitly modeling the evolutionary relationship between two sequences [12]. A standard pair-HMM for alignment contains three primary states: match states (emitting aligned residues from both sequences), insert states (emitting a residue from one sequence only), and delete states (equivalent to inserts but for the other sequence) [12]. The Viterbi algorithm applied to this model finds the most likely alignment path, while the forward algorithm computes the probability that the sequences are related.

Unlike ad-hoc scoring schemes, pair-HMMs naturally incorporate evolutionary distances through their transition and emission probabilities. The log-odds score of the sequences being related versus unrelated provides a statistically sound measure of alignment significance [12]. This framework extends to multiple sequence alignment and has influenced the development of probabilistic alignment methods that handle indels and substitutions in a evolutionarily consistent manner.

Context-Sensitive HMMs for Structural RNA

Context-sensitive HMMs (csHMMs) represent a significant advancement for modeling RNA sequences with secondary structure constraints [12]. Unlike standard HMMs where states are independent, csHMMs introduce long-range dependencies between non-adjacent positions that form base pairs in RNA secondary structures. This allows simultaneous modeling of sequence conservation and structural requirements in RNA families.

Profile-csHMMs combine the advantages of profile HMMs with the structural modeling capabilities of csHMMs, making them particularly effective for identifying non-coding RNA genes and predicting their secondary structures [12]. These models have been successfully applied to problems including RNA structural alignment, acceleration of RNA folding, and fast noncoding RNA annotation [12].

Advanced Applications in Comparative Genomics

Detecting Introgression with PhyloNet-HMM

A sophisticated application of HMMs in evolutionary genomics is the PhyloNet-HMM framework, which detects introgression (the transfer of genetic material between species through hybridization) by combining phylogenetic networks with HMMs [4] [16]. This approach addresses the key challenge of distinguishing true introgression from spurious signals caused by incomplete lineage sorting (ILS) - a major confounding factor in comparative genomics [4].

The PhyloNet-HMM framework defines a set of random variables representing the parental species tree at each genomic site and uses the HMM to model transitions between these trees along the chromosome [4]. By incorporating both the phylogenetic relationships among species and the dependencies between adjacent sites, the method can identify genomic regions of introgressive origin with high confidence. Application to mouse genome data successfully detected a previously reported adaptive introgression event involving the rodent poison resistance gene Vkorc1, along with numerous novel introgressed regions [4] [16].

Figure 1: Workflow for detecting introgression using PhyloNet-HMM framework

Pseudogene Identification with Profile HMMs

Profile HMMs provide an effective method for distinguishing functional genes from pseudogenes in genomic sequences [15]. This application leverages the fact that pseudogenes (non-functional copies of protein-coding genes) typically follow different molecular evolutionary paths compared to functional genes, accumulating mutations that disrupt the coding capacity while retaining sequence similarity [15].

In practice, pseudogene identification uses open reading frame length combined with sequence bit scores from HMM profile analysis [15]. Functional genes typically maintain longer open reading frames and yield higher bit scores when evaluated against profile HMMs of their protein family. This approach is particularly valuable for DNA barcoding and metabarcoding studies, where inadvertent amplification of pseudogenes can lead to overestimates of species diversity or misidentification [15].

Table 2: Characteristics distinguishing functional genes from pseudogenes in HMM analysis

Feature	Functional Genes	Pseudogenes	Analysis Method
Open Reading Frame	Full-length, no stop codons	Truncated, internal stop codons	ORF length analysis
Codon Usage	Species-specific bias	Randomized	Likelihood comparison
Selection Signal	Purifying selection (dN/dS < 1)	Neutral evolution (dN/dS ≈ 1)	Evolutionary rate analysis
HMM Bit Score	High	Low	Profile HMM analysis
GC Content	Typical genomic distribution	Often atypical	Composition analysis

Experimental Protocols

Protocol 1: Building a Profile HMM for Protein Family Annotation

Purpose: To construct a profile HMM for identifying members of a protein family in genomic sequences.

Materials and Reagents:

• Set of related protein sequences (minimum 5-10 representative members)
• Multiple sequence alignment software (e.g., MAFFT, Clustal Omega)
• HMM construction software (HMMER package)
• Target genomic or protein database for searching

Procedure:

• Curate seed alignment: Collect representative sequences of the protein family, ensuring diversity while maintaining reliable homology
• Create multiple sequence alignment: Use alignment software with parameters appropriate for the evolutionary distance (typically default parameters)
• Build profile HMM: Run hmmbuild from HMMER package on the alignment to create the profile HMM
• Calibrate the model: Use hmmpress to prepare the model for searching (this step calculates E-value parameters)
• Search databases: Use hmmsearch to identify homologs in target databases
• Evaluate results: Filter hits based on E-value (typically < 0.01) and bit score thresholds

Troubleshooting Tips:

• If the model retrieves too many false positives, increase the cutoff thresholds or refine the seed alignment
• If the model misses known family members, include more diverse representatives in the seed alignment
• For domain architecture analysis, use the --domtblout option to get domain-level annotations

Protocol 2: Detecting Introgression Using PhyloNet-HMM

Purpose: To identify genomic regions of introgressive origin in comparative genomic data.

Materials and Reagents:

• Whole-genome sequencing data from multiple closely related species
• Reference genome for alignment
• PhyloNet software package [4]
• Pre-defined phylogenetic network hypothesis

Procedure:

• Data preparation: Align whole-genome sequences from target species and outgroups to a reference genome
• Variant calling: Identify single nucleotide polymorphisms (SNPs) across all genomes
• Window selection: Divide the genome into windows of appropriate size (typically 1-10 kb depending on recombination rate)
• Gene tree estimation: Infer phylogenetic trees for each window using maximum likelihood or Bayesian methods
• Configure PhyloNet-HMM: Specify the set of possible parental species trees based on evolutionary hypotheses [4]
• Run analysis: Execute PhyloNet-HMM to compute posterior probabilities of each parental tree along the genome
• Identify introgressed regions: Define introgressed intervals as contiguous windows where the posterior probability for the introgressive tree exceeds a threshold (typically > 0.95)

Validation Methods:

• Simulate data under alternative evolutionary scenarios to validate statistical power
• Compare with independent methods (e.g., D-statistics) for consistency
• Examine functional annotations of introgressed regions for biological relevance

Figure 2: State architecture of a profile HMM for sequence analysis

The Scientist's Toolkit

Essential Software and Databases

Table 3: Key research reagents and computational tools for HMM-based genomic analysis

Resource Name	Type	Function	Application Context
HMMER	Software Suite	Building and searching with profile HMMs	Protein family annotation, remote homology detection
PhyloNet	Software Package	Phylogenetic network analysis	Introgression detection, reticulate evolution
Pfam	Database	Curated collection of protein family HMMs	Functional annotation of novel sequences
TMB	R Package	Accelerated parameter estimation for HMMs	Efficient model fitting for large datasets [17]
Custom HMM Scripts	Code Repository	Implementation of specialized HMM architectures	Novel applications beyond standard tools

Parameter Estimation and Computational Considerations

Effective application of HMMs in genomics requires careful attention to parameter estimation and computational efficiency. The standard approach for estimating HMM parameters is maximum likelihood estimation, typically implemented using the Baum-Welch algorithm (a special case of the Expectation-Maximization algorithm) [17]. For models with complex state architectures or large datasets, direct numerical maximization of the likelihood may be preferable [17].

The Template Model Builder (TMB) R package provides significant computational advantages for HMM parameter estimation, particularly for large genomic datasets [17]. TMB uses automatic differentiation to compute exact derivatives of the likelihood function, enabling efficient gradient-based optimization and calculation of confidence intervals [17]. This approach can dramatically reduce computation time compared to standard implementation while providing robust uncertainty estimates for parameter values.

When working with HMMs for genomic applications, several practical considerations are essential:

• Initial parameter values: Poor initialization can lead to convergence to local maxima; multiple random restarts are recommended
• Sequence length effects: Very long sequences may cause numerical underflow; use log-space computations to maintain precision
• Model complexity: Balance model flexibility with parameter count to avoid overfitting, using model selection criteria like AIC or BIC
• Computational scaling: For genome-scale applications, consider distributed computing approaches for tractable runtime

Hidden Markov Models have established themselves as indispensable tools in genomic sequence analysis due to their principled probabilistic foundation, flexibility in modeling diverse biological features, and proven performance across numerous applications. From basic gene prediction to sophisticated detection of evolutionary events like introgression, HMMs provide a unified framework for inferring hidden biological structure from observable sequence data.

The future of HMMs in genomics lies in developing more complex architectures that incorporate additional biological realities, such as context-sensitive dependencies in RNA structure [12] and phylogenetic networks for reticulate evolution [4] [3]. As genomic datasets continue growing in size and complexity, computational efficiency will remain a critical focus, with approaches like TMB-based acceleration playing increasingly important roles [17]. The integration of HMMs with emerging machine learning methods represents another promising direction, potentially combining the interpretability and statistical foundation of HMMs with the pattern recognition power of deep learning.

For researchers investigating comparative genomic introgression, HMM-based approaches like PhyloNet-HMM offer a powerful framework for distinguishing true introgression from confounding signals like incomplete lineage sorting [4] [16]. As these methods continue maturing and incorporating more realistic evolutionary models, they will further illuminate the complex history of hybridization and genetic exchange that has shaped eukaryotic genomes.

Within the framework of a thesis on hidden Markov models for comparative genomic introgression research, a precise understanding of three core parameters is fundamental. A Hidden Markov Model (HMM) is a statistical model used to represent systems where an underlying, unobservable process (the hidden states) generates a sequence of observable data [13] [18]. The model is defined by its state space, transition probabilities, and emission probabilities [19] [20]. These parameters allow researchers to infer hidden evolutionary events, such as introgression—the transfer of genetic material between species—from observed genomic sequences [4] [16]. Proper specification of these parameters is critical for developing accurate models that can distinguish introgression from other evolutionary forces like incomplete lineage sorting [4].

Parameter Definitions and Mathematical Formalisms

State Space (Q)

The state space, denoted as ( Q = q1, q2, ..., q_N ), is the set of all possible hidden states in the model [20]. In the context of genomic introgression, these states typically represent different evolutionary histories or genealogical trees that a genomic region may have inherited from [4]. For example, in a simple three-species scenario, hidden states could represent whether a genomic locus has evolved through the species tree ((A,(B,C))) or through a reticulate history involving introgression between species B and C ((A,(C,B))) [4]. The variable ( N ) represents the total number of possible hidden states. The system is always in one of these states at any given time, but the actual state is not directly observable [19].

Transition Probability Matrix (A)

The transition probability matrix, ( A = a{11} ...a{ij}...a{NN} ), is an ( N \times N ) matrix where each element ( a{ij} ) represents the probability of transitioning from state ( i ) to state ( j )—that is, ( a{ij} = P(\text{next state is } j \mid \text{current state is } i) ) [13] [20]. The matrix dictates the dynamics of the hidden Markov process. A critical mathematical constraint is that the probabilities of transitioning from any given state to all possible subsequent states must sum to unity: ( \sum{j=1}^{N} a_{ij} = 1 \quad \forall i ) [20]. In genomic applications, these probabilities are often related to the recombination rate between loci; loci that are physically closer on a chromosome have a lower probability of transitioning between different genealogical histories due to lower recombination rates [4].

Emission Probabilities (B)

The emission probabilities, ( B = bi(ot) ), also known as observation likelihoods, define the probability of generating a specific observation ( ot ) given that the system is in hidden state ( i ): ( bi(ot) = P(ot \mid qi) ) [20]. The observations ( O = o1, o2, ..., oT ) are drawn from a vocabulary ( V = v1, v2, ..., vV ) [20]. For each possible hidden state ( qi ), the emission probabilities over all possible observations must form a valid probability distribution, summing to 1: ( \sum{k=1}^{V} bi(v_k) = 1 ) [20]. In comparative genomics, an "observation" is typically a column in a multiple sequence alignment, and the emission probability is calculated based on a phylogenetic substitution model that describes how DNA sequences evolve along the branches of the genealogy specified by the hidden state [4].

Table 1: Summary of Core HMM Parameters and Their Properties

Parameter	Symbol	Description	Mathematical Constraint	Biological Interpretation in Introgression
State Space	( Q )	Set of all possible hidden states.	( N ) states.	Possible local genealogical trees/network parental species trees.
Transition Matrix	( A )	Probabilities of moving between hidden states.	Rows sum to 1: ( \sumj a{ij} = 1 ).	Governed by recombination rate between genomic loci.
Emission Probabilities	( B )	Probabilities of observations given a hidden state.	For each state ( i ), ( \sumk bi(v_k) = 1 ).	Likelihood of aligned DNA data under a phylogenetic model for a given genealogy.

Parameterization for Genomic Introgression

The PhyloNet-HMM Framework

The PhyloNet-HMM framework is a specific implementation designed for detecting introgression in eukaryotes [4] [16]. This framework integrates a phylogenetic network, which models reticulate evolutionary events like hybridization, with an HMM that captures dependencies along the genome [4]. In this model:

• The state space is defined by a set of parental species trees. For example, in a three-taxon scenario, one state may represent the primary species tree, while another represents a tree reflecting an introgression event [4].
• Transition probabilities between these states model the probability that a recombination event has occurred between two adjacent genomic sites, switching the local genealogy [4].
• Eission probabilities for each state are computed using standard phylogenetic likelihood methods applied to the corresponding species tree or network, calculating the probability of observing a particular column of the multiple sequence alignment given that tree and an evolutionary model [4].

Workflow for Introgression Detection

The following diagram illustrates the logical workflow of the PhyloNet-HMM framework for detecting introgression, from input data to the final analysis of genomic regions.

Figure 1: PhyloNet-HMM Introgression Detection Workflow

Experimental Protocols for Parameter Inference

The Baum-Welch Algorithm

The Baum-Welch algorithm, a special case of the Expectation-Maximization (EM) algorithm, is the standard method for learning HMM parameters from data when the true states are unknown [19]. It iteratively refines the parameters ( A ) and ( B ) to maximize the likelihood of the observed data [18].

Protocol:

• Initialization: Initialize the transition matrix ( A ) and emission matrix ( B ) with random or heuristic values that satisfy the probability constraints. The initial state distribution ( \pi ) is also initialized [18].
• Expectation Step (E-Step): For the current parameter estimates, use the Forward-Backward algorithm to compute two key probabilities for each position in the observation sequence:
  • The probability of being in state ( i ) at time ( t ), given the entire observation sequence (( \gammat(i) )).
  • The probability of being in state ( i ) at time ( t ) and state ( j ) at time ( t+1 ), given the entire observation sequence (( \xit(i, j) )).
• Maximization Step (M-Step): Use the counts from the E-step to update the parameter estimates:
  • Update transition probabilities: ( a{ij} = \frac{\text{Expected number of transitions from state i to state j}}{\text{Expected number of transitions from state i}} ).
  • Update emission probabilities: ( bi(vk) = \frac{\text{Expected number of times in state i and observing } vk}{\text{Expected number of times in state i}} ).
• Iteration: Repeat steps 2 and 3 until the parameters converge (i.e., the change in log-likelihood between iterations falls below a predefined threshold).

Handling Optimization Challenges

A known limitation of the Baum-Welch algorithm is its susceptibility to converge to local optima [21] [22]. Recent research has explored strategies to mitigate this, especially when scaling HMMs:

• Neural Reparameterization and Initialization: Using neural networks to initialize HMM parameters or to reparameterize them during training can lead to more effective optimization and better local optima [21]. Studies have found that even linear reparameterizations can be as effective as non-linear ones, and that these strategies are complementary [21].
• Global Optimization Heuristics: For complex models, metaheuristic algorithms like Particle Swarm Optimization (PSO) and Genetic Algorithms (GA) have been formulated to address the HMM parameter estimation problem, helping to avoid poor local minima [22].

Table 2: Key Reagent Solutions for HMM-based Introgression Research

Research Reagent / Tool	Function / Role	Example Use Case
Multiple Sequence Alignment	Provides the observation sequence (O). Each column is a single observation.	Input genomic data for PhyloNet-HMM analysis [4].
Phylogenetic Network Model	Defines the state space (Q) of possible evolutionary histories.	Specifying the set of parental species trees for the HMM [4].
Substitution Model (e.g., GTR)	Used to calculate the emission probability ( bi(ot) ) for a DNA observation ( o_t ) given a genealogy.	Computing the likelihood of an aligned site given a hidden state (genealogy) [4].
Baum-Welch / EM Algorithm	Infers optimal HMM parameters (A, B) from unlabeled observation sequences.	Training the HMM on genomic data to learn transition and emission parameters [19].
Viterbi Algorithm	Finds the single most likely sequence of hidden states given observations and model parameters.	Decoding the most probable genealogical history path along a chromosome [4] [19].

Advanced Parameterization and Visualization

Integrating Population Genetic Parameters

In advanced models like PhyloNet-HMM, the core HMM parameters are not standalone but are directly informed by underlying population genetic parameters. The transition probabilities are a function of population genetic parameters such as recombination rates, divergence times, and effective population sizes, which can be explicitly incorporated into the model [4]. Similarly, the emission probabilities depend on mutation rates and the demographic parameters embedded within the phylogenetic network. This integration allows the model to simultaneously account for multiple evolutionary processes, including point mutations, recombination, incomplete lineage sorting (ILS), and introgression [4] [16].

Visualizing the Integrated HMM Architecture

The following diagram illustrates the integrated architecture of an HMM, showing the relationship between hidden states, observations, and the core parameters, contextualized for a genomic application.

Figure 2: HMM Architecture for Genomic Analysis

Hidden Markov Models (HMMs) are powerful statistical tools for modeling doubly-embedded stochastic processes, consisting of an invisible process of hidden states and a visible process of observable symbols [12]. In molecular biology, HMMs have been extensively applied to problems ranging from gene prediction and protein family profiling to multiple sequence alignment [12]. Their ability to capture correlations between adjacent symbols or domains makes them particularly suitable for genomic analyses where functional regions display distinct statistical properties.

In comparative genomic introgression research, HMMs provide a framework for detecting genetic material transferred between species through hybridization and subsequent back-crossing [4]. The PhyloNet-HMM framework exemplifies this application by combining phylogenetic networks with HMMs to identify introgressed genomic regions while accounting for confounding factors like incomplete lineage sorting (ILS) [4]. This approach enables researchers to scan aligned genomes, inspecting local genealogies for incongruence that signals introgression.

The Three Canonical HMM Problems

Problem 1: Evaluation - Computing the Observation Probability

The evaluation problem addresses computing the probability that an observed sequence was generated by a given HMM. For a set of aligned genomic sequences, this involves calculating P(x|Θ), the probability of the observed alignment x given the model parameters Θ. The Forward Algorithm efficiently solves this problem using dynamic programming to recursively compute probabilities while avoiding the combinatorial explosion of considering all possible state paths.

In genomic introgression studies, the evaluation problem enables computing the likelihood of observed genetic variation data under different evolutionary scenarios, including those with and without introgression. For example, when analyzing chromosome 7 in mouse genomes, researchers can compute probabilities under different phylogenetic network models to identify which regions show higher likelihood under introgression hypotheses [4].

Table 1: Variables for the Forward Algorithm

Variable	Description	Biological Interpretation in Introgression Research
x = x₁...xₗ	Observed symbol sequence	Aligned genomic sequences from multiple species
y = y₁...yₗ	Hidden state sequence	Evolutionary history at each genomic position (e.g., species tree vs. introgression)
Θ	HMM parameters	Evolutionary model including mutation rates, population sizes, introgression probabilities
αₜ(i)	Forward variable: probability of partial observations and state i at time t	Joint probability of observed genetic data up to position t and specific evolutionary history at t
t(i,j)	Transition probability from state i to j	Probability of different evolutionary histories at adjacent genomic loci, influenced by recombination
e(x	i)	Emission probability of symbol x in state i	Probability of observed genetic variants given a specific evolutionary history

Protocol: Forward Algorithm for Genomic Evaluation

Purpose: Compute the probability of observed genomic sequences under a phylogenetic network model with introgression.

Inputs:

• Multiple sequence alignment of genomes from relevant species
• Phylogenetic network model with proposed introgression events
• Estimated parameters: transition probabilities between evolutionary states, emission probabilities for genetic variants

Procedure:

• Initialization: For each state i in the state space S, compute: α₁(i) = π(i) * e(x₁|i) where π(i) represents the prior probability of evolutionary history i

• Recursion: For each genomic position t from 2 to L and each state j in S: αₜ(j) = [Σᵢ αₜ₋₁(i) * t(i,j)] * e(xₜ|j) This step incorporates the transition probabilities between evolutionary histories at adjacent positions

• Termination: Compute the total observation probability: P(x|Θ) = Σᵢ αₗ(i)

• Interpretation: Compare probabilities across different phylogenetic network models to identify regions better explained by introgression

Problem 2: Decoding - Finding the Optimal State Path

The decoding problem involves determining the most likely sequence of hidden states given the observations and model parameters. The Viterbi algorithm solves this problem using dynamic programming to efficiently find the single best path through the state space. In genomic applications, this identifies the most likely evolutionary history at each genomic position.

For introgression research, decoding reveals which genomic regions originate from introgressive descent versus those following the species tree. When applied to mouse genomic data, this approach successfully identified the Vkorc1 gene region as introgressed, along with approximately 9% of sites on chromosome 7 covering about 13 Mbp and over 300 genes [4].

Table 2: Viterbi Algorithm Variables for Introgression Detection

Variable	Description	Role in Introgression Detection
Vₜ(i)	Maximum probability of state path ending in state i at time t	Best score for evolutionary history i at genomic position t
ψₜ(i)	Backpointer to best previous state	Tracks most likely evolutionary history path across genome
P*	Probability of optimal state path	Maximum joint probability of evolutionary histories across entire genomic region
y* = y₁...yₗ	Optimal state path	Inferred evolutionary history (species tree vs. introgression) at each genomic position

Protocol: Viterbi Algorithm for Introgression Mapping

Purpose: Identify genomic regions most likely originating from introgression events.

Inputs:

• Observation probabilities from Forward Algorithm
• Transition probabilities between evolutionary states
• Phylogenetic network models representing alternative evolutionary histories

Procedure:

• Initialization: For each state i: V₁(i) = π(i) * e(x₁|i) ψ₁(i) = 0

• Recursion: For t = 2 to L and each state j: Vₜ(j) = maxᵢ [Vₜ₋₁(i) * t(i,j)] * e(xₜ|j) ψₜ(j) = argmaxᵢ [Vₜ₋₁(i) * t(i,j)]

• Termination: P* = maxᵢ Vₗ(i) yₗ* = argmaxᵢ Vₗ(i)

• Path Backtrace: For t = L-1 to 1: yₜ* = ψₜ₊₁(yₜ₊₁*)

• Interpretation: The resulting state path y* indicates at each genomic position whether the evolutionary history follows the species tree or shows evidence of introgression

Viterbi Algorithm for Introgression Detection

Problem 3: Learning - Estimating Model Parameters

The learning problem addresses estimating optimal HMM parameters from observed data. The Baum-Welch algorithm, an expectation-maximization (EM) procedure, iteratively refines parameter estimates to maximize the likelihood of the observations. For each iteration, it computes expected state counts using the Forward-Backward algorithm, then updates parameters to maximize these expectations.

In comparative genomics, parameter learning enables estimation of evolutionary parameters such as introgression timing, population sizes, and recombination rates. When analyzing eukaryotic genes, HMMs can learn distinct emission probabilities for different functional regions, such as exons versus introns, based on their characteristic sequence composition [12].

Protocol: Baum-Welch Algorithm for Evolutionary Parameter Estimation

Purpose: Estimate evolutionary parameters from genomic data without prior knowledge of introgression locations.

Inputs:

• Multiple sequence alignment from relevant taxa
• Initial estimates of transition and emission probabilities
• Phylogenetic network topology

Procedure:

• Initialization: Set initial parameters Θ⁽⁰⁾ = (π, t, e) arbitrarily or based on biological priors

• Expectation Step: Compute using Forward-Backward algorithm:

  • γₜ(i) = P(yₜ = i | x, Θ), the probability of state i at position t given observations
  • ξₜ(i,j) = P(yₜ = i, yₜ₊₁ = j | x, Θ), the probability of transition from i to j at t

• Maximization Step: Update parameters:

  • π(i) = γ₁(i)
  • t(i,j) = Σₜ ξₜ(i,j) / Σₜ γₜ(i)
  • e(x|i) = Σₜ γₜ(i) * δ(xₜ = x) / Σₜ γₜ(i)

• Iteration: Repeat steps 2-3 until convergence of the log-likelihood log P(x|Θ)

• Validation: Apply learned model to independent test data or use cross-validation

Baum-Welch Parameter Estimation Process

Application Notes for Introgression Research

PhyloNet-HMM Framework for Comparative Genomics

The PhyloNet-HMM framework represents a significant advancement for detecting introgression in genomes by combining phylogenetic networks with HMMs [4]. This approach simultaneously captures potentially reticulate evolutionary history and dependencies within genomes while accounting for incomplete lineage sorting. The method scans multiple aligned genomes for signatures of introgression by incorporating both phylogenetic networks and hidden Markov models, allowing researchers to distinguish true introgression signals from spurious ones arising from population effects.

Application of this model to variation data from chromosome 7 in the mouse (Mus musculus domesticus) genome successfully detected a recently reported adaptive introgression event involving the rodent poison resistance gene Vkorc1, along with other newly detected introgressed genomic regions [4]. The analysis estimated that approximately 9% of sites within chromosome 7 were of introgressive origin, covering about 13 Mbp and over 300 genes.

Workflow for Genomic Introgression Analysis

Genomic Introgression Analysis Workflow

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Research Reagents and Computational Tools for HMM-Based Introgression Studies

Category	Item	Function	Example Applications
Computational Frameworks	PhyloNet-HMM	Detects introgression using phylogenetic networks combined with HMMs	Scanning aligned genomes for signatures of introgression while accounting for ILS [4]
	PhyloNet	Open-source package for phylogenetic network analysis	Evolutionary hypothesis testing, network inference [4]
Biological Data Resources	Multi-species genome alignments	Reference data for comparative genomic analysis	Establishing evolutionary relationships, identifying conserved regions
	Population genomic variation data	Polymorphism patterns within and between species	Distinguishing introgression from ancestral polymorphism
Algorithmic Components	Forward-Backward algorithm	Computes posterior state probabilities	Estimating evolutionary history probabilities at each genomic position
	Viterbi algorithm	Finds optimal state path	Identifying most likely evolutionary history across genomic regions
	Baum-Welch algorithm	Estimates HMM parameters from data	Learning evolutionary parameters without labeled training data
Model Parameters	Transition probabilities	Govern switches between evolutionary histories	Modeling recombination between genomic regions with different histories
	Emission probabilities	Define symbol distributions for each state	Representing mutation patterns under different evolutionary scenarios
	Initial state probabilities	Specify starting state distributions	Incorporating prior knowledge about evolutionary histories

Advanced Methodological Considerations

Handling Incomplete Lineage Sorting

A critical challenge in introgression research involves distinguishing true introgression from incomplete lineage sorting (ILS), where ancestral polymorphisms persist through speciation events [4]. The PhyloNet-HMM framework addresses this by explicitly modeling both processes, enabling more accurate detection of introgression. The method incorporates the multispecies coalescent model to account for ILS, preventing false positive introgression signals that might otherwise arise from ancestral population variation.

Modeling Genomic Dependencies

Biological sequences exhibit complex dependencies due to recombination, selection, and functional constraints. HMMs capture local dependencies through their Markov structure, but additional modeling techniques may be required for more complex dependency patterns. Context-sensitive HMMs and related extensions provide frameworks for capturing longer-range dependencies in genomic sequences [12], which can be particularly important for modeling evolutionary processes that operate across different genomic scales.

Validation and Benchmarking Strategies

Robust validation is essential for reliable introgression detection. Recommended approaches include:

• Simulation studies: Generating synthetic genomic data under known evolutionary scenarios to evaluate method performance [4]
• Negative controls: Applying methods to datasets where introgression is not expected [4]
• Biological validation: Testing predictions against independent experimental evidence
• Cross-validation: Assessing parameter stability across different genomic regions or subsets of data

These strategies help ensure that detected introgression signals reflect true biological processes rather than methodological artifacts or random variation.

From Theory to Practice: Implementing HMMs for Introgression Detection

Hidden Markov Models (HMMs) have emerged as powerful computational frameworks for detecting introgression—the transfer of genetic material between species through hybridization. This application note details two significant HMM-based frameworks: PhyloNet-HMM, designed for detecting introgression in eukaryotes by combining phylogenetic networks with HMMs, and Ancestry_HMM-S, which builds upon local ancestry inference methods to detect and quantify adaptive introgression. These tools address the critical challenge of distinguishing true introgression signals from spurious ones caused by evolutionary processes like incomplete lineage sorting (ILS) [4] [23]. The development of these methods represents a substantial advancement in comparative genomics, enabling researchers to systematically analyze genome-wide data to uncover the role of hybridization in evolution, adaptation, and species diversification.

The table below summarizes the core characteristics, applications, and requirements of PhyloNet-HMM and Ancestry_HMM-S.

Table 1: Comparative Overview of PhyloNet-HMM and Ancestry_HMM-S

Feature	PhyloNet-HMM	Ancestry_HMM-S
Primary Function	Detects introgressed genomic regions in eukaryotes [4]	Infers adaptive introgression and quantifies selection strength [24] [25]
Core Innovation	Integrates phylogenetic networks with HMMs [4]	Extends a local ancestry inference framework to model selection [24]
Key Advantage	Accounts for incomplete lineage sorting (ILS) and dependence across loci [4] [26]	Works on standard population genomic data sets and models selection coefficients [24]
Typical Input	Multiple aligned genomes from multiple species [4]	Genomic data (e.g., VCF) and a ploidy file from an admixed population [24]
Required Parameters	Set of parental species trees [4]	Effective population size (--ne), time and fraction of introgression (-p), and analysis mode [24]
Software Availability	Available as a Jar file or tarball [27]	Available via Bioconda or GitHub [24]

PhyloNet-HMM: Protocol for Detecting Introgression

Workflow and Interpretation

PhyloNet-HMM operates by scanning aligned genomes to calculate the probability that each site evolved under a specific parental species tree, which is represented within a phylogenetic network that models both vertical descent and hybridization [4]. The following diagram illustrates the core logical workflow of the PhyloNet-HMM framework.

Application Notes and Key Findings

In its foundational study, PhyloNet-HMM was applied to variation data from chromosome 7 of the house mouse (Mus musculus domesticus) [4] [23] [26]. The analysis successfully detected a previously known adaptive introgression event involving the Vkorc1 gene, which confers resistance to rodent poison [23] [26]. Beyond this validation, the analysis revealed that approximately 9% of sites on chromosome 7 (covering about 13 Mbp and over 300 genes) were of introgressive origin, a finding that significantly expanded the understood scope of introgression in this model organism [4] [26]. The framework's robustness was confirmed by the absence of false positive introgression signals in a negative control data set [4] [23].

Ancestry_HMM-S: Protocol for Detecting Adaptive Introgression

Workflow and Analysis Modes

Ancestry_HMM-S is designed to identify loci undergoing adaptive introgression by analyzing the local ancestry patterns in admixed populations. Its workflow requires specific population genetic parameters and offers different analysis modes.

Required Parameters and Usage

A typical command for running Ancestry_HMM-S involves specifying several key parameters [24]:

• Population Genetic Parameters: The effective population size of the admixed population is set with --ne.
• Introgression Pulse: The time and fraction of introgression are specified using the -p flag. For example, -p 1 100000 0.9 -p 0 100 0.1 models a 10% introgression pulse from population 0 into population 1 that occurred 100 generations ago.
• Analysis Mode: The --gss (Golden Section Search) mode is generally recommended for inferring adaptive introgression from real data, as it finds the selection coefficient at each site that yields the highest likelihood ratio.
• Performance Flags: Using flags like --traj 4 (for a faster 4-point approximation method) and --window p 10 (to limit the Markov chain to 10% of the chromosome length on each side of the focal site) is recommended to improve runtime without significant loss of accuracy [24].

Application Notes and Key Findings

Ancestry_HMM-S has been extensively validated and applied to real data sets. When applied to an admixed Drosophila melanogaster population from South Africa, the method identified 17 loci with signatures of adaptive introgression. Notably, four of these loci were previously known to confer resistance to insecticides, confirming the method's power to uncover biologically relevant adaptive events [25].

The Scientist's Toolkit

Table 2: Essential Research Reagents and Computational Tools

Item Name	Function/Description	Relevance in HMM-Based Introgression Studies
PhyloNet Software Package	A platform for evolutionary analysis and phylogenetic network inference [4].	PhyloNet-HMM is distributed as part of the open-source PhyloNet package, making it the essential software environment for running the analyses [4] [27].
Ancestry_HMM-S Software	A specialized program for inferring adaptive introgression [24].	The core analytical tool itself, available via Bioconda or by compilation from source code on GitHub [24].
Armadillo C++ Library	A high-performance linear algebra library [24].	A critical dependency that must be installed on the system for compiling and running Ancestry_HMM-S from source [24].
Multiple Sequence Alignment	A data structure representing aligned genomic sequences from multiple individuals/species [4].	The primary input format for PhyloNet-HMM, providing the raw comparative genomic data for analysis [4].
VCF File	Variant Call Format file storing genetic variation data [24].	A standard input format for population genomic tools like Ancestry_HMM-S, containing genotype information for the admixed population and its ancestors [24].
Parental Species Trees	A set of trees representing the non-reticulate evolutionary histories of the species involved [4].	A required input for PhyloNet-HMM, constraining the possible evolutionary histories during the search for introgressed regions [4].

Modeling Evolutionary Histories with Phylogenetic Networks

Evolutionary biology has traditionally relied on phylogenetic trees to represent the divergence history of species. However, a growing body of evidence indicates that reticulate evolutionary processes—including hybridization, introgression, and horizontal gene transfer—are pervasive across the Tree of Life [28]. Phylogenetic networks provide a powerful framework for modeling these complex histories, moving beyond the limitations of strictly tree-like representations. When modeling introgression—the integration of genetic material from one species into another through hybridization—researchers must account for confounding signals from processes like incomplete lineage sorting (ILS), where different genomic regions have conflicting genealogies due to ancestral polymorphism [4].

The integration of hidden Markov models (HMMs) with phylogenetic networks has emerged as a particularly powerful approach for detecting introgression from genomic data. These methods can distinguish true introgression signals from spurious ones that arise due to population effects, while simultaneously accounting for dependencies across genomic loci [4]. This application note details protocols for applying these computational frameworks to identify introgressed genomic regions, with specific applications spanning from pesticide resistance in mice to archaic hominin admixture in humans.

Computational Frameworks for Introgression Detection

Key Methodological Approaches

Table 1: Comparative Overview of Phylogenetic Network Methods for Introgression Detection

Method Name	Statistical Framework	Key Features	Targeted Evolutionary Processes	Reference Applications
PhyloNet-HMM	Phylogenetic networks + HMM	Accounts for ILS and dependency across loci; models local genealogies	Hybridization, introgression, recombination	Mouse chromosome 7 analysis (Vkorc1 gene) [4]
Ancestry_HMM-S	Hidden Markov Model	Quantifies selection strength on introgressed regions; models ancestry tract lengths	Adaptive introgression, selection in admixed populations	Drosophila pesticide resistance [29]
ArchaicSeeker 2.0	HMM + likelihood framework	Models multiple-wave admixture; detects unknown archaic lineages	Archaic hominin introgression, multiple admixture waves	Eurasian and Oceanian populations [30]
NetPlacer	Phylogenetic network placement	Places queries into phylogenetic networks rather than trees	Non-tree-like evolution, microbial evolution	Genome and metagenome placement [31]

Table 2: Key Computational Tools and Resources for Phylogenetic Network Analysis

Tool/Resource	Type	Primary Function	Application Context
PhyloNet	Software package	Inference and analysis of phylogenetic networks	Implementation of PhyloNet-HMM for genome-wide introgression scanning [4]
Ancestry_HMM-S	Software tool	Local ancestry inference with selection modeling	Detection and quantification of adaptive introgression [29]
ArchaicSeeker 2.0	Computational method	Archaic sequence detection and admixture history modeling	Multiple-wave archaic hominin introgression analysis [30]
Normal networks	Network class	Biologically relevant network representation with mathematical tractability	Framework for eukaryotic species evolution reconstruction [32]
Tree-child networks	Network class	Phylogenetic networks where every internal node has a tree child	Foundation for normal networks; ensures vertex visibility [32]

Experimental Protocols and Workflows

Protocol 1: Genome-Wide Introgression Scanning with PhyloNet-HMM

Purpose: To detect introgressed genomic regions while accounting for incomplete lineage sorting and dependencies across loci.

Input Requirements:

• A set of aligned genomes, each of length L
• A set of parental species trees representing potential evolutionary histories
• Genomic data in aligned format with appropriate annotation

Methodological Steps:

• Data Preparation and Alignment

  • Obtain whole-genome sequences from target species and outgroups
  • Perform multiple sequence alignment using preferred aligner (e.g., MAFFT, MUSCLE)
  • Verify alignment quality and address potential artifacts

• Parental Species Tree Specification

  • Define possible evolutionary histories based on established phylogenies
  • For the mouse protocol example [4], specify two parental trees corresponding to alternate evolutionary scenarios
  • Incorporate known speciation and potential hybridization events

• PhyloNet-HMM Configuration

  • Define hidden states corresponding to different local genealogies
  • Set transition probabilities between states based on recombination rates
  • Configure emission probabilities based on evolutionary substitution models

• Model Training and Execution

  • Implement dynamic programming algorithms for efficient computation [4]
  • Employ multivariate optimization heuristics for parameter estimation
  • Execute the HMM to compute probabilities for each site: P(Ψ_i = ψ | X) for every site i and parental species tree ψ [4]

• Result Interpretation

  • Identify genomic regions with significant probability of introgression
  • Calculate the proportion of sites of introgressive origin
  • Annotate genes within introgressed regions for functional analysis

Validation Approach:

• Apply method to negative control data sets where no introgression is expected
• Validate known introgressed regions (e.g., Vkorc1 in mice [4])
• Implement simulation studies under coalescent models with recombination, isolation, and migration

Figure 1: PhyloNet-HMM workflow for genome-wide introgression detection.

Protocol 2: Detecting Adaptive Introgression with Ancestry_HMM-S

Purpose: To identify genes undergoing adaptive introgression and quantify the strength of selection acting on them.

Input Requirements:

• Genomic data from an admixed focal population
• Two unadmixed, ancestral reference populations
• Optionally, low-coverage data from unphased diploid samples or pooled sequencing data

Methodological Steps:

• Data Collection and Preprocessing

  • Sequence or obtain genomic data from admixed population and reference populations
  • For the Drosophila example [29], collect data from South African admixed populations and ancestral references
  • Perform standard variant calling and quality filtering

• Neutral Admixture Modeling

  • Implement HMM framework for local ancestry inference
  • Infer timing of multiple admixture pulses based on ancestry patterns
  • Establish baseline introgression fraction across the genome

• Selection Detection Parameters

  • Identify loci exceeding baseline introgression fractions significantly
  • Analyze ancestry tract lengths - longer than expected tracts suggest selection
  • Model the impacts of natural selection during admixture

• Selection Strength Quantification

  • Estimate selective coefficients for candidate loci
  • Calculate statistical significance of adaptive introgression signals
  • Implement false discovery rate controls for multiple testing

• Functional Annotation

  • Annotate genes showing signatures of adaptive introgression
  • For the Drosophila analysis [29], focus on pesticide resistance genes (Cyp6 genes, Ace)
  • Perform enrichment analysis for biological pathways

Validation Approach:

• Extensive forward simulations under various selection scenarios
• Comparison with known adaptive loci (e.g., insecticide resistance genes)
• Replication in independent population samples

Figure 2: Ancestry_HMM-S workflow for detecting adaptive introgression.

Protocol 3: Multiple-Wave Archaic Introgression Detection with ArchaicSeeker 2.0

Purpose: To identify introgressed hominin sequences and model multiple-wave admixture from archaic hominins into modern humans.

Input Requirements:

• Modern human genomic data from target populations
• Reference archaic hominin genomes (Neanderthal, Denisovan)
• African genomes as non-introgressed outgroups

Methodological Steps:

• Data Preparation and Quality Control

  • Obtain modern human genomic data from worldwide populations
  • Integrate high-quality archaic hominin reference genomes
  • Implement stringent filtration procedures to remove potential artifacts

• HMM Configuration for Archaic Sequence Detection

  • Implement HMM to describe genomes of mixed ancestries
  • Define states representing archaic vs. modern human ancestry
  • Set transition probabilities based on recombination maps

• Ancestral Source Determination

  • Apply likelihood-based approach to trace ancestral sources
  • Classify introgressed sequences into different archaic lineages
  • Identify sequences from potentially unknown archaic hominins

• Multiple-Wave Admixture Modeling

  • Implement general discrete admixture model (MultiWaver methodology)
  • Infer timing and magnitude of distinct admixture waves
  • Estimate divergence times between modern humans and archaic sources

• Functional and Phenotypic Analysis

  • Annotate introgressed regions with functional genomic elements
  • Test for associations with phenotypic traits (e.g., immune function, UV response)
  • Identify "archaic deserts" - regions depleted of archaic ancestry

Performance Metrics:

• Precision: 93.0% (95% CI, 89.4-95.9%) based on length-based comparison [30]
• True-positive rate: 90.4% (95% CI, 84.1-94.1%) [30]
• False-positive rate: 0.14% (95% CI, 0.07-0.22%) [30]

Figure 3: ArchaicSeeker 2.0 workflow for archaic introgression analysis.

Key Applications and Case Studies

Rodenticide Resistance in House Mice (PhyloNet-HMM Application)

Biological Context: The analysis of chromosome 7 in the house mouse (Mus musculus domesticus) revealed a previously reported adaptive introgression event involving the rodent poison resistance gene Vkorc1 [4]. This introgression event likely contributed to the spread of resistance against warfarin and other rodenticides.

Methodological Implementation:

• Applied PhyloNet-HMM to variation data from chromosome 7
• Scanned for regions with signatures of introgression
• Accounted for incomplete lineage sorting and dependencies across loci

Key Findings:

• Approximately 9% of sites within chromosome 7 were of introgressive origin
• This covered about 13 Mbp of chromosome 7 and over 300 genes
• Successfully detected the Vkorc1 region as introgressed
• No introgression detected in negative control data sets, demonstrating specificity

Pesticide Resistance in Drosophila (Ancestry_HMM-S Application)

Biological Context: Analysis of an admixed Drosophila melanogaster population from South Africa with ancestry from both ancestral and cosmopolitan populations [29]. This population has experienced recent admixture (~30 years prior to sampling), providing an opportunity to study adaptive introgression.

Methodological Implementation:

• Implemented Ancestry_HMM-S to detect adaptive introgression
• Modeled selection strength on introgressed haplotypes
• Analyzed ancestry tract length distributions

Key Findings:

• Identified 17 loci showing signatures of adaptive introgression
• Four loci contained genes previously shown to confer insecticide resistance
• Cosmopolitan haplotypes carrying insecticide resistance alleles showed elevated frequencies
• Demonstrated that selection has shaped ancestry patterns in recently admixed populations

Archaic Hominin Introgression in Modern Humans (ArchaicSeeker 2.0 Application)

Biological Context: Modern human populations outside of Africa contain 1-3% archaic hominin ancestry from Neanderthals and Denisovans [30]. The timing, number, and locations of these admixture events remain active areas of research.

Methodological Implementation:

• Applied ArchaicSeeker 2.0 to worldwide modern human populations
• Modeled multiple-wave admixture from different archaic sources
• Detected introgressed sequences without stringent length filters

Key Findings:

• Identified two waves of introgression from both Denisovan-like and Neanderthal-like hominins in Eurasian populations
• Estimated early Denisovan-like introgression occurred ~118.8-94.0 kya
• Found only one episode of Denisovan-like admixture in indigenous peoples east of the Wallace Line
• Suggested early dispersal of modern humans throughout Asia before the Toba eruption (74 kya)
• Identified archaic sequences involved in immune function, cardiovascular traits, UV response, and metabolism

Technical Considerations and Best Practices

Network Class Selection for Biological Relevance

When applying phylogenetic networks to evolutionary questions, the selection of appropriate network classes is crucial for biological interpretability:

• Normal networks are emerging as a leading class for eukaryotic evolution modeling, as they provide mathematical tractability while maintaining biological relevance [32]
• Tree-child networks ensure that every internal vertex has at least one child that is a tree vertex or leaf, providing the important property of "visibility" where evolutionary signals from vertices are preserved in leaf data [32]
• Binary hybridization networks, which are particularly relevant for modeling introgression in eukaryotic species, have been shown to be normal networks [32]

Handling Methodological Challenges

Incomplete Lineage Sorting: PhyloNet-HMM and related methods specifically address the challenge of distinguishing true introgression from spurious signals caused by ILS [4]. This is particularly important for closely related species or recent divergences.

Multiple-Wave Admixture: ArchaicSeeker 2.0 addresses the limitation of methods that can only detect single-wave introgression by implementing a general discrete admixture model that can infer complex admixture histories [30].

Computational Scalability: Methods like NetPlacer address the computational challenges of moving from tree-based to network-based phylogenetic placement, though further scalability improvements are needed [31].

Validation and Performance Assessment

Robust validation is essential for introgression detection methods:

• Simulation studies: Methods should be validated using forward simulations under realistic demographic scenarios [29] [30]
• Negative controls: Application to data sets where no introgression is expected provides specificity measures [4]
• Positive controls: Known introgressed regions (e.g., Vkorc1 in mice) provide validation benchmarks [4]
• Performance metrics: Precision, true-positive rates, and false-positive rates should be reported using length-based, SNP-based, and segment-based comparisons [30]

Integrating HMMs to Account for Incomplete Lineage Sorting (ILS)

Incomplete lineage sorting (ILS) is a prevalent evolutionary phenomenon that occurs when ancestral genetic polymorphisms persist through rapid speciation events, leading to incongruences between gene trees and the overarching species tree. This process generates widespread genomic discordance, which can obscure phylogenetic relationships and complicate the identification of true evolutionary histories. Hidden Markov Models (HMMs) offer a powerful statistical framework for addressing this complexity by modeling the variation in genealogical histories across genomic alignments. Their application allows researchers to distinguish between genuine signatures of introgression and spurious signals resulting from ILS, thereby providing a more accurate reconstruction of evolutionary pathways. This protocol details the integration of HMMs for detecting and accounting for ILS in comparative genomic studies, with a specific focus on applications in introgression research.

Background and Quantitative Significance of ILS

Incomplete lineage sorting is not a rare occurrence but a major force in eukaryotic evolution. Recent empirical studies across diverse taxa have quantified its substantial impact on genomes.

• Hominids: Over 30% of the human genome supports conflicting phylogenetic trees due to ILS, affecting numerous genes with known morphological functions [33].
• Marsupials: Phylogenomic analyses reveal that over 50% of marsupial genomes are affected by ILS. In the South American monito del monte, a key species in marsupial evolution, over 31% of its genome is closer to the Diprotodontia order due to ILS during ancient radiation, rather than reflecting the true species tree [34].

The following table summarizes the quantitative evidence for ILS across different clades, highlighting its pervasive nature:

Table 1: Empirical Evidence of Genome-Wide Incomplete Lineage Sorting

Clade	Key Study Species	Genomic Region Analyzed	Percentage of Genome Affected by ILS	Key Reference
Hominids (Great Apes)	Homo sapiens	Genome-wide	> 30%	[33]
Marsupials	Dromiciops (Monito del Monte)	Genome-wide	> 50% (with >31% showing specific discordance)	[34]
Mouse	Mus musculus domesticus	Chromosome 7	~9% of sites (covering ~13 Mbp)	[4]

The phenotypic consequences of ILS, often misinterpreted as convergent evolution, are significant. Incongruent traits are frequently found in the craniofacial and appendicular skeletons, underscoring the importance of accurately modeling ILS to avoid misinterpreting evolutionary history [33].

Core Protocol: PhyloNet-HMM for ILS and Introgression Analysis

This protocol centers on PhyloNet-HMM, a comparative genomic framework that integrates phylogenetic networks with Hidden Markov Models to detect introgression while explicitly accounting for ILS [4]. Its primary strength lies in its ability to simultaneously model the reticulate evolutionary history of genomes and the dependencies between loci.

Input Data Requirements and Preparation

• Genomic Data: A multiple sequence alignment (MSA) of genomes from the species of interest. The alignment should be in a standard format (e.g., FASTA, MAF).
• Species Tree Hypothesis: A rooted species tree (or a set of candidate trees) representing the presumed evolutionary relationships without considering introgression.
• Phylogenetic Network Hypothesis (Optional but Recommended): A hypothesized phylogenetic network that includes proposed introgression events. This guides the model but is not strictly required, as the method can infer these events.

Step-by-Step Workflow

Table 2: PhyloNet-HMM Experimental Protocol

Step	Procedure	Parameters & Notes
1. Data Curation	Compile and align whole-genome sequences for the taxa of interest. Ensure data quality and format consistency.	Use tools like LASTZ for alignment. Curate data to minimize missing data and sequencing errors [34].
2. Model Training	Train the PhyloNet-HMM model on the aligned genomes and the provided species tree/network.	The model uses dynamic programming and multivariate optimization to learn the transition and emission probabilities [4].
3. Genome Scanning	The trained model scans the genome alignment site-by-site (or in windows).	For each site (i), it calculates the probability (P(\Psii \mid Alignment)), where (\Psii) is the parental species tree for that site [4].
4. Result Interpretation	Identify genomic regions with strong support for a phylogenetic history that differs from the species tree.	Regions where the probability for a discordant tree (e.g., due to introgression) is high are flagged as candidate introgressed regions.
5. Validation	Perform functional experiments or independent phylogenetic analyses on candidate regions to validate findings.	As done in marsupial studies, functional tests can confirm the phenotypic effects suggested by ILS [34].

The logical workflow of the PhyloNet-HMM analysis is outlined below:

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Resources for HMM-based ILS and Introgression Studies

Category	Item / Software	Function and Application
Computational Tools	PhyloNet-HMM [4]	The core framework for detecting introgression in the presence of ILS using HMMs.
	PhyloNet Software Package [4]	A broader suite of tools for phylogenetic network analysis, which includes PhyloNet-HMM.
	CoalHMM [34]	A related HMM-based method used for estimating population genetic parameters from genomic data.
Data Resources	NCBI Sequence Read Archive (SRA) [34]	Primary repository for raw genomic sequencing data.
	CNSA of CNGBdb [34]	Alternative repository for archiving and accessing genomic data sets.
Biological Materials	Reference Genomes	High-quality genomes for the species clade of interest (e.g., Hominids, Marsupials). Essential for accurate alignment and analysis.
	Voucher Specimens	Well-identified tissue or cell lines from which genomic data is derived, ensuring taxonomic accuracy.

Visualization of Evolutionary Relationships and Model Architecture

A core challenge is distinguishing between the signals of ILS and introgression, as both can produce similar patterns of genealogical discordance. The following diagram illustrates the evolutionary scenarios that HMM-based methods are designed to decipher.

The PhyloNet-HMM model architecture is designed to parse these complex signals. It functions as a classic HMM where the hidden states represent different parental species trees (or genealogies), and the observed states are the aligned genomic sequences.

Application Notes and Future Directions

Integrating HMMs to account for ILS is a cornerstone of modern comparative genomics. The PhyloNet-HMM framework provides a statistically robust method to dissect the intertwined effects of deep ancestral polymorphism and reticulate evolution. Future advancements will likely involve tighter integration with machine learning approaches, such as the variational autoencoders used in multi-omics data integration [35], to handle ever-larger and more complex genomic datasets. Furthermore, as demonstrated in hominid and marsupial research, bridging genomic findings with phenotypic data through functional experiments is the critical next step for moving beyond correlation to establish causative links in evolutionary biology [33] [34]. This protocol establishes a foundation for such integrative research, empowering scientists to unravel the complex network of life.

The genomic landscapes of admixed populations provide a powerful record of evolutionary selection pressures, including those exerted by insecticides. Adaptive introgression, the process by which advantageous genetic variants are incorporated into a population's gene pool through hybridization and back-crossing, is a major source of adaptive genetic variation [36]. In the context of pest species like Drosophila melanogaster, detecting signatures of selection in admixed genomes can reveal loci underlying resistance to commonly used insecticides such as organophosphates and pyrethroids [37].

This application note details the use of Ancestry_HMM-S (AHMMS), a computational tool designed to infer adaptive introgression from population genomic data by modeling local ancestry as a hidden Markov model (HMM) [24]. We demonstrate its application in detecting insecticide resistance loci in Drosophila melanogaster, providing a validated protocol for researchers investigating genomic adaptations in pest populations.

Background

Adaptive Introgression and Selection

Admixture can simultaneously introduce multiple selected alleles into a population. Unlike single mutation events, this process can create circumstances where linked selected sites significantly affect each other's evolutionary dynamics through genetic hitchhiking or Hill-Robertson interference [38]. Traditional methods for detecting selection in admixed populations often only account for selection at single sites, potentially overlooking the effects of linkage among multiple selected alleles and leading to overestimation of both the number of selected sites and their selection coefficients [38].

Insecticide Resistance Mechanisms in Drosophila

Resistance to insecticides in Drosophila melanogaster involves two primary mechanisms: target-site resistance and metabolic resistance [37].

• Target-site resistance involves mutations in the protein targeted by the insecticide, reducing binding affinity. For organophosphates (OPs), the target is acetylcholinesterase (Ace), while for pyrethroids, it is the voltage-gated sodium channel (Vssc) [37].
• Metabolic resistance involves increased production or activity of detoxification enzymes that break down insecticides before they reach their targets, primarily cytochrome P450s (Cyp genes), carboxylesterases, and glutathione S-transferases [37] [39].

Table 1: Key Insecticide Resistance Genes in Drosophila melanogaster

Insecticide Class	Target Site Genes	Metabolic Genes	Primary Resistance Mechanism
Organophosphates (e.g., parathion)	Ace (Acetylcholinesterase)	Cyp6a17, Cyp6a23	Target-site mutations [37]
Pyrethroids (e.g., deltamethrin)	Vssc (Voltage-gated sodium channel)	Cyp6a23, Cyp6a17, CG7627	Metabolic resistance via P450s [37]

AHMM-S Methodology and Workflow

AHMM-S extends the Ancestry_HMM framework by explicitly modeling selection during admixture in sequence alignment data. The method calculates the expected local ancestry landscape in an admixed population for a given selection model, then maximizes the likelihood of that model given the observed data [38] [24]. The hidden states of the HMM represent the local ancestries at ancestry-informative positions, while the observed states are the alleles in aligned reads at these positions [38].

Key innovations in AHMM-S include:

• Modeling the effects of selection on ancestry tract lengths, which tend to be longer around selected sites due to hitchhiking effects
• Introducing new techniques to infer expected transition probabilities between ancestry types at adjacent sites under generalized models of selection
• Using a numerical method to track the expected distribution of haplotypes after t generations of admixture [38]

Analytical Workflow

The following diagram illustrates the complete AHMM-S analysis workflow for detecting selective sweeps associated with insecticide resistance:

Diagram 1: AHMM-S analytical workflow (76 characters)

Protocol: Detecting Insecticide Resistance Loci

Software Installation and Dependencies

AHMMS Installation via Bioconda (Recommended)

Dependencies

• C++ linear algebra library Armadillo (libarmadillo-dev)
• Google performance tools (google-perftools) for decreased runtime [24]

Input File Preparation

AHMMS requires three primary inputs:

• Genomic Data File (-i filename): Contains variant information across genomes
• Ploidy File (-s filename): Specifies ploidy of all individuals in the datafile
• Population Parameters: Effective population size (--ne) and admixture history (-p) [24]

Table 2: Required Input Parameters for AHMM-S Analysis

Parameter	Format Example	Description	Note
Population size	--ne 1000	Effective population size of admixed population	Critical for accurate modeling
Admixture pulse	-p 1 100000 0.9 -p 0 100 0.1	Donor (0) and recipient (1) populations with time and proportion	Specify population 1 BEFORE population 0 [24]
Analysis window	--chr_win 1 5000	Limit region on chromosome for analysis	Optional for focused analysis
Markov chain window	--window p 10	Extend Markov chain 10% of chromosome length each side of selected site	Recommended for performance [24]

Command Execution and Analysis Modes

AHMMS offers three analysis modes, with Golden Section Search (GSS) recommended for inferring adaptive introgression from real data:

Table 3: AHMM-S Analysis Mode Comparison

Mode	Command Flag	Parameters	Use Case
Golden Section Search	--gss	WINDOWSTART, WINDOWEND, STEP, MINS, MAXS	Recommended for real data analysis [24]
Grid Search	--grid	START, STOP, STEP, SSTART, SSTOP, S_STEP	Comprehensive likelihood landscape
Single Site	--site	POSITION, S_VALUE	Testing specific site/coefficient

Critical Parameter Notes:

• For GSS, set MINS and MAXS to 0.001 and 0.15 as recommended starting values [24]
• Use --traj 4 to implement the fast and accurate 4-point approximation method
• The --window p 10 flag speeds computation by limiting the Markov chain to 10% of chromosome length each direction from the focal site [24]

Interpretation of Output

AHMMS generates likelihood ratios indicating sites with signatures of selection. Key outputs include:

• Genomic Positions with significantly elevated likelihood ratios
• Estimated selection coefficients (s) for each candidate locus
• Local ancestry patterns surrounding candidate loci

Candidate resistance loci should be cross-referenced with known resistance genes (Ace, Vssc, Cyp genes) and validated experimentally.

Research Reagent Solutions

Table 4: Essential Research Reagents and Resources for Drosophila Insecticide Resistance Studies

Reagent/Resource	Function/Application	Example Use in Resistance Studies
AHMM-S Software	Inference of adaptive introgression from genomic data	Detecting selected loci in admixed Drosophila populations [24]
DGRP (Drosophila Genetic Reference Panel)	GWAS panel of fully sequenced, mostly homozygous lines	Identifying resistance variants via genome-wide association [37]
Organophosphate Insecticides (e.g., parathion)	Selecting for resistance alleles	Bioassays to quantify resistance levels in Drosophila lines [37]
Pyrethroid Insecticides (e.g., deltamethrin)	Selecting for resistance alleles	Contact bioassays in coated vials to assess resistance [37]
UAS-RNAi Lines	Gene knockdown functional validation	Testing candidate gene involvement in resistance (e.g., Cyp6a23) [37]
Gal4 Drivers (e.g., Actin5c-Gal4)	Targeted gene expression	Conditional expression of resistance candidate genes [37]

Validation and Experimental Design

Complementary Wet-Lab Validation

Bioinformatic predictions from AHMM-S require experimental validation:

Insecticide Bioassays

• Use the residual contact application method with insecticides dissolved in acetone
• Coat scintillation vials with 0.5ml insecticide solution (e.g., 1.5 µg/ml parathion, 0.7 µg/ml deltamethrin)
• Expose 5-8 day old adult males (n=20 per line) to treated vials [37]

Functional Validation via Gene Knockdown

• Cross UAS-RNAi lines with ubiquitous drivers (e.g., Actin5c-Gal4)
• Assess resistance phenotypes in F1 progeny
• Confirm involvement of candidate genes like Cyp6a23 in pyrethroid resistance [37]

Case Study: Resistance Loci Detection

A GWAS study using the Drosophila Genetic Reference Panel revealed:

• Parathion (OP) resistance: Primarily associated with mutations in the target gene Ace
• Deltamethrin (pyrethroid) resistance: Strongly associated with Cyp6a23, a detoxification enzyme [37]

These findings demonstrate how AHMM-S can pinpoint both target-site and metabolic resistance mechanisms, providing a comprehensive view of resistance genetics in Drosophila populations.

Ancestry_HMM-S provides an powerful methodological framework for detecting genomic signatures of insecticide resistance in admixed Drosophila populations. By modeling local ancestry patterns under selection, researchers can identify both known and novel resistance loci, enabling more comprehensive monitoring of resistance evolution in pest populations. The integration of computational predictions with experimental validation creates a robust pipeline for characterizing the genetic architecture of insecticide resistance.

Introgression, the stable incorporation of genetic material from one species into the gene pool of another through hybridization and repeated back-crossing, is a recognized evolutionary force with significant implications for adaptation and speciation [23]. Detecting these introgressed regions in genomic data is computationally challenging, as their signatures can be confounded by other evolutionary processes, most notably Incomplete Lineage Sorting (ILS), which arises from ancestral population polymorphism and recombination [23].

This application note details a case study utilizing the PhyloNet-HMM framework to detect introgression in the mouse genome (Mus musculus domesticus), with a focus on the adaptive introgression of the Vkorc1 locus, a gene conferring resistance to rodenticides [23]. We provide a detailed protocol for applying this comparative genomic method, which synergistically combines phylogenetic networks with hidden Markov models (HMMs) to distinguish introgressive origins from ILS while accounting for dependencies along the genome [23].

Background

The PhyloNet-HMM Framework

PhyloNet-HMM is designed to infer the evolutionary history of genomes where hybridization and ILS are both at play. Its core innovation lies in modeling a walk along an alignment of multiple genomes, where the hidden states are local genealogies that can change due to recombination or introgression [23]. The model observes the genomic sequence data and infers the most probable sequence of underlying genealogies, some of which may signal an introgression event.

The framework is an extension of the Coalescent Hidden Markov Model (CoalHMM) concept [40]. In a standard CoalHMM, the hidden states represent different coalescent histories, and transitions between states are governed by population parameters and the recombination rate [40]. PhyloNet-HMM enhances this by expanding the set of possible hidden states to include genealogies that are consistent with a given phylogenetic network (which represents hybridization) rather than just a species tree [23].

The Vkorc1 Gene and Adaptive Introgression

The Vkorc1 gene encodes the vitamin K epoxide reductase complex subunit 1, the molecular target of warfarin and other anticoagulant compounds [41] [42]. In humans, polymorphisms in VKORC1 are a major genetic determinant of warfarin dosing sensitivity [43]. In mice, adaptive introgression has been reported where alleles of Vkorc1 from a donor species, conferring resistance to rodenticides, introgressed into M. m. domesticus populations, providing a selective advantage in environments where poisons were deployed [23]. This makes it an ideal candidate for testing and demonstrating the PhyloNet-HMM method.

Methodology/Experimental Protocol

Computational Workflow

The following diagram illustrates the logical workflow for applying PhyloNet-HMM to detect introgression, as implemented in this case study.

Step-by-Step Protocol

Step 1: Input Data Preparation

• Objective: Generate a multiple whole-genome alignment for the species or populations under study.
• Procedure:
  • Select Taxa: For the mouse case study, this included:
    • The focal species (M. m. domesticus)
    • The putative donor species
    • One or more closely related outgroup species.
  • Sequence Alignment: Assemble whole-genome sequencing data for all taxa. Perform a reference-based or de novo multiple sequence alignment, for example, using tools like MAFFT or MUSCLE. The output should be a single, column-wise alignment file (e.g., in FASTA or MAF format).
• Notes: Data quality is critical. Ensure high coverage and perform standard variant calling and filtering to minimize artifacts.

Step 2: Phylogenetic Network Hypothesis

• Objective: Define the putative evolutionary history, including a hybridization event, to be tested by the model.
• Procedure:
  • Infer a Species Tree: Use a standard phylogenomic method on non-introgressed, non-linked regions to establish the primary species tree.
  • Postulate Hybridization: Based on prior biological knowledge (e.g., from population statistics like D-statistics), specify the hybridizing lineages and the direction of introgression. For the mouse study, the network involved M. m. domesticus and a closely related species as the donor [23].

Step 3: PhyloNet-HMM Configuration

• Objective: Set up the HMM with states that reflect the possible genealogies under the network.
• Procedure:
  • State Definition: The hidden states are the set of possible local genealogies. Under a network, this includes trees that follow the primary species tree and trees that reflect the introgressive history [23].
  • Parameter Initialization:
    • Emission Probabilities: Calculate the probability of observing a column in the alignment given a specific local genealogy. This requires a nucleotide substitution model (e.g., HKY or GTR) and associated parameters [40].
    • Transition Probabilities: Define the probabilities of switching from one genealogy to another along the alignment. These are a function of the population parameters (ancestral effective population sizes, speciation times) and the per-site recombination rate (ρ) [23] [40].

Step 4: Model Training and Likelihood Computation

• Objective: Find the model parameters that best explain the observed genome alignment.
• Procedure:
  • Use an Expectation-Maximization (EM) algorithm or a similar optimization technique to estimate the model parameters (e.g., population sizes, recombination rate) that maximize the likelihood of the data.
  • The Forward-Backward algorithm is used to compute the likelihood of the alignment given the current parameter estimates [40].

Step 5: Decoding and Path Prediction

• Objective: Infer the most likely sequence of genealogical histories along the chromosomes.
• Procedure:
  • Apply the Viterbi algorithm to find the single most likely path of hidden states (genealogies) through the genome.
  • Alternatively, compute posterior probabilities for each state at every site using the Forward-Backward algorithm. A site can be assigned to the genealogy with the highest posterior probability [40].

Step 6: Identification of Introgressed Regions

• Objective: Translate the predicted state path into biologically interpretable regions of introgression.
• Procedure:
  • Annotation: Genomic intervals where the predicted state corresponds to the genealogy indicative of introgression are flagged as candidate introgressed regions.
  • Summary: Calculate summary statistics such as the total number and combined length of introgressed segments, and their genomic density.

Results and Data Analysis

Key Findings from the Mouse Chromosome 7 Analysis

Application of PhyloNet-HMM to chromosome 7 of M. m. domesticus successfully recapitulated and expanded upon known introgression events.

Table 1: Summary of Introgression Detection Results on Mouse Chromosome 7

Metric	Result
Total introgressed sites	~9% of all sites [23]
Total physical length	~13 Megabase pairs (Mbp) [23]
Number of genes covered	>300 genes [23]
Key identified locus	Vkorc1 gene [23]
Negative control result	No introgression detected [23]

The analysis confirmed the Vkorc1 locus as a clear signal of adaptive introgression, consistent with its known role in rodenticide resistance [23]. The gene is located on chromosome 7 at position 127,492,235-127,494,789 (GRCm39 assembly) [44]. The model's specificity was validated by its performance on a negative control dataset, where it correctly detected no introgression [23].

Performance on Simulated Data

The accuracy and robustness of PhyloNet-HMM were quantified using synthetic data simulated under a coalescent model with recombination, isolation, and migration.

Table 2: Performance Metrics from Coalescent Simulations

Performance Metric	Outcome
Introgression detection accuracy	Accurate detection and inference [23]
ILS distinction	Effective teasing apart from introgression signals [23]
Parameter estimation	Accurate inference of population genetic parameters [23]

The Scientist's Toolkit

This section catalogues the essential research reagents and computational tools central to implementing the PhyloNet-HMM framework as described in this case study.

Table 3: Essential Research Reagents and Tools for PhyloNet-HMM Analysis

Item Name	Type	Function/Description
Whole-Genome Sequencing Data	Data	High-coverage genome sequences from the focal, donor, and outgroup species. The primary input for alignment [23].
Multiple Genome Alignment	Data/Resource	A column-wise alignment of all input genomes, representing homologous regions across species [23].
PhyloNet-HMM Software	Software Tool	The core implementation of the HMM-based comparative genomic framework for introgression detection [23].
Phylogenetic Network	Model	A mathematical graph representing the hypothesized evolutionary relationships, including hybridization events [23].
Population Genetic Parameters	Parameters	Estimates of ancestral effective population sizes, speciation times, and recombination rates that parameterize the HMM [40].
Vkorc1 Genomic Coordinates	Reference	Specific location of the gene of interest for validation (Chr7:127492235-127494789 in GRCm39) [44].

Discussion

The successful identification of the introgressed Vkorc1 locus underscores the power of PhyloNet-HMM to uncover biologically relevant adaptive introgression events that may be missed by methods confounded by ILS. The estimate that 9% of chromosome 7 is of introgressive origin highlights the potential scale of this phenomenon and its impact on genome evolution [23].

This framework is part of a broader trend in evolutionary genomics that moves beyond summary statistics to fully model-based probabilistic approaches, which explicitly incorporate the coalescent process and recombination to provide fine-scale insights [3]. Future directions include extending the models to handle more complex demographic scenarios and integrating machine learning approaches, such as treating introgression detection as a semantic segmentation task [3].

For the field of drug development, understanding the evolutionary history of loci like Vkorc1 is more than an academic exercise. It provides crucial context for pharmacogenetics; the same gene that underlies warfarin resistance in rodents also explains a significant portion of warfarin dosing variability in humans [41] [43]. Recognizing that key drug targets or metabolizing enzymes may have mosaic evolutionary histories influenced by introgression can inform studies of population-specific drug responses and the prediction of resistance mechanisms.

Navigating Challenges and Enhancing HMM Performance

Hidden Markov Models (HMMs) provide a powerful statistical framework for analyzing sequential data where the underlying system states are not directly observable. In comparative genomic introgression research, HMMs have proven invaluable for detecting regions of the genome that have been transferred between species through hybridization, a process with significant implications for evolution, adaptation, and drug target discovery [4] [16]. The effectiveness of HMMs in these applications rests fundamentally on two core assumptions: the Markov property, which governs the dynamics of the hidden states, and observation independence, which defines the relationship between states and their emitted symbols [45] [12]. This article explores these foundational assumptions within the context of genomic introgression research, providing experimental protocols and analytical frameworks that enable researchers to properly leverage HMMs while understanding their theoretical limitations.

Theoretical Foundation of HMM Assumptions

The Markov Property in Hidden States

The Markov property, also known as the homogeneous Markov assumption, stipulates that the probability of transitioning to a future state depends exclusively on the current state, independent of all previous states [13] [45]. Formally, for a hidden state sequence (X = (X1, X2, ..., XT)) and observation sequence (O = (O1, O2, ..., OT)), this property is defined as:

[P(X{t+1} = j | Xt = i, X{t-1} = i{t-1}, ..., X1 = i1) = P(X{t+1} = j | Xt = i) = a_{ij}]

where (a_{ij}) represents the transition probability from state (i) to state (j) [45] [12]. This first-order Markov assumption implies that the hidden state process has limited memory, with each state encapsulating all relevant information from history for predicting future states. In genomic introgression research, this translates to modeling how genomic regions transition between different evolutionary histories (e.g., species tree vs. introgression tree) along a chromosome, where each position's evolutionary pathway depends primarily on the immediately preceding region due to recombination breakpoints [4].

Observation Independence Assumption

The observation independence assumption posits that each observation (Ot) depends solely on the current hidden state (Xt), independent of all other states and observations [45]. Mathematically, this is expressed as:

[P(Ot = vk | Xt = j, X{t-1}, O{t-1}, ..., X1, O1) = P(Ot = vk | Xt = j) = b_j(k)]

where (bj(k)) represents the emission probability of symbol (vk) from state (j) [45] [12]. In biological sequence analysis, this means that the observed nucleotide or amino acid at a specific genomic position is determined exclusively by the underlying functional or evolutionary state at that position, without influence from adjacent positions. While this assumption facilitates computational efficiency, it represents a simplification of biological reality where dependencies often exist between adjacent sites in molecular sequences [12].

Table 1: Core Parameters Defining a Hidden Markov Model

Parameter	Symbol	Description	Biological Interpretation in Introgression Research
State Space	(Q = {q1, ..., qN})	Set of all possible hidden states	Possible evolutionary histories (e.g., species tree, introgression tree)
Observation Space	(V = {v1, ..., vM})	Set of all possible observable symbols	Molecular sequence alphabet (e.g., {A, C, G, T} for DNA)
Transition Probability	(a{ij} = P(X{t+1} = j | X_t = i))	Probability of moving from state i to state j	Likelihood of change in evolutionary history between adjacent genomic regions
Emission Probability	(bj(k) = P(Ot = vk | Xt = j))	Probability of emitting symbol k from state j	Likelihood of observing specific genetic variants under a given evolutionary history
Initial State Distribution	(\pii = P(X1 = i))	Probability of starting in state i	Prior expectation for evolutionary history at chromosome start

HMM Architecture and Computational Framework

The following diagram illustrates the general architecture of an HMM, highlighting the relationship between hidden states and observations, governed by the Markov property and observation independence assumption:

Diagram 1: HMM Architecture. The diagram visualizes the dual stochastic process of an HMM. Hidden states (yellow circles) follow a Markov process where each state depends only on its immediate predecessor (red arrows). Observations (blue squares) are generated by emission probabilities (green arrows) dependent only on the current hidden state.

Algorithmic Implementation for Genomic Introgression

The forward algorithm enables efficient computation of observation sequence probabilities given model parameters through dynamic programming. For genomic introgression detection, this allows evaluation of how well observed genetic variation data fits different evolutionary scenarios [45]. The algorithm proceeds through three stages:

• Initialization: [ \alpha1(i) = \pii \cdot bi(O1), \quad i = 1,2,...,N ]

• Recursion: [ \alpha{t+1}(j) = \left(\sum{i=1}^N \alphat(i) \cdot a{ij}\right) \cdot bj(O{t+1}), \quad t = 1,2,...,T-1 ]

• Termination: [ P(O|\lambda) = \sum{i=1}^N \alphaT(i) ]

where (\alphat(i) = P(O1, O2, ..., Ot, X_t = i | \lambda)) represents the probability of the partial observation sequence up to time (t) and being in state (i) at time (t) [45].

For decoding the most probable hidden state sequence (e.g., identifying introgressed vs. non-introgressed genomic regions), the Viterbi algorithm employs a similar dynamic programming approach but maximizes rather than sums probabilities:

• Initialization: [ \delta1(i) = \pii \cdot bi(O1), \quad \psi_1(i) = 0 ]

• Recursion: [ \deltat(j) = \max{1 \leq i \leq N} [\delta{t-1}(i) \cdot a{ij}] \cdot bj(Ot), \quad \psit(j) = \arg\max{1 \leq i \leq N} [\delta{t-1}(i) \cdot a{ij}] ]

• Termination and Backtracking: [ P^* = \max{1 \leq i \leq N} \deltaT(i), \quad XT^* = \arg\max{1 \leq i \leq N} \deltaT(i), \quad Xt^* = \psi{t+1}(X{t+1}^*) ]

where (\deltat(i)) represents the highest probability along a single path ending in state (i) at time (t), and (\psit(i)) tracks the optimal previous state [45].

Experimental Protocol for Introgression Detection

PhyloNet-HMM Framework for Eukaryotic Genomes

The PhyloNet-HMM framework represents an advanced implementation specifically designed for detecting introgression in eukaryotic genomes while accounting for incomplete lineage sorting (ILS) [4] [16]. The methodology combines phylogenetic networks with HMMs to simultaneously capture reticulate evolutionary relationships between species and dependencies within genomes.

Table 2: PhyloNet-HMM Experimental Workflow

Step	Protocol Description	Key Parameters	Output
1. Data Preparation	Obtain multiple sequence alignments from at least three closely related species with suspected historical hybridization	Alignment length (>10,000 bp recommended), sequence quality filters, reference genome coordinates	Processed multiple sequence alignment in FASTA or VCF format
2. Model Specification	Define phylogenetic network topology incorporating putative introgression events	Reticulation nodes, inheritance probabilities, population size parameters	Species network model with annotated introgression branches
3. Parameter Estimation	Apply expectation-maximization (Baum-Welch) to estimate transition and emission probabilities	Convergence threshold (ε=1e-6), maximum iterations (100+), regularization parameters	Trained HMM with optimized transition and emission matrices
4. Decoding	Implement Viterbi algorithm to identify most likely evolutionary history for each genomic region	Window size (1kb-10kb), state posterior probability threshold (>0.95)	Annotation of genomic regions by evolutionary history (species tree vs. introgression)
5. Validation	Perform statistical testing using parametric bootstrap or cross-validation	False discovery rate control (e.g., FDR < 0.05), significance thresholds	Statistically supported introgressed regions with confidence estimates

Computational Implementation Protocol

Materials and Software Requirements:

• Whole-genome sequencing data from multiple individuals per species
• High-performance computing cluster with ≥16GB RAM
• PhyloNet software package (v3.8 or higher)
• Python with hmmlearn, NumPy, and SciPy libraries
• Multiple sequence alignment software (MAFFT, MUSCLE)

Step-by-Step Procedure:

• Data Preprocessing:

  • Quality filter raw sequencing reads using FastQC and Trimmomatic
  • Map reads to reference genome using BWA-MEM or Bowtie2
  • Call variants using GATK HaplotypeCaller following best practices
  • Generate multiple sequence alignment for target genomic regions
  • Convert alignment to HMM-readable format (e.g., binary matrix)

• Model Training:

• Genome Scanning:

  • Divide genome into non-overlapping windows (1-10kb depending on recombination rate)
  • For each window, compute posterior probabilities for each state using forward-backward algorithm
  • Apply Viterbi decoding to identify optimal state path
  • Merge adjacent windows with same state assignment

• Statistical Significance Testing:

  • Generate null distribution of test statistics using simulations without introgression
  • Compute likelihood ratio test statistics for each candidate introgressed region
  • Apply multiple testing correction (Benjamini-Hochberg FDR control)
  • Validate against known introgressed loci (e.g., Vkorc1 in mice)

Case Study: Adaptive Introgression in Mouse Genomes

Application of the PhyloNet-HMM framework to chromosome 7 variation data from Mus musculus domesticus successfully detected the previously reported adaptive introgression event involving the rodent poison resistance gene Vkorc1 [4] [16]. The analysis revealed that approximately 9-12% of sites within chromosome 7 were of introgressive origin, covering about 13-18 Mbp and over 300 genes [4] [16]. This finding demonstrates the power of HMM-based approaches to identify both known and novel introgressed regions at genomic scale.

The following diagram illustrates the evolutionary scenario modeled in this case study, where hybridization occurs between species B and C after their initial divergence:

Diagram 2: Reticulate Evolution with Introgression. The model depicts a hybridization event between species B and C followed by back-crossing and selection, resulting in genomic introgression. HMMs help distinguish regions with different evolutionary histories (species tree vs. introgression) along the genome.

Table 3: Quantitative Results from Mouse Chromosome 7 Introgression Analysis

Genomic Feature	Value	Biological Significance
Total introgressed sites	9-12% of chromosome 7	Extensive historical gene flow between mouse species
Physical coverage	13-18 Mbp	Large genomic regions affected by introgression
Genes in introgressed regions	>300 genes	Potential for adaptive trait transfer
Vkorc1 region	Statistically significant detection	Validation of known adaptive introgression
Negative control data sets	No introgression detected	Specificity of the method confirmed

Research Reagent Solutions

Table 4: Essential Research Materials and Computational Tools for HMM-based Introgression Studies

Resource Type	Specific Tool/Resource	Application in Introgression Research
Software Packages	PhyloNet [4]	Phylogenetic network inference and HMM implementation
	hmmlearn (Python) [18]	General-purpose HMM implementation and training
	BEAGLE [4]	High-performance likelihood computation
	BWA/GATK [4]	Sequence alignment and variant calling
Biological Data Resources	1000 Eukaryotic Genomes Project	Comparative genomic data for multiple species
	VCF file format	Standardized variant call data exchange
	UCSC Genome Browser [4]	Visualization and annotation of introgressed regions
Computational Resources	High-performance computing cluster	Handling genome-scale data and computations
	Minimum 16GB RAM	Processing whole-genome sequence data
	Python/R statistical environment	Data analysis, visualization, and custom algorithms

The Markov property and observation independence assumptions provide the mathematical foundation that enables efficient inference in hidden Markov models for comparative genomic introgression research. While these assumptions represent simplifications of biological reality, frameworks like PhyloNet-HMM demonstrate their remarkable utility in detecting introgression signals while accounting for confounding factors like incomplete lineage sorting [4] [16]. The experimental protocols and analytical frameworks presented here offer researchers a structured approach to implementing HMM-based analyses while maintaining awareness of theoretical assumptions and their biological implications. As genomic datasets continue to expand across diverse taxa, HMMs will remain indispensable tools for deciphering the complex evolutionary histories written in genomes, ultimately advancing our understanding of adaptation, speciation, and the genetic basis of traits relevant to human health and disease.

Selecting the Optimal Number of Hidden States and Model Complexity

In the field of comparative genomics, accurately identifying introgressed genomic regions—where genetic material has been transferred between species through hybridization—is essential for understanding evolutionary processes and adaptive traits. Hidden Markov Models (HMMs) have emerged as a powerful statistical framework for detecting these introgression events by modeling the underlying genealogical transitions along genomic sequences. The PhyloNet-HMM framework represents a significant advancement in this domain, combining phylogenetic networks with HMMs to distinguish true introgression signatures from confounding signals like Incomplete Lineage Sorting (ILS) while accounting for dependencies across genomic loci [4]. A critical challenge in applying these models lies in selecting the optimal number of hidden states, which directly impacts model accuracy, biological interpretability, and the reliability of downstream analyses. This article provides detailed application notes and protocols for determining this crucial parameter, with specific emphasis on genomic introgression research.

The State Selection Problem in Genomic HMMs

Theoretical Foundations and Challenges

In an HMM, hidden states represent unobserved conditions that generate the observed data through state-specific probability distributions. For genomic introgression analysis, these states typically correspond to different evolutionary histories or genealogical relationships [4] [13]. The number of states (K) fundamentally determines model complexity; too few states can oversimplify the biological reality, while too many can lead to overfitting, where the model captures noise rather than signal [46].

The core challenge arises because likelihood-based metrics alone cannot reliably determine the optimal K. As model complexity increases, the likelihood of the observed data typically improves, creating a misleading impression that more states always yield better models [46]. This overfitting phenomenon is particularly problematic in genomic applications where the goal is biological interpretation rather than mere data fitting. An overfit HMM may identify spurious introgressed regions or fail to detect genuine evolutionary events, ultimately compromising research conclusions.

Biological Considerations for State Number Selection

In comparative genomics, the number of hidden states should reflect plausible biological scenarios. For example, when analyzing three species with potential introgression between two of them, the PhyloNet-HMM framework might utilize states representing different parental species trees [4]. The biological context provides critical constraints—the number of states should correspond to meaningful evolutionary categories rather than serving as mere statistical abstractions.

Quantitative Framework and Methodologies

Model Selection Criteria

The following table summarizes the primary metrics and approaches for determining the optimal number of hidden states in HMMs:

Table 1: Methods for Selecting the Optimal Number of Hidden States in HMMs

Method	Principle	Advantages	Limitations
Cross-Validation	Splits data into training/test sets; evaluates generalization performance [46]	Directly measures predictive accuracy; reduces overfitting	Computationally intensive; results may vary with different splits
Information Criteria (AIC/BIC)	Adds penalty terms for model complexity to the log-likelihood [46]	Balanced trade-off between fit and complexity; computationally efficient	Assumes large sample sizes; may not align with biological interpretability
Ordered HMM (oHMMed)	Restricts transitions to neighboring states based on emission parameter ordering [47]	Reduces estimable parameters; enables meaningful state labeling; intuitive diagnostics	Requires autocorrelated data; specific to ordered emission distributions
Bayesian Methods	Uses MCMC sampling to obtain posterior distributions of model parameters [47]	Fully probabilistic uncertainty quantification; avoids overfitting through priors	Computationally demanding; complex implementation and diagnostics

The oHMMed Framework for Genomic Applications

The oHMMed (ordered HMM with emission densities) approach offers particular advantages for genomic data where states naturally exhibit ordering [47]. By parameterizing emission densities as convex functions and restricting transitions to neighboring states within an imposed order, oHMMed implements a tridiagonal transition matrix that reflects the autocorrelation patterns common in genomic landscapes. This structured approach significantly reduces the number of estimable parameters compared to conventional HMMs, thereby mitigating overfitting while providing intuitive diagnostic criteria for state number selection.

In practice, oHMMed has been successfully applied to genomic features such as GC content and chromatin accessibility data, where it identifies statistically distinguishable regions with different mean emissions [47]. The ordering constraint not only prevents the "label switching" problem common in MCMC methods but also enables meaningful comparison of states across different runs or algorithms.

Experimental Protocols and Implementation

Protocol 1: Standardized Pipeline for State Number Selection

Objective: Determine the optimal number of hidden states (K) for a PhyloNet-HMM analysis of genomic introgression.

Input Requirements:

• Multiple sequence alignments from the species of interest
• A set of candidate parental species trees
• Computational resources for HMM training and evaluation

Procedure:

• Data Preparation

  • Format aligned genomic sequences as required by PhyloNet-HMM [4]
  • Define the set of potential parental species trees based on evolutionary hypotheses

• Model Training

  • For each candidate K (typically K = 2 to 5 for initial exploration):
    • Initialize PhyloNet-HMM parameters with K hidden states
    • Implement multiple random restarts (e.g., 10) to avoid local optima [48]
    • Execute the Baum-Welch algorithm until convergence (e.g., EM precision threshold of 1e-6)
    • Record the maximum log-likelihood achieved across restarts

• Model Evaluation

  • Calculate AIC/BIC values for each K:
    • AIC = -2logLik + 2P
    • BIC = -2logLik + Plog(N) where P is the number of parameters and N is the sequence length
  • Alternatively, implement cross-validation:
    • Partition data into training (e.g., 66%) and test sets [49]
    • Train on the training set and evaluate log-likelihood on the test set
    • Repeat with multiple random partitions to assess stability

• Biological Validation

  • For the optimal K, examine the identified introgressed regions for known adaptive genes (e.g., Vkorc1 in mice) [4]
  • Compare with negative control datasets where no introgression is expected [4]

Output Interpretation:

• The optimal K typically corresponds to the model with the lowest AIC/BIC or highest cross-validated likelihood
• Ensure the resulting states have biological meaning within the comparative genomic context

Protocol 2: oHMMed-Specific Implementation for Genomic Landscapes

Objective: Apply oHMMed to segment genomes into compositionally distinct regions with ordered states.

Input Requirements:

• Genomic feature data (e.g., GC content, gene density, epigenetic markers) in non-overlapping windows
• Pre-specification of emission distribution type (Normal mixture for continuous data, Poisson-Gamma for count data)

Procedure:

• Model Specification

  • Select the number of states K based on prior biological knowledge or preliminary analyses
  • Choose appropriate emission distributions for the genomic features:
    • Normal distributions for GC proportions (shared standard deviation, state-specific means) [47]
    • Gamma distributions for gene counts (shared shape parameter, state-specific rates) [47]

• Parameter Estimation

  • Implement MCMC sampling (e.g., Gibbs sampler) to obtain posterior distributions
  • Utilize the ordering constraint to ensure state-specific parameters are identifiable
  • Run multiple chains to assess convergence

• State Number Diagnostics

  • Apply oHMMed's diagnostic criteria for selecting K [47]
  • Compare models with different K using posterior predictive checks
  • Evaluate the biological coherence of the resulting genomic segmentation

• Annotation and Interpretation

  • Use the Viterbi algorithm to determine the most likely state sequence [48]
  • Characterize the genomic features associated with each state
  • Examine transition patterns between states along chromosomes

Visualization and Computational Framework

The following diagram illustrates the core architecture of an HMM for genomic introgression analysis, highlighting the relationship between hidden states (evolutionary histories) and observed genomic data:

Diagram 1: HMM structure for genomic analysis. Hidden states (X) represent evolutionary histories, while observations (Y) represent genomic sequence data. The Markov property dictates that each state depends only on its immediate predecessor.

The Scientist's Toolkit: Research Reagent Solutions

Table 2: Essential Computational Tools for HMM-Based Introgression Analysis

Tool/Resource	Function	Application Context	Implementation Considerations
PhyloNet	Infers phylogenetic networks and detects introgression [4]	Comparative genomic analysis of eukaryotic species	Open-source package; implements PhyloNet-HMM for genome-wide scans
oHMMed R Package	Ordered HMM with emission densities for genomic segmentation [47]	Characterizing autocorrelated genomic landscapes (GC content, epigenetic markers)	Available on CRAN; models continuous emission distributions
HMMTool	Graphical interface for HMM analysis of sequence data [48]	Rainfall modeling, ecological applications, and general time series analysis	GUI implementation; supports multiple spatial models and distribution types
Baum-Welch Algorithm	Expectation-Maximization approach for HMM parameter estimation [4] [48]	Training HMMs on observed sequence data	Multiple random restarts recommended to avoid local optima
Viterbi Algorithm	Finds the most probable sequence of hidden states [48]	Annotating genomic regions with evolutionary histories	Efficient dynamic programming solution; provides discrete state assignments
Markov Chain Monte Carlo (MCMC)	Bayesian parameter estimation for HMMs [47]	oHMMed implementation and other Bayesian HMM variants	Provides full posterior distributions; requires convergence diagnostics

Selecting the optimal number of hidden states represents a critical step in applying HMMs to comparative genomic introgression research. Neither purely statistical metrics nor biological intuition alone suffices; rather, an integrated approach that combines cross-validation, information criteria, and biological validation provides the most robust framework. The emergence of specialized methods like oHMMed, which incorporates domain knowledge through state ordering and emission parameter constraints, offers promising avenues for enhancing biological interpretability while maintaining statistical rigor. As genomic datasets continue to grow in size and complexity, these methodological considerations will become increasingly vital for extracting meaningful evolutionary insights from the genomic mosaics shaped by introgression events throughout the tree—or network—of life.

Mitigating Overfitting and Underfitting in Model Training

In the specialized field of comparative genomic introgression research, the reliability of analytical models is paramount. Hidden Markov Models (HMMs) have emerged as a powerful probabilistic framework for detecting introgressive regions by distinguishing true hybridization signals from spurious patterns caused by evolutionary processes like incomplete lineage sorting (ILS) [4] [16]. However, like all machine learning models, HMMs are susceptible to the dual pitfalls of overfitting and underfitting, which can severely compromise the interpretation of genomic landscapes [50]. Overfitting occurs when a model learns the training data too closely, including its noise and idiosyncrasies, thereby failing to generalize to new, unseen data [51] [50]. Conversely, underfitting arises when a model is too simplistic to capture the underlying patterns in the data, leading to poor performance even on training data [51] [52]. For researchers investigating questions of adaptive introgression—such as the transfer of the rodenticide resistance gene Vkorc1 in mouse populations [4] [16]—a model's generalization capability directly impacts the biological validity of its findings. This application note provides detailed protocols and frameworks to diagnose, prevent, and mitigate these critical issues within the context of HMM-based genomic analysis.

The following table summarizes the core characteristics and quantitative indicators of overfitting and underfitting, which are crucial for diagnosing model performance in genomic analyses.

Table 1: Diagnostic Characteristics of Overfitting and Underfitting

Aspect	Overfitting	Underfitting
Core Problem	Model is overly complex, learning noise and spurious patterns from the training data [51] [50].	Model is overly simplistic, failing to learn the underlying data patterns [51] [52].
Performance on Training Data	Performance is typically very high (e.g., low loss, high accuracy) [50] [53].	Performance is poor (e.g., high loss, low accuracy) [50] [53].
Performance on Validation/Test Data	Performance is significantly worse than on training data; model fails to generalize [51] [50].	Performance is poor and similar to that on training data; model cannot generalize effectively [50].
Key Visual Cue in Learning Curves	A growing gap between training and validation loss curves as training progresses [54].	Training and validation loss curves converge at a high value, showing minimal decrease [54].
Primary Error Type	High variance [51].	High bias [51].

Experimental Protocols for Mitigation

This section outlines detailed, actionable protocols for mitigating overfitting and underfitting, with specific adaptations for HMM-based genomic research.

Protocol 1: Mitigating Overfitting in Model Architecture and Training

Application: This protocol is essential when training HMMs or other models (e.g., neural networks) on genomic sequence alignments to ensure robust detection of introgressed loci.

Materials:

• Aligned genomic sequences (e.g., VCF files, multiple sequence alignments).
• Computational framework for HMM training (e.g., PhyloNet-HMM [4] [16] or custom implementations).
• Training and validation data partitions.

Procedure:

• Apply Regularization:
  • Introduce a complexity penalty to the model's objective function. In the context of HMMs, this can be achieved by imposing priors on the transition and emission probabilities during the Baum-Welch algorithm to prevent them from becoming over-specialized to the training data [51] [53].
  • For other model types, such as linear models used in preliminary feature selection, employ L1 (Lasso) or L2 (Ridge) regularization. The loss functions are defined as:
    • L1 Regularization: Loss = MSE + λ ∑|β| [53]
    • L2 Regularization: Loss = MSE + λ ∑(β)² [53]

• Implement Early Stopping:

  • Continuously monitor the model's loss on a held-out validation set throughout the training process.
  • Terminate training when the validation loss fails to improve for a pre-defined number of epochs (the "patience" parameter) [54] [53].
  • Benchmark: The OverfitGuard approach, which uses a history-based classifier, has been shown to stop training at least 32% earlier than standard early stopping while achieving the same or better model selection rate [54].

• Utilize History-Based Detection:

  • Employ a trained time-series classifier (e.g., OverfitGuard) on the validation loss history to automatically flag an overfit model post-training. This method has achieved an F1 score of 0.91 in detection tasks, outperforming other non-intrusive methods [54].

• Conduct Feature Selection:

  • Prior to model training, perform rigorous exploratory data analysis (EDA) to remove features (e.g., specific genomic loci) that are unrelated to the evolutionary signal of interest, as they can contribute to learning noise [53].
  • Treat outliers in the training data that may cause the model to compute aberrantly high coefficients [53].

Protocol 2: Mitigating Underfitting in Model Design

Application: Use this protocol when a model demonstrates consistently poor performance on both training and validation data, indicating a failure to capture the core evolutionary signals.

Materials:

• The same aligned genomic sequences and computational frameworks as in Protocol 1.

Procedure:

• Increase Model Complexity:
  • For HMMs, consider increasing the number of hidden states to better represent the complexity of the evolutionary process, such as multiple introgression scenarios [18] [55].
  • If using tree-based models for auxiliary analysis, switch from a simple Decision Tree to a more complex ensemble method like Gradient Boosting [53].

• Enhance Feature Engineering:

  • Ensure the input data contains the necessary informative features. In genomics, this may involve incorporating additional phylogenetic invariants or summary statistics that are known to be sensitive to introgression [3] [53].
  • Create new, more informative features from existing data (e.g., polynomial features) to help the model capture complex, non-linear relationships [51].

• Reduce Regularization:

  • If a model is underfitting, review and potentially decrease the strength of the regularization parameters (e.g., the λ value in L1/L2 regularization) that may be overly constraining the model's capacity to learn [52] [53].

• Select the Appropriate Algorithm:

  • Confirm that the chosen model is inherently capable of capturing the patterns under investigation. For instance, a simple linear model may be fundamentally unsuitable for detecting the complex, reticulate signals of introgression that HMMs are designed to identify [4] [53].

Workflow Visualization

The following diagram illustrates the integrated decision-making process for diagnosing and addressing model fit issues within a genomic analysis pipeline.

Diagram 1: Model Fit Diagnosis Workflow

The Scientist's Toolkit: Research Reagent Solutions

Table 2: Essential Computational Tools for Introgression Analysis with HMMs

Research Reagent	Function/Description	Application in Mitigation
PhyloNet-HMM	A computational framework that integrates phylogenetic networks with HMMs to detect introgression while accounting for ILS and dependencies across genomic loci [4] [16].	Core model for analysis; its parameters (state number, probabilities) are tuned to prevent under/overfitting.
Training History Data	The byproduct of model training, recording metrics like loss over each epoch [54].	Used for visual diagnosis from learning curves and as input for automated overfitting detection tools like OverfitGuard [54].
Regularization Parameters (λ)	Scalar values that control the strength of the penalty applied to model coefficients to discourage complexity [51] [53].	Increasing λ helps mitigate overfitting; decreasing λ can help mitigate underfitting.
Validation Dataset	A subset of the genomic data not used during model training, reserved for evaluating generalization performance [50].	Critical for implementing early stopping, evaluating model fit, and tuning hyperparameters.
Baum-Welch Algorithm	An Expectation-Maximization (EM) algorithm used for training HMMs by iteratively updating transition and emission probabilities [18] [55].	The core training procedure where regularization priors can be incorporated to prevent overfitting.

Effectively managing model complexity by mitigating overfitting and underfitting is not merely a technical exercise but a foundational requirement for producing biologically meaningful results in comparative genomic introgression research. By implementing the detailed protocols and diagnostic workflows outlined in this document—ranging from architectural strategies like regularization and early stopping to systematic performance evaluation—researchers can significantly enhance the reliability and interpretability of their HMM-based analyses. This rigorous approach ensures that detected signatures of introgression, such as those linked to adaptive traits, reflect true evolutionary history rather than artifacts of statistical over-optimization or under-performance.

Handling Data Sparsity, Noise, and Missing Values in Genomic Data

Genomic data analysis is fundamentally challenged by technical artifacts including data sparsity, technical noise, and missing values. These issues are particularly critical in comparative genomic introgression research, where accurately identifying ancestral genomic regions requires robust statistical models that can distinguish true biological signal from artifact. Hidden Markov Models (HMMs) provide a powerful framework for modeling sequential genomic data, but their performance depends heavily on how these data quality challenges are addressed. This application note provides detailed protocols for handling these issues within HMM-based genomic analyses, enabling more reliable detection of introgressed regions and improving the accuracy of evolutionary inferences.

Understanding the Data Challenges

Characterization of Data Imperfections

Genomic datasets exhibit several types of data quality issues that can significantly impact downstream analysis:

• Data Sparsity: In single-cell Hi-C data, which measures chromatin interactions, a high percentage of potential contacts remain undetected, creating sparse matrices that obscure true biological structures [56]. Similarly, in genomic windows analysis, many regions may lack sufficient coverage for reliable variant calling.

• Technical Noise: This includes non-biological fluctuations introduced during experimental processes. In single-cell RNA-seq data, technical noise manifests as "dropout" events where transcripts fail to be detected despite being present in the cell [56]. Background noise in single-cell experiments can constitute 3-35% of total counts per cell, directly impacting marker gene detectability [57].

• Missing Data Mechanisms: Missing values in genomic datasets fall into three primary categories [58]:

  • Missing Completely at Random (MCAR): Missingness independent of both observed and unobserved data
  • Missing at Random (MAR): Missingness depends on observed data but not unobserved values
  • Missing Not at Random (MNAR): Missingness depends on unobserved factors

Impact on Genomic Analysis

These data imperfections directly affect HMM performance in genomic introgression studies by obscuring the true state transitions and emission probabilities. Noisy or sparse data can lead to incorrect inference of hidden states (e.g., ancestral vs. introgressed regions), reduce statistical power to detect true introgressed segments, and increase false positive rates in identifying introgressed regions.

Computational Frameworks for HMM-Based Analysis

Advanced HMM Architectures for Genomic Data

Several specialized HMM frameworks have been developed to address genomic data challenges:

• Sparsely Correlated HMM (scHMM): This approach models correlations between multiple genomic datasets without the computational intractability of full multivariate HMMs. scHMM uses penalized regression to select a subset of correlated data series and estimate their effects on transition probabilities, effectively capturing biological interactions while maintaining computational feasibility [59].

• Ordered HMM with Emission Densities (oHMMed): This class of HMMs incorporates continuous emission distributions and orders hidden states by their emission means, restricting transitions to neighboring states. This architecture naturally models the autocorrelated nature of genomic landscapes and provides robust segmentation of genomes into homogeneous regions [47].

• Swarm Intelligence-Based Pliable HMM: This framework combines HMM with swarm optimization techniques (Particle Swarm Optimization, Whale Optimization Algorithm) to determine optimal hidden states, demonstrating particular utility for complex biomedical signal classification [60].

Table 1: Specialized HMM Frameworks for Genomic Data Challenges

HMM Framework	Primary Data Challenge Addressed	Key Mechanism	Genomic Applications
scHMM [59]	Multiple correlated datasets	Penalized regression for sparse correlations	Joint analysis of histone modifications
oHMMed [47]	Autocorrelated continuous data	Ordered states with restricted transitions	GC-content segmentation, chromatin accessibility
SI-HMM [60]	Complex signal noise	Swarm intelligence for state optimization	EEG classification (potential for genomic signals)
TRANSIT HMM [61]	Sparse insertion data	Essentiality state modeling	Genomic essentiality in transposon mutagenesis

Denoising Methods for Genomic Data

Effective denoising is crucial for preparing genomic data for HMM analysis:

• SAMDUDE: A denoising method specifically designed for aligned genomic data (SAM files) that employs a Discrete Universal Denoiser algorithm incorporating both read sequences and quality scores. This approach has demonstrated improved variant identification in whole genome sequencing data, identifying nearly 2,000 additional true variants while eliminating over 1,500 erroneous calls [62].

• RECODE/iRECODE: A high-dimensional statistics-based platform for technical noise reduction in single-cell data that simultaneously addresses technical noise and batch effects while preserving data dimensions. The method maps gene expression to an essential space using noise variance-stabilizing normalization and singular value decomposition [56].

• Background Noise Removal: For single-cell data, methods like CellBender, DecontX, and SoupX specifically target background noise from ambient RNA or barcode swapping. Evaluation studies indicate CellBender provides the most precise estimates of background noise levels [57].

Experimental Protocols

Protocol 1: Handling Missing Data in Genomic Features

Purpose: To impute missing values in genomic feature matrices prior to HMM analysis for introgression detection.

Materials:

• Genomic variant dataset with missing annotations
• R software environment with NAsImpute package [58]
• Computational resources (8GB+ RAM for genome-scale data)

Procedure:

• Data Preparation: Load genomic variant matrix with features (e.g., CADD, DANN, GERP, phyloP scores) and identify missing value patterns [58].
• Missing Mechanism Diagnosis: Use the NAsImpute package to assess whether missingness follows MCAR, MAR, or MNAR patterns.
• Method Selection: Based on missing data mechanism:
  • For MCAR/MAR: Implement Random Forest (missForest) or k-Nearest Neighbors (DMwR) imputation
  • For MNAR: Apply Multivariate Imputation by Chained Equations (MICE) or Amelia
• Imputation Execution:

• Validation: Calculate imputation accuracy using root mean square error (RMSE) and mean absolute error (MAE) on held-out complete cases [58].

Troubleshooting:

• If imputation produces biologically implausible values, constrain imputed ranges based on known biological limits
• For large datasets, use subsetting to manage computational load
• Validate imputation stability through multiple runs with different random seeds

Protocol 2: Denoising Aligned Genomic Data for Variant Calling

Purpose: To reduce noise in aligned sequencing data to improve variant calling accuracy for HMM-based introgression analysis.

Materials:

• Aligned genomic data in SAM/BAM format
• SAMDUDE software (https://github.com/ihwang/SAMDUDE) [62]
• Reference genome sequence

Procedure:

• Data Preparation: Sort and index BAM files using samtools.
• SAMDUDE Execution:

• Parameter Optimization: Adjust context window size based on read length and error profiles.
• Quality Score Recalibration: SAMDUDE simultaneously updates base calls and quality scores based on alignment context.
• Variant Calling: Perform variant calling on both original and denoised BAM files using standard tools (e.g., GATK, bcftools).
• Validation: Compare variant call sets using precision-recall metrics with validated variant databases.

Notes: SAMDUDE specifically targets substitution errors predominant in Illumina sequencing technologies and demonstrates improved variant calling performance compared to other denoisers like Musket and RACER [62].

Protocol 3: Joint Analysis of Multiple Genomic Datasets with scHMM

Purpose: To perform simultaneous HMM inference across multiple correlated genomic datasets while maintaining computational tractability.

Materials:

• Multiple genomic data series (e.g., different histone modifications, transcription factor binding)
• scHMM package (http://sourceforge.net/p/schmm/) [59]
• Linux computing environment

Procedure:

• Data Formatting: Format each genomic data series as a separate wiggle file with consistent genomic coordinates.
• Initial State Estimation: Run independent HMMs on each series to obtain initial hidden state estimates.
• Correlation Assessment: Identify significantly correlated series using penalized regression:

• Transition Probability Adjustment: For each series, adjust transition probabilities based on hidden states of correlated series using logistic regression models [59].
• Forward-Backward Algorithm: Apply the forward-backward algorithm with adjusted transition probabilities to infer final hidden states.
• Result Integration: Generate genome annotation based on joint hidden state predictions across all series.

Validation: Compare scHMM results with independent HMM analyses and validated genomic annotations. scHMM has demonstrated improved recovery of characterized histone modifications compared to independent HMMs [59].

Visualization and Workflow Diagrams

scHMM Framework for Correlated Genomic Series

Comprehensive Genomic Data Preprocessing Pipeline

The Scientist's Toolkit

Table 2: Essential Research Reagents and Computational Tools

Tool/Resource	Primary Function	Application Context	Key Reference
scHMM Package	Sparse correlated HMM inference	Joint analysis of multiple genomic datasets	[59]
SAMDUDE	Denoising aligned genomic data	Improving variant calling accuracy	[62]
RECODE/iRECODE	Technical noise reduction	Single-cell RNA-seq, Hi-C, spatial transcriptomics	[56]
NAsImpute R Package	Evaluation of imputation methods	Genomic feature matrices with missing values	[58]
oHMMed R Package	Ordered HMM with continuous emissions	Genomic landscape segmentation	[47]
CellBender	Background noise removal	Single-cell RNA-seq data	[57]
TRANSIT HMM	Genomic essentiality profiling	Transposon insertion sequencing (Tn-Seq)	[61]

Performance Evaluation and Benchmarking

Quantitative Assessment of Methods

Table 3: Performance Comparison of Imputation Methods for Genomic Data

Imputation Method	MCAR Scenario	MAR Scenario	MNAR Scenario	Computational Efficiency	Recommended Use
Random Forest (missForest)	Best performance (Lowest RMSE)	Best performance (Lowest RMSE)	Good performance	Computationally intensive	Small to medium datasets
k-Nearest Neighbors	Good performance	Good performance	Moderate performance	Moderate	Large datasets
MICE	Moderate performance	Good performance	Good performance	Moderate	MNAR scenarios
Mean/Median Imputation	Poor performance	Poor performance	Poor performance	High	Baseline only
Amelia	Moderate performance	Moderate performance	Good performance	Moderate	MNAR scenarios

Note: Performance assessment based on evaluation of 31,245 genomic variants with 13 genome-wide features [58].

Effective handling of data sparsity, noise, and missing values is essential for robust HMM-based analysis of genomic introgression. The protocols and frameworks presented here provide comprehensive solutions for these data challenges, enabling more accurate detection of introgressed regions and advancing our understanding of evolutionary genomics. By implementing specialized HMM architectures like scHMM and oHMMed alongside rigorous data preprocessing protocols, researchers can significantly improve the reliability of their genomic inferences. As genomic technologies continue to evolve, further development of integrated analysis frameworks will remain crucial for extracting maximal biological insight from increasingly complex and noisy genomic data.

Hidden Markov Models (HMMs) have become indispensable statistical frameworks in bioinformatics for modeling sequential data with hidden, unobservable states. In comparative genomics, they are particularly valuable for detecting complex evolutionary patterns such as introgression—the integration of genetic material from one species into another through hybridization. The power of HMMs lies in their capacity to capture dependencies between adjacent genomic loci while distinguishing true introgression signals from spurious patterns that arise due to other evolutionary forces like incomplete lineage sorting (ILS). The core HMM framework consists of a hidden state sequence, representing unknown evolutionary histories, and an observed sequence, representing the genomic data itself. These models are parameterized by: (i) a state transition probability matrix (A), governing transitions between hidden states; (ii) an emission probability matrix (B), defining the probability of observations given hidden states; and (iii) an initial state distribution (π). The parameter set λ = (A, B, π) provides a complete description of the model for computational analysis [45].

The application of HMMs to genomic introgression represents a significant methodological advancement. Traditional phylogenetic methods assume a strictly tree-like evolutionary history, which fails to account for reticulate processes such as hybridization and introgression. The PhyloNet-HMM framework, for instance, addresses this limitation by combining phylogenetic networks with HMMs to simultaneously capture potentially reticulate evolutionary histories and dependencies within genomes. This approach can scan aligned genomes for signatures of introgression while accounting for confounding factors like ILS and recombination. In one notable application to chromosome 7 data from house mice (Mus musculus domesticus), this method detected a previously reported adaptive introgression event involving the rodent poison resistance gene Vkorc1, along with other newly identified introgressed regions. The analysis estimated that approximately 9% of sites on chromosome 7 were of introgressive origin, covering about 13 Mbp and over 300 genes [4].

Core Algorithmic Framework

The Three Canonical Problems of HMMs

The practical application of Hidden Markov Models revolves around solving three fundamental problems [45]:

• Evaluation Problem: Given an HMM λ = (A, B, π) and an observation sequence O = {o₁, o₂, ..., oₜ}, calculate the likelihood P(O | λ). This probability represents how well the model explains the observed data. The Forward-Backward Algorithm provides an efficient solution using dynamic programming to compute this likelihood without enumerating all possible hidden state sequences.

• Decoding Problem: Given λ and O, determine the most probable hidden state sequence X = {x₁, x₂, ..., xₜ} that produced the observations. This is typically solved using the Viterbi Algorithm, which employs dynamic programming to find the single best path (state sequence) that maximizes P(X | O, λ).

• Learning Problem: Given an observation sequence O, adjust the model parameters λ = (A, B, π) to maximize P(O | λ). The Baum-Welch Algorithm, a special case of the Expectation-Maximization (EM) algorithm, is the standard method for this problem when training sequences are unlabeled.

The Viterbi Algorithm for Path Decoding

The Viterbi algorithm is a dynamic programming algorithm that identifies the most likely sequence of hidden states given a sequence of observations and a fully parameterized HMM. The algorithm operates by recursively computing the probability of the most probable path ending in each possible state at each time step.

The algorithm utilizes four key steps [45]:

• Initialization: For each state i, initialize the path probability δ₁(i) = πᵢ · bᵢ(o₁) and the backtracking pointer ψ₁(i) = 0.
• Recursion: For each observation position t from 2 to T and for each state j, compute: δₜ(j) = max₁≤ᵢ≤N [δₜ₋₁(i) · aᵢⱼ] · bⱼ(oₜ) and ψₜ(j) = argmax₁≤ᵢ≤N [δₜ₋₁(i) · aᵢⱼ].
• Termination: Identify the probability of the most likely path P* = max₁≤ᵢ≤N δT(i) and the starting state for backtracking xₜ* = argmax₁≤ᵢ≤N δT(i).
• Backtracking: For t from T-1 down to 1, reconstruct the optimal path: xₜ* = ψₜ₊₁(xₜ₊₁*).

In genomic introgression analysis, the Viterbi algorithm is deployed to annotate genomes by identifying the most likely evolutionary history (e.g., species tree vs. introgression history) at each genomic position. For example, in the PhyloNet-HMM framework, the hidden states represent different local genealogical trees embedded within a phylogenetic network. The Viterbi path thus provides a detailed map of which genomic regions evolved under which parental species tree, directly pinpointing introgressed segments [4].

The Baum-Welch Algorithm for Parameter Learning

The Baum-Welch algorithm addresses the challenge of learning HMM parameters from unlabeled observation sequences. As an expectation-maximization algorithm, it iteratively refines model parameters to maximize the likelihood of the observed data.

The algorithm proceeds as follows [45] [63]:

• Initialization: Start with an initial estimate for model parameters λ = (A, B, π).
• Expectation Step (E-step): Using the current parameter estimates, compute the expected number of transitions between states and the expected number of emissions from each state. This involves calculating:
  • ξₜ(i, j) = P(xₜ = i, xₜ₊₁ = j | O, λ): The probability of being in state i at time t and state j at time t+1, given the observations.
  • γₜ(i) = P(xₜ = i | O, λ): The probability of being in state i at time t, given the observations.
• Maximization Step (M-step): Update the model parameters using the expectations computed in the E-step:
  • aᵢⱼ = (Expected number of transitions from state i to j) / (Expected total transitions from state i)
  • bⱼ(k) = (Expected number of times in state j emitting symbol k) / (Expected total times in state j)
  • πᵢ = γ₁(i): The expected frequency of being in state i at the first time step.

Advanced variants of this algorithm have been developed for specific genomic applications. The constrained Baum-Welch (cBW) algorithm extends the standard approach to handle situations with partially labeled training sequences, which is common when limited experimental validation data is available. This method constrains the possible hidden state paths during the E-step to only those consistent with the available partial labels, leading to more accurate parameter estimation when full supervision is not feasible [63] [64].

Table 1: Key Algorithms for HMM Optimization in Genomics

Algorithm	Core Function	Genomic Application Example	Key Advantage
Viterbi	Finds optimal hidden state sequence	Decoding evolutionary histories in PhyloNet-HMM [4]	Provides precise genomic localization
Standard Baum-Welch	Trains parameters from unlabeled data	General sequence modeling and annotation [45]	Works without supervised training data
Constrained Baum-Welch (cBW)	Trains parameters with partial labels	Detecting plasmodesmata signals with sparse experimental data [63]	Leverages limited experimental validation

Application Notes for Introgression Research

Experimental Design and Workflow

The successful application of HMMs to detect introgression requires careful experimental design and execution. The following workflow outlines the key steps from data preparation to biological interpretation:

• Data Collection and Alignment: Collect whole-genome sequencing data from the focal species and closely related reference species. Perform multiple sequence alignment to obtain positionally homologous sites across all genomes. For mouse introgression studies, for example, this would involve sequencing and aligning genomes from Mus musculus domesticus and its close relatives [4].

• Model Configuration: Define the phylogenetic network topology based on known or hypothesized evolutionary relationships, including potential introgression events. Determine the set of possible hidden states, where each state typically corresponds to a different genealogical history within the network.

• Parameter Training: Use the Baum-Welch algorithm to estimate transition and emission probabilities from the aligned genomic data. If partial biological knowledge is available (e.g., certain genomic regions confirmed through experimental studies to be introgressed), implement the constrained Baum-Welch algorithm to incorporate these constraints during training [63].

• Genome Scanning and Decoding: Apply the Viterbi algorithm to scan the aligned genomes and decode the most likely hidden state at each genomic position. This generates a detailed annotation of the genome into regions with different evolutionary histories.

• Statistical Validation: Assess the significance of identified introgressed regions through comparison with negative control datasets and simulation studies. In the PhyloNet-HMM analysis, no introgression was detected in negative control data, confirming method specificity [4].

• Biological Interpretation: Annotate genes and functional elements within predicted introgressed regions. Perform enrichment analyses to identify potential adaptive advantages conferred by introgressed alleles, such as the Vkorc1 poison resistance variant in mice [4].

Diagram 1: HMM Analysis Workflow for Introgression Detection. This workflow outlines the key steps from data preparation to biological interpretation in genomic introgression studies.

Protocol: Detecting Adaptive Introgression with PhyloNet-HMM

This protocol provides a detailed methodology for detecting adaptive introgression using the PhyloNet-HMM framework, based on the successful application to mouse genomic data [4].

Materials and Software Requirements:

• Whole-genome sequencing data from at least three closely related species/populations
• High-performance computing cluster with sufficient memory for genomic analyses
• Multiple sequence alignment software (e.g., MAFFT, MUSCLE)
• PhyloNet software package (includes PhyloNet-HMM implementation)
• Programming environment (Python/R) for downstream analysis

Procedure:

• Data Preprocessing:

  • Perform quality control on raw sequencing reads using FastQC.
  • Map reads to a reference genome using BWA-MEM or similar aligner.
  • Call variants using GATK or bcftools, generating a VCF file.
  • Extract sequences from focal chromosomal regions and create multiple sequence alignments.

• Phylogenetic Network Specification:

  • Define the set of parental species trees based on known phylogeny.
  • For the mouse chromosome 7 analysis, the network included trees representing: (i) the consensus species relationship, and (ii) an alternative history with introgression from a donor species [4].
  • Specify the model structure, with hidden states corresponding to genealogies evolving within different parental trees.

• Model Training:

  • Initialize transition probabilities with higher values for self-transitions (to model genomic block structure).
  • Initialize emission probabilities based on evolutionary substitution models.
  • Run the Baum-Welch algorithm to train model parameters on the aligned genomic sequences.
  • Check for convergence by monitoring log-likelihood changes between iterations.

• Introgression Detection:

  • Execute the Viterbi algorithm on the trained model to decode the most likely state path across the genome.
  • Identify contiguous genomic regions assigned to the introgression state.
  • Apply a minimum length threshold (e.g., 1 kb) to filter spurious short predictions.

• Validation and Analysis:

  • Compare predictions with known introgressed regions (if available) to calculate accuracy metrics.
  • Annotate genes within predicted introgressed regions using genome annotation files (GTF/GFF).
  • Perform functional enrichment analysis using tools like g:Profiler or DAVID.

Troubleshooting Tips:

• If the model fails to converge, try different parameter initializations or increase the number of Baum-Welch iterations.
• If predictions show excessive fragmentation, adjust the transition probability matrix to favor longer state segments.
• Validate findings with independent methods, such as D-statistics or phylogenetic incongruence tests.

Advanced Methodological Adaptations

Handling Partial Labels with Constrained Training

Biological sequence analysis often faces the challenge of limited experimental validation, where only partial information is available about the true state of certain sequence regions. To address this, the constrained Baum-Welch (cBW) algorithm was developed to incorporate partial label information during model training [63] [64].

The key innovation of cBW lies in modifying the E-step of the standard Baum-Welch algorithm. While the traditional approach sums over all possible hidden state paths, the constrained version restricts this summation to only those paths that are consistent with the available partial labels. This incorporates prior biological knowledge directly into the parameter estimation process, leading to more accurate models when fully labeled training data is unavailable. In practice, this method has shown significant improvements in training accuracy compared to approaches that only consider partial labels at individual positions without constraining the entire path [63].

Application areas for cBW in genomics include:

• Refining predictions of plasmodesmata targeting signals in plant proteins using limited experimental data [63]
• Improving gene structure annotation when only partial experimental evidence is available
• Enhancing detection of functional domains in proteins with partially characterized motifs

Modeling Autocorrelation with Ordered HMMs

Genomic data frequently exhibits substantial autocorrelation, where observations at nearby loci are more similar than those further apart. The oHMMed (ordered HMM with emission densities) framework addresses this property by incorporating two key features [47]:

• Continuous Emission Densities: Instead of discrete emission probabilities, oHMMed models observations as realizations of continuous probability distributions (e.g., normal or gamma-Poisson mixtures) with state-specific parameters.

• Ordered State Transitions: The hidden states are ordered by their emission means, and transitions are restricted to neighboring states in this ordering. This creates a tridiagonal transition matrix that naturally captures the autocorrelation pattern commonly observed in genomic landscapes.

This approach has been successfully applied to segment genomes based on GC content, gene density, and chromatin accessibility data. In human genomic studies, oHMMed identified compositionally distinct regions while modeling the gradual transitions between different genomic domains, providing a more nuanced understanding of genomic architecture than traditional segmentation methods [47].

Table 2: Advanced HMM Variants for Genomic Analysis

HMM Variant	Methodological Innovation	Genomic Application	Advantage Over Standard HMM
Constrained Baum-Welch (cBW)	Incorporates partial labels during training	Signal peptide detection with sparse experimental data [63]	More accurate training with limited biological validation
oHMMed	Ordered states with continuous emissions	GC content segmentation, chromatin landscape analysis [47]	Better models genomic autocorrelation patterns
PhyloNet-HMM	Combines phylogenetic networks with HMM	Genome-wide introgression scanning [4]	Teases apart introgression from incomplete lineage sorting

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Tools for HMM-based Genomic Analysis

Resource Type	Specific Tool/Resource	Function in HMM Analysis	Application Context
Software Packages	PhyloNet[cite

Assessing Accuracy and Comparative Analysis of HMM Frameworks

Within comparative genomic introgression research, establishing the statistical validity of findings is paramount, particularly when using sophisticated models like hidden Markov models (HMMs). The analysis of genomic data to detect introgressed regions—where genetic material has transferred between species—is susceptible to both false positives and false negatives. Forward simulation provides a critical framework for benchmarking the performance of these analytical methods by creating synthetic datasets where the true evolutionary history is known. This allows researchers to rigorously quantify a method's statistical power to detect true introgression and its ability to control the false discovery rate (FDR). This protocol details the application of forward simulations for validating HMM-based methods in introgression research, providing a standardized approach for researchers and drug development scientists to assess the reliability of their genomic analyses [65] [4].

The need for such validation is underscored by the complex challenges in comparative genomics. Evolutionary processes like incomplete lineage sorting (ILS) can generate phylogenetic patterns that mimic introgression, leading to spurious conclusions if not properly accounted for [4]. Furthermore, the analysis of high-throughput genomic data constitutes a massive multiple testing problem, necessitating robust FDR control to ensure that reported introgressed regions are reliable [66]. The framework described herein addresses these challenges by leveraging the known ground truth of simulated data to evaluate how well HMM-based approaches can distinguish introgression from other confounding evolutionary signals.

Core Concepts and Principles

The Role of Forward Simulation in Validation

In the context of method validation, forward simulations are not merely a test but a foundational component of a rigorous evaluation cycle. They allow for the generation of synthetic genomic sequences under a specified evolutionary model that includes parameters for mutation, recombination, population divergence, and introgression. By applying an HMM-based method to these simulated datasets and comparing its inferences to the known simulated truth, one can obtain unbiased estimates of its performance characteristics [65] [4].

This process is intrinsically linked to the concept of predictive capability, which extends beyond simple agreement with empirical data. It involves using a computational model to forecast system responses under conditions for which empirical data may not be available, a common scenario in evolutionary studies of past introgression events [65]. The validation framework thus assesses whether the HMM possesses the fidelity required for its intended application, whether that is the discovery of adaptive introgression or the comprehensive scanning of a genome for introgressed loci [65].

Key Performance Metrics: Power and FDR

Two metrics are central to the benchmarking process:

• Statistical Power is defined as the probability that the method correctly identifies a genuinely introgressed genomic region. Also known as the true positive rate or sensitivity, it is crucial for ensuring that biologically relevant introgression events are not missed.
• False Discovery Rate (FDR) is the expected proportion of falsely reported introgressed regions among all regions declared significant. Controlling the FDR, for instance at a 5% threshold, is essential for the credibility of reported findings [66].

The interplay between these metrics is often visualized using a power-FDR curve, which illustrates the trade-off between sensitivity and specificity as the method's detection threshold is varied. An optimal method maintains high power while keeping the FDR below a desired level.

Computational Protocols

Workflow for Simulation-Based Validation

The following workflow outlines the primary steps for conducting a validation study, from simulation to performance assessment. The process is designed to be iterative, allowing researchers to refine their HMM models and analysis parameters based on initial benchmarking results [65].

Diagram 1: The forward simulation validation workflow for benchmarking HMM performance.

Experimental Setup and Data Simulation

The first phase involves generating the synthetic datasets that will serve as the benchmark.

Protocol 1: Simulating Genomic Data with Introgression

• Define the Evolutionary Model: Specify the phylogenetic network or population history, including speciation times, effective population sizes, and introgression events (timing, direction, and proportion of genetic material transferred). For example, a simple three-species model (Species A, B, C) with a hybridization event between B and C after their initial divergence is a common starting point [4].
• Set Genomic Parameters: Determine the total genome length, recombination rate, and mutation rate. The genome should be sufficiently long (e.g., simulating multiple chromosomes or a whole genome) to obtain reliable estimates of power and FDR.
• Incorporate Evolutionary Complexities: To create a realistic and challenging benchmark, simulate data that includes confounding factors such as:
  • Incomplete Lineage Sorting (ILS): Model ancestral population polymorphism by setting appropriate effective population sizes and divergence times [4].
  • Variable Recombination Rates: Use a recombination map instead of a uniform rate.
  • Selection: Introduce sites under positive selection within introgressed regions to model adaptive introgression.
• Execute Simulations: Use a forward simulation tool such as SLiM (for a more detailed, individual-based model) or ms/msprime (for a coalescent-based approach, which is computationally faster for neutral scenarios). The output is a set of ancestral recombination graphs (ARGs) or genotype data for all simulated individuals.
• Generate Sequence Alignments: Convert the simulated genealogies or genotypes into multiple sequence alignments (MSAs) for the extant taxa, which will serve as the direct input for the HMM-based introgression detection method [4].

HMM Analysis and Performance Benchmarking

The second phase involves applying the HMM method to the simulated data and evaluating its performance.

Protocol 2: Benchmarking the HMM Method

• Method Application: Run the HMM-based introgression detection method (e.g., PhyloNet-HMM) on the simulated MSAs. Provide the method with the correct set of possible parental species trees but not the true history of introgression [4].
• Output Parsing: For each genomic site or window, record the HMM's inference. This is typically the posterior probability of each possible parental species tree or a binary call (introgressed/not introgressed) based on a predefined probability threshold [66] [4].
• Comparison with Ground Truth: For each site/window, compare the HMM's inference with the known evolutionary history from the simulation. Classify the result into one of four categories:
  • True Positive (TP): Introgression was both simulated and detected.
  • False Positive (FP): Introgression was detected but not simulated (e.g., a region affected by ILS was mistaken for introgression).
  • True Negative (TN): No introgression was simulated and none was detected.
  • False Negative (FN): Introgression was simulated but not detected.

• Performance Calculation: Aggregate the counts from step 3 across the entire simulated genome to compute the following metrics:

  Power = TP / (TP + FN) FDR = FP / (TP + FP)

• Threshold Analysis: Repeat the classification and calculation (steps 3-4) across a range of probability thresholds used to call introgression. This generates the data needed to plot a power-FDR curve, which is essential for understanding the operational characteristics of the method.

Table 1: Example of Benchmarking Results for an HMM Method Under Different Evolutionary Scenarios

Simulation Scenario	Introgression Proportion	Power (%)	FDR (%)	Key Observation
Baseline (No ILS)	5%	98.5	1.2	High performance in simple case.
With High ILS	5%	85.3	8.7	ILS leads to increased false positives.
Ancient Introgression	2%	65.1	15.4	Power drops for older events.
Recent Introgression	10%	96.8	3.1	Method excels for recent hybridization.
Low Sequencing Coverage	5%	72.6	12.9	Data quality significantly impacts performance.

The Scientist's Toolkit

A successful validation study relies on a suite of computational tools and reagents. The following table outlines essential components for implementing the described protocols.

Table 2: Research Reagent Solutions for Simulation-Based Validation

Category & Item	Function in Validation	Examples & Notes
Simulation Software
ms / msprime	Coalescent-based simulator; fast for neutral evolution.	Ideal for large-scale benchmarking; models complex demography and ILS [4].
SLiM	Forward, individual-based simulator; detailed modelling of selection.	Better for scenarios involving adaptive introgression and complex phenotypes.
HMM Analysis Tools
PhyloNet-HMM	Detects introgression from MSAs using phylogenetic network HMMs.	Specifically designed to tease apart introgression from ILS [4].
Custom HMMs (e.g., in R)	Implements HMMs with local FDR control for ordered genomic data.	Provides flexibility for model customization as described in two-stage designs [66].
Analysis & Visualization
R / Python (Biopython)	Data analysis, statistics, and generation of power-FDR curves.	Biopython provides parsers for bioinformatics formats and integration with analysis pipelines [67].
GenomeDiagram (Biopython)	Visualizes annotated genomic features, such as predicted introgressed regions.	Useful for creating publication-quality figures of results along chromosomes [68].

Advanced Application: A Two-Stage HMM Design

A powerful extension in biomarker discovery is the two-stage HMM design, which optimizes cost-efficiency without sacrificing statistical rigor. This approach is highly relevant for drug development where validation budgets are constrained.

Protocol 3: Implementing a Two-Stage HMM with FDR Control

• Stage 1 - Discovery:

  • Apply the HMM to a large set of genomic features (e.g., all genes or taxa) from an initial, potentially expensive assay (e.g., whole-genome sequencing).
  • Calculate the local FDR (LIS) for each feature, which is the probability of the null hypothesis (no introgression/difference) given the data [66].
  • Select the top-ranking features for follow-up based on an optimized LIS threshold that ensures a high proportion of truly introgressed features is carried forward.

• Stage 2 - Validation:

  • Test the selected features in a new, independent sample set using a different, often more targeted, assay (e.g., PCR) [66].
  • Integrated Analysis: Do not discard the Stage 1 data. Instead, perform a joint analysis by combining the LIS values or test statistics from both stages. This integrated approach is more powerful than analyzing the second stage alone [66].

• Overall FDR Control: The threshold for significance in the joint analysis is determined to control the overall FDR across both stages at the desired level (e.g., 5%). This ensures the reliability of the final list of validated introgressed regions [66].

Diagram 2: Two-stage HMM design workflow for cost-effective validation.

Validation through forward simulations is an indispensable strategy for establishing the credibility of HMM-based inferences in comparative genomics. The protocols outlined herein provide a clear roadmap for benchmarking the power and false discovery rates of methods designed to detect introgression. By simulating data under realistic evolutionary scenarios that include challenges like ILS, and by rigorously quantifying performance, researchers and drug developers can have greater confidence in the genomic signals they report. The adoption of these validation practices, including advanced designs like the two-stage approach, is fundamental to advancing the field of comparative genomics and ensuring its findings are robust and reproducible.

The identification of introgressed genomic regions—genetic material transferred between species through hybridization—is crucial for understanding evolutionary processes, adaptation, and speciation. Accurately detecting these regions is methodologically challenging, as signals of introgression can be confounded by other evolutionary processes such as incomplete lineage sorting (ILS). The field has witnessed the development of diverse computational approaches to address this challenge, with Hidden Markov Models (HMMs) establishing themselves as a foundational methodology. More recently, neural networks, particularly Convolutional Neural Networks (CNNs), have emerged as powerful alternatives. This article provides a comparative analysis of these distinct methodological paradigms, framing them within the context of genomic introgression research and providing detailed protocols for their application.

Methodological Frameworks at a Glance

The table below summarizes the core attributes, strengths, and limitations of HMM-based and CNN-based approaches for introgression detection.

Table 1: Comparative Overview of HMM and Neural Network Approaches for Introgression Detection

Feature	HMM-Based Approaches	Neural Network (CNN) Approaches
Core Principle	Combines phylogenetic models with a probabilistic state transition system to model the genome as a sequence of hidden states (e.g., introgressed vs. not) [4].	Uses deep learning with multiple convolution layers to automatically learn informative spatial patterns from genotype matrices for classification [69].
Key Example Tools	PhyloNet-HMM [4] [26], hmmix [70]	genomatnn [69]
Handling of ILS	Explicitly accounts for incomplete lineage sorting within the phylogenetic network model [4].	Learns to distinguish patterns of introgression from ILS indirectly through training on simulated data that includes both processes [69].
Data Input	Multiple sequence alignments and a set of putative parental species trees [4]; or a single ingroup genome and an outgroup population [70].	A genotype matrix derived from donor, recipient, and outgroup populations, formatted as a 2D image-like structure [69].
Output	Probability of introgression for each genomic site and the evolutionary history of each site [4].	A probability score that a genomic region has undergone adaptive introgression [69].
Key Strengths	Model parameters are biologically interpretable (e.g., relate to coalescence times) [70]. Provides a coherent probabilistic framework that jointly models evolution and dependency.	High accuracy (>95% on simulated data) [69]. Does not require an explicit analytical model of evolution. Effective even with unphased genomic data [69].
Primary Limitations	Relies on a predefined set of parental species trees or a well-defined outgroup [4] [70]. Can be computationally intensive for complex models.	Function as a "black box," with limited direct biological interpretability [69]. Performance is heavily dependent on the quality and scope of training simulations.

Experimental Protocols and Workflows

Protocol for HMM-Based Introgression Detection (usinghmmix)

hmmix is a reference-free HMM method designed to detect introgressed segments without an archaic reference genome, instead leveraging an introgression-free outgroup population [70]. The following workflow details its application.

Workflow Diagram: HMM-based Introgression Detection with hmmix

Step-by-Step Protocol:

• Input Data Preparation

  • Ingroup and Outgroup Variants: Provide a VCF file containing genomic variants for the target (ingroup) individual(s) and the individuals from the introgression-free outgroup population [70].
  • Callability Mask: Provide a BED file indicating regions of the reference genome that can be confidently called, which is crucial for accurate parameter estimation [70].
  • Reference and Ancestral Genomes: Provide the species reference genome in FASTA format. An inferred ancestral genome sequence is optional but recommended [70].

• Identify Derived Variants in the Outgroup: Process the outgroup VCF to identify derived (non-ancestral) alleles. This step catalogs genetic variation that arose after the divergence of the ingroup and outgroup [70].

• Estimate Local Mutation Rate: Calculate a genome-wide mutation rate. hmmix can also incorporate local variation in mutation rates for more refined modeling [70].

• Find Private Variants in the Ingroup: Process the ingroup genome's VCF to identify "private" derived variants—those not found in the outgroup catalog. The density of these private variants is the key observational signal for the HMM. In non-introgressed regions, this density is low due to recent coalescence with the outgroup. In introgressed regions, which derive from an archaic population that diverged much earlier, the density of private variants is significantly higher [70].

• Train the HMM Parameters: Use the Expectation-Maximization (Baum-Welch) algorithm to train the HMM on the observed private variant data. This step learns the emission probabilities (which distinguish the "human" and "archaic" states based on SNP density) and the transition probabilities (which govern the switching between states along the chromosome) [70]. The trained parameters have direct biological interpretations: emissions relate to coalescence times (Tarchaic vs. Tingroup), and transitions relate to admixture proportion and time [70].

• Decode to Identify Archaic Segments: Use the Viterbi or posterior decoding algorithm on the trained HMM to identify the most likely sequence of hidden states (archaic or human) across the ingroup genome. The output is typically a BED file listing introgressed genomic segments, each annotated with statistics like posterior probability for downstream filtering (e.g., applying a posterior probability cutoff of 0.9) [70].

Protocol for CNN-Based Adaptive Introgression Detection (usinggenomatnn)

genomatnn is a deep learning method that frames introgression detection as an image recognition problem, using CNNs to scan genotype data for patterns indicative of adaptive introgression (AI) [69].

Workflow Diagram: CNN-based Adaptive Introgression Detection

Step-by-Step Protocol:

• Input Data and Genotype Matrix Construction: Gather sequence data from three populations: the donor (source of introgression), the recipient (population that received the genetic material), and an unadmixed outgroup (sister group to the recipient) [69]. For a given genomic region, construct a 2D genotype matrix where rows represent haplotypes (or diploid genotypes) and columns represent genomic positions (binned). The matrix entries contain counts of minor alleles.

• Genome Partitioning and Matrix Assembly: Partition the genome into sequential windows (e.g., 100 kbp in size). For each window, create a separate genotype matrix [69].

• Data Sorting and Concatenation: Within each population, sort the haplotypes/genotypes by their similarity to the donor population. Then, concatenate the sorted matrices from the donor, recipient, and outgroup populations into a single, composite genotype matrix. This creates the "image" that the CNN will analyze [69].

• CNN Training on Simulated Data: Train the CNN using a vast dataset of simulated genotype matrices. These simulations must encompass a wide range of evolutionary scenarios, including neutral evolution, selective sweeps, and adaptive introgression with various selection coefficients and onset times. The CNN architecture typically uses a series of convolutional layers with a 2x2 step size to efficiently extract high-level features, followed by fully connected layers that output the probability of adaptive introgression [69].

• Application to Empirical Data and Visualization: Apply the trained CNN to the concatenated genotype matrices from real genomic data. The model outputs a probability score for adaptive introgression in each genomic window. To aid interpretation, genomatnn can generate saliency maps, which highlight the regions of the input genotype matrix that most influenced the CNN's prediction, offering a glimpse into the "black box" [69].

Successful application of these bioinformatic methods requires a curated set of data resources and software tools.

Table 2: Essential Research Reagents and Resources for Introgression Detection

Resource Type	Specific Example(s)	Function and Application Notes
Software Packages	PhyloNet-HMM [4], hmmix [70], genomatnn [69]	Core analytical tools for implementing the HMM or CNN-based detection protocols.
Pre-computed Genomic Resources	Outgroup variation files, callability masks, ancestral sequences (e.g., provided for hg19/hg38 by hmmix [70])	Critical for running reference-free methods like hmmix without reprocessing raw data, ensuring reproducibility and saving computational time.
Reference Genomes	Species-specific reference assemblies (e.g., GRCh38 for human, GRCm39 for mouse)	Provides the coordinate system for alignment and variant calling; essential for all analyses.
Population Genomic Datasets	1000 Genomes Project, Human Genome Diversity Project (HGDP) [70]	Sources of ingroup and outgroup variation data for analysis, particularly in human evolutionary studies.
Simulation Frameworks	stdpopsim (with selection modules used by genomatnn) [69]	Generating realistic training data for CNN approaches and for benchmarking method performance under controlled evolutionary scenarios.

The choice between HMMs and neural networks for detecting introgression is not a matter of one being universally superior, but rather a strategic decision based on the biological question and available data.

HMMs offer a principled, model-based approach. Their parameters are directly tied to evolutionary quantities like coalescence times and admixture proportions, providing not just detection but also biological interpretation [70]. Methods like PhyloNet-HMM are powerful for testing specific phylogenetic hypotheses where the relationships between parental lineages are known or can be proposed [4]. However, their performance can degrade if the underlying model is misspecified.

In contrast, CNNs represent a powerful, data-driven paradigm. They excel in complex discrimination tasks, achieving high accuracy by learning directly from data without requiring an exact analytical model of the evolutionary process [69]. This makes them particularly suited for scanning genomes to identify candidate regions for adaptive introgression, especially when the selective pressure is unknown or complex. Their main drawback is their "black box" nature, though techniques like saliency maps are helping to improve interpretability.

In conclusion, HMMs and neural networks are complementary tools in the evolutionary genomicist's toolkit. HMMs provide a transparent, interpretable framework for testing well-defined hypotheses, while CNNs offer robust, high-throughput detection capabilities for exploratory analysis. As the field progresses, the integration of these approaches—for instance, using CNN-based pre-screening to identify candidate regions for subsequent fine-scale HMM analysis—promises to further decode the complex genomic landscapes of introgression.

Within the broader context of hidden Markov models (HMMs) for comparative genomic research, the empirical validation of predicted introgressed loci represents a critical step in confirming evolutionary history. The expansion of genomic datasets across diverse taxa, coupled with advances in methodological development, has created new opportunities to investigate the impact of introgression along individual genomes [3]. This application note details protocols for confirming known introgressed regions in model organisms, focusing on the integration of HMM-based frameworks with empirical data to validate evolutionary inferences. We present detailed methodologies, quantitative validation data, and practical implementation tools to support researchers in confirming introgression signals, with a particular emphasis on the well-characterized Vkorc1 locus in house mice (Mus musculus domesticus) [4].

HMM-Based Genomic Frameworks for Introgression Detection

Hidden Markov Models provide a powerful computational framework for detecting introgressed genomic regions by modeling the underlying genealogical transitions along chromosomal sequences. These methods effectively distinguish true introgression signatures from spurious signals arising from convergent evolutionary processes such as incomplete lineage sorting (ILS) [4].

The PhyloNet-HMM framework represents a significant methodological advance that combines phylogenetic networks with HMMs to simultaneously capture potentially reticulate evolutionary history and genomic dependencies [4]. This approach explicitly accounts for ILS and dependence across loci, enabling more accurate detection of introgressed segments. Similarly, Ancestry_HMM-S extends this foundation to specifically identify adaptive introgression by quantifying the strength of selection acting on introgressed loci [29].

Table 1: Key HMM-Based Methods for Introgression Detection

Method Name	Computational Approach	Key Features	Applicable Scenarios
PhyloNet-HMM	Phylogenetic networks + HMM	Accounts for ILS and dependencies across loci; provides fine-scale insights	General introgression detection across diverse species [4]
Ancestry_HMM-S	HMM with selection modeling	Quantifies selection strength on introgressed alleles; models ancestry transitions	Adaptive introgression in admixed populations [29]
IntroUNET	Semantic segmentation with deep learning	Identifies introgressed alleles in individual genomes; handles ghost introgression	Pixel-level classification of introgressed alleles in two-population alignments [71]

Figure 1: HMM-Based Computational Workflow for Introgression Detection. The pipeline begins with aligned genomic sequences, processes them through specialized HMM frameworks, and outputs introgressed regions with probability scores for empirical validation.

Experimental Protocol: Validating Introgressed Loci Using HMM Frameworks

Protocol 1: PhyloNet-HMM Analysis for Genome-Wide Introgression Scanning

Purpose: To systematically detect and validate introgressed regions across entire chromosomes, distinguishing true introgression from ILS.

Input Data Requirements:

• Whole-genome alignment data from at least three populations/species (donor, recipient, and outgroup)
• Phased haplotypes for optimal resolution
• Pre-defined set of parental species trees representing potential evolutionary histories

Methodology:

• Data Preparation: Format multiple sequence alignments in FASTA or VCF format, ensuring consistent coordinate systems across genomes.
• Model Configuration: Specify the set of possible parental species trees based on established phylogenetic relationships.
• Parameter Training: Use dynamic programming algorithms paired with multivariate optimization heuristics to train the model on genomic data [4].
• Genome Scanning: Execute the PhyloNet-HMM algorithm to compute for each genomic site the probability that it evolved under a specific parental species tree (Equation 1) [4].
• Region Identification: Identify contiguous genomic regions with significantly elevated probabilities for the introgressive parental tree.
• Statistical Validation: Compare detected regions against known introgressed loci (e.g., Vkorc1 in mice) to confirm methodological accuracy.

Output Interpretation:

• Genomic coordinates of introgressed regions with posterior probabilities
• Estimated percentage of sites of introgressive origin
• Gene annotations within detected regions for functional assessment

Protocol 2: Ancestry_HMM-S for Detecting Adaptive Introgression

Purpose: To identify and quantify selection on introgressed loci, distinguishing neutral from adaptive introgression.

Input Data Requirements:

• Genomic samples from an admixed focal population
• Two unadmixed, ancestral reference populations
• Genetic map or recombination rate estimates

Methodology:

• Local Ancestry Inference: Apply Ancestry_HMM to infer local ancestry patterns across genomes of admixed individuals [29].
• Selection Modeling: Expand the framework to explicitly model impacts of natural selection during admixture (Ancestry_HMM-S) [29].
• Parameter Estimation: Infer selection coefficients for individual loci based on deviations from neutral ancestry tract length distributions.
• Outlier Detection: Identify loci with significantly elevated proportions of introgressing ancestry relative to genome-wide baseline.
• Convergence Assessment: Validate findings through forward simulations under plausible selection scenarios.

Output Interpretation:

• Estimated selection coefficients for candidate regions
• Posterior probabilities for adaptive introgression
• Significance assessment through false discovery rate control

Empirical Validation Data from Model Organisms

Empirical validation of HMM-based introgression detection methods has been successfully demonstrated in several model organisms, with the house mouse system providing particularly compelling evidence.

Table 2: Empirical Validation of Introgressed Loci in Model Organisms

Organism	Introgressed Locus	HMM Method Applied	Validation Outcome	Functional Significance
House mouse (Mus musculus domesticus)	Vkorc1 (Chromosome 7)	PhyloNet-HMM	Accurately detected previously reported adaptive introgression event [4]	Rodent poison resistance [4]
House mouse (Mus musculus domesticus)	Chromosome 7 (genome-wide)	PhyloNet-HMM	9% of sites (13 Mbp, >300 genes) showed introgressive origin [4]	Potential adaptive traits beyond Vkorc1
Drosophila melanogaster (South Africa population)	17 identified loci	Ancestry_HMM-S	Detected signatures of adaptive introgression; 4 loci confer insecticide resistance [29]	Cyp6 P450 genes (DDT resistance); Ace (organophosphate resistance)

Application of PhyloNet-HMM to chromosome 7 variation data in house mice demonstrated the method's validation capabilities, successfully recovering the previously reported adaptive introgression event involving the Vkorc1 gene, which confers resistance to rodenticides [4]. Beyond this known locus, the genome-wide scan revealed approximately 9% of sites within chromosome 7 (covering about 13 Mbp and over 300 genes) showed evidence of introgressive origin, suggesting a more extensive history of hybridization than previously recognized [4].

In Drosophila melanogaster, Ancestry_HMM-S applied to an admixed population from South Africa identified 17 loci with signatures of adaptive introgression, four of which corresponded to previously known insecticide resistance genes including Cyp6 P450 genes (conferring DDT resistance) and Ace (conferring organophosphate resistance) [29]. This validation against known resistance loci confirms the method's accuracy in detecting biologically relevant adaptive introgression.

Figure 2: Empirical Validation Feedback Loop. Known introgressed loci provide ground truth for validating HMM methods, which in turn detect new regions, creating an expanding knowledge base for understanding genomic introgression.

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Research Reagents and Computational Tools for Introgression Validation

Resource Type	Specific Tool/Reagent	Function in Validation Protocol	Implementation Notes
Computational Software	PhyloNet-HMM	Detects introgression in genomes while accounting for ILS	Open-source; part of PhyloNet distribution [4]
Computational Software	Ancestry_HMM-S	Identifies adaptive introgression and quantifies selection	Available at https://github.com/jesvedberg/Ancestry_HMM-S/ [29]
Reference Data	Parental species trees	Defines possible evolutionary histories for HMM	Must be based on established phylogenetic relationships [4]
Reference Data	Annotated reference genome	Provides genomic coordinates and gene annotations	Enables functional interpretation of detected regions [4]
Validation Standard	Known introgressed loci (e.g., Vkorc1)	Positive control for method validation	Confirms methodological accuracy in empirical tests [4]

The empirical validation of introgressed loci using HMM-based frameworks represents a robust methodology for confirming evolutionary gene flow events. Through the application of specialized tools like PhyloNet-HMM and Ancestry_HMM-S to model organism systems, researchers can confidently identify introgressed regions while distinguishing true introgression from confounding signals like ILS. The successful detection of known adaptive loci such as Vkorc1 in mice and insecticide resistance genes in Drosophila provides strong validation of these approaches and enables more accurate characterization of the genomic landscapes of introgression across diverse taxa.

Evaluating Performance in Distinguishing Introgression from ILS

The precise identification of genomic regions affected by introgression—the transfer of genetic material between species through hybridization—is a cornerstone of modern evolutionary genomics. However, this task is significantly complicated by incomplete lineage sorting (ILS), a phenomenon where the genealogical histories of different genes conflict with the species tree due to the retention of ancestral polymorphisms [72]. Distinguishing between the genomic signatures of ILS and introgression is methodologically challenging, as both processes can produce similar patterns of shared genetic variation [73]. This application note examines the performance of various computational approaches for discriminating between these evolutionary forces, with a specific focus on the growing role of hidden Markov models (HMMs) and related probabilistic frameworks within comparative genomics research.

The challenge is particularly pronounced in closely related species and those with large effective population sizes, where ILS is more common. As genomic datasets expand across diverse taxa, researchers have developed increasingly sophisticated methods to tease apart these confounding signals. Understanding the relative strengths and limitations of these approaches is essential for making accurate inferences about evolutionary history, including the identification of adaptively introgressed loci [3].

Performance Evaluation of Methodological Approaches

Current methods for detecting introgression can be broadly classified into three categories, each with distinct theoretical foundations and performance characteristics [3]. Summary statistics approaches calculate metrics based on sequence divergence and polymorphism to identify outliers suggestive of introgression. Probabilistic modeling methods, including HMMs, explicitly incorporate evolutionary processes to infer introgression probabilities across genomes. Supervised learning represents an emerging paradigm that frames introgression detection as a classification problem, often achieving high accuracy when sufficient training data exists.

Quantitative Performance Comparison

The table below summarizes the key performance characteristics of representative methods across these categories:

Table 1: Performance Comparison of Methods for Detecting Introgression

Method Category	Specific Method/Statistic	Key Performance Strengths	Key Limitations	Power/Accuracy Notes
Summary Statistics	dXY, FST [74]	Robust to some selection effects; simple computation.	Low sensitivity to recent or low-frequency migration.	High power only when migration is recent and strong.
Summary Statistics	dmin [74]	High power when assumptions are met; sensitive to rare lineages.	Assumes no mutation rate variation; sensitive to parameter violations.	Powerful for detecting even rare introgressed lineages.
Summary Statistics	RNDmin [74]	Robust to mutation rate variation; reliable with inaccurate divergence times.	Requires phased haplotypes and outgroup.	Modest power increase over related tests; remains reliable.
Probabilistic Modeling	PhyloNet-HMM [4]	Accounts for ILS and dependence across loci; provides fine-scale insights.	Computationally intensive; requires specification of parental species trees.	Accurately detects introgression in simulated and empirical data (e.g., mouse Vkorc1).
Supervised Learning	Semantic Segmentation [3]	Great potential for complex scenarios; can integrate diverse features.	Dependent on training data quality and representativeness.	Emerging approach with promising but not yet fully quantified performance.

Impact of Genomic Context on Performance

The performance of all methods is context-dependent, influenced by genomic landscape features. Research in primates shows that ILS proportions increase with recombination rate and are reduced in exons compared to intergenic regions (approximately 31% reduction in local effective population size) [75]. This variation is primarily driven by selection at linked sites, meaning methods must account for heterogeneous genomic landscapes to avoid false positives or negatives. Regions under strong selective constraint or with low recombination rates can mimic signatures of introgression by reducing genetic diversity, thereby confounding methods that do not explicitly model these forces [75].

Detailed Experimental Protocols

The RNDmin statistic offers a robust approach for detecting introgressed regions between sister species while accounting for mutation rate variation [74].

Research Reagent Solutions:

• Genomic Data: Phased haplotype data from two sister species and an outgroup species.
• Software Tools: Population genomic analysis pipelines (e.g., custom Python/R scripts).
• Computing Resources: Standard computational workstation sufficient for summary statistic calculations.

Procedure:

• Data Preparation: Obtain whole-genome, phased haplotype data for Population X, Population Y (sister species), and an outgroup species O. The outgroup should not have experienced introgression with X or Y.
• Calculate Minimum Distance: For each genomic window, compute dmin = min{d(x,y)}, the minimum number of sequence differences between any haplotype in X and any haplotype in Y.
• Calculate Outgroup Divergence: Compute the average distance from X and Y to the outgroup O: d_out = (d(X,O) + d(Y,O)) / 2.
• Compute RNDmin: Calculate the ratio RNDmin = dmin / d_out.
• Identify Outliers: Scan the genome for regions where the RNDmin value falls significantly below the genomic background. These regions represent candidates for recent introgression.
• Validation: Compare candidates with other lines of evidence, such as phylogenetic incongruence or functional annotation.

Protocol for HMM-Based Analysis (PhyloNet-HMM)

PhyloNet-HMM provides a powerful probabilistic framework for detecting introgression while explicitly accounting for ILS and dependencies between genomic sites [4].

Research Reagent Solutions:

• Genomic Data: A multiple sequence alignment from multiple individuals across at least three species.
• Evolutionary Hypothesis: A set of putative parental species trees (including possible introgressive histories).
• Software Tools: PhyloNet software package (available from its distribution site).
• Computing Resources: High-performance computing cluster may be necessary for genome-scale analyses.

Procedure:

• Input Preparation: Compile a whole-genome multiple sequence alignment (O(1)...O(L)) for genomes from the focal and related species. Define the set of possible parental species trees (S1...SK) based on the evolutionary hypotheses to be tested.
• Model Training: Use the PhyloNet-HMM algorithm to train the model parameters. This involves employing a dynamic programming approach paired with an optimization heuristic to estimate the transition and emission probabilities of the HMM.
• Probability Calculation: For each site i in the alignment, compute the posterior probability P(H(i) = Sk | O) that the site evolved under the parental species tree Sk. This is achieved using the forward-backward algorithm for HMMs [45].
• State Decoding: Use the Viterbi algorithm to find the most likely sequence of hidden states (i.e., the most probable parental tree at each position) across the genome [45].
• Interpretation: Identify genomic regions where the most probable history corresponds to an introgressive species tree. The distribution of lengths of these regions can provide insights into the timing and selection acting on the introgressed haplotypes.

The logical workflow and dependencies of the PhyloNet-HMM method are illustrated below:

Figure 1: PhyloNet-HMM Analysis Workflow. The diagram outlines the key steps in applying PhyloNet-HMM to distinguish introgression from ILS.

Protocol for Simulation-Based Performance Assessment

Simulation studies are crucial for evaluating the statistical power and false positive rates of different methods under controlled conditions with known evolutionary parameters [72] [75].

Research Reagent Solutions:

• Software Tools: Coalescent simulation software with recombination and introgression capabilities (e.g., ms or msprime).
• Analysis Pipelines: Implementation of the introgression detection methods to be evaluated.
• Computing Resources: High-performance computing cluster for large-scale simulations.

Procedure:

• Parameterization: Define realistic demographic parameters (effective population size, divergence times) and introgression parameters (time, rate, and direction of gene flow).
• Data Simulation: Use coalescent simulator to generate hundreds of replicate genomic sequences under models with both ILS and introgression. Also generate control datasets under a null model with ILS only.
• Method Application: Run each introgression detection method (e.g., summary statistics, PhyloNet-HMM) on all simulated datasets.
• Power Calculation: Calculate the proportion of true introgressed regions correctly identified by each method (True Positive Rate).
• False Positive Assessment: Calculate the proportion of regions falsely identified as introgressed in the null (ILS-only) simulations (False Positive Rate).
• Performance Comparison: Compare methods using metrics like Receiver Operating Characteristic (ROC) curves, considering factors like introgression timing, rate, and selection strength.

The Scientist's Toolkit

Table 2: Essential Research Reagents and Tools for Introgression Analysis

Item Name	Function/Application	Example Use Case
Phased Haplotype Data	Required for methods like dmin and RNDmin that rely on haplotype comparisons.	Analyzing population-level data to detect rare introgressed haplotypes.
Outgroup Genome	Provides an independent estimate of divergence to normalize for mutation rate variation.	Rooting the phylogenetic analysis and calculating statistics like RNDmin.
PhyloNet Software	Implements phylogenetic network models, including the PhyloNet-HMM framework.	Detecting introgression from genome alignments while accounting for ILS.
Coalescent Simulator	Generates synthetic genomic data under evolutionary models with ILS and introgression.	Benchmarking method performance and evaluating statistical power.
Reference Genome Annotation	Provides functional context for identified genomic regions (e.g., genes, regulatory elements).	Interpreting the potential phenotypic impact of candidate introgressed regions.

Integrated Analysis Workflow

To robustly distinguish introgression from ILS in empirical datasets, an integrated, multi-method approach is recommended. The conceptual relationship between different evolutionary processes and their genomic signatures is complex, as visualized below:

Figure 2: Relationship Between Evolutionary Processes and Detection Methods. The diagram shows how different processes create shared genetic variation and which methods are designed to disentangle them.

A robust analytical workflow should:

• Apply Multiple Methods: Use both summary statistics (e.g., RNDmin) and model-based approaches (e.g., PhyloNet-HMM) on the same dataset. Consistent signals across methods provide stronger evidence for introgression.
• Compare Allopatric and Parapatric Populations: When geographic data is available, compare genetic differentiation between species in zones of sympathy/parapatry versus allopatry. Significantly lower differentiation in parapatry suggests ongoing gene flow rather than ILS alone [73].
• Leverage Multi-omic Data: Integrate information from differently inherited markers (e.g., biparental nuclear, maternal mitochondrial, paternal chloroplast). Incongruent patterns of sharing can help distinguish ancient ILS from sex-biased introgression [73].
• Perform Coalescent Simulations: Simulate data under both ILS-only and ILS-with-introgression models specific to the study system to quantify the false discovery rate and ensure the chosen methods have power under the inferred demographic history.

Distinguishing introgression from ILS remains a central challenge in evolutionary genomics. Summary statistics like RNDmin provide a robust and computationally efficient first pass for detecting introgressed regions, particularly when migration is recent and strong. However, probabilistic models, especially HMM-based frameworks like PhyloNet-HMM, offer superior performance for fine-scale inference by simultaneously accounting for ILS, mutation, recombination, and ancestral polymorphism. The emerging field of supervised learning shows great promise for handling increasingly complex evolutionary scenarios. Ultimately, the most reliable inferences will come from integrating multiple methodological approaches, leveraging complementary strengths to illuminate the complex history of gene flow that has shaped the genomes of species across the tree—and network—of life.

Strengths and Limitations of Current HMM Tools in Genomics

Hidden Markov Models (HMMs) have become fundamental tools in modern genomic analysis, often described as the "workhorse" of biological sequence investigation [47]. Their application is particularly crucial in the evolving field of comparative genomic introgression research, which seeks to identify genetic material transferred between species through hybridization and back-crossing. In this context, HMMs provide a powerful statistical framework for deciphering the complex evolutionary histories embedded within genomic sequences. The ability of HMMs to model dependencies between adjacent sites in a genome makes them exceptionally suited for detecting the mosaic patterns characteristic of introgressed regions, which appear as contiguous segments with distinct phylogenetic signatures [4] [47].

As genomic datasets expand across diverse taxa, the demand for sophisticated analytical methods has intensified. Current research aims to quantify variation along chromosomes, account for it during inference procedures, and ultimately determine the causal evolutionary processes behind it [47]. HMM-based approaches are at the forefront of this endeavor, enabling researchers to distinguish true introgression signals from confounding factors such as incomplete lineage sorting (ILS), recombination, and ancestral polymorphism [4]. This Application Note examines the current capabilities and constraints of HMM tools in genomics, with a specific focus on their application in comparative genomic introgression research, providing detailed protocols and analytical frameworks for researchers and drug development professionals.

Current HMM Tool Landscape in Genomics

The landscape of HMM-based genomic tools has diversified significantly, with implementations tailored to specific evolutionary questions and data types. Within introgression research, two primary categories of HMM approaches have emerged: those incorporating explicit phylogenetic network models and those designed for genomic landscape characterization.

Phylogenetic Network-Integrated HMMs

PhyloNet-HMM represents a sophisticated framework that combines phylogenetic networks with HMMs to detect introgression while simultaneously accounting for ILS and dependence across loci [4]. This method scans multiple aligned genomes for signatures of introgression by modeling the evolutionary history of genomic regions as a hidden state process, where each state corresponds to a different phylogenetic history, including those involving reticulate evolution. The model can identify regions of introgressive descent by calculating the probability that a given genomic site evolved under a specific phylogenetic network, enabling the detection of even ancient introgression events [4].

Application studies using PhyloNet-HMM on chromosome 7 variation data from house mice (Mus musculus domesticus) successfully detected a previously reported adaptive introgression event involving the rodent poison resistance gene Vkorc1, along with numerous newly identified introgressed regions [4]. The analysis estimated that approximately 9% of sites within chromosome 7 (covering about 13 Mbp and over 300 genes) were of introgressive origin, demonstrating the method's ability to provide genome-wide quantification of introgression.

Genomic Landscape Characterization HMMs

For broader genomic landscape analysis, tools like oHMMed (ordered HMM with emission densities) have been developed to segment genomes into comparably homogeneous regions based on autocorrelated observed sequences [47]. This approach models genomic features such as base composition, recombination rates, and gene density as emissions from hidden states with ordered means, effectively capturing the inherent autocorrelation patterns along chromosomes while remaining agnostic to specific biological mechanisms.

oHMMed employs a restricted, tridiagonal transition matrix that only allows transitions between neighboring states in the ordered sequence, effectively modeling the autocorrelation patterns observed in genomic features like GC content [47]. This design avoids the "label switching" problem common in MCMC methods and reduces the number of estimable parameters, improving algorithm behavior and guarding against over-fitting.

Table 1: Key HMM Tools for Genomic Introgression Research

Tool Name	Methodological Approach	Primary Application	Key Features
PhyloNet-HMM	Phylogenetic networks + HMMs	Introgression detection in comparative genomics	Accounts for ILS and dependence across loci; Models phylogenetic network history
oHMMed	Ordered HMM with emission densities	Genomic landscape segmentation	Orders hidden states by emission means; Restricted transitions between neighboring states; Models autocorrelation
Saguaro	HMMs with artificial neural networks	Local phylogenetic incongruence	Categorizes local genealogies; Combines multiple analytical approaches

Strengths of HMM Approaches

HMM-based methods offer several distinct advantages for genomic analysis, particularly in the complex domain of introgression research.

Handling Dependencies and Autocorrelation

A fundamental strength of HMMs is their ability to model dependencies between adjacent sites in genomic sequences [4] [47]. Unlike methods that assume independence across loci, HMMs explicitly account for the autocorrelation inherent in genomic data due to factors such as linkage and variable recombination rates. This capability is crucial for accurately identifying contiguous introgressed regions and distinguishing them from sporadic phylogenetic incongruences caused by other evolutionary forces.

The oHMMed framework specifically leverages this autocorrelation to identify biologically meaningful segments in genomic landscapes, with applications demonstrating effective segmentation of GC content, gene density, and epigenetic markers along chromosomes [47]. By modeling the sequential nature of genomic data, HMMs can more accurately reconstruct the boundaries of introgressed segments and provide better estimates of their lengths and distributions.

Distinguishing Complex Evolutionary Processes

Modern HMM implementations excel at teasing apart confounding signals from different evolutionary processes. PhyloNet-HMM specifically addresses the challenge of distinguishing introgression from incomplete lineage sorting, both of which can produce similar patterns of phylogenetic incongruence [4]. By incorporating a model that explicitly includes both processes, the method can attribute genomic patterns to their most likely evolutionary cause, significantly reducing false positive introgression calls.

This discriminatory power was validated through analysis of negative control datasets, where PhyloNet-HMP correctly detected no introgression, and through accurate performance on synthetic datasets simulated under the coalescent model with recombination, isolation, and migration [4]. Such validation demonstrates the robustness of well-designed HMM frameworks for making accurate evolutionary inferences.

Flexibility and Modeling Power

The modular structure of HMMs allows for adaptation to diverse evolutionary scenarios and data types. The emission distributions can be customized for different types of genomic data, from discrete character matrices to continuous measurements of genomic features [47]. Similarly, the state space can be designed to represent various evolutionary histories, including those involving reticulation events.

The ordered HMM approach exemplified by oHMMed demonstrates how incorporating biological knowledge into model structure (through state ordering and restricted transitions) can improve performance and interpretability while reducing parameter space [47]. This flexibility ensures that HMM approaches can continue to evolve alongside genomic technologies and research questions.

Table 2: Quantitative Performance of HMM Tools in Genomic Studies

Tool/Study	Data Analyzed	Detection Accuracy	Key Findings
PhyloNet-HMM	Mouse chromosome 7	Detected known Vkorc1 introgression; 9% of sites introgressed	13 Mbp introgressed covering 300+ genes; No false positives in negative control
oHMMed	Human, mouse, fruit fly genomes	Effective segmentation of GC content and gene density	Identified compositionally distinct regions; Characterized continuous variation patterns
General HMM applications	Various eukaryotic taxa	Accurate for many clades	Revealed loci linked to immunity, reproduction, and environmental adaptation

Limitations and Challenges

Despite their considerable strengths, current HMM approaches face several important limitations that researchers must consider when designing genomic studies.

Model Complexity and Computational Demands

The sophisticated statistical models underlying tools like PhyloNet-HMM come with substantial computational requirements [4]. The need to integrate over possible phylogenetic networks and account for multiple evolutionary processes simultaneously creates significant computational challenges, particularly as the number of taxa or genomic scale increases. This can limit application to very large genomic datasets or complex phylogenetic scenarios with many potential hybridization events.

For oHMMed, while the restricted state transitions reduce parameter space, the MCMC sampling approach still requires substantial computational resources for genome-wide analyses [47]. These computational constraints often force practical compromises in analysis, such as analyzing chromosomes separately or reducing sampling intensity, which may affect result quality.

Dependence on Model Specification

HMM performance is highly dependent on correct model specification, including appropriate state space design, emission distributions, and transition structures. Misspecification at any of these levels can lead to biased results or failure to detect true biological signals. For example, in introgression detection, failure to include the correct phylogenetic networks in the model can prevent identification of true introgressed regions or generate false positives [4].

The problem of selecting the appropriate number of hidden states is particularly challenging, with oHMMed developers noting that "development of model selection criteria for HMMs is complex and no consensus criteria exist" [47]. While oHMMed proposes intuitive diagnostic criteria for state number selection, this remains an area of methodological difficulty across HMM applications.

Limitations in Specific Evolutionary Contexts

Current HMM tools have specific limitations in certain evolutionary contexts. Methods like PhyloNet-HMM primarily focus on distinguishing introgression from ILS but may not adequately account for other confounding processes such as gene duplication and loss, or selection [4]. Additionally, most methods assume simple hybridization scenarios and may struggle with complex histories involving multiple introgression events, ghost lineages, or high levels of subsequent recombination.

The assumption of a fixed set of parental species trees in PhyloNet-HMM may also limit applicability in scenarios where the parental lineages are poorly defined or where extensive ghost introgression has occurred [4] [3]. Similarly, the ordered state assumption in oHMMed, while beneficial for modeling autocorrelated genomic features, may not be appropriate for all types of genomic variation.

Experimental Protocols

Protocol 1: Detecting Introgression with PhyloNet-HMM

Application: Identifying genomic regions of introgressive descent in comparative genomic data.

Principle: The method combines phylogenetic networks with HMMs to compute the probability that each genomic site evolved under a specific evolutionary history, including those involving hybridization events [4].

Materials and Reagents

Table 3: Research Reagent Solutions for PhyloNet-HMM Analysis

Reagent/Resource	Specification	Function/Purpose
Genomic sequences	Multiple aligned genomes from target and reference species	Primary data for introgression detection
PhyloNet software package	Version including PhyloNet-HMM implementation	Core analytical tool for network-based HMM analysis
Parental species trees	User-specified set of possible species trees	Defines constrained evolutionary scenarios for the HMM
Computational resources	High-performance computing cluster recommended	Handles computationally intensive phylogenetic calculations

Procedure

• Data Preparation and Alignment

  • Obtain whole-genome sequences for all studied individuals/species
  • Perform multiple sequence alignment using appropriate tools (e.g., MAFFT, MUSCLE)
  • Format alignment data in compatible format (e.g., FASTA, NEXUS)

• Model Specification

  • Define set of possible parental species trees based on established phylogenetic relationships
  • Specify phylogenetic network models to be considered, including potential hybridization events
  • Set HMM parameters including state space (possible local genealogies) and transition probabilities

• Parameter Estimation and Training

  • Use dynamic programming algorithms to compute probabilities of hidden states given observed data [4]
  • Apply multivariate optimization heuristic to train model parameters on genomic data
  • Iterate until parameter convergence or maximum iterations reached

• Decoding and Interpretation

  • Compute for each site the probability of each parental species tree using forward-backward algorithm [4]
  • Identify genomic regions with high probability of introgressive origin (e.g., posterior probability > 0.95)
  • Annotate introgressed regions with gene content and functional information

• Validation and Downstream Analysis

  • Perform negative control analysis using datasets without expected introgression
  • Compare results with known introgressed loci if available
  • Investigate functional implications of identified introgressed regions

PhyloNet-HMM Workflow for Introgression Detection

Protocol 2: Genomic Landscape Segmentation with oHMMed

Application: Partitioning genomes into compositionally distinct regions based on autocorrelated features.

Principle: oHMMed models observed genomic features (e.g., GC content, gene density) as emissions from hidden states with ordered means, using a restricted transition matrix that only permits moves between neighboring states [47].

Materials and Reagents

Table 4: Research Reagent Solutions for oHMMed Analysis

Reagent/Resource	Specification	Function/Purpose
Genomic sequence or feature data	Reference genome or measured genomic features	Input data for segmentation analysis
oHMMed R package	Available on CRAN	Implementation of ordered HMM with emission densities
Computational environment	R statistical computing platform	Execution environment for oHMMed algorithms
Window definition parameters	Size and step for sliding window analysis	Defines resolution of genomic landscape analysis

Procedure

• Data Preparation and Windowing

  • Obtain genomic sequence or feature data for analysis
  • Divide genome into consecutive windows (size depends on analysis goals)
  • Calculate feature of interest (GC content, gene count, epigenetic marks) for each window

• Model Initialization

  • Specify number of hidden states K (can use diagnostic criteria for selection)
  • Define emission distributions (normal mixture for continuous data, Poisson-gamma for counts)
  • Initialize state-specific emission parameters with ordered means
  • Set up tridiagonal transition matrix allowing only neighboring state transitions

• Parameter Inference via MCMC

  • Run Markov Chain Monte Carlo sampling to obtain posterior distributions
  • Iterate until convergence (monitor likelihood and parameter stability)
  • Use label matching to maintain consistent state ordering across iterations

• Genome Annotation and Segmentation

  • Compute marginal probabilities of hidden states at each genomic position
  • Annotate each window with its most probable hidden state
  • Identify boundaries between different genomic segments

• Downstream Analysis and Interpretation

  • Characterize segmented regions by their emission parameters
  • Compare segmentation across different genomic features or species
  • Investigate biological correlates of distinct genomic segments

oHMMed Model Structure with Ordered States

Future Directions

The field of HMM development for genomic analysis continues to evolve rapidly, with several promising directions emerging. Integration with supervised learning approaches represents a frontier where HMMs could be combined with neural networks or other machine learning techniques to improve detection accuracy [3]. As noted in recent reviews, "Supervised learning is an emerging approach with great potential, particularly when the detection of introgressed loci is framed as a semantic segmentation task" [3].

Methodological improvements focusing on computational efficiency will be crucial for handling the increasing scale of genomic datasets. Development of approximate inference techniques, distributed computing implementations, and cloud-native solutions will make sophisticated HMM analyses accessible to more research groups. Additionally, extending HMM frameworks to better handle complex evolutionary scenarios involving multiple introgression events, ghost lineages, and varying selection pressures will enhance their biological realism and applicability.

The integration of HMMs with functional genomic data represents another promising direction, allowing evolutionary inferences to be connected with regulatory annotations, chromatin states, and phenotypic associations. Such integration could ultimately support drug development professionals in identifying functionally relevant introgressed regions that may contribute to disease susceptibility or treatment response.

HMM-based tools have established themselves as indispensable for genomic introgression research, providing powerful frameworks for detecting evolutionary signals in complex genomic data. Tools like PhyloNet-HMM and oHMMed demonstrate the versatility of HMM approaches, from specific introgression detection to broader genomic landscape characterization. While challenges remain in computational efficiency, model specification, and handling extreme evolutionary complexity, ongoing methodological developments continue to expand the capabilities of these approaches. As genomic datasets grow in size and diversity, HMM-based methods will remain essential for deciphering the intricate evolutionary histories written in genomes, with important implications for understanding biodiversity, adaptation, and the genetic basis of phenotypic variation.

Conclusion

Hidden Markov Models have firmly established themselves as a powerful and flexible statistical framework for detecting introgression in comparative genomics. By effectively teasing apart the confounding effects of incomplete lineage sorting and capturing the mosaic structure of admixed genomes, tools like Ancestry_HMM-S and PhyloNet-HMM enable the precise identification of adaptively introgressed regions, as demonstrated in cases of pesticide resistance and pathogen immunity. Future directions will involve refining these models to handle more complex demographic histories, integrating with deep learning approaches for enhanced power, and expanding applications into clinical and pharmaceutical domains. The ability to pinpoint adaptive introgression events holds profound implications for understanding evolutionary mechanisms, identifying genetic variants underlying complex diseases, and discovering new drug targets, ultimately bridging evolutionary biology with precision medicine.

References

• 1. Adaptive Introgression as an Evolutionary Force: A Meta ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12181397/]
• 2. Adaptive Introgression: An Untapped Evolutionary ... [https://www.frontiersin.org/journals/plant-science/articles/10.3389/fpls.2019.00004/full]
• 3. Review Decoding genomic landscapes of introgression [https://www.sciencedirect.com/science/article/pii/S0168952525001660]
• 4. An HMM-Based Comparative Genomic Framework for ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC4055573/]
• 5. A Hidden Markov Model Approach for Simultaneously ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5242547/]
• 6. Enhancement of hidden Markov model analyses ... [https://carta.anthropogeny.org/libraries/bibliography/enhancement-hidden-markov-model-analyses-improved-inference-archaic]
• 7. Performance evaluation of adaptive introgression ... [https://peercommunityjournal.org/articles/10.24072/pcjournal.617/]
• 8. Introgression shapes fruit color convergence in invasive ... [https://elifesciences.org/articles/64165]
• 9. Hybrid introgression as a mechanism of rapid evolution ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12331986/]
• 10. Introgression impacts the evolution of bacteria, but species ... [https://www.nature.com/articles/s41467-025-64947-1]
• 11. Recent Applications of Hidden Markov Models in ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5172443/]
• 12. Hidden Markov Models and their Applications in Biological ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2766791/]
• 13. Hidden Markov model [https://en.wikipedia.org/wiki/Hidden_Markov_model]
• 14. A sequence-profile-based HMM for predicting and ... [https://pubmed.ncbi.nlm.nih.gov/12169530/]
• 15. Profile hidden Markov model sequence analysis can help ... [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-04180-x]
• 16. [1310.7989] An HMM-based Comparative Genomic ... [https://arxiv.org/abs/1310.7989]
• 17. A gentle tutorial on accelerated parameter and confidence ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC9796807/]
• 18. Hidden Markov Model in Machine learning [https://www.geeksforgeeks.org/machine-learning/hidden-markov-model-in-machine-learning/]
• 19. Hidden Markov Models Explained [https://www.digitalocean.com/community/tutorials/hidden-markov-models]
• 20. Understanding emission probability in HMM definition [https://stats.stackexchange.com/questions/539038/understanding-emission-probability-in-hmm-definition]
• 21. Optimizing Hidden Markov Language Models: An Empirical ... [https://aclanthology.org/2025.findings-naacl.429/]
• 22. Optimal determination of hidden Markov model parameters ... [https://www.sciencedirect.com/science/article/pii/S2468227622000825]
• 23. An HMM-Based Comparative Genomic Framework for ... [https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003649]
• 24. jesvedberg/Ancestry_HMM-S [https://github.com/jesvedberg/Ancestry_HMM-S]
• 25. Inferring Adaptive Introgression Using Hidden Markov ... [https://researchprofiles.ku.dk/en/publications/inferring-adaptive-introgression-using-hidden-markov-models]
• 26. An HMM-based comparative genomic framework for detecting ... [https://pubmed.ncbi.nlm.nih.gov/24922281/]
• 27. PhyloNet-HMM [https://phylogenomics.rice.edu/html/phyloHMM.html]
• 28. Phylogenetic Networks And Evolutionary History [https://www.nature.com/research-intelligence/nri-topic-summaries/phylogenetic-networks-and-evolutionary-history-micro-25647]
• 29. Inferring Adaptive Introgression Using Hidden Markov ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8097282/]
• 30. Refining models of archaic admixture in Eurasia with ... [https://www.nature.com/articles/s41467-021-26503-5]
• 31. Phylogenetic Placement of Aligned Genomes and ... [https://pubmed.ncbi.nlm.nih.gov/40811323/]
• 32. “Normal” phylogenetic networks may be emerging as the ... [https://arxiv.org/html/2107.10414v3]
• 33. Phenotypic signatures of incomplete lineage sorting in ... [https://pubmed.ncbi.nlm.nih.gov/41266166/]
• 34. Incomplete lineage sorting and phenotypic evolution in ... [https://www.sciencedirect.com/science/article/pii/S0092867422003440]
• 35. A Guide to Multi-omics Integration Strategies [https://frontlinegenomics.com/a-guide-to-multi-omics-integration-strategies/]
• 36. Adaptive Introgression: An Untapped Evolutionary Mechanism ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6367218/]
• 37. Signatures of Insecticide Selection in the Genome of ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6222576/]
• 38. Inferring multi-locus selection in admixed populations - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC10707604/]
• 39. Copy number variants underlie major selective sweeps in ... [https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3002898]
• 40. Ancestral Population Genomics: The Coalescent Hidden ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2746150/]
• 41. Warfarin Therapy and VKORC1 and CYP Genotype - NCBI [https://www.ncbi.nlm.nih.gov/books/NBK84174/]
• 42. Thirteen novel VKORC1 mutations associated with oral ... [https://www.sciencedirect.com/science/article/pii/S1538783622083167]
• 43. Pharmacogenetic testing of CYP2C9 and VKORC1 alleles ... [https://www.nature.com/articles/gim200821]
• 44. MGI:106442 - vitamin K epoxide reductase complex, subunit 1 [https://www.informatics.jax.org/marker/MGI:106442]
• 45. The hidden Markov model and its applications in ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12505677/]
• 46. Number of states in HMM - Cross Validated [https://stats.stackexchange.com/questions/492005/number-of-states-in-hmm]
• 47. Inference of genomic landscapes using ordered Hidden ... [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-024-05751-4]
• 48. Hidden Markov Models Tool Tutorial [https://iri.columbia.edu/climate/forecast/stochasticTools/HMM_Overview.html]
• 49. Hidden Markov Models methods for selecting optimal ... [https://stackoverflow.com/questions/23233219/hidden-markov-models-methods-for-selecting-optimal-number-of-states]
• 50. Overfitting, Underfitting and General Model Overconfidence ... [https://www.ncbi.nlm.nih.gov/books/NBK610560/]
• 51. Mastering Model Complexity: Avoiding Underfitting and ... [https://www.pecan.ai/blog/machine-learning-model-underfitting-and-overfitting/]
• 52. How to Solve Underfitting and Overfitting Data Models [https://allcloud.io/blog/how-to-solve-underfitting-and-overfitting-data-models/]
• 53. Methods to Control Underfitting and Overfitting - Rhydham Gupta [https://rhydhamgupta.medium.com/methods-to-control-underfitting-and-overfitting-72663901d069]
• 54. A History-Based Approach to Mitigate Overfitting [https://arxiv.org/html/2401.10359v1]
• 55. Master Natural Language Processing — Part 26 A Deep ... [https://medium.com/@conniezhou678/master-natural-language-processing-part-26-a-deep-dive-for-dynamic-programming-in-hidden-markov-d2997013098f]
• 56. Comprehensive noise reduction in single-cell data with the ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12570323/]
• 57. The effect of background noise and its removal on the analysis ... [https://genomebiology.biomedcentral.com/articles/10.1186/s13059-023-02978-x]
• 58. Evaluation of different approaches for missing data imputation ... [https://biodatamining.biomedcentral.com/articles/10.1186/s13040-021-00274-7]
• 59. Sparsely correlated hidden Markov models with application ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3582268/]
• 60. Sparse measures with swarm-based pliable hidden ... [https://www.frontiersin.org/journals/computational-neuroscience/articles/10.3389/fncom.2022.1016516/full]
• 61. HMM — TRANSIT v3.3.13 documentation [https://transit.readthedocs.io/en/latest/method_HMM.html]
• 62. Denoising of Aligned Genomic Data | Scientific Reports [https://www.nature.com/articles/s41598-019-51418-z]
• 63. A new algorithm to train hidden Markov models for biological ... [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-04080-0]
• 64. A new algorithm to train hidden Markov models for ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC7995745/]
• 65. Computer Simulation | springerprofessional.de Validation [https://www.springerprofessional.de/computer-simulation-validation/16644808]
• 66. A two-stage hidden Markov model design for biomarker ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6116560/]
• 67. Biopython Tutorial and Cookbook [https://biopython.org/DIST/docs/tutorial/Tutorial-1.83.html]
• 68. Graphics including GenomeDiagram [https://biopython.org/docs/dev/Tutorial/chapter_graphics.html]
• 69. Detecting adaptive introgression in human evolution using ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8192126/]
• 70. Best practices and pitfalls in using hmmix for reference-free ... [https://www.pivotscipub.com/hpgg/5/4/0007/html]
• 71. IntroUNET: Identifying introgressed alleles via semantic ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC10906877/]
• 72. Determining the probability of hemiplasy in the presence ... [https://elifesciences.org/articles/63753]
• 73. Importance of incomplete lineage sorting and introgression ... [https://www.nature.com/articles/hdy201672]
• 74. Powerful methods for detecting introgressed regions from ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC4899106/]
• 75. Day #10 — ILS x Introgression Analysis [1/3] [https://medium.com/@learningthoughts/day-10-ils-x-introgression-analysis-1-3-a4ce3725b43]