---
title: "Phenotype-guided haplotype refinement"
author: "Shuchao Ren"
date: "2026-08-31"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{12 Phenotype-guided haplotype refinement}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
gtr_pkgdown <- tolower(Sys.getenv("GENETRACKR_PKGDOWN", "false")) %in%
  c("1", "true", "yes")
gtr_vignette_dpi <- if (gtr_pkgdown) 300 else 144

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 12.5,
  fig.height = 8,
  dpi = gtr_vignette_dpi,
  fig.retina = 1,
  out.width = "100%",
  fig.align = "center"
)
```

## Introduction

```{r load-package}
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`:

```{r installation, eval=FALSE}
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:

```{r haplotype-refinement-01}
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:

```{r haplotype-refinement-02}
hap_gene <- hap_gene_variant(
  vcf,
  annotation = gp,
  gene_id = "GeneA",
  genotype_mode = "string"
)
hap_gene
hap_gene$haplotypes
```

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:

```{r haplotype-refinement-03}
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
```

The deterministic GeneA design is hierarchical:

```text
          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:

```{r haplotype-refinement-04}
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
protein_before$pvalue
```

The expected phenotype means are:

```text
Hap1 ≈ 38
Hap2 ≈ 38
Hap3 ≈ 44
Hap4 ≈ 44
```

This correspondence is intentional:

```text
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:

```{r haplotype-refinement-05}
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
refined_protein$refined_haplotypes
```

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

```text
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:

```{r haplotype-refinement-06}
names(refined_protein)
```

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:

```{r haplotype-refinement-07}
refined_protein$haplotype_map
refined_protein$sample_refined_haplotypes
```

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

```{r haplotype-refinement-08}
refined_protein$refined_hap
refined_protein$refined_hap$haplotypes
```

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:

```{r haplotype-refinement-09}
refined_protein$pairwise_test[, c(
  "trait",
  "group1",
  "group2",
  "mean_group1",
  "mean_group2",
  "abs_mean_diff",
  "p_value",
  "p_adj",
  "can_merge"
)]
```

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:

```{r haplotype-refinement-10}
refined_protein$trait_group_map
```

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

```{r haplotype-refinement-11}
refined_protein$haplotype_map
```

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:

```{r haplotype-refinement-12}
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:

```{r haplotype-refinement-13}
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:

```{r haplotype-refinement-14}
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
refined_protein_pheno$pvalue
```

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:

```{r haplotype-refinement-15}
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
```

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:

```{r haplotype-refinement-16}
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
```

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:

```{r haplotype-refinement-17}
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
```

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:

```{r haplotype-refinement-18}
refined_protein$original_hap$haplotypes
refined_protein$haplotype_map
refined_protein$refined_haplotypes
```

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:

```{r haplotype-refinement-19}
# 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

```{r session-info}
sessionInfo()
```