Phenotype-guided haplotype refinement

Introduction

library(GeneTrackR)

refine_haplotype() performs phenotype-guided grouping of existing haplotypes. It does not discover new variants and it does not rebuild the original genotype patterns. Instead, it starts from an existing HapVariant, compares the phenotype distributions of the original haplotypes, and groups haplotypes that are sufficiently similar for the selected trait or traits.

This module uses the same GeneA haplotypes created in the haplotype workflow and the same phenotype table used in the phenotype workflow. The demo is deliberately structured so that the refinement behavior can be checked against known expectations:

  • protein_content: the genotype-nearest pair Hap1/Hap2 shares the low phenotype class and the genotype-nearest pair Hap3/Hap4 shares the high phenotype class, providing the main 4 -> 2 refinement example;
  • seed_weight: all four haplotypes retain distinct means, but genotype-nearest haplotypes remain phenotypically closer than cross-cluster haplotypes;
  • flowering_time: no GeneA haplotype effect was designed, so the four haplotypes can collapse into one phenotype group under the default significance criterion.

The workflow follows nine steps:

  1. prepare the original haplotypes and phenotype data;
  2. inspect genotype similarity and phenotype correspondence before refinement;
  3. refine GeneA haplotypes using protein_content;
  4. inspect the HapRefined result object;
  5. inspect pairwise tests and the original-to-refined mapping;
  6. plot the grouped and collapsed refined haplotype structures;
  7. plot phenotype distributions for refined groups;
  8. evaluate trait choice, multiple traits, and refinement thresholds;
  9. interpret refined groups and pass them downstream.

Package documentation: https://renscq.github.io/GeneTrackR/.

Installation

After Bioconductor acceptance, install GeneTrackR with BiocManager:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::install("GeneTrackR")
library(GeneTrackR)

The development version remains available from the GeneTrackR GitHub repository.

Step 1. Prepare the original haplotypes and phenotype data

Locate and read the same demo files used in the previous workflows:

gp_file <- system.file(
  "extdata",
  "gtr_demo.genePredExt",
  package = "GeneTrackR"
)
vcf_file <- system.file(
  "extdata",
  "gtr_demo_variants.vcf",
  package = "GeneTrackR"
)
pheno_file <- system.file(
  "extdata",
  "gtr_demo_pheno.tsv",
  package = "GeneTrackR"
)
gp <- read_genepred(
  gp_file,
  format = "genePredExt",
  verbose = FALSE,
  progress = FALSE
)
vcf <- read_vcf(
  vcf_file,
  mode = "memory",
  keep_genotype = TRUE,
  verbose = FALSE,
  progress = FALSE
)
pheno <- read_pheno(
  pheno_file,
  verbose = FALSE,
  progress = FALSE
)

Rebuild the standard GeneA gene-body haplotypes using the same definition as in the haplotype article:

hap_gene <- hap_gene_variant(
  vcf,
  annotation = gp,
  gene_id = "GeneA",
  genotype_mode = "string"
)
#> [GeneTrackR] Retrieved variants: 11.
hap_gene
#> <HapVariant>
#>   region    : chr1:12340001-12352000
#>   variants  : 11
#>   samples   : 36
#>   haplotypes: 4
hap_gene$haplotypes
#> Key: <hap_id>
#>    hap_id sample_n                             samples varA01 varA02 varA03
#>    <char>    <int>                              <char> <char> <char> <char>
#> 1:   Hap1        9 S10;S11;S12;S13;S14;S15;S16;S17;S18      A      C      C
#> 2:   Hap2        9 S01;S02;S03;S04;S05;S06;S07;S08;S09      A      T      C
#> 3:   Hap3        9 S19;S20;S21;S22;S23;S24;S25;S26;S27      G      C      G
#> 4:   Hap4        9 S28;S29;S30;S31;S32;S33;S34;S35;S36      G      T      G
#>    varLD01 varLD02 varLD03 varLD04 varLD05 varLD06 varA04 varA05
#>     <char>  <char>  <char>  <char>  <char>  <char> <char> <char>
#> 1:       G       A       C       G       T       A      A     i3
#> 2:       G       A       C       G       T       A      A     i3
#> 3:       A       T       T       C       G       C     i3     i3
#> 4:       A       T       T       C       G       C     i3      A

For the bundled demo, the expected starting point is:

  • 11 GeneA variants;
  • 36 samples;
  • 4 original haplotypes;
  • 9 samples per haplotype.

