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

## Introduction

GeneTrackR uses `BwgTrack` as the common object for continuous genomic signals from **bedGraph**, **wig**, and **bigWig** files. Signal reading, regional querying, and writing are implemented in pure R. The normal signal workflow is:

1. read signal metadata or records into a `BwgTrack`;
2. validate and inspect sample metadata;
3. retrieve only the genomic interval needed for analysis;
4. optionally normalize, merge, or bin tracks;
5. write standardized signal files when required;
6. visualize signal at gene, transcript, or regional scale;
7. combine signal with the gene model and other genomic tracks.

The deterministic demo provides separate plus- and minus-strand RNA-seq and Ribo-seq bedGraph files.

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

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

rnaseq_files <- system.file(
  "extdata",
  c(
    "gtr_demo_rnaseq_plus.bedgraph",
    "gtr_demo_rnaseq_minus.bedgraph"
  ),
  package = "GeneTrackR"
)

riboseq_files <- system.file(
  "extdata",
  c(
    "gtr_demo_riboseq_plus.bedgraph",
    "gtr_demo_riboseq_minus.bedgraph"
  ),
  package = "GeneTrackR"
)
```

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.

## Read signal tracks

For small inputs or repeated local analysis, `mode = "memory"` loads all signal records immediately.

```{r read-signal}
rnaseq <- read_bwg(
  rnaseq_files,
  format = "bedgraph",
  sample_names = c("RNA_seq_plus", "RNA_seq_minus"),
  strand = c("+", "-"),
  mode = "memory",
  verbose = FALSE
)

riboseq <- read_bwg(
  riboseq_files,
  format = "bedgraph",
  sample_names = c("Ribo_seq_plus", "Ribo_seq_minus"),
  strand = c("+", "-"),
  mode = "memory",
  verbose = FALSE
)

summary_bwg(rnaseq)
summary_bwg(riboseq)
```

A `BwgTrack` keeps sample-level metadata separately from signal records. Inspect `$samples` when strand, format, normalization, or backing-file information matters.

```{r inspect-samples}
rnaseq$samples[, c("sample_id", "format", "strand"), with = FALSE]
```

## Lazy versus memory mode

For large real-world bigWig or indexed bedGraph files, prefer `mode = "lazy"`. Lazy mode keeps file metadata and retrieves only the requested interval instead of loading the whole signal track.

```{r lazy-signal, eval=FALSE}
signal_lazy <- read_bwg(
  files = c("sample1.bigwig", "sample2.bigwig"),
  format = "bigwig",
  sample_names = c("sample1", "sample2"),
  mode = "lazy"
)
```

Use memory mode when the track is small enough to fit comfortably in memory and will be queried repeatedly over many loci. Full-memory bigWig loading can be expensive; GeneTrackR therefore warns when a bigWig is loaded in memory mode.

## Strand semantics

Standard bigWig and wig files do not encode strand at the record level. GeneTrackR therefore treats them as unstranded signals. For paired plus/minus **bedGraph** files, an explicit `strand = c("+", "-")` records the source-file strand in the `BwgTrack` sample table and enables strand-aware filtering.

With gene- or transcript-based retrieval, `strand = "auto"` can use the strand of the annotation target.

## Retrieve a genomic interval

`retrieve_bwg()` supports direct genomic coordinates and annotation-aware gene/transcript retrieval. This is the preferred interface for lazy signal objects.

```{r retrieve-gene}
rnaseq_genea <- retrieve_bwg(
  rnaseq,
  annotation = gp,
  gene_id = "GeneA",
  strand = "auto",
  as = "data.table"
)

head(rnaseq_genea)
```

The same operation can return a new in-memory `BwgTrack` when a reusable regional object is more convenient:

```{r retrieve-object}
rnaseq_region <- retrieve_bwg(
  rnaseq,
  chrom = "chr1",
  start = 12339001,
  end = 12352000,
  strand = "+",
  as = "BwgTrack"
)

