fix: add hs1/hs1_chrR TxDb support to R annotation scripts (ccbr_annotate_bed.R, ccbr_annotate_peaks.R)

[kopardev]

Summary

ccbr_annotate_bed.R and ccbr_annotate_peaks.R do not handle hs1 or hs1_chrR genomes. Any rule that calls these scripts crashes with object 'tdb' not found because neither genome is present in the genome→TxDb dispatch chain.

Discovered during the Biowulf test run for PR #120.

---

Error

Error in is(TxDb, "GRanges") : object 'tdb' not found

Affected rules: atac_genrich_peakcalling, atac_annotate_fixed_width_consensus_peaks, annotate_roi
Scripts: workflow/scripts/ccbr_annotate_bed.R, workflow/scripts/ccbr_annotate_peaks.R

---

Root cause

Both scripts dispatch on args$genome to assign tdb (TxDb annotation database) and adb (OrgDb). The current if-chain covers: hg19, hg38, mm9, mm10, mmul10, bosTau9 — but not hs1 or hs1_chrR. When neither matches, tdb is never assigned and downstream usage crashes.

---

Modern viable solutions (in order of preference)

1. Use a real T2T TxDb — a pre-built one now exists ✅ (preferred)

The BiocT2T package provides a TxDb built from the actual CHM13v2.0 annotation. Install via:

BiocT2T::install_early_t2t_txdb()
library(TxDb.Hsapiens.NCBI.CHM13v2)

This gives correct T2T coordinates rather than an hg38 stand-in.

2. Build the TxDb from the NCBI GFF (officially documented method)

makeTxDbFromUCSC() does not support hs1, so the recommended route is the GFF file provided by NCBI. Download GCF_009914755.1_T2T-CHM13v2.0_genomic.gff.gz from NCBI, then:

library(rtracklayer); library(txdbmaker)
gff <- import("GCF_009914755.1_T2T-CHM13v2.0_genomic.gff.gz")
seqinfo(gff) <- Seqinfo(genome="T2T-CHM13v2.0")
txdb <- makeTxDbFromGRanges(gff, taxonomyId=9606)

This produces a proper TxDb with ~188,205 transcripts on the T2T-CHM13v2.0 genome. The resulting .sqlite file could be shipped as a cached resource in ASPEN's resource bundle so it isn't rebuilt each run.

3. BSgenome / sequence side

BSgenome.Hsapiens.NCBI.T2T.CHM13v2.0 wraps the T2T-CHM13v2.0 assembly (GCA_009914755.4) and is installable via BiocManager. There's also BSgenome.Hsapiens.UCSC.hs1 if UCSC-style seqnames are needed.

OrgDb (adb): org.Hs.eg.db is correct as-is

An OrgDb holds no positional data — it maps gene identifiers and pathways only. Reusing org.Hs.eg.db for hs1/hs1_chrR is correct; only the TxDb needs a real T2T source.

---

⚠️ Seqname convention caveat

T2T RefSeq/Liftoff annotations use NCBI-style seqnames (1, 2, ..., X, Y, MT) while ASPEN's hs1 index uses UCSC-style names (chr1, chr2, ..., chrX, chrY). Any TxDb built from the NCBI GFF must be harmonized before use, or annotation will silently produce no overlaps:

seqlevelsStyle(tdb) <- "UCSC"

This must be applied to tdb before passing it to annotatePeak().

---

Recommended fix for the PR

Add the hs1 / hs1_chrR case to both scripts, assigning:

• tdb = T2T TxDb (option 1 via BiocT2T, or option 2 pre-cached .sqlite), not hg38
• adb = "org.Hs.eg.db" (already correct)
• Apply seqlevelsStyle(tdb) <- "UCSC" after loading

---

Blocked by / Dependencies

• feat(ccbr_atacseq): add BiocT2T / TxDb.Hsapiens.NCBI.CHM13v2 for T2T annotation support (Dockerfile.v12) Dockers2#422 — build ccbr_atacseq:v12 with TxDb.Hsapiens.NCBI.CHM13v2

Related

• Genome added by issue feat: add hs1_chrR as a supported reference genome #117 / PR feat: add hs1_chrR (T2T-CHM13 + rDNA unit) as a supported genome #120
• Companion TSS ZeroDivisionError fixed in issue fix: ZeroDivisionError in _ccbr_counts2density.py when flanking TSS bins have zero reads #122 (commit 19dcbf5)