Refinement is therefore performed on a fixed biological definition of the original haplotypes. Do not silently change the gene region, variant filter, missing-data rule, or sample set immediately before refinement, because those changes would redefine the groups being compared.

Step 2. Inspect genotype similarity and phenotype correspondence before refinement

Refinement is easiest to interpret when the phenotype grouping is compatible with the underlying genotype geometry. First calculate the number of allele-state differences between the four original GeneA haplotypes:

hap_variant_ids <- as.character(hap_gene$variants$variant_id)
hap_alleles <- as.matrix(
  hap_gene$haplotypes[, ..hap_variant_ids]
)
hap_ids <- as.character(hap_gene$haplotypes$hap_id)
hap_distance <- outer(
  seq_along(hap_ids),
  seq_along(hap_ids),
  Vectorize(function(i, j) {
    sum(hap_alleles[i, ] != hap_alleles[j, ])
  })
)
dimnames(hap_distance) <- list(hap_ids, hap_ids)
hap_distance
#>      Hap1 Hap2 Hap3 Hap4
#> Hap1    0    1    9   11
#> Hap2    1    0   10   10
#> Hap3    9   10    0    2
#> Hap4   11   10    2    0

The deterministic GeneA design is hierarchical:

          Hap1  Hap2  Hap3  Hap4
Hap1         0     1     9    11
Hap2         1     0    10    10
Hap3         9    10     0     2
Hap4        11    10     2     0

Therefore:

  • Hap1 and Hap2 are the closest genotype pair and differ at only one of the 11 GeneA variants;
  • Hap3 and Hap4 are the second closest pair and differ at only two variants;
  • comparisons between these two clusters differ at 9-11 variants.

The phenotype design follows the same structure. Inspect protein_content before refinement:

protein_before <- plot_hap_pheno(
  hap_gene,
  phenotype = pheno,
  traits = "protein_content",
  min_hap_samples = 3,
  test_method = "t.test",
  p_adjust = "BH",
  p_value_type = "adjusted",
  show_points = TRUE,
  fill_palette = "Paired",
  text_size = 9
)
protein_before$figure

protein_before$summary
#>              trait hap_id sample_n median_value mean_value min_value max_value
#>             <char> <fctr>    <int>        <num>      <num>     <num>     <num>
#> 1: protein_content   Hap2        9           38         38     37.58     38.42
#> 2: protein_content   Hap1        9           38         38     37.58     38.42
#> 3: protein_content   Hap3        9           44         44     43.58     44.42
#> 4: protein_content   Hap4        9           44         44     43.58     44.42
protein_before$pvalue
#>              trait group1 group2      p_value        p_adj method
#>             <char> <char> <char>        <num>        <num> <char>
#> 1: protein_content   Hap2   Hap1 1.000000e+00 1.000000e+00 t.test
#> 2: protein_content   Hap2   Hap3 2.707879e-18 4.061818e-18 t.test
#> 3: protein_content   Hap2   Hap4 2.707879e-18 4.061818e-18 t.test
#> 4: protein_content   Hap1   Hap3 2.707879e-18 4.061818e-18 t.test
#> 5: protein_content   Hap1   Hap4 2.707879e-18 4.061818e-18 t.test
#> 6: protein_content   Hap3   Hap4 1.000000e+00 1.000000e+00 t.test

The expected phenotype means are:

Hap1 ≈ 38
Hap2 ≈ 38
Hap3 ≈ 44
Hap4 ≈ 44

This correspondence is intentional:

genotype-nearest pair        phenotype class
Hap1 + Hap2                  low protein_content
Hap3 + Hap4                  high protein_content

varA03 follows the same cluster split: Hap1/Hap2 carry the REF class and Hap3/Hap4 carry the ALT class. The main refinement example therefore groups haplotypes that are similar in both multi-variant genotype structure and phenotype, rather than joining genetically distant haplotypes only because their phenotype means happen to match.

Step 3. Refine haplotypes using protein_content

Run phenotype-guided refinement:

