Genomic signal tracks

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.

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:

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.

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)
#>        sample_id sum_signal mean_signal max_signal covered_bases
#>           <char>      <num>       <num>      <num>         <int>
#> 1:  RNA_seq_plus   632305.6    13.71596     33.264         46100
#> 2: RNA_seq_minus   361448.2    12.46373     28.512         29000
summary_bwg(riboseq)
#>         sample_id sum_signal mean_signal max_signal covered_bases
#>            <char>      <num>       <num>      <num>         <int>
#> 1:  Ribo_seq_plus      83318    4.124858         21         20199
#> 2: Ribo_seq_minus      45395    3.977133         19         11414

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

rnaseq$samples[, c("sample_id", "format", "strand"), with = FALSE]
#>        sample_id   format strand
#>           <char>   <char> <char>
#> 1:  RNA_seq_plus bedgraph      +
#> 2: RNA_seq_minus bedgraph      -

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.

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.

rnaseq_genea <- retrieve_bwg(
  rnaseq,
  annotation = gp,
  gene_id = "GeneA",
  strand = "auto",
  as = "data.table"
)

head(rnaseq_genea)
#>       sample_id  chrom    start      end value strand
#>          <char> <char>    <int>    <int> <num> <char>
#> 1: RNA_seq_plus   chr1 12340001 12340100 25.76      +
#> 2: RNA_seq_plus   chr1 12340101 12340200 28.00      +
#> 3: RNA_seq_plus   chr1 12340201 12340300 30.24      +
#> 4: RNA_seq_plus   chr1 12340301 12340400 26.88      +
#> 5: RNA_seq_plus   chr1 12340401 12340500 29.12      +
#> 6: RNA_seq_plus   chr1 12340501 12340600 25.76      +

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

rnaseq_region <- retrieve_bwg(
  rnaseq,
  chrom = "chr1",
  start = 12339001,
  end = 12352000,
  strand = "+",
  as = "BwgTrack"
)

rnaseq_region
#> <BwgTrack>
#>   samples: 2
#>   mode   : memory
#>   backend: GeneTrackR-native-R
#>   coordinate: 1-based closed
#>   schema : 2
#>   seqinfo: 2 sample/chromosome record(s)
#>   records: 57

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

rnaseq_gr <- as_granges(rnaseq_region)
rnaseq_gr
#> GRanges object with 57 ranges and 2 metadata columns:
#>        seqnames            ranges strand |    sample_id     score
#>           <Rle>         <IRanges>  <Rle> |  <character> <numeric>
#>    [1]     chr1 12340001-12340100      + | RNA_seq_plus     25.76
#>    [2]     chr1 12340101-12340200      + | RNA_seq_plus     28.00
#>    [3]     chr1 12340201-12340300      + | RNA_seq_plus     30.24
#>    [4]     chr1 12340301-12340400      + | RNA_seq_plus     26.88
#>    [5]     chr1 12340401-12340500      + | RNA_seq_plus     29.12
#>    ...      ...               ...    ... .          ...       ...
#>   [53]     chr1 12351501-12351600      + | RNA_seq_plus    24.472
#>   [54]     chr1 12351601-12351700      + | RNA_seq_plus    26.600
#>   [55]     chr1 12351701-12351800      + | RNA_seq_plus    28.728
#>   [56]     chr1 12351801-12351900      + | RNA_seq_plus    25.536
#>   [57]     chr1 12351901-12352000      + | RNA_seq_plus    27.664
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

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

rnaseq_gene_gr <- retrieve_bwg(
  rnaseq,
  annotation = gp,
  gene_id = "GeneA",
  strand = "auto",
  as = "GRanges"
)

rnaseq_gene_gr
#> GRanges object with 57 ranges and 2 metadata columns:
#>        seqnames            ranges strand |    sample_id     score
#>           <Rle>         <IRanges>  <Rle> |  <character> <numeric>
#>    [1]     chr1 12340001-12340100      + | RNA_seq_plus     25.76
#>    [2]     chr1 12340101-12340200      + | RNA_seq_plus     28.00
#>    [3]     chr1 12340201-12340300      + | RNA_seq_plus     30.24
#>    [4]     chr1 12340301-12340400      + | RNA_seq_plus     26.88
#>    [5]     chr1 12340401-12340500      + | RNA_seq_plus     29.12
#>    ...      ...               ...    ... .          ...       ...
#>   [53]     chr1 12351501-12351600      + | RNA_seq_plus    24.472
#>   [54]     chr1 12351601-12351700      + | RNA_seq_plus    26.600
#>   [55]     chr1 12351701-12351800      + | RNA_seq_plus    28.728
#>   [56]     chr1 12351801-12351900      + | RNA_seq_plus    25.536
#>   [57]     chr1 12351901-12352000      + | RNA_seq_plus    27.664
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

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.

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:

signal_all <- merge_bwg(rnaseq, riboseq)
summary_bwg(signal_all)
#>         sample_id sum_signal mean_signal max_signal covered_bases
#>            <char>      <num>       <num>      <num>         <int>
#> 1:   RNA_seq_plus   632305.6   13.715957     33.264         46100
#> 2:  RNA_seq_minus   361448.2   12.463731     28.512         29000
#> 3:  Ribo_seq_plus    83318.0    4.124858     21.000         20199
#> 4: Ribo_seq_minus    45395.0    3.977133     19.000         11414

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.

genea_binned <- bin_bwg(
  rnaseq_genea,
  bin_size = 50
)

head(genea_binned)
#>       sample_id  chrom   bin    start      end value
#>          <char> <char> <num>    <int>    <int> <num>
#> 1: RNA_seq_plus   chr1     1 12340001 12340100 25.76
#> 2: RNA_seq_plus   chr1     3 12340101 12340200 28.00
#> 3: RNA_seq_plus   chr1     5 12340201 12340300 30.24
#> 4: RNA_seq_plus   chr1     7 12340301 12340400 26.88
#> 5: RNA_seq_plus   chr1     9 12340401 12340500 29.12
#> 6: RNA_seq_plus   chr1    11 12340501 12340600 25.76

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.

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
#>        sample_id                                                        file
#>           <char>                                                      <char>
#> 1:  RNA_seq_plus  /tmp/RtmpCT2eho/GeneTrackR_demo_bigwig/RNA_seq_plus.bigwig
#> 2: RNA_seq_minus /tmp/RtmpCT2eho/GeneTrackR_demo_bigwig/RNA_seq_minus.bigwig
#>    format
#>    <char>
#> 1: bigwig
#> 2: 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.

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:

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.

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.

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:

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

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