Skip to content

feat: add hs1_chrR (T2T-CHM13 + rDNA unit) as a supported genome - #120

Open
kopardev wants to merge 23 commits into
mainfrom
issue_117
Open

feat: add hs1_chrR (T2T-CHM13 + rDNA unit) as a supported genome#120
kopardev wants to merge 23 commits into
mainfrom
issue_117

Conversation

@kopardev

@kopardev kopardev commented Aug 3, 2026

Copy link
Copy Markdown
Contributor

Changes

Add hs1_chrR as a new supported genome in ASPEN — the T2T-CHM13 assembly with the rDNA repeat unit (chrR) appended.

  • Add hs1_chrR.blacklist.fa.gz, hs1.promoters.bed.gz, frip beds, and tssBed tarballs via Git LFS
  • Add .gitattributes LFS tracking entries for the new resource files
  • Add hs1_chrR genome stanza to config/config.yaml with user-facing comment
  • Add hs1_chrR to SUPPORTED_GENOMES in the aspen CLI wrapper
  • Document hs1_chrR in docs/deployment.md and docs/limitations.md (with rDNA scientific rationale tip box)
  • Add resources/chroms/hs1_chrR.chroms listing chr1–chr22, chrX, chrY, chrR
  • Add CHANGELOG.md entry

Issues

Closes #117

PR Checklist

  • This comment contains a description of changes with justifications, with any relevant issues linked.
  • Update docs if there are any API changes.
  • Update CHANGELOG.md with a short description of any user-facing changes and reference the PR number. Guidelines: https://keepachangelog.com/en/1.1.0/
  • Test run completes successfully on biowulf.

…flanking bins are zero

When all 20 flanking TSS bins have zero reads (e.g. genes on chrR with no
ATAC background), flankavg=0 caused a ZeroDivisionError. Now emits 0.0 for
all density values and reports TSS enrichment as 0.0 rather than crashing.

Closes #122

⚡ generated using AI ⚡
@kopardev

kopardev commented Aug 4, 2026

Copy link
Copy Markdown
Contributor Author

Applied fix for issue #122 in commit 19dcbf5:

Bug fixed: ZeroDivisionError in workflow/scripts/_ccbr_counts2density.py (rule atac_tss) when all 20 flanking TSS bins have zero reads.

Change: Added a zero-guard after computing flankavg. If flankavg == 0, the script now emits 0.0 for all density values and reports # TSS enrichment: 0.0000 rather than crashing. This is the correct biological interpretation — a gene with no background flanking signal has undefined (effectively zero) TSS enrichment.

Issue #123 (missing hs1/hs1_chrR cases in the R annotation scripts) will be addressed in a follow-up commit on this same branch.

kopardev added 10 commits August 4, 2026 20:04
…iner to v13-feat

- Add hs1/hs1_chrR branch to ccbr_annotate_bed.R and ccbr_annotate_peaks.R
  genome dispatch chains; loads pre-built T2T-CHM13v2.0 TxDb from
  /opt2/annotation/TxDb.Hsapiens.NCBI.T2T.CHM13v2.0.sqlite via
  AnnotationDbi::loadDb(); adb remains org.Hs.eg.db
- Bump masterdocker from nciccbr/ccbr_atacseq:v11-feat to v13-feat
  (ships the pre-built TxDb sqlite; no seqlevelsStyle<- needed)
- Update --genome help strings in both scripts
- Update CHANGELOG

Closes #123

⚡ generated using AI ⚡
…pWriter bug

The monolithic ccbr_atacseq:v13-feat container ships a cutadapt version
that triggers 'AttributeError: _ThreadedGzipWriter object has no attribute
name' when running with -j >1 on gzip output. Pin the trim rule to the
dedicated nciccbr/ccbr_cutadapt_1.18:v3 container (cutadaptdocker), which
ships cutadapt 1.18 — predates the _ThreadedGzipWriter class entirely.

- Add cutadaptdocker key to config/config.yaml
- Change trim rule container from masterdocker to cutadaptdocker
- Update CHANGELOG

⚡ generated using AI ⚡
… cutadapt bug in v13-feat

ccbr_cutadapt_1.18:v3 lacks /opt2/TruSeq_and_nextera_adapters.consolidated.fa.
Use v11-feat for the trim rule instead — has working cutadapt and the
required adapter FASTA. v13-feat remains masterdocker for all other rules
(annotation TxDb support).

⚡ generated using AI ⚡
org.Bt.eg.db, org.Mmu.eg.db, TxDb.Btaurus.UCSC.bosTau9.refGene, and
TxDb.Mmulatta.UCSC.rheMac10.refGene were unconditionally loaded at the top
of ccbr_annotate_bed.R and ccbr_annotate_peaks.R. The v13-feat container
is missing org.Bt.eg.db, causing an immediate crash for all genomes
(including hs1_chrR) even though those packages are never needed.

Move TxDb loads into their respective dispatch blocks; OrgDb packages
(org.Bt.eg.db, org.Mmu.eg.db) are dropped from explicit library() calls
entirely — ChIPseeker loads them dynamically via annoDb=.

⚡ generated using AI ⚡
…cripts