refined_protein <- refine_haplotype(
  hap_gene,
  phenotype = pheno,
  traits = "protein_content",
  min_hap_samples = 3,
  test_method = "t.test",
  p_adjust = "BH",
  alpha = 0.05,
  effect_threshold = 0.5
)
refined_protein
#> <HapRefined>
#>   original haplotypes: 4
#>   refined haplotypes : 2
#>   traits             : protein_content
refined_protein$refined_haplotypes
#>    hap_id original_hap_n original_hap_ids sample_n
#>    <char>          <int>           <char>    <int>
#> 1:  RHap1              2        Hap1;Hap2       18
#> 2:  RHap2              2        Hap3;Hap4       18
#>                                                                    samples
#>                                                                     <char>
#> 1: S10;S11;S12;S13;S14;S15;S16;S17;S18;S01;S02;S03;S04;S05;S06;S07;S08;S09
#> 2: S19;S20;S21;S22;S23;S24;S25;S26;S27;S28;S29;S30;S31;S32;S33;S34;S35;S36
#>    varA01 varA02 varA03 varLD01 varLD02 varLD03 varLD04 varLD05 varLD06 varA04
#>    <char> <char> <char>  <char>  <char>  <char>  <char>  <char>  <char> <char>
#> 1:      A  mixed      C       G       A       C       G       T       A      A
#> 2:      G  mixed      G       A       T       T       C       G       C     i3
#>    varA05
#>    <char>
#> 1:     i3
#> 2:  mixed

The expected demo result is 4 original haplotypes -> 2 refined haplotypes, with 18 samples in each refined group:

Hap1 + Hap2  -> one refined group
Hap3 + Hap4  -> one refined group

This is the key truth condition for the example: refinement follows the two nearest genotype clusters rather than crossing between them.

The two refinement criteria have different roles:

  • alpha = 0.05: an original haplotype pair can be merged only when its adjusted pairwise P value is greater than the cutoff;
  • effect_threshold = 0.5: the absolute difference between the two haplotype phenotype means must also be no larger than 0.5 protein-content units.

Using an effect threshold is recommended when the biological interpretation depends on similarity. A non-significant P value alone means that a difference was not detected; it is not proof that the two groups are biologically equivalent.

Step 4. Inspect the HapRefined object

refine_haplotype() returns a HapRefined object rather than only a plotting result:

names(refined_protein)
#>  [1] "original_hap"              "refined_hap"              
#>  [3] "refined_haplotypes"        "sample_refined_haplotypes"
#>  [5] "haplotype_map"             "trait_group_map"          
#>  [7] "phenotype_summary"         "pairwise_test"            
#>  [9] "plot_data"                 "parameters"

The main components are:

  • original_hap: the original HapVariant supplied to refine_haplotype();
  • refined_hap: a refined HapVariant compatible with the standard haplotype plotting/phenotype functions;
  • refined_haplotypes: one row per refined haplotype group;
  • sample_refined_haplotypes: sample-level refined group assignments while retaining the original haplotype ID;
  • haplotype_map: mapping from each original haplotype to its refined haplotype;
  • trait_group_map: trait-specific phenotype grouping before multiple traits are combined;
  • phenotype_summary: phenotype summary statistics for the original haplotypes;
  • pairwise_test: pairwise phenotype tests and merge decisions;
  • plot_data: matched sample-level phenotype data used during refinement;
  • parameters: the refinement settings used to generate the object.

Inspect the original-to-refined mapping directly:

refined_protein$haplotype_map
#>    hap_id refined_hap_id refined_key protein_content
#>    <char>         <char>      <char>          <char>
#> 1:   Hap1          RHap1          G1              G1
#> 2:   Hap2          RHap1          G1              G1
#> 3:   Hap3          RHap2          G2              G2
#> 4:   Hap4          RHap2          G2              G2
refined_protein$sample_refined_haplotypes
#>     sample_id hap_id original_hap_id             hap_pattern
#>        <char> <char>          <char>                  <char>
#>  1:       S10  RHap1            Hap1  A|C|C|G|A|C|G|T|A|A|i3
#>  2:       S11  RHap1            Hap1  A|C|C|G|A|C|G|T|A|A|i3
#>  3:       S12  RHap1            Hap1  A|C|C|G|A|C|G|T|A|A|i3
#>  4:       S13  RHap1            Hap1  A|C|C|G|A|C|G|T|A|A|i3
#>  5:       S14  RHap1            Hap1  A|C|C|G|A|C|G|T|A|A|i3
#>  6:       S15  RHap1            Hap1  A|C|C|G|A|C|G|T|A|A|i3
#>  7:       S16  RHap1            Hap1  A|C|C|G|A|C|G|T|A|A|i3
#>  8:       S17  RHap1            Hap1  A|C|C|G|A|C|G|T|A|A|i3
#>  9:       S18  RHap1            Hap1  A|C|C|G|A|C|G|T|A|A|i3
#> 10:       S01  RHap1            Hap2  A|T|C|G|A|C|G|T|A|A|i3
#> 11:       S02  RHap1            Hap2  A|T|C|G|A|C|G|T|A|A|i3
#> 12:       S03  RHap1            Hap2  A|T|C|G|A|C|G|T|A|A|i3
#> 13:       S04  RHap1            Hap2  A|T|C|G|A|C|G|T|A|A|i3
#> 14:       S05  RHap1            Hap2  A|T|C|G|A|C|G|T|A|A|i3
#> 15:       S06  RHap1            Hap2  A|T|C|G|A|C|G|T|A|A|i3
#> 16:       S07  RHap1            Hap2  A|T|C|G|A|C|G|T|A|A|i3
#> 17:       S08  RHap1            Hap2  A|T|C|G|A|C|G|T|A|A|i3
#> 18:       S09  RHap1            Hap2  A|T|C|G|A|C|G|T|A|A|i3
#> 19:       S19  RHap2            Hap3 G|C|G|A|T|T|C|G|C|i3|i3
#> 20:       S20  RHap2            Hap3 G|C|G|A|T|T|C|G|C|i3|i3
#> 21:       S21  RHap2            Hap3 G|C|G|A|T|T|C|G|C|i3|i3
#> 22:       S22  RHap2            Hap3 G|C|G|A|T|T|C|G|C|i3|i3
#> 23:       S23  RHap2            Hap3 G|C|G|A|T|T|C|G|C|i3|i3
#> 24:       S24  RHap2            Hap3 G|C|G|A|T|T|C|G|C|i3|i3
#> 25:       S25  RHap2            Hap3 G|C|G|A|T|T|C|G|C|i3|i3
#> 26:       S26  RHap2            Hap3 G|C|G|A|T|T|C|G|C|i3|i3
#> 27:       S27  RHap2            Hap3 G|C|G|A|T|T|C|G|C|i3|i3
#> 28:       S28  RHap2            Hap4  G|T|G|A|T|T|C|G|C|i3|A
#> 29:       S29  RHap2            Hap4  G|T|G|A|T|T|C|G|C|i3|A
#> 30:       S30  RHap2            Hap4  G|T|G|A|T|T|C|G|C|i3|A
#> 31:       S31  RHap2            Hap4  G|T|G|A|T|T|C|G|C|i3|A
#> 32:       S32  RHap2            Hap4  G|T|G|A|T|T|C|G|C|i3|A
#> 33:       S33  RHap2            Hap4  G|T|G|A|T|T|C|G|C|i3|A
#> 34:       S34  RHap2            Hap4  G|T|G|A|T|T|C|G|C|i3|A
#> 35:       S35  RHap2            Hap4  G|T|G|A|T|T|C|G|C|i3|A
#> 36:       S36  RHap2            Hap4  G|T|G|A|T|T|C|G|C|i3|A
#>     sample_id hap_id original_hap_id             hap_pattern
#>        <char> <char>          <char>                  <char>
#>     non_missing_variant_n
#>                     <int>
#>  1:                    11
#>  2:                    11
#>  3:                    11
#>  4:                    11
#>  5:                    11
#>  6:                    11
#>  7:                    11
#>  8:                    11
#>  9:                    11
#> 10:                    11
#> 11:                    11
#> 12:                    11
#> 13:                    11
#> 14:                    11
#> 15:                    11
#> 16:                    11
#> 17:                    11
#> 18:                    11
#> 19:                    11
#> 20:                    11
#> 21:                    11
#> 22:                    11
#> 23:                    11
#> 24:                    11
#> 25:                    11
#> 26:                    11
#> 27:                    11
#> 28:                    11
#> 29:                    11
#> 30:                    11
#> 31:                    11
#> 32:                    11
#> 33:                    11
#> 34:                    11
#> 35:                    11
#> 36:                    11
#>     non_missing_variant_n
#>                     <int>

refined_protein$refined_hap is intentionally a normal HapVariant-compatible object:

refined_protein$refined_hap
#> <HapVariant>
#>   region    : chr1:12340001-12352000
#>   variants  : 11
#>   samples   : 36
#>   haplotypes: 2
refined_protein$refined_hap$haplotypes
#>    hap_id original_hap_n original_hap_ids sample_n
#>    <char>          <int>           <char>    <int>
#> 1:  RHap1              2        Hap1;Hap2       18
#> 2:  RHap2              2        Hap3;Hap4       18
#>                                                                    samples
#>                                                                     <char>
#> 1: S10;S11;S12;S13;S14;S15;S16;S17;S18;S01;S02;S03;S04;S05;S06;S07;S08;S09
#> 2: S19;S20;S21;S22;S23;S24;S25;S26;S27;S28;S29;S30;S31;S32;S33;S34;S35;S36
#>    varA01 varA02 varA03 varLD01 varLD02 varLD03 varLD04 varLD05 varLD06 varA04
#>    <char> <char> <char>  <char>  <char>  <char>  <char>  <char>  <char> <char>
#> 1:      A  mixed      C       G       A       C       G       T       A      A
#> 2:      G  mixed      G       A       T       T       C       G       C     i3
#>    varA05
#>    <char>
#> 1:     i3
#> 2:  mixed

This makes it possible to reuse the standard haplotype plotting and phenotype-analysis infrastructure without rebuilding variants or genotypes.

Step 5. Inspect pairwise tests and merge decisions

Do not treat refinement as a black box. Inspect the pairwise evidence used to generate the refined groups:

refined_protein$pairwise_test[, c(
  "trait",
  "group1",
  "group2",
  "mean_group1",
  "mean_group2",
  "abs_mean_diff",
  "p_value",
  "p_adj",
  "can_merge"
)]
#>              trait group1 group2 mean_group1 mean_group2 abs_mean_diff
#>             <char> <char> <char>       <num>       <num>         <num>
#> 1: protein_content   Hap2   Hap1          38          38             0
#> 2: protein_content   Hap2   Hap3          38          44             6
#> 3: protein_content   Hap2   Hap4          38          44             6
#> 4: protein_content   Hap1   Hap3          38          44             6
#> 5: protein_content   Hap1   Hap4          38          44             6
#> 6: protein_content   Hap3   Hap4          44          44             0
#>         p_value        p_adj can_merge
#>           <num>        <num>    <lgcl>
#> 1: 1.000000e+00 1.000000e+00      TRUE
#> 2: 2.707879e-18 4.061818e-18     FALSE
#> 3: 2.707879e-18 4.061818e-18     FALSE
#> 4: 2.707879e-18 4.061818e-18     FALSE
#> 5: 2.707879e-18 4.061818e-18     FALSE
#> 6: 1.000000e+00 1.000000e+00      TRUE

For the protein_content demo, the two mergeable comparisons are expected to be Hap1 versus Hap2 and Hap3 versus Hap4. Both pairs have:

  • nearly identical means;
  • adjusted P values greater than 0.05;
  • abs_mean_diff <= 0.5;
  • can_merge = TRUE.

All comparisons between the Hap1/Hap2 and Hap3/Hap4 genotype clusters have a large phenotype difference and should not be merged.

The trait-specific grouping can be inspected separately:

refined_protein$trait_group_map
#>              trait hap_id trait_group mean_value median_value sample_n
#>             <char> <char>      <char>      <num>        <num>    <int>
#> 1: protein_content   Hap1          G1         38           38        9
#> 2: protein_content   Hap2          G1         38           38        9
#> 3: protein_content   Hap3          G2         44           44        9
#> 4: protein_content   Hap4          G2         44           44        9

The final original-to-refined relationship is stored in:

refined_protein$haplotype_map
#>    hap_id refined_hap_id refined_key protein_content
#>    <char>         <char>      <char>          <char>
#> 1:   Hap1          RHap1          G1              G1
#> 2:   Hap2          RHap1          G1              G1
#> 3:   Hap3          RHap2          G2              G2
#> 4:   Hap4          RHap2          G2              G2

This separation is useful when multiple traits are used, because each trait is grouped independently before GeneTrackR combines the trait-specific group labels into the final refined haplotypes.

Step 6. Plot the refined haplotype structure

plot_refined_hap_variant() provides two complementary views.

Keep the original haplotypes inside refined blocks

The default collapse_refined = FALSE keeps the original haplotype rows and arranges them into refined groups:

refined_variant_grouped <- plot_refined_hap_variant(
  refined_protein,
  annotation = gp,
  min_hap_samples = 3,
  collapse_refined = FALSE,
  show_reference_row = TRUE,
  variant_label = "pos",
  gene_pos_x_angle = 90,
  gene_track_legend_position = "top",
  gene_palette = "Paired",
  table_palette = "Paired",
  variant_palette = "Paired",
  text_size = 9
)
refined_variant_grouped

