Export and reproducibility

Introduction

GeneTrackR outputs are intentionally heterogeneous. Genomic tracks should be written with the corresponding GeneTrackR writer, plotting functions may return a figure directly or store it inside a structured result object, statistical analyses contain several tables with different biological meanings, and complete S3 objects preserve the state needed to resume an analysis. Saving only the final image or flattening an entire result class into one table loses important provenance.

library(GeneTrackR)

This article defines a practical export contract around the deterministic GeneA demo. A reproducible project separates four output layers:

  1. genomic track files for use by GeneTrackR or genome browsers;
  2. PDF/PNG figures for reports and manuscripts;
  3. tabular results for inspection and downstream statistics;
  4. complete R objects for exact reloading without recalculation.

This module follows eight steps:

  1. prepare a structured output directory and representative GeneTrackR objects;
  2. export annotation, signal, and variant tracks;
  3. save direct and class-contained figures as PDF/PNG;
  4. export haplotype, phenotype, LD, refinement, and variant-effect tables;
  5. export matrix-shaped results such as the LD matrix;
  6. optionally collect related tables into an Excel workbook;
  7. save complete GeneTrackR objects with saveRDS();
  8. record a manifest and session information and verify that saved objects reload correctly.

Step 1. Prepare output directories and representative objects

Create separate directories for track files, figures, tables, and serialized R objects:

export_root <- file.path(tempdir(), "GeneTrackR_export")
export_dirs <- file.path(
  export_root,
  c("tracks", "figures", "tables", "objects")
)

invisible(lapply(
  export_dirs,
  dir.create,
  recursive = TRUE,
  showWarnings = FALSE
))

track_dir <- export_dirs[1]
figure_dir <- export_dirs[2]
table_dir <- export_dirs[3]
object_dir <- export_dirs[4]

Read the demo inputs and build representative downstream objects. These are the same object classes used in the preceding workflow modules:

gp <- read_genepred(
  system.file("extdata", "gtr_demo.genePredExt", package = "GeneTrackR"),
  format = "genePredExt",
  verbose = FALSE,
  progress = FALSE
)

vcf <- read_vcf(
  system.file("extdata", "gtr_demo_variants.vcf", package = "GeneTrackR"),
  mode = "memory",
  keep_genotype = TRUE,
  verbose = FALSE,
  progress = FALSE
)

pheno <- read_pheno(
  system.file("extdata", "gtr_demo_pheno.tsv", package = "GeneTrackR"),
  verbose = FALSE,
  progress = FALSE
)

hap_gene <- hap_gene_variant(
  vcf,
  annotation = gp,
  gene_id = "GeneA",
  genotype_mode = "string"
)
#> [GeneTrackR] Retrieved variants: 11.

hap_effect <- hap_gene_variant(
  vcf,
  annotation = gp,
  gene_id = "GeneA",
  genotype_mode = "code"
)
#> [GeneTrackR] Retrieved variants: 11.

hap_pheno <- plot_hap_pheno(
  hap_gene,
  phenotype = pheno,
  traits = "seed_weight",
  min_hap_samples = 3,
  p_value_type = "adjusted",
  text_size = 10
)

ld_result <- compute_ld_block(
  vcf,
  chrom = "chr1",
  start = 12342620,
  end = 12355500,
  variant_type = "snp",
  method = "r2",
  verbose = FALSE
)

ld_result <- plot_ld_block(
  ld_result,
  annotation = gp,
  show_region = TRUE,
  show_variant_labels = FALSE,
  font = 10
)

refined <- refine_haplotype(
  hap_gene,
  phenotype = pheno,
  traits = "protein_content",
  min_hap_samples = 3,
  effect_threshold = 0.5
)

refined_variant_figure <- plot_refined_hap_variant(
  refined,
  annotation = gp,
  min_hap_samples = 3,
  text_size = 10
)

refined_pheno <- plot_refined_hap_pheno(
  refined,
  phenotype = pheno,
  traits = "protein_content",
  min_hap_samples = 3,
  p_value_type = "adjusted",
  text_size = 10
)

effect_result <- plot_variant_effect(
  hap_effect,
  phenotype = pheno,
  traits = "protein_content",
  min_group_samples = 3,
  effect_type = "absolute",
  x_axis = "position",
  text_size = 10
)

