Introduction to HAPTRACE

Shweta Sahoo, Sara de las Heras-Saldana, Julius H. J. van der Werf, Mohammad H. Ferdosi

2026-07-22

Introduction

In livestock breeding programs, genetically superior animals from multiple breeds are often selected and mated to enhance genetic merit and genetic diversity, resulting in an admixed population. The genomic prediction of an animal’s genetic merit requires the use of pedigree and genotype information. Tracking breed-specific genomic segments from different populations is essential to characterise population admixture and assess the feasibility of the animal breeding program.

The HAPTRACE package allows haplotype-based tracking in a simulated admixed population for the accurate estimation of breed composition (Sahoo et al., 2025). This package provides a wide range of applications, offering the capability to simulate complex animal breeding programs. Output files are provided in the format required for PLINK (Purcell et al., 2007) and SNeP software (Barbato et al., 2015), facilitating population diversity analysis. Additionally, this package allows the visualization of haplotypes, QTL effects, recombination and admixture to better describe the population.

This vignette provides a small example of how to use HAPTRACE for simulating pure and admixed populations, as well as visualize the outputs.

library(HAPTRACE)
set.seed(1)

Input data for simulating the population

HAPTRACE supports three simulation modes, each differing in the source of the base population: 1) an internally generated base population; 2) a base population from QMSim outputs; or 3) use real genotype data as the starting point. However, to obtain the haplotype information, the user should ensure that the genotypes are phased. This vignette describes the internal simulation of a base population and its subsequent use in generating an admixed population.

First, a MAP and QTL effects should be simulated using the functions generateMAP and QTLeffects, respectively. The MAP file includes information on marker ID, chromosome number and position, while the QTL effects file describes QTLID, chromosome number, QTL position and QTL effects. The mapQTLfile function must be used to combine QTL effects and markers to obtain the total number of markers. In real data analysis, SNP markers are often observed to be around the QTL. However, for simulation purposes in this package, the QTL position is included along with the map file to get a combined marker file. The output from the combined file will be used as input to generate a historical population.

Creating MAP and QTL file

MAPfile <- generateMAP(nChr = 10, n_markers = 1000, len_Chr = 100)
QTLfile <- QTLeffects(n_QTL_Chr = 10, nChr = 10, len_Chr = 100, effect_distribution = "gamma")

MAP_QTL_Pop <- mapQTLfile(map = MAPfile, QTL = QTLfile)
Effect <- as.vector(as.numeric(MAP_QTL_Pop$Effect))

Creating a historical population

After generating a combined MAP and QTL file, the function generateHP will generate haplotype information for the historical population across multiple generations, as described below. The historical population provides haplotype information for the base population. Users can define the number of sires and dams to mate and simulate the next generation of progenies, as well as the selection types (e.g. random, based on phenotypes or true breeding values). In this case, we opted for random selection for demonstration.

Historical_Population <- generateHP(n_generations = 5, n_animals = 2000, 
                                   map = MAP_QTL_Pop, nSire = 200, nDam = 500, 
                                   nProgeny = 1000, mutationRate = 2.5 * 10^-5,
                                   SelType = "random", Effect = Effect, nChr = 10, 
                                   h2 = 0.3, trait_mean = 40, VarE = 0.6,recL = 10, 
                                   minLength = 200)
#> Working on Generation: 1
#> Working on Generation: 2
#> Working on Generation: 3
#> Working on Generation: 4
#> Working on Generation: 5

Sim_HP <- Historical_Population$Haplotype
str(Sim_HP)
#>  num [1:10000, 1:10100] 0 1 0 1 0 0 0 0 0 0 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : chr [1:10000] "H_1" "H_1" "H_2" "H_2" ...
#>   ..$ : NULL
class(Sim_HP)
#> [1] "matrix" "array"
Sim_HP[1:10, 1:10]
#>     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> H_1    0    0    0    1    1    0    0    0    0     0
#> H_1    1    1    0    0    1    0    1    1    0     0
#> H_2    0    1    1    0    0    0    0    0    0     0
#> H_2    1    1    0    0    1    0    1    1    0     0
#> H_3    0    1    1    0    1    0    0    1    0     1
#> H_3    0    1    1    1    1    0    1    1    0     1
#> H_4    0    1    1    0    0    0    0    0    0     0
#> H_4    0    1    1    1    1    0    1    1    0     1
#> H_5    0    1    1    0    1    0    0    1    0     1
#> H_5    0    0    1    1    0    0    0    0    1     0

Next, the historical population is divided into populations A and B for further simulation.

After that, we can rescale the QTL effect according to the genomic information from the animals and the phenotypic standard deviation to obtain the desirable variance components for the simulated trait.

Population_A <- Sim_HP[1:(nrow(Sim_HP)/2),]
Population_B <- Sim_HP[(nrow(Sim_HP)/2 + 1):nrow(Sim_HP), ]