rnaseq_region
```

`slice_bwg()` provides a compact coordinate-only interface when annotation-aware resolution is not required.

## Bioconductor interoperability

In-memory `BwgTrack` objects can be converted directly to `GenomicRanges::GRanges`. Genomic coordinates and strand are retained, while `sample_id` and the signal value are represented as metadata (`score`):

```{r signal-granges}
rnaseq_gr <- as_granges(rnaseq_region)
rnaseq_gr
```

Lazy tracks should remain lazy until a bounded interval is needed. `retrieve_bwg(..., as = "GRanges")` provides that interval directly in a standard Bioconductor container.

```{r retrieve-granges}
rnaseq_gene_gr <- retrieve_bwg(
  rnaseq,
  annotation = gp,
  gene_id = "GeneA",
  strand = "auto",
  as = "GRanges"
)

rnaseq_gene_gr
```

Bioconductor's `rtracklayer` package provides broad standard genomic-file import/export infrastructure. GeneTrackR retains its native pure-R bedGraph/WIG/BigWig layer because the `BwgTrack` workflow combines multi-sample metadata, lazy/memory access, strand policy, normalization state, and a single API across the three signal formats. The native backend is not intended to replace `rtracklayer`; `GRanges` is the exchange boundary when GeneTrackR results need to participate in the wider Bioconductor ecosystem.

## Normalize and merge signals

`norm_bwg()` records and applies supported sample-level normalization strategies. For RPM/CPM normalization, supply appropriate library sizes; for externally defined scaling, use named scale factors.

```{r normalize-signal}
rnaseq_scaled <- norm_bwg(
  rnaseq,
  method = "scale",
  scale_factor = c(
    RNA_seq_plus = 1,
    RNA_seq_minus = 1
  )
)
```

RNA-seq and Ribo-seq tracks have distinct sample IDs in the demo and can be merged directly:

```{r merge-signal}
signal_all <- merge_bwg(rnaseq, riboseq)
summary_bwg(signal_all)
```

When duplicated sample IDs exist, resolve them deliberately through `sample_conflict`; do not silently merge biologically distinct samples.

## Bin retrieved signal

For broad regions, fixed-width binning can reduce plotting density or produce a compact summary table.

```{r bin-signal}
genea_binned <- bin_bwg(
  rnaseq_genea,
  bin_size = 50
)

head(genea_binned)
```

Binning changes the resolution of the signal representation. It should therefore be chosen according to the biological scale of the analysis rather than only for visual convenience.

## Write bigWig output

`write_bwg()` writes bedGraph, wig, or bigWig from an in-memory `BwgTrack`. All three formats are written in pure R. BigWig output requires chromosome sizes, but it does not require a compiled library or an external `bedGraphToBigWig` executable.

```{r write-bigwig}
chrom_sizes_file <- system.file(
  "extdata",
  "gtr_demo.chrom.sizes",
  package = "GeneTrackR"
)

bigwig_dir <- file.path(tempdir(), "GeneTrackR_demo_bigwig")
dir.create(bigwig_dir, recursive = TRUE, showWarnings = FALSE)

rnaseq_bigwig <- write_bwg(
  rnaseq,
  outdir = bigwig_dir,
  format = "bigwig",
  chrom_sizes = chrom_sizes_file,
  overwrite = TRUE
)

rnaseq_bigwig
```

GeneTrackR stores in-memory intervals as 1-based closed coordinates and handles the file-format conversion during writing. Do not manually shift coordinates before `write_bwg()`.

## Plot signal over a gene

`plot_signal_gene()` resolves the genomic span from the annotation. For strand-specific bedGraph inputs, `strand = "auto"` selects the track matching the target gene strand.

```{r plot-rnaseq-gene}
rnaseq_gene_figure <- plot_signal_gene(
  signal = rnaseq,
  annotation = gp,
  gene_id = "GeneA",
  plot_type = "bar",
  strand = "auto",
  signal_palette = "Blues",
  signal_palette_direction = -1,
  signal_y_scale = "fixed",
  signal_track_height = 4,
  gene_track_height = 1,
  direction_mode = "end",
  plot_theme = "classic",
  text_size = 10
)