Keeping these objects separate is useful because their public contracts differ. plot_hap_variant() returns a figure directly, whereas plot_hap_pheno() and plot_variant_effect() return result objects containing $figure; plot_ld_block() returns an updated LDTrack with its figure stored in $figure.

Step 2. Export genomic track files

Annotation tracks

Use write_feature() for annotation conversion. Annotation coordinate conventions and format-specific details are covered in the annotation module; the export module only centralizes representative calls:

write_feature(
  gp,
  file = file.path(track_dir, "GeneA_demo.genePredExt"),
  format = "genepredext",
  overwrite = TRUE
)

write_feature(
  gp,
  file = file.path(track_dir, "GeneA_demo.bed12"),
  format = "bed12",
  overwrite = TRUE
)

Signal tracks

Read an in-memory signal track and export it with write_bwg():

rnaseq_plus <- read_bwg(
  system.file(
    "extdata",
    "gtr_demo_rnaseq_plus.bedgraph",
    package = "GeneTrackR"
  ),
  format = "bedgraph",
  sample_names = "RNA_seq_plus",
  strand = "+",
  mode = "memory",
  verbose = FALSE
)

bedgraph_files <- write_bwg(
  rnaseq_plus,
  outdir = track_dir,
  format = "bedgraph",
  overwrite = TRUE
)

bedgraph_files
#>       sample_id                                                           file
#>          <char>                                                         <char>
#> 1: RNA_seq_plus /tmp/RtmpCT2eho/GeneTrackR_export/tracks/RNA_seq_plus.bedgraph
#>      format
#>      <char>
#> 1: bedgraph

bigWig output uses GeneTrackR’s built-in pure-R writer and requires chromosome sizes:

chrom_sizes_file <- system.file(
  "extdata",
  "gtr_demo.chrom.sizes",
  package = "GeneTrackR"
)

bigwig_files <- write_bwg(
  rnaseq_plus,
  outdir = track_dir,
  format = "bigwig",
  chrom_sizes = chrom_sizes_file,
  overwrite = TRUE
)

bigwig_files
#>       sample_id                                                         file
#>          <char>                                                       <char>
#> 1: RNA_seq_plus /tmp/RtmpCT2eho/GeneTrackR_export/tracks/RNA_seq_plus.bigwig
#>    format
#>    <char>
#> 1: bigwig

Variant tracks

Use write_vcf() when a standardized site-level VCF is required:

write_vcf(
  vcf,
  file = file.path(track_dir, "gtr_demo.sites.vcf"),
  overwrite = TRUE
)

The current write_vcf() writes the standard site columns CHROM, POS, ID, REF, ALT, QUAL, FILTER, and INFO. It does not guarantee round-trip preservation of FORMAT and sample genotype columns. Keep the original genotype-containing VCF as the archival source whenever genotypes must be preserved exactly.

Step 3. Save figures as PDF and PNG

Use the actual figure object rather than passing an entire GeneTrackR result class to ggsave().

Define a small local helper when both PDF and PNG are required. The helper keeps a publication-oriented default of 500 dpi, while the executable vignette calls below pass the build-aware gtr_vignette_dpi value so package checks do not spend time rendering manuscript-scale raster files. Set dpi = 300 to 500 explicitly for final publication export.

save_gtr_figure <- function(
  plot,
  stem,
  outdir,
  width = 8,
  height = 6,
  dpi = 500
) {
  ggplot2::ggsave(
    filename = file.path(outdir, paste0(stem, ".pdf")),
    plot = plot,
    width = width,
    height = height
  )

  ggplot2::ggsave(
    filename = file.path(outdir, paste0(stem, ".png")),
    plot = plot,
    width = width,
    height = height,
    dpi = dpi
  )
}

A direct-return figure can be saved immediately:

hap_variant_figure <- plot_hap_variant(
  hap_gene,
  annotation = gp,
  min_hap_samples = 3,
  text_size = 10
)

save_gtr_figure(
  hap_variant_figure,
  stem = "GeneA_haplotype_variant",
  outdir = figure_dir,
  width = 10,
  height = 7,
  dpi = gtr_vignette_dpi
)

For structured result classes, extract $figure explicitly:

save_gtr_figure(
  hap_pheno$figure,
  stem = "GeneA_seed_weight_haplotype",
  outdir = figure_dir,
  width = 7,
  height = 5,
  dpi = gtr_vignette_dpi
)

save_gtr_figure(
  ld_result$figure,
  stem = "GeneA_LD_gradient",
  outdir = figure_dir,
  width = 10,
  height = 8,
  dpi = gtr_vignette_dpi
)

save_gtr_figure(
  refined_variant_figure,
  stem = "GeneA_refined_haplotype_variant",
  outdir = figure_dir,
  width = 10,
  height = 7,
  dpi = gtr_vignette_dpi
)

save_gtr_figure(
  refined_pheno$figure,
  stem = "GeneA_refined_protein_content",
  outdir = figure_dir,
  width = 7,
  height = 5,
  dpi = gtr_vignette_dpi
)

save_gtr_figure(
  effect_result$figure,
  stem = "GeneA_protein_variant_effect",
  outdir = figure_dir,
  width = 9,
  height = 5,
  dpi = gtr_vignette_dpi
)

The same rule applies to plot_variant_pheno() and plot_refined_hap_pheno(): save their $figure component. plot_refined_hap_variant() returns its figure directly.

Step 4. Export analysis tables

A HapVariant contains several biologically different tables. Export them separately instead of flattening the complete class:

data.table::fwrite(
  hap_gene$variants,
  file.path(table_dir, "GeneA.variants.tsv"),
  sep = "\t"
)

data.table::fwrite(
  hap_gene$haplotypes,
  file.path(table_dir, "GeneA.haplotypes.tsv"),
  sep = "\t"
)

data.table::fwrite(
  hap_gene$sample_haplotypes,
  file.path(table_dir, "GeneA.sample_haplotypes.tsv"),
  sep = "\t"
)

data.table::fwrite(
  hap_gene$genotype_long,
  file.path(table_dir, "GeneA.genotype_long.tsv"),
  sep = "\t"
)

For phenotype association results, retain the statistical tables used to create the figure:

data.table::fwrite(
  hap_pheno$pvalue,
  file.path(table_dir, "GeneA.seed_weight.pvalue.tsv"),
  sep = "\t"
)

data.table::fwrite(
  hap_pheno$summary,
  file.path(table_dir, "GeneA.seed_weight.summary.tsv"),
  sep = "\t"
)

data.table::fwrite(
  hap_pheno$plot_data,
  file.path(table_dir, "GeneA.seed_weight.plot_data.tsv"),
  sep = "\t"
)

Export pairwise LD rather than only the heatmap image:

data.table::fwrite(
  ld_result$data,
  file.path(table_dir, "GeneA.LD_pairs.tsv"),
  sep = "\t"
)

data.table::fwrite(
  ld_result$variants,
  file.path(table_dir, "GeneA.LD_variants.tsv"),
  sep = "\t"
)

For phenotype-guided refinement, the mapping and decision tables are essential for explaining why original haplotypes were merged:

data.table::fwrite(
  refined$haplotype_map,
  file.path(table_dir, "GeneA.refined_haplotype_map.tsv"),
  sep = "\t"
)

data.table::fwrite(
  refined$refined_haplotypes,
  file.path(table_dir, "GeneA.refined_haplotypes.tsv"),
  sep = "\t"
)

data.table::fwrite(
  refined$sample_refined_haplotypes,
  file.path(table_dir, "GeneA.sample_refined_haplotypes.tsv"),
  sep = "\t"
)

data.table::fwrite(
  refined$pairwise_test,
  file.path(table_dir, "GeneA.refinement_pairwise_test.tsv"),
  sep = "\t"
)

data.table::fwrite(
  refined$phenotype_summary,
  file.path(table_dir, "GeneA.refinement_pheno_summary.tsv"),
  sep = "\t"
)

data.table::fwrite(
  refined_pheno$pvalue,
  file.path(table_dir, "GeneA.refined_protein.pvalue.tsv"),
  sep = "\t"
)

data.table::fwrite(
  refined_pheno$summary,
  file.path(table_dir, "GeneA.refined_protein.summary.tsv"),
  sep = "\t"
)

For regional variant-effect prioritization, retain both the summarized effects and the sample-level plotting table:

data.table::fwrite(
  effect_result$effect,
  file.path(table_dir, "GeneA.protein_variant_effect.tsv"),
  sep = "\t"
)

data.table::fwrite(
  effect_result$plot_data,
  file.path(table_dir, "GeneA.protein_variant_effect_plot_data.tsv"),
  sep = "\t"
)

Step 5. Export the LD matrix cleanly

A matrix should be converted to a rectangular table with an explicit row identifier before writing it:

ld_matrix_table <- data.table::as.data.table(
  ld_result$matrix,
  keep.rownames = "variant_id"
)

data.table::fwrite(
  ld_matrix_table,
  file.path(table_dir, "GeneA.LD_matrix.tsv"),
  sep = "\t"
)

This preserves both row and column variant IDs and is easier to reopen in R, Python, spreadsheet software, or downstream matrix tools than a matrix written without row names.

Step 6. Optionally collect tables in an Excel workbook

For collaborators who prefer one multi-sheet workbook, the same result tables can be collected with the optional openxlsx package. This is an interoperability step rather than a GeneTrackR dependency:

if (requireNamespace("openxlsx", quietly = TRUE)) {
  openxlsx::write.xlsx(
    list(
      variants = as.data.frame(hap_gene$variants),
      haplotypes = as.data.frame(hap_gene$haplotypes),
      sample_haps = as.data.frame(hap_gene$sample_haplotypes),
      hap_pvalue = as.data.frame(hap_pheno$pvalue),
      hap_summary = as.data.frame(hap_pheno$summary),
      LD_pairs = as.data.frame(ld_result$data),
      LD_matrix = as.data.frame(ld_matrix_table),
      refined_map = as.data.frame(refined$haplotype_map),
      refined_test = as.data.frame(refined$pairwise_test),
      variant_effect = as.data.frame(effect_result$effect)
    ),
    file = file.path(table_dir, "GeneA.analysis.xlsx"),
    overwrite = TRUE
  )
}

TSV remains the dependency-free default and is preferable for very large tables. Excel is most useful for compact result summaries intended for manual inspection.

Step 7. Save complete GeneTrackR objects with RDS

Tabular exports are easy to inspect, but they do not preserve the complete object state. Save the corresponding S3 objects when an analysis may need to be reopened without recomputation:

saveRDS(
  hap_gene,
  file.path(object_dir, "GeneA.HapVariant.rds")
)

saveRDS(
  hap_pheno,
  file.path(object_dir, "GeneA.seed_weight.GeneTrackRPhenoPlot.rds")
)

saveRDS(
  ld_result,
  file.path(object_dir, "GeneA.LDTrack.rds")
)

saveRDS(
  refined,
  file.path(object_dir, "GeneA.HapRefined.rds")
)

saveRDS(
  effect_result,
  file.path(object_dir, "GeneA.VariantEffect.rds")
)

The RDS file preserves class information, metadata, parameters, intermediate tables, and stored figures. It complements rather than replaces human-readable tables and standard genomic track files.

Step 8. Record provenance and verify the export

Record the package/session environment together with the result directory:

writeLines(
  capture.output(sessionInfo()),
  con = file.path(export_root, "sessionInfo.txt")
)

export_manifest <- data.table::data.table(
  file = list.files(export_root, recursive = TRUE)
)

data.table::fwrite(
  export_manifest,
  file.path(export_root, "manifest.tsv"),
  sep = "\t"
)

