binGenome() assigns exon bins to wrong genes due to table() alphabetical sort

[aromanowski]

Using the O tauri genome I mentioned in the previous issue, I noticed a plastidial gene' exons were incorrectly being assigned to a different nuclear chromosome:

OstapCp10 in ASpli bins (featuresb):
                   locus    seqnames start   end strand feature    event
OstapCp10:I001 OstapCp10 NC_008289_1  9773 13508      -       I        -
OstapCp10:E001 OstapCp10 NC_014426_2 13009 13125      +       E external
OstapCp10:E002 OstapCp10 NC_014426_2 13494 15890      +       E external

OstapCp10 in GTF (all attributes):
GTF columns containing 'OstapCp10': gene_id locus_tag 

         seqnames        type start   end strand   gene_id gene            transcript_id
50201 NC_008289.1        gene  8852 14039      - OstapCp10 atpB                         
50202 NC_008289.1  transcript  8852 14039      - OstapCp10 atpB unassigned_transcript_92
50203 NC_008289.1        exon 13509 14039      - OstapCp10 atpB unassigned_transcript_92
50204 NC_008289.1        exon  8852  9772      - OstapCp10 atpB unassigned_transcript_92
50205 NC_008289.1         CDS 13509 14039      - OstapCp10 atpB unassigned_transcript_92
50206 NC_008289.1         CDS  8855  9772      - OstapCp10 atpB unassigned_transcript_92
50207 NC_008289.1 start_codon 14037 14039      - OstapCp10 atpB unassigned_transcript_92
50208 NC_008289.1  stop_codon  8852  8854      - OstapCp10 atpB unassigned_transcript_92

OstapCp10 in TxDb (transcriptsBy):
GRanges object with 1 range and 2 metadata columns:
         seqnames     ranges strand |     tx_id                tx_name
            <Rle>  <IRanges>  <Rle> | <integer>            <character>
  [1] NC_008289.1 8852-14039      - |        72 unassigned_transcrip..
  -------
  seqinfo: 22 sequences from an unspecified genome; no seqlengths

ASpli featuresb locus for NC_008289 bins:
                                 locus    seqnames start   end feature
OstapCp10:I001               OstapCp10 NC_008289_1  9773 13508       I
OT_ostta20g00850:E001 OT_ostta20g00850 NC_008289_1  8852  9772       E
OT_ostta20g00850:E002 OT_ostta20g00850 NC_008289_1 13509 14039       E

I first thought makeTxDbFromGFF might be incorrectly linking exons to genes for organellar chromosomes. However, The TxDb is completely correct — both OstapCp10 and OT_ostta01g00090 have exactly the right exons on the right chromosomes. The bug is entirely within ASpli's .createGRangesExons.getExonBinNames() in binGenome_functions.R. table(names(aGRanges)) sorts alphabetically, but the GRanges object it's counting is in TxDb tx_id order. When a gene from an organellar chromosome (OstapCp10, tx_id 72) appears earlier in the TxDb than nuclear genes but lands at a different position alphabetically relative to them, the sorted names vector diverges from the ordered coordinate vector — every gene gets the next gene's exons.

In binGenome_functions.R, .createGRangesExons.getExonBinNames(), lines 136–144:

nExonsByGene  <- table( names( aGRanges ) )            # sorts alphabetically
exon.bins.num <- unlist( lapply( nExonsByGene, seq ) )
exonBinIds    <- sprintf('E%03d', exon.bins.num)
exonGeneNames <- rep( names(nExonsByGene), nExonsByGene)  # alphabetical order
exonBinNames  <- paste(exonGeneNames, exonBinIds, sep=":")
return( list( exonBinNames, exonGeneNames, exonBinIds) )

Here, table() returns its result sorted alphabetically by name. However, aGRanges (the unlisted disjoint exon GRanges) is in the order returned by exonsBy(txdb, by="gene") — TxDb internal order by tx_id — which does not match alphabetical order when genes from different chromosomes or annotation sources are present. The gene name vector exonGeneNames is therefore in a different order from the coordinate vector aGRanges, causing every gene to be assigned the coordinates of a different gene.

Environment
R 4.5.2 (2025-10-31 ucrt)
Bioconductor3.22
ASpli2.20.0

Proposed fix
Replace table() with rle() (run-length encoding), which counts consecutive identical names while preserving their original order:

# BEFORE (buggy):
nExonsByGene  <- table( names( aGRanges ) )
exon.bins.num <- unlist( lapply( nExonsByGene, seq ) )
exonBinIds    <- sprintf('E%03d', exon.bins.num)
exonGeneNames <- rep( names(nExonsByGene), nExonsByGene)
exonBinNames  <- paste(exonGeneNames, exonBinIds, sep=":")
return( list( exonBinNames, exonGeneNames, exonBinIds) )

# AFTER (fixed):
gene_rle      <- rle( names( aGRanges ) )
nExonsByGene  <- gene_rle$lengths
exonGeneNames <- rep( gene_rle$values, gene_rle$lengths )
exon.bins.num <- unlist( lapply( nExonsByGene, seq ) )
exonBinIds    <- sprintf( 'E%03d', exon.bins.num )
exonBinNames  <- paste( exonGeneNames, exonBinIds, sep = ":" )
return( list( exonBinNames, exonGeneNames, exonBinIds ) )