rnaseq_gene_figure
```

`signal_track_height` and `gene_track_height` control only panel proportions; they do not rescale the biological signal.

## Plot signal over a transcript

`plot_signal_transcript()` can use genomic coordinates or spliced transcript coordinates. For RNA-seq coverage, a bar plot in transcript coordinates removes intronic distance:

```{r plot-rnaseq-transcript}
plot_signal_transcript(
  signal = rnaseq,
  annotation = gp,
  transcript_id = "TxA1",
  coordinate = "transcript",
  plot_type = "bar",
  strand = "auto",
  signal_palette = "Blues",
  signal_palette_direction = -1,
  signal_y_scale = "fixed",
  direction_mode = "end",
  plot_theme = "classic",
  text_size = 10
)
```

## Ribo-seq reading-frame view

For Ribo-seq P-site-like signals, `plot_type = "frame"` maps counts to transcript reading frame. This view is only available for transcript-level plotting because the frame is defined relative to a transcript CDS.

```{r plot-riboseq-frame}
plot_signal_transcript(
  signal = riboseq,
  annotation = gp,
  transcript_id = "TxA1",
  coordinate = "transcript",
  plot_type = "frame",
  strand = "auto",
  frame_palette = "Set1",
  signal_track_height = 4,
  gene_track_height = 1,
  direction_mode = "end",
  plot_theme = "classic",
  text_size = 10
)
```

The bundled Ribo-seq demo intentionally contains a dominant frame-0 signal with lower irregular frame-1/frame-2 counts. It is designed to exercise frame-aware plotting rather than to represent a statistical periodicity test.

## Plot a genomic region

`plot_signal_region()` is appropriate for intervals containing several genes or for regions not naturally defined by a single gene/transcript.

```{r plot-rnaseq-region}
plot_signal_region(
  signal = rnaseq,
  annotation = gp,
  chrom = "chr1",
  start = 12339001,
  end = 12374500,
  strand = "both",
  plot_type = "bar",
  signal_palette = "Blues",
  signal_palette_direction = -1,
  signal_y_scale = "free",
  signal_y_ticks = "pretty",
  signal_track_height = 4,
  gene_track_height = 1,
  direction_mode = "end",
  plot_theme = "classic",
  text_size = 10
)
```

For comparisons between samples, choose `signal_y_scale` deliberately. `"fixed"` preserves a common y scale; `"free"` emphasizes within-track shape but should not be used to visually compare absolute signal magnitude between panels.

## Signal summaries and visual transformations

The plotting functions support `signal_summary` for grouped summaries and `signal_transform` for visual transformation (`none`, `log2`, `log10`, or `sqrt`). These controls affect the displayed representation and should not be confused with biological normalization performed by `norm_bwg()`.

Similarly, `bin_size` and heatmap binning parameters alter plotting resolution. Preserve unmodified signal objects when the downstream analysis requires original base-resolution values.

## Integrated track figure

`plot_tracks()` combines the gene model with signal and, in later workflows, feature or variant tracks. For a positive-strand GeneA view, select the plus-strand RNA-seq and Ribo-seq samples explicitly:

```{r integrated-tracks}
plot_tracks(
  annotation = gp,
  signal = signal_all,
  gene_id = "GeneA",
  samples = c("RNA_seq_plus", "Ribo_seq_plus"),
  strand = "+",
  signal_type = "bar",
  signal_palette = "Dark2",
  signal_y_scale = "free",
  signal_y_ticks = "pretty",
  heights = c(
    signal = 4,
    gene = 1,
    feature = 0.8,
    variant = 0.7
  ),
  direction_mode = "end",
  plot_theme = "classic",
  text_size = 10
)
```

The named `heights` vector controls the relative vertical allocation of each panel. Missing feature or variant inputs do not need separate placeholder tracks.

## Practical guidance for large signal data

For large sequencing datasets:

- prefer `mode = "lazy"` for bigWig and indexed bedGraph;
- retrieve only the required locus before materializing signal records;
- retain explicit sample IDs and strand metadata;
- normalize with documented library sizes or scale factors rather than modifying values ad hoc;
- use binning for visualization only when base-resolution detail is unnecessary;
- keep the original `BwgTrack` unchanged when exploring alternative transformations or plotting settings.

This keeps the annotation coordinate frame, signal values, and downstream variant/haplotype analysis separable and reproducible.

## Session information

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