Skip to content

Historical GRCh38 RefSeq alignments - #51 - #124

Merged
davmlaw merged 1 commit into
mainfrom
issue-51-historical-grch38
Aug 13, 2026
Merged

Historical GRCh38 RefSeq alignments - #51#124
davmlaw merged 1 commit into
mainfrom
issue-51-historical-grch38

Conversation

@davmlaw

@davmlaw davmlaw commented Aug 13, 2026

Copy link
Copy Markdown
Contributor

🤖 Written by Claude

Closes #51 (code side; the data lands with the next data release).

Adds NCBI's historical GRCh38 transcript alignments (RefSeq_historical_alignments, RS_2024_08 set from ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/historical/GRCh38/) as a refseq/GRCh38 source. Compared against the current cdot-0.2.33.refseq.GRCh38.json.gz this contributes 40,885 new transcript versions (39,566 NM_, 1,319 NR_), mostly old clinical-era versions (eg NM_000014.3, NM_000016.2) replaced before any annotation release cdot ingests.

Pipeline changes

  • cdot_transcripts.yaml: new source placed after UTA and before the official releases, so anything from an official annotation overwrites it (the ordering decided in Historical GRCh38 refseq #51 back in 2023).
  • Snakefile: generic multi-URL source support. A urls: list is downloaded as parts and concatenated in list order (annotation first, then alignments; concatenated gzip streams are valid gzip), and a per-source cdot_args: key passes extra flags to cdot_json.py (here --skip-missing-parents). Verified with a snakemake dry-run.

Why this stalled in 2023, and why it is fine now

  • The duplicated exons were the old shell script cat-ing the annotation file twice, not a parser bug.
  • The validation failures ("cDNA match starts at 3 not 1", exon sums shorter than the sequence) are genuine partial alignments in NCBI's data: eg NM_000016.3 is 2454 nt but only bases 1..2423 align (unaligned poly-A tail), and NM_002521.1's alignment starts at base 2. Since the start_codon/stop_codon: documented 1-based, and undefined across an inter-exon transcript gap #123 coordinate work the parser represents these faithfully; only positions inside an unaligned region cannot be projected.

Tests

New fixture refseq_test.historical_RS_2024_08.gff cut from the real NCBI files (both ACADM historical versions, gapped alignments on both strands, a partial-start case, and an alignment-only orphan) plus test_refseq_gff3_historical covering no-duplicate-exons, Gap passthrough, CDS+alignment combination, partial starts and skip_missing_parents behaviour. Full suite: 381 passed.

🤖 Generated with Claude Code

https://claude.ai/code/session_0152Ksi6d76CwsAooS7dBXVp

Add NCBI's RefSeq_historical_alignments (RS_2024_08 set) as a GRCh38
refseq source: 40,885 transcript versions not in the current data
(39,566 NM_, 1,319 NR_), mostly old clinical-era versions replaced
before any annotation release cdot ingests.

The source is two files (annotation + alignments) concatenated before
parsing, so the Snakefile gains generic multi-URL source support
(urls: list in cdot_transcripts.yaml) and per-source cdot_args
(here --skip-missing-parents, for alignments whose transcript has no
annotation entry). Merged after UTA and before all official releases,
so official annotation always wins.

The 2023 blockers are resolved: the duplicated exons came from the old
shell script cat-ing the annotation twice, and the short/offset cDNA
coordinates are genuine partial alignments (eg NM_000016.3 is 2454 nt
but only 1..2423 aligns; NM_002521.1 starts at base 2), which the
parser now represents faithfully since the #123 coordinate work.
@davmlaw davmlaw mentioned this pull request Aug 13, 2026
@davmlaw
davmlaw merged commit aa3cd01 into main Aug 13, 2026
5 checks passed
@davmlaw
davmlaw deleted the issue-51-historical-grch38 branch August 13, 2026 11:43
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Historical GRCh38 refseq

1 participant