rle() encodes consecutive identical values without sorting, preserving the original order of names as they appear in aGRanges. This should ensure that gene names and exon coordinates remain aligned.

Impact
This bug affects every multi-exonic gene in every dataset. It is not specific to organellar genes — any dataset where the alphabetical order of gene IDs differs from the TxDb internal order will produce mis-labelled exon bins. The bug was identified by comparing exonsBy(txdb, by="gene") coordinates (correct) against featuresb(features) locus assignments (shifted by one position throughout the genome).

Note
The intron bins produced by .createGRangesIntrons() are not affected because that function uses names(introns.by.gene) directly via rep() without passing through table().

Alternative outside fix
make the TxDb's internal order match alphabetical order by rebuilding it from a GTF where gene groups are sorted alphabetically by gene_id. Then table() and exonsBy will always agree.

[aromanowski]

For the 'alternative outside fix', I also had to rename the 'Ostap' gene IDs to match the nuclear gene IDs prefix. That seemed to do the trick:

# Rename 'Ostap' Genes so that they match nuclear name prefixes
# We use across() to apply the rename to gene_id, transcript_id, and any 
# other associated ID columns so your TxDb builds cleanly without orphaned exons.
gtf_df_full <- gtf_df_full %>%
  dplyr::mutate(
    dplyr::across(
      dplyr::any_of(c("gene_id", "transcript_id", "protein_id", "locus_tag")),
      ~ stringr::str_replace_all(., 
                                 c("^Ostap" = "OT_ostta"))
    )
  )

# Sort: first alphabetically by gene_id, then by original row order within group
# to preserve gene → transcript → exon feature ordering
gtf_df_sorted <- gtf_df_full %>%
  dplyr::mutate(orig_order = dplyr::row_number()) %>%
  dplyr::arrange(gene_id, orig_order) %>%
  dplyr::select(-orig_order)

# Are they identical? If TRUE, sorting had no effect (already alphabetical)
cat("\nOrder identical to original:", 
    identical(unique(gtf_df_full$gene_id), 
              unique(gtf_df_sorted$gene_id)), "\n")

# Convert back to GRanges and export sorted GTF
gtf_sorted_gr <- makeGRangesFromDataFrame(
  gtf_df_sorted,
  keep.extra.columns = TRUE,
  seqnames.field     = "seqnames",
  start.field        = "start",
  end.field          = "end",
  strand.field       = "strand"
)

sorted_gtf_path <- file.path(dirname(params$gtf_file), "genome_alphabetical.gtf")
export(gtf_sorted_gr, sorted_gtf_path, format = "gtf")
cat("Sorted GTF written to:", sorted_gtf_path, "\n")

# Rebuild TxDb from alphabetically sorted GTF
txdb_sorted_cache <- file.path(dirname(params$gtf_file), "otauri_txdb_sorted.sqlite")

if (!file.exists(txdb_sorted_cache)) {
  message("Building alphabetically-ordered TxDb...")
  txdb_sorted <- txdbmaker::makeTxDbFromGFF(
    file       = sorted_gtf_path,
    format     = "gtf",
    organism   = "Ostreococcus tauri",
    taxonomyId = 70448L
  )
  saveDb(txdb_sorted, txdb_sorted_cache)
} else {
  txdb_sorted <- loadDb(txdb_sorted_cache)
}

After applying this fix I got from binGenome(txdb_sorted, logTo = features_sorted_log):

> features_sorted <- binGenome(txdb_sorted, logTo = features_sorted_log)

|=======================================================================================================================================================| 100%, Elapsed 01:40
binGenome log (sorted TxDb):
* Number of extracted Genes = 7833 
* Number of extracted Exon Bins = 3439 
* Number of extracted intron bins = 2046 
* Number of extracted trascripts = 7834 
* Number of extracted junctions = 2046 
* Number of AS bins (not include external) = 0
* Number of AS bins (include external) = 0
* Classified as: 
  ES bins = 0	(NaN%)
  IR bins = 0	(NaN%)
  Alt5'ss bins = 0	(NaN%)
  Alt3'ss bins = 0	(NaN%)
  Multiple AS bins = 0	(NaN%)
  classified as:
  		ES bins = 0	(NaN%)
  		IR bins = 0	(NaN%)
  		Alt5'ss bins = 0	(NaN%)
  		Alt3'ss bins = 0	(NaN%)

I verified the proper assignments:

OT_osttaCp10 bins (sorted TxDb):
                         locus    seqnames start   end strand feature
OT_osttaCp10:E001 OT_osttaCp10 NC_008289_1  8852  9772      -       E
OT_osttaCp10:I001 OT_osttaCp10 NC_008289_1  9773 13508      -       I
OT_osttaCp10:E002 OT_osttaCp10 NC_008289_1 13509 14039      -       E

OT_ostta01g00090 bins (sorted TxDb):
                                 locus    seqnames start   end strand feature
OT_ostta01g00090:E001 OT_ostta01g00090 NC_014426_2 13009 13125      +       E
OT_ostta01g00090:I001 OT_ostta01g00090 NC_014426_2 13126 13493      +       I
OT_ostta01g00090:E002 OT_ostta01g00090 NC_014426_2 13494 15890      +       E