EffectA <- Rescale(Haplotype = Population_A, Effect = Effect, Phenosd = 3, h2 = 0.3)

Assigning of population code and haplotype recording

To track the pedigree, we need to assign distinct codes to populations A and B, which will be reflected in the animals’ ID. Similarly, the haplotypes of these two populations should be recoded to facilitate tracking them. In the following example, the resulting haplotype information was recoded from (0,1) to (-1,1) and (-2,2) for populations A and B, respectively.

# Assignment of population code
Population_A <- PopulationID(Population = Population_A, PopCode = "P1")
Population_B <- PopulationID(Population = Population_B, PopCode = "P2")

# Recode haplotype of populations A and B with 1 and 2, respectively
Population_A <- Recode_Haplotype(df = Population_A, popBase = 1)
Population_B <- Recode_Haplotype(df = Population_B, popBase = 2)

Modelling the breeding program

Populations A and B are simulated for 3 generations with different numbers of sires and dams assigned for random mating using GenOnePopulation. This function provides information on the pedigree, phenotypes, haplotypes coded according to the breed of origin, the original haplotype information, true breeding value and recombination points from sires and dams. The parameter recL defines the average number of recombinations per animal, which can be specified in the function. In addition to random selection, users can apply selection based on either phenotype or true breeding values.

# Population A
popA <- GenOnePopulation(ngenerations = 1, map = MAP_QTL_Pop,
                         populationSim = Population_A, nSire = 20, nDam = 100,
                         mutationRate = 2.5 * 10^-5, SelType = "random", 
                         h2 = 0.3, trait_mean = 10, VarE = 0.6,Effect = Effect, 
                         recL = 10, nChr = 10, minLength = 200, IndStartVal=1,
                         prefixID = "A", nProgeny = 300)
#> Working on Generation: 1

# Population B
popB <- GenOnePopulation(ngenerations = 1, map = MAP_QTL_Pop,
                         populationSim = Population_B, nSire = 10, nDam = 100,
                         mutationRate = 2.5 * 10^-5, SelType = "random", 
                         h2 = 0.3, trait_mean = 30, VarE = 0.6 ,Effect = Effect,
                         recL = 10,nChr = 10, minLength = 200, IndStartVal=1,
                         prefixID = "B", nProgeny = 300)
#> Working on Generation: 1

# Simulation of a crossbred population derived from populations A and B

The animals from the last generation in populations A and B are crossed using the function TwoPopCross to get a crossbred population named F1Cross. The animals in this population were inter crossed further for the next 2 generations using GenOnePopulation function.

# First cross between populations A and B
F1Cross <- TwoPopCross(map = MAP_QTL_Pop, populationSim_A = popA$Haplotype,
                       populationSim_B = popB$Haplotype, nSireA = 20, nDamA = 100,
                       nSireB = 10, nDamB = 100, Sire_From = "A",
                       mutationRate = 2.5 * 10^-5, SelType = "random", h2 = 0.3,
                       trait_mean = 20,VarE_PopA = 0.6, VarE_PopB = 0.6, 
                       Effect_PopA = Effect, Effect_PopB = Effect,
                       IndStartVal = 1, prefixID = "F1C", nProgeny = 300, recL = 10, 
                       nChr = 10, minLength = 200)
#> Working on two populations cross
# Simulation of an admixed population
F3 <- GenOnePopulation(ngenerations = 2, map = MAP_QTL_Pop,
                       populationSim = F1Cross$Haplotype,
                       nSire = 20,nDam = 100,mutationRate = 2.5 * 10^-5,
                       SelType = "random",h2 = 0.3, trait_mean = 20,
                       VarE = 0.6, Effect = Effect, recL = 10, 
                       nChr = 10, minLength = 200, IndStartVal = 1,
                       prefixID = "C2", nProgeny = 300)
#> Working on Generation: 1
#> Working on Generation: 2

Phenotype data can be extracted using Ext_Pheno, while Ped_Pheno provides both pedigree and phenotype information.

# Extracting pedigree and phenotype
F3Pheno = Ext_Pheno(df = F3$population)
head(F3Pheno)
#>       AnimID      TBV phenotype
#> F1C_1  F1C_1 11.07333  30.22665
#> F1C_2  F1C_2 15.21855  35.35315
#> F1C_3  F1C_3 14.51038  34.20899
#> F1C_4  F1C_4 15.45068  35.61448
#> F1C_5  F1C_5 14.69087  34.19344
#> F1C_6  F1C_6 19.73718  38.28955
F3Ped_Pheno = Ped_Pheno(Pedigree = F3$parents, Phenotype = F3Pheno)
head(F3Ped_Pheno)
#>     AnimID SireID   DamID Sex       TBV phenotype
#> 1     C2_1 F1C_60 F1C_299   1  8.056027  27.68618
#> 112   C2_2 F1C_60 F1C_299   1 12.029614  31.48033
#> 223   C2_3 F1C_60 F1C_299   1  8.247797  27.64909
#> 334   C2_4 F1C_60 F1C_299   1 15.172398  35.28659
#> 445   C2_5 F1C_60 F1C_223   1 14.923840  34.66154
#> 556   C2_6 F1C_60 F1C_223   1 10.083562  28.35290