This is usually the most informative refinement plot because it shows which original genotype patterns were grouped together while retaining the original variant states.

Collapse each refined group to one row

Set collapse_refined = TRUE when a compact summary is required:

refined_variant_collapsed <- plot_refined_hap_variant(
  refined_protein,
  annotation = gp,
  min_hap_samples = 3,
  collapse_refined = TRUE,
  show_reference_row = TRUE,
  variant_label = "pos",
  gene_palette = "Paired",
  table_palette = "Paired",
  variant_palette = "Paired",
  text_size = 9
)
refined_variant_collapsed

When original haplotypes within one refined group carry different alleles at a variant, the collapsed refined row uses the configured mixed_label (default "mixed"). A mixed state does not mean that an individual sample is heterozygous. It means that the refined phenotype group contains more than one original haplotype state at that variant.

For mechanistic interpretation, keep the grouped view available even when the collapsed view is used in a summary figure.

Step 7. Plot phenotype distributions after refinement

Use plot_refined_hap_pheno() to visualize the phenotype distributions of the refined groups:

refined_protein_pheno <- plot_refined_hap_pheno(
  refined_protein,
  phenotype = pheno,
  traits = "protein_content",
  min_hap_samples = 3,
  test_method = "t.test",
  p_adjust = "BH",
  p_value_type = "adjusted",
  show_points = TRUE,
  fill_palette = "Paired",
  text_size = 9
)
refined_protein_pheno$figure

refined_protein_pheno$summary
#>              trait hap_id sample_n median_value mean_value min_value max_value
#>             <char> <fctr>    <int>        <num>      <num>     <num>     <num>
#> 1: protein_content  RHap1       18           38         38     37.58     38.42
#> 2: protein_content  RHap2       18           44         44     43.58     44.42
refined_protein_pheno$pvalue
#>              trait group1 group2      p_value        p_adj method
#>             <char> <char> <char>        <num>        <num> <char>
#> 1: protein_content  RHap1  RHap2 2.008652e-37 2.008652e-37 t.test

The return value is the same GeneTrackRPhenoPlot structure used by plot_hap_pheno():

  • $figure;
  • $pvalue;
  • $summary;
  • $bracket;
  • $plot_data.

This is useful for exporting the numerical evidence together with the refined-group figure rather than saving only the plot.

Step 8. Evaluate trait choice and refinement thresholds

Refinement is trait dependent. The same original haplotypes can produce very different refined structures when a different phenotype is selected.

A strongly haplotype-associated trait may preserve all groups

seed_weight was designed with four distinct means while preserving the same genotype hierarchy: Hap1/Hap2 remain closer to one another than to Hap3/Hap4, and the same is true within the second pair. The within-pair mean differences are nevertheless larger than effect_threshold = 0.5, so refinement should retain all four groups:

refined_seed_weight <- refine_haplotype(
  hap_gene,
  phenotype = pheno,
  traits = "seed_weight",
  min_hap_samples = 3,
  test_method = "t.test",
  p_adjust = "BH",
  alpha = 0.05,
  effect_threshold = 0.5
)
refined_seed_weight$refined_haplotypes
#>    hap_id original_hap_n original_hap_ids sample_n
#>    <char>          <int>           <char>    <int>
#> 1:  RHap1              1             Hap1        9
#> 2:  RHap2              1             Hap2        9
#> 3:  RHap3              1             Hap3        9
#> 4:  RHap4              1             Hap4        9
#>                                samples varA01 varA02 varA03 varLD01 varLD02
#>                                 <char> <char> <char> <char>  <char>  <char>
#> 1: S10;S11;S12;S13;S14;S15;S16;S17;S18      A      C      C       G       A
#> 2: S01;S02;S03;S04;S05;S06;S07;S08;S09      A      T      C       G       A
#> 3: S19;S20;S21;S22;S23;S24;S25;S26;S27      G      C      G       A       T
#> 4: S28;S29;S30;S31;S32;S33;S34;S35;S36      G      T      G       A       T
#>    varLD03 varLD04 varLD05 varLD06 varA04 varA05
#>     <char>  <char>  <char>  <char> <char> <char>
#> 1:       C       G       T       A      A     i3
#> 2:       C       G       T       A      A     i3
#> 3:       T       C       G       C     i3     i3
#> 4:       T       C       G       C     i3      A

