---
title: "Variant effect analysis"
author: "Shuchao Ren"
date: "2026-08-31"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{13 Variant effect analysis}
  %\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 = 10.5,
  fig.height = 6.5,
  dpi = gtr_vignette_dpi,
  fig.retina = 1,
  out.width = "100%",
  fig.align = "center"
)
```

## Introduction

```{r load-package}
library(GeneTrackR)
```

`plot_variant_effect()` evaluates every natural variant retained in a `HapVariant` object against one or more numeric phenotype traits. It is intended for **regional prioritization**: estimate genotype-associated phenotype effects across all variants, rank the strongest candidates, and then interpret those candidates together with LD, gene annotation, single-variant phenotype plots, and haplotype refinement.

The GeneA demo is deliberately useful for this workflow. `varA03` is the designed `protein_content` effect variant, but it shares the same `p13` genotype partition with several neighboring variants. Consequently, multiple linked variants have the same large phenotype effect. This demonstrates an important principle: a variant-effect scan can identify an associated variant set, but effect size alone cannot distinguish a causal variant from tightly linked passenger variants.

This module follows nine steps:

1. prepare a binary-coded GeneA `HapVariant` and phenotype table;
2. understand the genotype coding used for effect direction;
3. calculate the primary `protein_content` variant-effect scan;
4. inspect the `GeneTrackRVariantEffectPlot` result object;
5. rank and summarize candidate variants;
6. interpret tied effects in the context of shared genotype patterns and LD;
7. inspect signed effects and their interpretation limits;
8. validate the lead example with a direct single-variant phenotype analysis;
9. compare multiple traits and connect variant effects to refinement and downstream validation.

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 GeneA variants and phenotype data

Locate and read the same demo VCF, annotation, and phenotype files used in the haplotype and refinement workflows:

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

For effect screening, construct a second GeneA `HapVariant` with `genotype_mode = "code"`:

```{r variant-effect-02}
hap_effect <- hap_gene_variant(
  vcf,
  annotation = gp,
  gene_id = "GeneA",
  genotype_mode = "code"
)
hap_effect
```

The region is the same canonical GeneA gene body used in the haplotype workflow: 11 variants and 36 samples. Only the genotype representation changes.

## Step 2. Understand genotype coding before interpreting effect direction

With `genotype_mode = "code"`, GeneTrackR uses a binary presence representation:

- `0`: reference genotype with no ALT allele present;
- `1`: at least one ALT allele is present.

This is useful for variant-effect screening because a two-group signed effect becomes interpretable as:

```text
mean(ALT-present group) - mean(reference group)
```

However, this coding is **not an additive 0/1/2 dosage model**. Heterozygous and homozygous-ALT genotypes are both represented as `1`. If heterozygous and homozygous-ALT effects must be separated, inspect the variant directly with `plot_variant_pheno(..., genotype_mode = "string")` rather than interpreting the binary screen as an additive genetic model.

The original string-coded `hap_gene` from the haplotype workflow remains preferable for visualizing exact allele patterns. The binary-coded `hap_effect` is created specifically to make the direction of two-group variant effects easier to interpret.

## Step 3. Calculate the primary protein-content effect scan

The main demo trait is `protein_content`. Run one effect calculation for every GeneA variant:

```{r variant-effect-03}
effect_protein <- plot_variant_effect(
  hap_effect,
  phenotype = pheno,
  traits = "protein_content",
  min_group_samples = 3,
  test_method = "t.test",
  p_adjust = "BH",
  effect_type = "absolute",
  top_n = 10,
  variant_label = "variant_id",
  x_axis = "position",
  text_size = 10
)
effect_protein$figure
```

The default screening view uses absolute effect size on the y-axis. Point size also reflects absolute effect size, point transparency reflects `-log10(adjusted P)`, and point color retains the signed effect direction.

Using `x_axis = "position"` places variants at their genomic coordinates. Use `x_axis = "index"` when evenly spaced variant labels are easier to read than physical distance.

## Step 4. Inspect the result object

`plot_variant_effect()` returns a `GeneTrackRVariantEffectPlot` object rather than a bare ggplot:

```{r variant-effect-04}
class(effect_protein)
names(effect_protein)
```

The three main components are:

| Component | Content | Typical use |
| --- | --- | --- |
| `figure` | effect-size ggplot | visualization and export |
| `effect` | one row per trait × variant | candidate ranking and statistical interpretation |
| `plot_data` | long sample-level genotype/phenotype table | detailed QC or custom plotting |

Inspect the effect table directly:

```{r variant-effect-05}
effect_protein$effect[, .(
  variant_id,
  chrom,
  pos,
  group_n,
  sample_n,
  effect,
  abs_effect,
  low_group,
  high_group,
  low_group_mean,
  high_group_mean,
  p_value,
  p_adj
)]
```

For a binary `code`-mode variant, `effect` is the mean difference between genotype group `1` and genotype group `0`. `abs_effect` is its absolute magnitude. `p_adj` is adjusted within each trait using the method supplied through `p_adjust`.

## Step 5. Rank candidate variants

Sort variants by absolute phenotype effect:

```{r variant-effect-06}
protein_rank <- effect_protein$effect[
  order(-abs_effect, p_adj),
  .(
    variant_id,
    pos,
    effect,
    abs_effect,
    p_value,
    p_adj
  )
]
protein_rank
```

The deterministic GeneA demo has the following designed effect structure:

```text
varA01 / varA03 / varLD01-varLD06 / varA04 : absolute effect ~= 6
varA05                                       : absolute effect ~= 4
varA02                                       : absolute effect ~= 0
```

Retrieve the strongest class explicitly:

```{r variant-effect-07}
protein_top <- effect_protein$effect[
  abs_effect >= 5.9,
  .(variant_id, pos, effect, abs_effect, p_adj)
]
protein_top[order(pos)]
```

`varA03` is therefore correctly recovered as a high-effect variant, but it is **not uniquely ranked first**. That is expected and biologically informative rather than an error.

## Step 6. Interpret tied effects together with genotype structure and LD

The variants with an effect of approximately `+6` share the same `p13` genotype split in the deterministic demo: `Hap1/Hap2` are in the reference class and `Hap3/Hap4` are in the ALT-present class. They consequently partition the same 36 samples into the same two phenotype groups and receive the same effect estimate.

Inspect their variant metadata:

```{r variant-effect-08}
retrieve_vcf(
  vcf,
  variant_id = protein_top$variant_id,
  verbose = FALSE
)[, .(
  variant_id,
  chrom,
  pos,
  ref,
  alt,
  variant_type,
  info
)]
```

This is the key interpretation boundary of the effect scan:

```text
large phenotype effect
        !=
