---
title: "Export and reproducibility"
author: "Shuchao Ren"
date: "2026-08-31"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{15 Export and reproducibility}
  %\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

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.

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

```{r export-results-01}
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:

```{r export-results-02}
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"
)

hap_effect <- hap_gene_variant(
  vcf,
  annotation = gp,
  gene_id = "GeneA",
  genotype_mode = "code"
)

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:

```{r export-results-03}
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()`:

```{r export-results-04}
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
```

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

```{r export-results-05}
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
```

#### Variant tracks

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

```{r export-results-06}
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.

```{r export-results-07}
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:

```{r export-results-08}
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:

```{r export-results-09}
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:

```{r export-results-10}
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:

```{r export-results-11}
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:

```{r export-results-12}
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:

```{r export-results-13}
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:

```{r export-results-14}
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:

```{r export-results-15}
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:

```{r export-results-16}
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:

```{r export-results-17}
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:

```{r export-results-18}
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
```

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

```{r export-results-19}
ld_reloaded <- readRDS(
  file.path(object_dir, "GeneA.LDTrack.rds")
)

class(ld_reloaded)
head(ld_reloaded$data)
```

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

```{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.

## Session information

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