_{⚡ Generated using AI ⚡}

[kopardev]

Update: nciccbr/ccbr_atacseq:v13-feat is built and on Docker Hub — annotation approach revised

CCBR/Dockers2#422 is resolved. master_images/ccbr_atacseq/Dockerfile.v13 built successfully in CI (source@fb575c6) and is published as nciccbr/ccbr_atacseq:v13-feat. ASPEN should switch to this image.

A few things changed vs. the original plan in this issue's description — updating here so the R-script fix stays accurate:

BiocT2T does not exist

It was a hallucinated package name (404 on github.com/bioc/BiocT2T, absent from the full Bioconductor 3.23 package list). Do not reference BiocT2T::install_early_t2t_txdb() or TxDb.Hsapiens.NCBI.CHM13v2 — neither exists.

seqlevelsStyle(tdb) <- "UCSC" does not work for this genome either

T2T-CHM13v2.0's RefSeq accessions aren't covered by GenomeInfoDbData's rename tables, so this errors with "found no sequence renaming map...". Separately, even a manual renameSeqlevels() applied to an already-built TxDb does not survive AnnotationDbi::saveDb()/loadDb() — the rename has to happen on the source GRanges before txdbmaker::makeTxDbFromGRanges() is called, which only matters if we ever rebuild it, since:

The TxDb is now pre-built and shipped inside the image — no runtime build/rename needed

A TxDb was built from NCBI RefSeq GCF_009914755.1_T2T-CHM13v2.0_genomic.gff with seqlevels already renamed to UCSC style (chr1...chr22, chrX, chrY) and seqlengths set, then saved and gzip-compressed into the repo (master_images/ccbr_atacseq/resources/annotation/TxDb.Hsapiens.NCBI.T2T.CHM13v2.0.sqlite.gz, ~80 MB). The Dockerfile copies and gunzips it at build time, so it ships in the image at:

/opt2/annotation/TxDb.Hsapiens.NCBI.T2T.CHM13v2.0.sqlite

Loading it at runtime is just:

tdb <- AnnotationDbi::loadDb("/opt2/annotation/TxDb.Hsapiens.NCBI.T2T.CHM13v2.0.sqlite")
adb <- "org.Hs.eg.db"

No seqlevelsStyle<- call is needed — the seqnames are already chr1...chrY, matching ASPEN's hs1/hs1_chrR index naming. Verified end-to-end (fresh loadDb(), transcripts()/exons()/cds()/genes() counts, and ChIPseeker::annotatePeak() on test peaks all correct).

New requirement: same TxDb for both hs1 and hs1_chrR

Per discussion, ccbr_annotate_bed.R / ccbr_annotate_peaks.R should dispatch both hs1 and hs1_chrR to this same TxDb/OrgDb pair:

} else if (args$genome %in% c("hs1", "hs1_chrR")) {
    tdb <- AnnotationDbi::loadDb("/opt2/annotation/TxDb.Hsapiens.NCBI.T2T.CHM13v2.0.sqlite")
    adb <- "org.Hs.eg.db"

Caveat: the TxDb only has the 24 standard nuclear chromosomes (no chrM — it isn't in NCBI's T2T-CHM13v2.0 GFF3 either — and no chrR, the synthetic rDNA-model contig added in hs1_chrR). Peaks called on chrR (or chrM, if ever present) will simply get no gene/transcript overlap from annotatePeak() — same as any unannotated contig — which is expected and shouldn't need special-casing in the scripts.

Updated recommended fix for the PR

• Add one hs1/hs1_chrR branch (shared) to both scripts' genome dispatch, as above — no BiocT2T, no seqlevelsStyle<-, no runtime GFF download/build.
• adb = "org.Hs.eg.db" (unchanged, still correct).
• Base image bump to nciccbr/ccbr_atacseq:v13-feat (pending merge of the feature/atacseq_master_autobuild branch to dev/main in Dockers2).

Blocked by / Dependencies

• feat(ccbr_atacseq): add BiocT2T / TxDb.Hsapiens.NCBI.CHM13v2 for T2T annotation support (Dockerfile.v12) Dockers2#422 — done, see fb575c6 and nciccbr/ccbr_atacseq:v13-feat on Docker Hub.

_{⚡ Generated using AI ⚡}