tidyverse is absent from nciccbr/ccbr_atacseq:v13-feat. The three affected
scripts (fixed_width_peakSets_to_consensus_peakSet.R,
narrowPeak_normalize_pvalues.R, narrowPeak_to_fixed_width_peakSet.R) only
use arrange/mutate/filter, all of which are in dplyr — already installed.

⚡ generated using AI ⚡
When ChIPseeker's annotatePeak cannot map TxDb gene IDs through org.*.db
(e.g. T2T/hs1_chrR assemblies using a RefSeq-based SQLite TxDb), the
ENSEMBL, SYMBOL, and GENENAME columns are not added to the result data
frame. The subsequent column selection then crashes with 'undefined
columns selected'. Fill the missing columns with NA_character_ before
the merge so annotation still succeeds and emits NA for those fields.

⚡ generated using AI ⚡
The genome_is_known check in ccbr_atac_genrich_peak_calling.bash and
ccbr_atac_macs2_peak_calling.bash only recognised hg19/hg38/mm10, so
RUNCHIPSEEKER was forced to False for hs1/hs1_chrR. This caused the
per-replicate .annotated files never to be created, then the unconditional
mv in the FILTERPEAKS loop failed with 'No such file or directory'.

Both hs1 and hs1_chrR are now supported: nciccbr/ccbr_atacseq:v13-feat
ships a T2T CHM13v2 SQLite TxDb and the R annotation scripts have been
updated to dispatch to it.

⚡ generated using AI ⚡
Document the four fixes from this session that were committed but not
yet logged:
- org.Bt.eg.db on-demand loading (#124)
- tidyverse → dplyr in peak-width R scripts (#124)
- NA fallback for ENSEMBL/SYMBOL/GENENAME in annotation scripts (#123)
- hs1/hs1_chrR genome_is_known in peakcalling bash scripts (#117)

⚡ generated using AI ⚡
The following scripts are unused by any pipeline rule or helper script:
- picardcloud.jar: Picard is called as /opt2/picardcloud.jar inside the
  container; the local copy is never mounted or referenced.
- ccbr_atac_trim_align_pe.bash: superseded by separate trim + align rules.
- ccbr_bam_filter_by_mapq.py: no references anywhere in the codebase.
- atac_assign_multimappers.py: mentioned only in a comment in align.smk,
  never executed.
- script1.py: no references anywhere in the codebase.

⚡ generated using AI ⚡
The pieplot chunk in DESeq2.Rmd only built roiAnnodf for hg38/mm10,
so the diffatac rule crashed with "object 'roiAnnodf' not found" for
every peakcaller/method/contrast combination when genome is hs1 or
hs1_chrR. Add a matching branch that loads the T2T CHM13v2.0 TxDb,
mirroring the pattern already used in ccbr_annotate_bed.R/
ccbr_annotate_peaks.R (#123), and backfill ENSEMBL/SYMBOL/GENENAME
with NA when ChIPseeker can't map T2T gene IDs, since SYMBOL is
referenced later in the volcano plot chunk.

Fixes #125

_commit message is ai-generated_
…nge values

create_tn5bams failed at bedToBam with "Start Coordinate detected that
is < 0" for reads near chromosome boundaries, because Tn5 cutsite and
read-end shifts (+4/-5) were applied without bounds checking.

- ccbr_atac_bam2tn5bed.py: add parse_genome_sizes()/clamp_interval()/
  tn5_cutsite() helpers; require --genomefile to clamp every tn5 and
  read interval into [0, chrom_len]; skip unmapped/secondary/
  supplementary pairs and unknown chromosomes; emit per-read Tn5
  cutsite entries; report pairs_seen/written/skipped and
  clipped/skipped counters to stderr for QC.
- align.smk: pass --genomefile {params.genomefile} (already defined
  for the existing bedToBam calls) through to the script.

Fixes #124

_commit message is ai-generated_
… support

The diffatac/diffatac_aggregate rules ran in the baser container
(nciccbr/ccbr_baser:230531), which has all R packages DESeq2.Rmd needs
but lacks the T2T-CHM13v2.0 TxDb sqlite baked into masterdocker
(nciccbr/ccbr_atacseq:v13-feat). This caused hs1_chrR diffatac runs to
fail with "DB file '/opt2/annotation/TxDb.Hsapiens.NCBI.T2T.CHM13v2.0.sqlite'
not found" after the roiAnnodf fix in #125 started working.

masterdocker has been rebuilt with the previously-missing packages
(tidyverse, DT, pander, ggfortify, EnhancedVolcano, ComplexHeatmap; see
CCBR/Dockers2#426), so diffatac/diffatac_aggregate now use masterdocker
instead, and the now-unused baser container entry is removed from
config.yaml. Also drop two other container entries (fastqcdocker,
ucscdocker) confirmed unused by any rule.

Fixes #126

_commit message is ai-generated_
Use the latest ccbr_tools symlink (currently v0.7) instead of pinning
to v0.4, so the pipeline picks up ccbr_tools updates without requiring
a config change each time.

_commit message is ai-generated_
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

ASPEN RepoName

Projects

None yet

Development

Successfully merging this pull request may close these issues.

feat: add hs1_chrR as a supported reference genome

1 participant