The expected result is four refined groups, showing that refinement does not automatically force haplotypes to merge.

A negative-control trait can collapse all groups

The demo flowering_time trait contains no designed GeneA haplotype effect:

refined_flowering <- refine_haplotype(
  hap_gene,
  phenotype = pheno,
  traits = "flowering_time",
  min_hap_samples = 3,
  test_method = "t.test",
  p_adjust = "BH",
  alpha = 0.05,
  effect_threshold = 0.5
)
refined_flowering$refined_haplotypes
#>    hap_id original_hap_n    original_hap_ids sample_n
#>    <char>          <int>              <char>    <int>
#> 1:  RHap1              4 Hap1;Hap2;Hap3;Hap4       36
#>                                                                                                                                            samples
#>                                                                                                                                             <char>
#> 1: S10;S11;S12;S13;S14;S15;S16;S17;S18;S01;S02;S03;S04;S05;S06;S07;S08;S09;S19;S20;S21;S22;S23;S24;S25;S26;S27;S28;S29;S30;S31;S32;S33;S34;S35;S36
#>    varA01 varA02 varA03 varLD01 varLD02 varLD03 varLD04 varLD05 varLD06 varA04
#>    <char> <char> <char>  <char>  <char>  <char>  <char>  <char>  <char> <char>
#> 1:  mixed  mixed  mixed   mixed   mixed   mixed   mixed   mixed   mixed  mixed
#>    varA05
#>    <char>
#> 1:  mixed

For this deterministic demo, all four original haplotypes have the same mean flowering time and can collapse into one refined group.

This example is also a reminder that a refined group is conditional on the selected phenotype. It should not be interpreted as a universal genetic equivalence class.

Multiple traits make the refinement criterion stricter

Several numeric traits can be supplied together:

refined_multi_trait <- refine_haplotype(
  hap_gene,
  phenotype = pheno,
  traits = c("protein_content", "seed_weight"),
  min_hap_samples = 3,
  test_method = "t.test",
  p_adjust = "BH",
  alpha = 0.05,
  effect_threshold = 0.5
)
refined_multi_trait$haplotype_map
#>    hap_id refined_hap_id refined_key protein_content seed_weight
#>    <char>         <char>      <char>          <char>      <char>
#> 1:   Hap1          RHap1       G1|G2              G1          G2
#> 2:   Hap2          RHap2       G1|G1              G1          G1
#> 3:   Hap3          RHap3       G2|G4              G2          G4
#> 4:   Hap4          RHap4       G2|G3              G2          G3

Each trait is grouped independently. GeneTrackR then combines the trait-specific group labels, so two original haplotypes remain together in the final result only when their complete multi-trait grouping pattern is the same. In this demo, adding the strongly discriminating seed_weight trait prevents the 4 -> 2 collapse produced by protein_content alone.

For real data, evaluate the sensitivity of the result to:

  • min_hap_samples;
  • test_method;
  • p_adjust;
  • alpha;
  • effect_threshold;
  • the selected trait set.

Do not choose these parameters only to obtain a desired number of refined groups.

Step 9. Interpret and reuse refined haplotypes

Phenotype-guided refinement changes the group label, not the underlying VCF or the original haplotype definition. Keep both levels of information:

refined_protein$original_hap$haplotypes
#> Key: <hap_id>
#>    hap_id sample_n                             samples varA01 varA02 varA03
#>    <char>    <int>                              <char> <char> <char> <char>
#> 1:   Hap1        9 S10;S11;S12;S13;S14;S15;S16;S17;S18      A      C      C
#> 2:   Hap2        9 S01;S02;S03;S04;S05;S06;S07;S08;S09      A      T      C
#> 3:   Hap3        9 S19;S20;S21;S22;S23;S24;S25;S26;S27      G      C      G
#> 4:   Hap4        9 S28;S29;S30;S31;S32;S33;S34;S35;S36      G      T      G
#>    varLD01 varLD02 varLD03 varLD04 varLD05 varLD06 varA04 varA05
#>     <char>  <char>  <char>  <char>  <char>  <char> <char> <char>
#> 1:       G       A       C       G       T       A      A     i3
#> 2:       G       A       C       G       T       A      A     i3
#> 3:       A       T       T       C       G       C     i3     i3
#> 4:       A       T       T       C       G       C     i3      A
refined_protein$haplotype_map
#>    hap_id refined_hap_id refined_key protein_content
#>    <char>         <char>      <char>          <char>
#> 1:   Hap1          RHap1          G1              G1
#> 2:   Hap2          RHap1          G1              G1
#> 3:   Hap3          RHap2          G2              G2
#> 4:   Hap4          RHap2          G2              G2
refined_protein$refined_haplotypes
#>    hap_id original_hap_n original_hap_ids sample_n
#>    <char>          <int>           <char>    <int>
#> 1:  RHap1              2        Hap1;Hap2       18
#> 2:  RHap2              2        Hap3;Hap4       18
#>                                                                    samples
#>                                                                     <char>
#> 1: S10;S11;S12;S13;S14;S15;S16;S17;S18;S01;S02;S03;S04;S05;S06;S07;S08;S09
#> 2: S19;S20;S21;S22;S23;S24;S25;S26;S27;S28;S29;S30;S31;S32;S33;S34;S35;S36
#>    varA01 varA02 varA03 varLD01 varLD02 varLD03 varLD04 varLD05 varLD06 varA04
#>    <char> <char> <char>  <char>  <char>  <char>  <char>  <char>  <char> <char>
#> 1:      A  mixed      C       G       A       C       G       T       A      A
#> 2:      G  mixed      G       A       T       T       C       G       C     i3
#>    varA05
#>    <char>
#> 1:     i3
#> 2:  mixed

For the protein_content demo, the final grouping is expected to be Hap1/Hap2 versus Hap3/Hap4. The two groups retain different allele classes at varA03 and at the broader p13 cluster, including varA01, varLD01-varLD06, and varA04. Within the low group, varA02 is heterogeneous; within the high group, varA02 and varA05 are heterogeneous. These positions appear as mixed in the collapsed view. This is the intended refinement logic: genetically close haplotypes with the same phenotype class are combined, while the variants that consistently separate the two genotype/phenotype clusters remain visible.

The refined object can be passed directly to the refined plotting wrappers, while the original hap_gene remains the preferred input for variant-level effect prioritization:

# Refined-group visualization
# plot_refined_hap_variant(refined_protein, annotation = gp, text_size = 9)
# plot_refined_hap_pheno(refined_protein, phenotype = pheno, traits = "protein_content", text_size = 9)
# Variant-level prioritization keeps the original sample genotypes
# effect <- plot_variant_effect(
#   hap_gene,
#   phenotype = pheno,
#   traits = "protein_content",
#   min_group_samples = 3
# )

For reproducible analysis, retain and export at least:

  • haplotype_map;
  • pairwise_test;
  • phenotype_summary;
  • refined_haplotypes;
  • sample_refined_haplotypes;
  • the refinement parameters;
  • the grouped and/or collapsed figure.

This preserves the evidence needed to explain why each original haplotype was or was not merged.

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] GeneTrackR_0.99.2 rmarkdown_2.32   
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6         jsonlite_2.0.0       crayon_1.5.3        
#>  [4] compiler_4.6.1       zip_3.0.2            Rcpp_1.1.2          
#>  [7] Biostrings_2.81.9    bitops_1.1-0         Rsamtools_2.29.0    
#> [10] GenomicRanges_1.65.4 parallel_4.6.1       jquerylib_0.1.4     
#> [13] IRanges_2.47.5       Seqinfo_1.3.2        scales_1.4.0        
#> [16] BiocParallel_1.47.0  yaml_2.3.12          fastmap_1.2.0       
#> [19] XVector_0.53.0       ggplot2_4.0.3        R6_2.6.1            
#> [22] patchwork_1.3.2      labeling_0.4.3       generics_0.1.4      
#> [25] openxlsx_4.2.9       knitr_1.52           BiocGenerics_0.59.12
#> [28] maketools_1.3.2      bslib_0.12.0         RColorBrewer_1.1-3  
#> [31] rlang_1.3.0          stringi_1.8.9        cachem_1.1.0        
#> [34] xfun_0.61            sass_0.4.10          sys_3.4.3           
#> [37] S7_0.2.2             otel_0.2.0           cli_3.6.6           
#> [40] withr_3.0.3          digest_0.6.39        grid_4.6.1          
#> [43] lifecycle_1.0.5      S4Vectors_0.51.10    vctrs_0.7.3         
#> [46] evaluate_1.0.5       glue_1.8.1           data.table_1.18.6.1 
#> [49] farver_2.1.2         codetools_0.2-20     buildtools_1.0.0    
#> [52] stats4_4.6.1         tools_4.6.1          htmltools_0.5.9