proof of causal variation
```

Variants that are identical or nearly identical in genotype pattern cannot be separated by association effect alone. Combine the effect ranking with:

- the linkage-disequilibrium article;
- gene position and functional annotation;
- coding/regulatory consequence information;
- independent populations or recombinants that break the linkage;
- experimental validation when causal interpretation is required.

In the demo, `varA03` is marked as the designed `protein_effect` variant, while `varLD01`-`varLD06` are explicitly retained as linked LD markers. The equal phenotype effects demonstrate why biological annotation and recombination information remain necessary after statistical prioritization.

## Step 7. Inspect signed effects carefully

A signed-effect plot can be generated directly:

```{r variant-effect-09}
effect_protein_signed <- plot_variant_effect(
  hap_effect,
  phenotype = pheno,
  traits = "protein_content",
  min_group_samples = 3,
  test_method = "t.test",
  p_adjust = "BH",
  effect_type = "signed",
  top_n = 10,
  variant_label = "variant_id",
  x_axis = "position",
  text_size = 10
)
effect_protein_signed$figure
```

For the present binary `code`-mode workflow, positive values mean that genotype group `1` has a higher phenotype mean than genotype group `0`. For `varA03`, this corresponds to ALT presence being associated with higher `protein_content`.

Check the designed site directly:

```{r variant-effect-10}
effect_protein_signed$effect[
  variant_id == "varA03",
  .(
    variant_id,
    effect,
    abs_effect,
    low_group,
    high_group,
    low_group_mean,
    high_group_mean,
    p_adj
  )
]
```

The expected effect is approximately `+6`.

Do not generalize this interpretation blindly to every possible `HapVariant`. With string-coded genotypes, the two-group order can follow genotype-label ordering rather than a universal REF-to-ALT direction. With more than two genotype groups, `plot_variant_effect()` reports the maximum-minus-minimum mean difference, so the result is primarily an effect magnitude rather than a simple ALT-minus-REF contrast. For region-wide prioritization, `effect_type = "absolute"` is therefore the safer default.

## Step 8. Validate varA03 with a direct single-variant phenotype plot

The regional effect scan tells us that `varA03` belongs to the high-effect set. Confirm the genotype/phenotype relationship directly with `plot_variant_pheno()`:

```{r variant-effect-11}
varA03_pheno <- plot_variant_pheno(
  variant = vcf,
  phenotype = pheno,
  variant_id = "varA03",
  traits = "protein_content",
  genotype_mode = "code",
  min_group_samples = 3,
  p_adjust = "BH",
  p_value_type = "adjusted",
  text_size = 10
)
varA03_pheno$figure
varA03_pheno$summary
varA03_pheno$pvalue
```

Retrieve the variant itself to connect the binary genotype classes to its alleles:

```{r variant-effect-12}
retrieve_vcf(
  vcf,
  variant_id = "varA03",
  verbose = FALSE
)[, .(
  variant_id,
  chrom,
  pos,
  ref,
  alt,
  variant_type,
  info
)]
```

The demo truth is:

```text
varA03: C > G
code 0: REF class, Hap1/Hap2, protein_content ~= 38
code 1: ALT-present class, Hap3/Hap4, protein_content ~= 44
ALT-present minus REF effect ~= +6
```

This direct check is preferable to inferring allele direction from the regional ranking plot alone.

## Step 9. Compare multiple traits and connect downstream evidence

`plot_variant_effect()` can evaluate several numeric traits in one call:

```{r variant-effect-13}
effect_multi <- plot_variant_effect(
  hap_effect,
  phenotype = pheno,
  traits = c(
    "protein_content",
    "seed_weight",
    "plant_height",
    "flowering_time"
  ),
  min_group_samples = 3,
  test_method = "t.test",
  p_adjust = "BH",
  effect_type = "absolute",
  top_n = 5,
  variant_label = "variant_id",
  x_axis = "position",
  text_size = 10
)
effect_multi$figure
```

Numeric traits may be stored as either integer or double columns. `plot_variant_effect()` normalizes selected numeric traits internally before reshaping them, so mixed numeric storage modes can be analyzed together without `data.table::melt()` coercion warnings.

Inspect the same variant across traits:

```{r variant-effect-14}
effect_multi$effect[
  variant_id == "varA03",
  .(
    trait,
    effect,
    abs_effect,
    p_value,
    p_adj
  )
]
```

For the deterministic demo, `varA03` has approximately:

```text
protein_content : +6
seed_weight     : +8
plant_height    : +6
flowering_time  :  0
```

The first three effects arise because these traits follow the broader GeneA genotype hierarchy, whereas `flowering_time` is the designed negative control. P-value adjustment is performed separately within each trait.

The complete evidence chain is therefore:

```text
GeneA HapVariant
    -> region-wide variant-effect ranking
    -> identify high-effect genotype partitions
    -> compare with LD and functional annotation
    -> validate specific variants with plot_variant_pheno()
    -> compare with haplotype refinement
    -> prioritize variants for biological validation
```

For the demo, `protein_content` refinement groups `Hap1/Hap2` and `Hap3/Hap4`, and `varA03` separates exactly those two refined genotype/phenotype clusters. The surrounding `p13` variants receive the same statistical effect because they carry the same genotype partition. That agreement is useful evidence for the regional signal, but it should not be mistaken for proof that every tied variant is causal.

## Session information

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