Annotation I/O and conversion

Introduction

GeneTrackR provides one annotation workflow for GenePred, GenePredExt, GTF, GFF3, and BED files. The readers normalize file-level coordinate conventions into a common internal representation, and the downstream API works on compatible GenePred / FeatureTrack objects instead of requiring format-specific code at every step.

This article follows the normal analysis order:

read annotation files;
inspect and validate the returned objects;
retrieve genes, transcripts, exons, or genomic intervals;
convert between standardized representations;
merge compatible annotation subsets;
write native or cross-format outputs;
visualize annotation structure before downstream signal, variant, haplotype, or LD analysis.

library(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.

Supported formats and coordinate conventions

Format	Typical content	Reader	Main return	Input coordinate convention
GenePred	10-column transcript model	`read_genepred(..., format = "genePred")`	`GenePred`	0-based half-open by default
GenePredExt	15-column transcript model with gene/status/frame fields	`read_genepred(..., format = "genePredExt")`	`GenePred`	0-based half-open by default
GTF	9-column hierarchical annotation	`read_gtf()`	`FeatureTrack`	1-based closed
GFF3	9-column hierarchy using `ID` / `Parent`	`read_gff()`	`FeatureTrack`	1-based closed
BED	BED3-BED12 interval annotation	`read_bed()`	`FeatureTrack`	0-based half-open by default

GeneTrackR stores imported annotation intervals internally as 1-based closed coordinates. Coordinate conversion therefore belongs at the file boundary: do not manually add or subtract one after reading an annotation unless an external downstream format explicitly requires a different convention.

write_feature() converts GenePred- and BED-like outputs back to UCSC-style coordinates by default with coordinate = "ucsc".

Load the deterministic annotation files

The bundled demo files represent the same small reference genome in multiple annotation formats.

gp_file <- system.file(
  "extdata",
  "gtr_demo.genePredExt",
  package = "GeneTrackR"
)
gtf_file <- system.file(
  "extdata",
  "gtr_demo.gtf",
  package = "GeneTrackR"
)
gff_file <- system.file(
  "extdata",
  "gtr_demo.gff3",
  package = "GeneTrackR"
)
bed_file <- system.file(
  "extdata",
  "gtr_demo_features.bed",
  package = "GeneTrackR"
)

Read all four representations through their native readers:

gp <- read_genepred(
  gp_file,
  format = "genePredExt",
  verbose = FALSE,
  progress = FALSE
)

gtf <- read_gtf(
  gtf_file,
  verbose = FALSE,
  progress = FALSE
)

gff <- read_gff(
  gff_file,
  verbose = FALSE,
  progress = FALSE
)

features <- read_bed(
  bed_file,
  verbose = FALSE,
  progress = FALSE
)

GenePred/GenePredExt and hierarchical GTF/GFF3 inputs provide gene-model information. BED is primarily treated as an interval feature track and should not be assumed to reconstruct transcript/exon blocks from generic BED input.

Annotation object contract

GeneTrackR intentionally exposes a small common object model:

Input	Main object	Important contents
GenePred / GenePredExt	`GenePred`	`$genes`, `$transcripts`, `$exons`, `$data`
Gene-model GTF / GFF3	`FeatureTrack`	`$data` plus derived `$genes`, `$transcripts`, `$exons` when available
BED	`FeatureTrack`	primarily the standardized interval `$data` table

The objects can be inspected directly without immediately flattening them:

gp
#> <Feature>
#>   records    : 284
#>   genes      : 20
#>   transcripts: 24
#>   exons      : 76
#>   format     : genePredExt
#>   coordinate : 1-based closed
head(gp$genes)
#>    gene_id  chrom strand gene_start gene_end n_transcripts gene_type
#>     <char> <char> <char>      <int>    <int>         <int>    <char>
#> 1:   GeneD   chr1      +    1000001  1008000             1    coding
#> 2:   GeneE   chr1      +    3000001  3009000             1    coding
#> 3:   GeneF   chr1      -    5000001  5010000             1    coding
#> 4:   GeneG   chr1      +    7000001  7007000             1    coding
#> 5:   GeneH   chr1      -    9000001  9009000             1    coding
#> 6:   GeneA   chr1      +   12340001 12352000             2    coding
head(gp$transcripts)
#>    transcript_id gene_id  chrom strand tx_start   tx_end cds_start  cds_end
#>           <char>  <char> <char> <char>    <int>    <int>     <int>    <int>
#> 1:          TxD1   GeneD   chr1      +  1000001  1008000   1000201  1007800
#> 2:          TxE1   GeneE   chr1      +  3000001  3009000   3000201  3008800
#> 3:          TxF1   GeneF   chr1      -  5000001  5010000   5000201  5009800
#> 4:          TxG1   GeneG   chr1      +  7000001  7007000   7000201  7006800
#> 5:          TxH1   GeneH   chr1      -  9000001  9009000   9000201  9008800
#> 6:          TxA1   GeneA   chr1      + 12340001 12352000  12340501 12351500
#>    exon_count gene_type score cds_start_stat cds_end_stat  exon_frames
#>         <int>    <char> <num>         <char>       <char>       <char>
#> 1:          3    coding     0       complete     complete    -1,-1,-1,
#> 2:          3    coding     0       complete     complete    -1,-1,-1,
#> 3:          3    coding     0       complete     complete    -1,-1,-1,
#> 4:          3    coding     0       complete     complete    -1,-1,-1,
#> 5:          3    coding     0       complete     complete    -1,-1,-1,
#> 6:          4    coding     0       complete     complete -1,-1,-1,-1,
head(gp$exons)
#>    transcript_id gene_id  chrom strand exon_number exon_start exon_end
#>           <char>  <char> <char> <char>       <int>      <int>    <int>
#> 1:          TxD1   GeneD   chr1      +           1    1000001  1000900
#> 2:          TxD1   GeneD   chr1      +           2    1003551  1004450
#> 3:          TxD1   GeneD   chr1      +           3    1007001  1008000
#> 4:          TxE1   GeneE   chr1      +           1    3000001  3000900
#> 5:          TxE1   GeneE   chr1      +           2    3004051  3004950
#> 6:          TxE1   GeneE   chr1      +           3    3008001  3009000
#>    exon_frame
#>         <int>
#> 1:         -1
#> 2:         -1
#> 3:         -1
#> 4:         -1
#> 5:         -1
#> 6:         -1

head(gtf$data)
#>    feature_id   name  chrom   start     end           type      level score
#>        <char> <char> <char>   <int>   <int>         <char>     <char> <num>
#> 1:      GTF_4  GeneD   chr1 1000001 1000200 five_prime_UTR subfeature    NA
#> 2:      GTF_3  GeneD   chr1 1000001 1000900           exon       exon    NA
#> 3:      GeneD  GeneD   chr1 1000001 1008000           gene       gene    NA
#> 4:       TxD1  GeneD   chr1 1000001 1008000     transcript transcript    NA
#> 5:      GTF_5  GeneD   chr1 1000201 1000900            CDS subfeature    NA
#> 6:      GTF_6  GeneD   chr1 1003551 1004450           exon       exon    NA
#>    strand          source gene_id transcript_id parent_id      gene_type
#>    <char>          <char>  <char>        <char>    <char>         <char>
#> 1:      + GeneTrackR_demo   GeneD          TxD1      <NA> protein_coding
#> 2:      + GeneTrackR_demo   GeneD          TxD1      <NA> protein_coding
#> 3:      + GeneTrackR_demo   GeneD          <NA>      <NA> protein_coding
#> 4:      + GeneTrackR_demo   GeneD          TxD1      <NA> protein_coding
#> 5:      + GeneTrackR_demo   GeneD          TxD1      <NA> protein_coding
#> 6:      + GeneTrackR_demo   GeneD          TxD1      <NA> protein_coding
#>    exon_number  phase
#>          <int> <char>
#> 1:           1      .
#> 2:           1      .
#> 3:          NA      .
#> 4:          NA      .
#> 5:           1      0
#> 6:           2      .
#>                                                                                                                        attribute
#>                                                                                                                           <char>
#> 1:                     gene_id "GeneD"; transcript_id "TxD1"; gene_name "GeneD"; exon_number "1"; gene_biotype "protein_coding";
#> 2:                     gene_id "GeneD"; transcript_id "TxD1"; gene_name "GeneD"; exon_number "1"; gene_biotype "protein_coding";
#> 3:                                                            gene_id "GeneD"; gene_name "GeneD"; gene_biotype "protein_coding";
#> 4: gene_id "GeneD"; transcript_id "TxD1"; gene_name "GeneD"; gene_biotype "protein_coding"; transcript_biotype "protein_coding";
#> 5:                     gene_id "GeneD"; transcript_id "TxD1"; gene_name "GeneD"; exon_number "1"; gene_biotype "protein_coding";
#> 6:                     gene_id "GeneD"; transcript_id "TxD1"; gene_name "GeneD"; exon_number "2"; gene_biotype "protein_coding";
head(features$data)
#>    feature_id                                    name  chrom    start      end
#>        <char>                                  <char> <char>    <int>    <int>
#> 1:      BED_4                      SeedWeight_QTL|QTL   chr1 12330001 12380000
#> 2:      BED_2                 GeneA_enhancer|enhancer   chr1 12338201 12338750
#> 3:      BED_1                 GeneA_promoter|promoter   chr1 12339001 12340000
#> 4:      BED_3 GeneA_candidate_region|candidate_region   chr1 12339701 12352000
#> 5:      BED_5                     GeneA_repeat|repeat   chr1 12343201 12343800
#> 6:      BED_6        GeneA_conserved|conserved_region   chr1 12350401 12350900
#>      type   level score strand source gene_id transcript_id parent_id gene_type
#>    <char>  <char> <num> <char> <char>  <char>        <char>    <char>    <char>
#> 1:    BED feature     0      *    BED    <NA>          <NA>      <NA>      <NA>
#> 2:    BED feature     0      *    BED    <NA>          <NA>      <NA>      <NA>
#> 3:    BED feature     0      +    BED    <NA>          <NA>      <NA>      <NA>
#> 4:    BED feature     0      *    BED    <NA>          <NA>      <NA>      <NA>
#> 5:    BED feature     0      *    BED    <NA>          <NA>      <NA>      <NA>
#> 6:    BED feature     0      *    BED    <NA>          <NA>      <NA>      <NA>
#>    exon_number  phase attribute
#>          <int> <char>    <char>
#> 1:          NA   <NA>      <NA>
#> 2:          NA   <NA>      <NA>
#> 3:          NA   <NA>      <NA>
#> 4:          NA   <NA>      <NA>
#> 5:          NA   <NA>      <NA>
#> 6:          NA   <NA>      <NA>

This object contract is important downstream. Gene-model objects can be passed directly to gene-model, signal, variant, haplotype, and LD visualization functions, while interval-style FeatureTrack objects can be overlaid as annotation tracks.

Validate and summarize annotations

Use validate_feature() after importing external annotation data, especially before cross-format conversion or merging.

gp_validation <- validate_feature(gp)
gtf_validation <- validate_feature(gtf)
gff_validation <- validate_feature(gff)
bed_validation <- validate_feature(features)

gp_validation$invalid_summary
#> Empty data.table (0 rows and 2 cols): reason,N
bed_validation$invalid_summary
#> Empty data.table (0 rows and 2 cols): reason,N

For gene-model objects, validation checks both the standardized feature table and transcript/exon consistency. For interval-style BED tracks, validation focuses on the standardized feature records.

summary_feature() provides a compact QC view at different hierarchy levels:

summary_feature(gp, level = "gene")
#>     chrom strand n_genes n_coding_genes n_noncoding_genes median_gene_length
#>    <char> <char>   <int>          <int>             <int>              <num>
#> 1:   chr1      +       6              5                 1               7500
#> 2:   chr1      -       4              4                 0               9500
#> 3:   chr2      +       6              6                 0               8000
#> 4:   chr2      -       4              3                 1               9000
summary_feature(gp, level = "transcript")
#>     chrom strand n_transcripts n_coding_transcripts n_noncoding_transcripts
#>    <char> <char>         <int>                <int>                   <int>
#> 1:   chr1      +             7                    6                       1
#> 2:   chr1      -             5                    5                       0
#> 3:   chr2      +             7                    7                       0
#> 4:   chr2      -             5                    4                       1
#>    median_transcript_length median_exon_count
#>                       <num>             <num>
#> 1:                     8000                 3
#> 2:                    10000                 3
#> 3:                     8000                 3
#> 4:                     9000                 3
summary_feature(gp, level = "exon")
#>     chrom strand n_exons median_exon_length
#>    <char> <char>   <int>              <num>
#> 1:   chr1      +      23               1000
#> 2:   chr1      -      17               1000
#> 3:   chr2      +      19               1000
#> 4:   chr2      -      17               1000
summary_feature(features, level = "feature")
#>     chrom   type n_features median_feature_length
#>    <char> <char>      <int>                 <num>
#> 1:   chr1    BED          9                  1000
#> 2:   chr2    BED          6                  6500

A region can be summarized without first constructing a separate subset:

summary_feature(
  gp,
  chrom = "chr1",
  start = 12339001,
  end = 12374500,
  level = "gene"
)
#>     chrom strand n_genes n_coding_genes n_noncoding_genes median_gene_length
#>    <char> <char>   <int>          <int>             <int>              <num>
#> 1:   chr1      +       2              2                 0               8750
#> 2:   chr1      -       1              1                 0              10500

Retrieve genes, transcripts, and regions

retrieve_feature() is the common annotation query interface. By default it preserves a compatible annotation object; request as = "data.table" only when a flat table is required.

Exact gene or transcript retrieval

gp_gene_a <- retrieve_feature(
  gp,
  gene_id = "GeneA"
)

gp_tx_a1 <- retrieve_feature(
  gp,
  transcript_id = "TxA1"
)

gp_gene_a$genes
#>    gene_id  chrom strand gene_start gene_end n_transcripts gene_type
#>     <char> <char> <char>      <int>    <int>         <int>    <char>
#> 1:   GeneA   chr1      +   12340001 12352000             2    coding
gp_tx_a1$transcripts
#>    transcript_id gene_id  chrom strand tx_start   tx_end cds_start  cds_end
#>           <char>  <char> <char> <char>    <int>    <int>     <int>    <int>
#> 1:          TxA1   GeneA   chr1      + 12340001 12352000  12340501 12351500
#>    exon_count gene_type score cds_start_stat cds_end_stat  exon_frames
#>         <int>    <char> <num>         <char>       <char>       <char>
#> 1:          4    coding     0       complete     complete -1,-1,-1,-1,

Genomic-region retrieval

gp_region <- retrieve_feature(
  gp,
  chrom = "chr1",
  start = 12339001,
  end = 12374500,
  mode = "overlap"
)

The region-selection modes have distinct semantics:

`mode`	Meaning
`"overlap"`	Keep records overlapping the requested interval.
`"within"`	Keep records fully contained in the requested interval.
`"trim"`	Keep overlapping hierarchical gene-model records and clip model boundaries to the requested interval.

For generic interval-style BED tracks, trim behaves as overlap selection because there is no transcript/exon hierarchy to rebuild.

Retrieve a flat table

gp_gene_table <- retrieve_feature(
  gp,
  chrom = "chr1",
  start = 12339001,
  end = 12374500,
  level = "gene",
  as = "data.table"
)

gp_gene_table
#>    gene_id  chrom strand gene_start gene_end n_transcripts gene_type
#>     <char> <char> <char>      <int>    <int>         <int>    <char>
#> 1:   GeneA   chr1      +   12340001 12352000             2    coding
#> 2:   GeneB   chr1      -   12356001 12366500             2    coding
#> 3:   GeneC   chr1      +   12369001 12374500             1    coding

BED-derived features can be retrieved with the same regional interface or by text pattern:

bed_gene_a <- retrieve_feature(
  features,
  chrom = "chr1",
  start = 12338201,
  end = 12352000,
  mode = "overlap"
)

bed_promoters <- retrieve_feature(
  features,
  pattern = "promoter"
)

head(bed_promoters$data)
#>    feature_id                    name  chrom    start      end   type   level
#>        <char>                  <char> <char>    <int>    <int> <char>  <char>
#> 1:      BED_1 GeneA_promoter|promoter   chr1 12339001 12340000    BED feature
#> 2:      BED_7 GeneB_promoter|promoter   chr1 12366501 12367500    BED feature
#> 3:     BED_12 GeneL_promoter|promoter   chr2  8518001  8519300    BED feature
#> 4:     BED_15 GeneT_promoter|promoter   chr2 16994001 16995000    BED feature
#>    score strand source gene_id transcript_id parent_id gene_type exon_number
#>    <num> <char> <char>  <char>        <char>    <char>    <char>       <int>
#> 1:     0      +    BED    <NA>          <NA>      <NA>      <NA>          NA
#> 2:     0      -    BED    <NA>          <NA>      <NA>      <NA>          NA
#> 3:     0      -    BED    <NA>          <NA>      <NA>      <NA>          NA
#> 4:     0      +    BED    <NA>          <NA>      <NA>      <NA>          NA
#>     phase attribute
#>    <char>    <char>
#> 1:   <NA>      <NA>
#> 2:   <NA>      <NA>
#> 3:   <NA>      <NA>
#> 4:   <NA>      <NA>

Convert annotations and expose standardized tables

The unified object model allows annotation formats to be converted without reparsing their original text representation.

# Hierarchical GTF/GFF3 -> GenePred-compatible model
gtf_as_gp <- as_genepred(gtf)
gff_as_gp <- as_genepred(gff)

# Standardized hierarchy tables
gene_dt <- as_gene_table(gp)
tx_dt <- as_transcript_table(gp)
exon_dt <- as_exon_table(gp)
feature_dt <- as_feature_table(gp)

head(gene_dt)
#>    gene_id  chrom strand gene_start gene_end n_transcripts gene_type
#>     <char> <char> <char>      <int>    <int>         <int>    <char>
#> 1:   GeneD   chr1      +    1000001  1008000             1    coding
#> 2:   GeneE   chr1      +    3000001  3009000             1    coding
#> 3:   GeneF   chr1      -    5000001  5010000             1    coding
#> 4:   GeneG   chr1      +    7000001  7007000             1    coding
#> 5:   GeneH   chr1      -    9000001  9009000             1    coding
#> 6:   GeneA   chr1      +   12340001 12352000             2    coding
head(tx_dt)
#>    transcript_id gene_id  chrom strand tx_start   tx_end cds_start  cds_end
#>           <char>  <char> <char> <char>    <int>    <int>     <int>    <int>
#> 1:          TxD1   GeneD   chr1      +  1000001  1008000   1000201  1007800
#> 2:          TxE1   GeneE   chr1      +  3000001  3009000   3000201  3008800
#> 3:          TxF1   GeneF   chr1      -  5000001  5010000   5000201  5009800
#> 4:          TxG1   GeneG   chr1      +  7000001  7007000   7000201  7006800
#> 5:          TxH1   GeneH   chr1      -  9000001  9009000   9000201  9008800
#> 6:          TxA1   GeneA   chr1      + 12340001 12352000  12340501 12351500
#>    exon_count gene_type score cds_start_stat cds_end_stat  exon_frames
#>         <int>    <char> <num>         <char>       <char>       <char>
#> 1:          3    coding     0       complete     complete    -1,-1,-1,
#> 2:          3    coding     0       complete     complete    -1,-1,-1,
#> 3:          3    coding     0       complete     complete    -1,-1,-1,
#> 4:          3    coding     0       complete     complete    -1,-1,-1,
#> 5:          3    coding     0       complete     complete    -1,-1,-1,
#> 6:          4    coding     0       complete     complete -1,-1,-1,-1,

as_feature() returns the common feature-compatible object representation, while as_granges() provides the shared Bioconductor exchange boundary used by annotation, signal, and variant tracks:

gp_feature <- as_feature(gp)
gene_gr <- as_granges(gp, level = "gene")

gp_feature
#> <Feature>
#>   records    : 284
#>   genes      : 20
#>   transcripts: 24
#>   exons      : 76
#>   format     : genePredExt
#>   coordinate : 1-based closed
gene_gr
#> GRanges object with 20 ranges and 2 metadata columns:
#>        seqnames            ranges strand |     gene_id   gene_type
#>           <Rle>         <IRanges>  <Rle> | <character> <character>
#>    [1]     chr1   1000001-1008000      + |       GeneD      coding
#>    [2]     chr1   3000001-3009000      + |       GeneE      coding
#>    [3]     chr1   5000001-5010000      - |       GeneF      coding
#>    [4]     chr1   7000001-7007000      + |       GeneG      coding
#>    [5]     chr1   9000001-9009000      - |       GeneH      coding
#>    ...      ...               ...    ... .         ...         ...
#>   [16]     chr2   8524001-8530000      + |       GeneM      coding
#>   [17]     chr2 10000001-10009000      - |       GeneQ      coding
#>   [18]     chr2 12000001-12008000      + |       GeneR      coding
#>   [19]     chr2 14000001-14006000      - |       GeneS  non-coding
#>   [20]     chr2 16995001-17006000      + |       GeneT      coding
#>   -------
#>   seqinfo: 2 sequences from an unspecified genome; no seqlengths

For annotation objects, level selects gene, transcript, exon, or feature intervals. In-memory BwgTrack and VariantTrack objects can be converted with the same as_granges() function when downstream analysis is performed with standard Bioconductor interval infrastructure.

Use the hierarchy-specific table helpers when downstream code explicitly requires tabular data. Preserve the full object for GeneTrackR workflows so gene/transcript/exon relationships are not silently discarded.

Merge compatible annotations

merge_feature() accepts Feature, FeatureTrack, and GenePred objects, including lists of objects.

gp_gene_b <- retrieve_feature(
  gp,
  gene_id = "GeneB"
)

gp_merged <- merge_feature(
  gp_gene_a,
  gp_gene_b,
  source_names = c("GeneA", "GeneB"),
  conflict = "deduplicate"
)

gp_merged$genes
#>    gene_id  chrom strand gene_start gene_end n_transcripts gene_type
#>     <char> <char> <char>      <int>    <int>         <int>    <char>
#> 1:   GeneA   chr1      +   12340001 12352000             2    coding
#> 2:   GeneB   chr1      -   12356001 12366500             2    coding
#>    track_source
#>          <char>
#> 1:        GeneA
#> 2:        GeneB

Duplicate identifiers are handled explicitly:

`conflict`	Behavior
`"deduplicate"`	Keep the first compatible copy according to input order. Default.
`"rename"`	Rename conflicting IDs in later inputs and update hierarchy references.
`"keep_all"`	Retain conflicting records unchanged.
`"error"`	Stop if duplicate identifiers are detected.
`"keep_first"`	Backward-compatible alias of `"deduplicate"`.

When two formats describe the same locus but both representations must be retained, rename prevents identifier collisions:

gtf_gene_a <- retrieve_feature(
  gtf,
  gene_id = "GeneA"
)

cross_format <- merge_feature(
  list(gp_gene_a, gtf_gene_a),
  source_names = c("GenePred", "GTF"),
  conflict = "rename"
)
#> Warning: Duplicated annotation identifiers were detected across input objects:
#> gene_id=1 [GeneA]; transcript_id=2 [TxA1, TxA2]; feature_id=27
#> [TxA1:5UTR:12340001:12340500, TxA2:5UTR:12340001:12340500, TxA2.1, TxA1.1,
#> GeneA, ...]. Applying `conflict = "rename"`. With `conflict = "deduplicate"`,
#> input order determines which record is retained.

cross_format$genes
#>      gene_id  chrom strand gene_start gene_end n_transcripts gene_type
#>       <char> <char> <char>      <int>    <int>         <int>    <char>
#> 1: GTF_GeneA   chr1      +   12340001 12352000             2    coding
#> 2:     GeneA   chr1      +   12340001 12352000             2    coding
#>    track_source
#>          <char>
#> 1:          GTF
#> 2:     GenePred

If the inputs are alternative encodings of the same annotation and only one copy is required, use conflict = "deduplicate" instead.

Write annotations and perform round-trip checks

write_feature() is the main annotation writer. It supports genepred, genepredext, gff, gtf, bed6, and bed12 outputs.

GenePred / GenePredExt

gp_out <- file.path(tempdir(), "GeneA.genePredExt")

write_feature(
  gp_gene_a,
  gp_out,
  format = "genepredext",
  overwrite = TRUE
)

write_genepred() remains available as a compatibility wrapper, but new workflows should generally use write_feature() because it provides the unified cross-format interface.

GTF / GFF3

gtf_out <- file.path(tempdir(), "GeneA.gtf")
gff_out <- file.path(tempdir(), "GeneA.gff3")

write_feature(
  gtf_gene_a,
  gtf_out,
  format = "gtf",
  overwrite = TRUE
)

write_feature(
  gtf_gene_a,
  gff_out,
  format = "gff",
  overwrite = TRUE
)

BED6 / BED12

Generic BED-derived intervals should normally be written as BED6:

bed_out <- file.path(tempdir(), "GeneA_features.bed6")

write_feature(
  bed_gene_a,
  bed_out,
  format = "bed6",
  overwrite = TRUE
)

BED12 requires transcript/exon model information and is therefore appropriate for GenePred or hierarchical GTF/GFF3 models:

bed12_out <- file.path(tempdir(), "GeneA.bed12")

write_feature(
  gp_gene_a,
  bed12_out,
  format = "bed12",
  overwrite = TRUE
)

For reproducible conversion workflows, read the written file back and compare the biological records rather than comparing raw coordinates across formats that use different native conventions:

gp_roundtrip <- read_genepred(
  gp_out,
  format = "genePredExt",
  verbose = FALSE,
  progress = FALSE
)

as_gene_table(gp_gene_a)[, .(gene_id, chrom, gene_start, gene_end, strand)]
#>    gene_id  chrom gene_start gene_end strand
#>     <char> <char>      <int>    <int> <char>
#> 1:   GeneA   chr1   12340001 12352000      +
as_gene_table(gp_roundtrip)[, .(gene_id, chrom, gene_start, gene_end, strand)]
#>    gene_id  chrom gene_start gene_end strand
#>     <char> <char>      <int>    <int> <char>
#> 1:   GeneA   chr1   12340001 12352000      +

Visualize annotation structure

Annotation QC is easier when the tables and genomic structure are inspected together. GeneTrackR provides complementary views for gene models and interval feature tracks.

Gene model

plot_gene(
  gp,
  gene_id = "GeneA",
  label_position = "axis",
  direction_mode = "end",
  plot_theme = "classic",
  text_size = 10
)

Transcript structure

plot_transcript(
  gp,
  transcript_id = "TxA1",
  coordinate = "genomic",
  plot_theme = "classic",
  text_size = 10
)

Spliced transcript coordinates can be requested with coordinate = "transcript" when introns should be removed from the x-axis representation.

Genomic region

plot_region(
  gp,
  chrom = "chr1",
  start = 12339001,
  end = 12374500,
  mode = "overlap",
  label_by = "gene",
  plot_theme = "classic",
  text_size = 10
)

Interval feature track

plot_feature_track(
  features,
  chrom = "chr1",
  start = 12338201,
  end = 12374500,
  plot_theme = "classic",
  text_size = 10
)

Gene-length distribution

The distribution helper separates calculation from plotting so the derived table can be inspected or exported independently.

gene_length_dt <- get_gene_length_distribution_table(gp)
head(gene_length_dt)
#>    feature   unit feature_id gene_id transcript_id  chrom strand gene_type
#>     <fctr> <char>     <char>  <char>        <char> <char> <char>    <char>
#> 1:    gene   gene      GeneD   GeneD          <NA>   chr1      +    coding
#> 2:    gene   gene      GeneE   GeneE          <NA>   chr1      +    coding
#> 3:    gene   gene      GeneF   GeneF          <NA>   chr1      -    coding
#> 4:    gene   gene      GeneG   GeneG          <NA>   chr1      +    coding
#> 5:    gene   gene      GeneH   GeneH          <NA>   chr1      -    coding
#> 6:    gene   gene      GeneA   GeneA          <NA>   chr1      +    coding
#>       start      end length_bp
#>       <int>    <int>     <num>
#> 1:  1000001  1008000      8000
#> 2:  3000001  3009000      9000
#> 3:  5000001  5010000     10000
#> 4:  7000001  7007000      7000
#> 5:  9000001  9009000      9000
#> 6: 12340001 12352000     12000

plot_gene_length_distribution(gp)

Recommended annotation workflow

For external annotation projects, a robust sequence is:

Read with the native parser and let GeneTrackR handle source-format coordinates.
Validate immediately with validate_feature() before downstream filtering or conversion.
Inspect hierarchy and summary tables to confirm gene/transcript/exon counts and chromosome/strand coverage.
Retrieve sub-objects, not only plain tables, when the result will continue through GeneTrackR.
Convert only at explicit interfaces using as_*() or write_feature() rather than manually reconstructing columns.
Resolve merge conflicts explicitly instead of silently dropping duplicate identifiers.
Round-trip important exports and compare standardized biological coordinates/IDs.
Visualize representative loci before integrating signal, variants, haplotypes, or LD.

The annotation object produced here is the common structural input for the remaining GeneTrackR workflows. The dedicated visualization, signal, variant, haplotype, and association articles build on the same object contract rather than introducing a second annotation representation.

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       compiler_4.6.1      
#>  [4] GenomicRanges_1.65.4 jquerylib_0.1.4      IRanges_2.47.5      
#>  [7] Seqinfo_1.3.2        scales_1.4.0         yaml_2.3.12         
#> [10] fastmap_1.2.0        ggplot2_4.0.3        R6_2.6.1            
#> [13] labeling_0.4.3       generics_0.1.4       knitr_1.52          
#> [16] BiocGenerics_0.59.12 maketools_1.3.2      bslib_0.12.0        
#> [19] RColorBrewer_1.1-3   rlang_1.3.0          cachem_1.1.0        
#> [22] xfun_0.61            sass_0.4.10          sys_3.4.3           
#> [25] S7_0.2.2             otel_0.2.0           cli_3.6.6           
#> [28] withr_3.0.3          digest_0.6.39        grid_4.6.1          
#> [31] lifecycle_1.0.5      S4Vectors_0.51.10    vctrs_0.7.3         
#> [34] evaluate_1.0.5       glue_1.8.1           data.table_1.18.6.1 
#> [37] farver_2.1.2         buildtools_1.0.0     stats4_4.6.1        
#> [40] tools_4.6.1          htmltools_0.5.9