Variant effect analysis

Introduction

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:

prepare a binary-coded GeneA HapVariant and phenotype table;
understand the genotype coding used for effect direction;
calculate the primary protein_content variant-effect scan;
inspect the GeneTrackRVariantEffectPlot result object;
rank and summarize candidate variants;
interpret tied effects in the context of shared genotype patterns and LD;
inspect signed effects and their interpretation limits;
validate the lead example with a direct single-variant phenotype analysis;
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:

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:

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":

hap_effect <- hap_gene_variant(
  vcf,
  annotation = gp,
  gene_id = "GeneA",
  genotype_mode = "code"
)
#> [GeneTrackR] Retrieved variants: 11.
hap_effect
#> <HapVariant>
#>   region    : chr1:12340001-12352000
#>   variants  : 11
#>   samples   : 36
#>   haplotypes: 4

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:

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:

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:

class(effect_protein)
#> [1] "GeneTrackRVariantEffectPlot" "list"
names(effect_protein)
#> [1] "figure"    "effect"    "plot_data"

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:

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
)]
#>     variant_id  chrom      pos group_n sample_n effect abs_effect low_group
#>         <char> <char>    <int>   <int>    <int>  <num>      <num>    <char>
#>  1:     varA01   chr1 12340250       2       36      6          6         0
#>  2:     varA02   chr1 12340600       2       36      0          0         0
#>  3:     varA03   chr1 12342550       2       36      6          6         0
#>  4:    varLD01   chr1 12342620       2       36      6          6         0
#>  5:    varLD02   chr1 12342710       2       36      6          6         0
#>  6:    varLD03   chr1 12342805       2       36      6          6         0
#>  7:    varLD04   chr1 12342920       2       36      6          6         0
#>  8:    varLD05   chr1 12343040       2       36      6          6         0
#>  9:    varLD06   chr1 12343180       2       36      6          6         0
#> 10:     varA04   chr1 12344500       2       36      6          6         0
#> 11:     varA05   chr1 12351050       2       36      4          4         0
#>     high_group low_group_mean high_group_mean      p_value        p_adj
#>         <char>          <num>           <num>        <num>        <num>
#>  1:          1             38              44 2.008652e-37 2.455019e-37
#>  2:          1             41              41 1.000000e+00 1.000000e+00
#>  3:          1             38              44 2.008652e-37 2.455019e-37
#>  4:          1             38              44 2.008652e-37 2.455019e-37
#>  5:          1             38              44 2.008652e-37 2.455019e-37
#>  6:          1             38              44 2.008652e-37 2.455019e-37
#>  7:          1             38              44 2.008652e-37 2.455019e-37
#>  8:          1             38              44 2.008652e-37 2.455019e-37
#>  9:          1             38              44 2.008652e-37 2.455019e-37
#> 10:          1             38              44 2.008652e-37 2.455019e-37
#> 11:          1             40              44 1.190766e-07 1.309843e-07

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:

protein_rank <- effect_protein$effect[
  order(-abs_effect, p_adj),
  .(
    variant_id,
    pos,
    effect,
    abs_effect,
    p_value,
    p_adj
  )
]
protein_rank
#>     variant_id      pos effect abs_effect      p_value        p_adj
#>         <char>    <int>  <num>      <num>        <num>        <num>
#>  1:     varA01 12340250      6          6 2.008652e-37 2.455019e-37
#>  2:     varA03 12342550      6          6 2.008652e-37 2.455019e-37
#>  3:    varLD01 12342620      6          6 2.008652e-37 2.455019e-37
#>  4:    varLD02 12342710      6          6 2.008652e-37 2.455019e-37
#>  5:    varLD03 12342805      6          6 2.008652e-37 2.455019e-37
#>  6:    varLD04 12342920      6          6 2.008652e-37 2.455019e-37
#>  7:    varLD05 12343040      6          6 2.008652e-37 2.455019e-37
#>  8:    varLD06 12343180      6          6 2.008652e-37 2.455019e-37
#>  9:     varA04 12344500      6          6 2.008652e-37 2.455019e-37
#> 10:     varA05 12351050      4          4 1.190766e-07 1.309843e-07
#> 11:     varA02 12340600      0          0 1.000000e+00 1.000000e+00

The deterministic GeneA demo has the following designed effect structure:

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

Retrieve the strongest class explicitly:

protein_top <- effect_protein$effect[
  abs_effect >= 5.9,
  .(variant_id, pos, effect, abs_effect, p_adj)
]
protein_top[order(pos)]
#>    variant_id      pos effect abs_effect        p_adj
#>        <char>    <int>  <num>      <num>        <num>
#> 1:     varA01 12340250      6          6 2.455019e-37
#> 2:     varA03 12342550      6          6 2.455019e-37
#> 3:    varLD01 12342620      6          6 2.455019e-37
#> 4:    varLD02 12342710      6          6 2.455019e-37
#> 5:    varLD03 12342805      6          6 2.455019e-37
#> 6:    varLD04 12342920      6          6 2.455019e-37
#> 7:    varLD05 12343040      6          6 2.455019e-37
#> 8:    varLD06 12343180      6          6 2.455019e-37
#> 9:     varA04 12344500      6          6 2.455019e-37

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:

retrieve_vcf(
  vcf,
  variant_id = protein_top$variant_id,
  verbose = FALSE
)[, .(
  variant_id,
  chrom,
  pos,
  ref,
  alt,
  variant_type,
  info
)]
#>    variant_id  chrom      pos    ref    alt variant_type
#>        <char> <char>    <int> <char> <char>       <char>
#> 1:     varA01   chr1 12340250      A      G          SNP
#> 2:     varA03   chr1 12342550      C      G          SNP
#> 3:    varLD01   chr1 12342620      G      A          SNP
#> 4:    varLD02   chr1 12342710      A      T          SNP
#> 5:    varLD03   chr1 12342805      C      T          SNP
#> 6:    varLD04   chr1 12342920      G      C          SNP
#> 7:    varLD05   chr1 12343040      T      G          SNP
#> 8:    varLD06   chr1 12343180      A      C          SNP
#> 9:     varA04   chr1 12344500      A    ATG          INS
#>                                  info
#>                                <char>
#> 1:        ROLE=core_haplotype,high_ld
#> 2: ROLE=core_haplotype,protein_effect
#> 3:           ROLE=high_ld,ld_gradient
#> 4:           ROLE=high_ld,ld_gradient
#> 5:           ROLE=high_ld,ld_gradient
#> 6:           ROLE=high_ld,ld_gradient
#> 7:           ROLE=high_ld,ld_gradient
#> 8:           ROLE=high_ld,ld_gradient
#> 9:          ROLE=core_haplotype,indel

This is the key interpretation boundary of the effect scan:

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:

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:

effect_protein_signed$effect[
  variant_id == "varA03",
  .(
    variant_id,
    effect,
    abs_effect,
    low_group,
    high_group,
    low_group_mean,
    high_group_mean,
    p_adj
  )
]
#>    variant_id effect abs_effect low_group high_group low_group_mean
#>        <char>  <num>      <num>    <char>     <char>          <num>
#> 1:     varA03      6          6         0          1             38
#>    high_group_mean        p_adj
#>              <num>        <num>
#> 1:              44 2.455019e-37

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():

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
#>              trait hap_id sample_n median_value mean_value min_value max_value
#>             <char> <fctr>    <int>        <num>      <num>     <num>     <num>
#> 1: protein_content      0       18           38         38     37.58     38.42
#> 2: protein_content      1       18           44         44     43.58     44.42
varA03_pheno$pvalue
#>              trait group1 group2      p_value        p_adj method
#>             <char> <char> <char>        <num>        <num> <char>
#> 1: protein_content      0      1 2.008652e-37 2.008652e-37 t.test

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

retrieve_vcf(
  vcf,
  variant_id = "varA03",
  verbose = FALSE
)[, .(
  variant_id,
  chrom,
  pos,
  ref,
  alt,
  variant_type,
  info
)]
#>    variant_id  chrom      pos    ref    alt variant_type
#>        <char> <char>    <int> <char> <char>       <char>
#> 1:     varA03   chr1 12342550      C      G          SNP
#>                                  info
#>                                <char>
#> 1: ROLE=core_haplotype,protein_effect

The demo truth is:

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:

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:

effect_multi$effect[
  variant_id == "varA03",
  .(
    trait,
    effect,
    abs_effect,
    p_value,
    p_adj
  )
]
#>              trait effect abs_effect      p_value        p_adj
#>             <char>  <num>      <num>        <num>        <num>
#> 1: protein_content      6          6 2.008652e-37 2.455019e-37
#> 2:     seed_weight      8          8 2.262764e-19 2.765601e-19
#> 3:    plant_height      6          6 2.669514e-13 3.262739e-13
#> 4:  flowering_time      0          0 1.000000e+00 1.000000e+00

For the deterministic demo, varA03 has approximately:

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:

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

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