The function Coded_to_Haplo converts the coded haplotype into (0,1) format and MakeGeno converts haplotypes to genotype information in (0,1,2) format.

# Extracting genotype information for F1Cross and F3 populations
CBgenotype <- rbind(MakeGeno(Coded_to_Haplo(F1Cross$Haplotype)),
                    MakeGeno(Coded_to_Haplo(F3$Haplotype)))
CBgenotype[1:10, 1:10]
#>        [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> F1C_1     0    1    0    1    1    1    1    1    0     0
#> F1C_2     0    1    1    1    1    1    2    0    0     0
#> F1C_3     1    0    1    1    1    0    1    1    1     2
#> F1C_4     2    0    1    0    2    0    1    1    1     1
#> F1C_5     1    0    0    0    1    1    1    0    2     2
#> F1C_6     2    0    0    1    1    1    1    0    2     1
#> F1C_7     0    1    1    0    1    2    2    0    1     2
#> F1C_8     1    0    0    0    1    1    1    0    2     2
#> F1C_9     0    1    1    0    1    2    2    0    1     2
#> F1C_10    2    1    0    1    1    1    1    1    2     1
# Coded haplotype for F1Cross and F3 populations
CBHaplotype <- rbind(F1Cross$Haplotype,F3$Haplotype)
CBHaplotype[1:10, 1:10]
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> F1C_1   -1   -1   -1    1    1    1    1   -1   -1    -1
#> F1C_1   -2    2   -2   -2   -2   -2   -2    2   -2    -2
#> F1C_2   -1   -1   -1    1    1    1    1   -1   -1    -1
#> F1C_2   -2    2    2   -2   -2   -2    2   -2   -2    -2
#> F1C_3   -1   -1    1    1   -1   -1    1   -1   -1     1
#> F1C_3    2   -2   -2   -2    2   -2   -2    2    2     2
#> F1C_4    1   -1   -1   -1    1   -1   -1   -1    1     1
#> F1C_4    2   -2    2   -2    2   -2    2    2   -2    -2
#> F1C_5    1   -1   -1   -1    1   -1   -1   -1    1     1
#> F1C_5   -2   -2   -2   -2   -2    2    2   -2    2     2

Visualization of an admixed population

The coded haplotype information recorded in the F3 admixed population helps to determine the accurate breed composition. Various plot functions are included in this package to provide an overview of how genomic segments are inherited from different populations. The plotHaplo function offers a snippet of haplotype segments inherited from different populations, whereas Admixplot depicts the accurate breed composition of an admixed population, where x-axis shows the number of animals and y-axis shows the proportion from different populations in different colors.

# Haplotype plot
smallF3 <- F3$Haplotype[1:100, ]
F3pAplot <- plotHaplo(Haplotype = smallF3, 
                     colors = c("#CC79A7", "#0072B2"), 
                     map = MAP_QTL_Pop, nChr = 10, 
                     legendlabels = c("A", "B"))
Figure 1: Haplotype plot for crossbred animals obtained from plotHaplo function
Figure 1: Haplotype plot for crossbred animals obtained from plotHaplo function
# Admixture plot
TBC <- trueBreedComp(Haplotype = smallF3)

# Arrange the population
TBC_data <- TBC[c(2,1), ]

Admixplot(trueBC = TBC_data, 
          color = c("#CC79A7", "#0072B2"),
          legendlabels = c("A", "B"))
Figure 2: Admixture plot for crossbred animals obtained from Admixplot function
Figure 2: Admixture plot for crossbred animals obtained from Admixplot function

References

Barbato, M., Orozco-terWengel, P., Tapio, M., & Bruford, M. W. (2015). SNeP: A tool to estimate trends in recent effective population size trajectories using genome-wide SNP data. Frontiers in Genetics, 6. https://doi.org/10.3389/fgene.2015.00109
Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., Maller, J., Sklar, P., De Bakker, P. I. W., Daly, M. J., & Sham, P. C. (2007). PLINK: A tool set for whole-genome association and population-based linkage analyses. The American Journal of Human Genetics, 81(3), 559–575. https://doi.org/10.1086/519795
Sahoo, S., Ferdosi, M., Werf, J. van der, & De las Heras-Saldana, S. (2025). Ascertaining the breed composition of an admixed sheep population using genomic information: A simulation study. Proc. Assoc. Advmt. Anim. Breed. Genet, 26, 323–326.