export_manifest
#>                                                  file
#>                                                <char>
#>  1:               figures/GeneA_haplotype_variant.pdf
#>  2:               figures/GeneA_haplotype_variant.png
#>  3:                     figures/GeneA_LD_gradient.pdf
#>  4:                     figures/GeneA_LD_gradient.png
#>  5:          figures/GeneA_protein_variant_effect.pdf
#>  6:          figures/GeneA_protein_variant_effect.png
#>  7:       figures/GeneA_refined_haplotype_variant.pdf
#>  8:       figures/GeneA_refined_haplotype_variant.png
#>  9:         figures/GeneA_refined_protein_content.pdf
#> 10:         figures/GeneA_refined_protein_content.png
#> 11:           figures/GeneA_seed_weight_haplotype.pdf
#> 12:           figures/GeneA_seed_weight_haplotype.png
#> 13:                      objects/GeneA.HapRefined.rds
#> 14:                      objects/GeneA.HapVariant.rds
#> 15:                         objects/GeneA.LDTrack.rds
#> 16: objects/GeneA.seed_weight.GeneTrackRPhenoPlot.rds
#> 17:                   objects/GeneA.VariantEffect.rds
#> 18:                                   sessionInfo.txt
#> 19:                        tables/GeneA.analysis.xlsx
#> 20:                    tables/GeneA.genotype_long.tsv
#> 21:                       tables/GeneA.haplotypes.tsv
#> 22:                        tables/GeneA.LD_matrix.tsv
#> 23:                         tables/GeneA.LD_pairs.tsv
#> 24:                      tables/GeneA.LD_variants.tsv
#> 25: tables/GeneA.protein_variant_effect_plot_data.tsv
#> 26:           tables/GeneA.protein_variant_effect.tsv
#> 27:            tables/GeneA.refined_haplotype_map.tsv
#> 28:               tables/GeneA.refined_haplotypes.tsv
#> 29:           tables/GeneA.refined_protein.pvalue.tsv
#> 30:          tables/GeneA.refined_protein.summary.tsv
#> 31:         tables/GeneA.refinement_pairwise_test.tsv
#> 32:         tables/GeneA.refinement_pheno_summary.tsv
#> 33:                tables/GeneA.sample_haplotypes.tsv
#> 34:        tables/GeneA.sample_refined_haplotypes.tsv
#> 35:            tables/GeneA.seed_weight.plot_data.tsv
#> 36:               tables/GeneA.seed_weight.pvalue.tsv
#> 37:              tables/GeneA.seed_weight.summary.tsv
#> 38:                         tables/GeneA.variants.tsv
#> 39:                           tracks/GeneA_demo.bed12
#> 40:                     tracks/GeneA_demo.genePredExt
#> 41:                         tracks/gtr_demo.sites.vcf
#> 42:                      tracks/RNA_seq_plus.bedgraph
#> 43:                        tracks/RNA_seq_plus.bigwig
#>                                                  file
#>                                                <char>

Finally, verify that a serialized analysis object can be restored with its expected class and results:

ld_reloaded <- readRDS(
  file.path(object_dir, "GeneA.LDTrack.rds")
)

class(ld_reloaded)
#> [1] "LDTrack"
head(ld_reloaded$data)
#>    variant_i variant_j index_i index_j    pos_i    pos_j distance_bp n_samples
#>       <char>    <char>   <int>   <int>    <int>    <int>       <int>     <int>
#> 1:   varLD01   varLD02       1       2 12342620 12342710          90        36
#> 2:   varLD01   varLD03       1       3 12342620 12342805         185        36
#> 3:   varLD01   varLD04       1       4 12342620 12342920         300        36
#> 4:   varLD01   varLD05       1       5 12342620 12343040         420        36
#> 5:   varLD01   varLD06       1       6 12342620 12343180         560        36
#> 6:   varLD01   varLD07       1       7 12342620 12352500        9880        36
#>            r        r2         D Dprime_signed Dprime   p_i   p_j        ld
#>        <num>     <num>     <num>         <num>  <num> <num> <num>     <num>
#> 1: 1.0000000 1.0000000 0.5000000             1      1   0.5   0.5 1.0000000
#> 2: 1.0000000 1.0000000 0.5000000             1      1   0.5   0.5 1.0000000
#> 3: 1.0000000 1.0000000 0.5000000             1      1   0.5   0.5 1.0000000
#> 4: 1.0000000 1.0000000 0.5000000             1      1   0.5   0.5 1.0000000
#> 5: 1.0000000 1.0000000 0.5000000             1      1   0.5   0.5 1.0000000
#> 6: 0.8888889 0.7901235 0.4444444             1      1   0.5   0.5 0.7901235
#>    method
#>    <char>
#> 1:     r2
#> 2:     r2
#> 3:     r2
#> 4:     r2
#> 5:     r2
#> 6:     r2

For a complete project, keep the original inputs, standard-format track exports, PDF/PNG figures, tabular results, RDS analysis objects, and session information together. This makes the analysis both easy to inspect and reproducible at the object level.

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.

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