55e11ee5939bfbe8b3f7283b1abb337de5a55ae3 max Tue Sep 15 05:48:03 2026 -0700 hgTracks: the exon and codon mouseovers disagreed on out-of-frame transcripts, refs #38353 An exon's mouseover worked its codon numbers out as (c+2)/3 and the per-codon mouseover took them from the codon list, which builds them from exonFrames. Where a transcript's annotated CDS does not begin on a codon boundary the two part company: baseColorCodonsFromGenePred makes codon 1 short rather than shifting every number after it, so the arithmetic is off by the missing bases for the rest of the transcript. About 9,000 transcripts in hg38's wgEncodeGencodeCompV48 and 101 in ncbiRefSeqCurated are affected. hg38.knownGene has no exonFrames and never was. The exon note now reads the numbers off the codon list when there is one, and otherwise - zoomed out past the level that builds it - counts them with the same short first codon in mind. The per-codon note measures its c. range off the codon's own bases instead of taking it as 3*codonIndex-2..3*codonIndex, which was wrong for the same reason, and gathers both pieces of a codon split across an intron. Checked on hg38 ncbiRefSeqCurated NM_001372106.1, whose first coding exon has frame 1: exon 3 now reads c.369-602 (p.124-201) and the first per-codon box inside it c.369-371 (p.124), where before the exon said p.123 and the codon box said c.370-372. Also in this file, two things the review turned up while reading it: loadTxCdsBatch() capped its IN list at 2000 accessions and dropped the rest without saying so, which would leave the transcripts past the cap numbered along the genome with nothing to show for it. It queries in chunks now. txCodonIndexForCodon()'s comment claimed any of a codon's three bases gives the same answer. It does not when there is an indel inside the codon, which is the case the whole feature exists for; the code already uses the 5'-most base, so only the comment was wrong. genbankCds txCds is zeroed before baseColorTxAliForGenePred() may leave it untouched. diff --git src/hg/hgTracks/cds.c src/hg/hgTracks/cds.c index ab0a63ea221..1d165cde550 100644 --- src/hg/hgTracks/cds.c +++ src/hg/hgTracks/cds.c @@ -860,69 +860,80 @@ * on a page, which is not worth paying on every codon-level view. */ { struct hash *pslHash; /* accession -> struct psl list */ struct hash *cdsHash; /* accession -> struct genbankCds in transcript coords */ }; static void loadTxCdsBatch(char *db, char *table, struct hash *pslHash, struct hash *cdsHash) /* Fill cdsHash with the CDS, in its own transcript coordinates, of every accession in * pslHash. This has to come from the transcript's own annotation and not from the * alignment: for ~800 transcripts on hg38 the alignment does not even reach the start of * the CDS, and anchoring on the alignment would quietly renumber from the wrong base. */ { struct hashEl *el, *elList = hashElListHash(pslHash); if (elList == NULL) return; +struct sqlConnection *conn = hAllocConn(db); +/* Where the CDS comes from, decided once for the whole list. */ +boolean fromNcbiRefSeq = (sameString(table, "ncbiRefSeqPsl") && hTableExists(db, "ncbiRefSeqCds")); +/* refGene's versionless accessions get their CDS from the genbank tables, which usually + * live in hgFixed, so these names arrive database-qualified and only sqlTableExists can + * see them; hTableExists looks inside db and would always say no. */ +boolean fromGenbank = (sameString(table, "refSeqAli") && + sqlTableExists(conn, gbCdnaInfoTable) && sqlTableExists(conn, cdsTable)); +if (fromNcbiRefSeq || fromGenbank) + { + /* In chunks, rather than one query with every accession in it: the window can hold + * more transcripts than one IN list should carry, and a cap that silently dropped the + * rest would leave the transcripts past it numbered along the genome with no sign + * that anything had been left out. */ + el = elList; + while (el != NULL) + { struct dyString *accs = dyStringNew(1024); -int n = 0; -for (el = elList; el != NULL && n < 2000; el = el->next, n++) + int n; + for (n = 0; el != NULL && n < 2000; el = el->next, n++) { if (n > 0) sqlDyStringPrintf(accs, ","); sqlDyStringPrintf(accs, "'%s'", el->name); } -struct sqlConnection *conn = hAllocConn(db); struct dyString *query = NULL; -if (sameString(table, "ncbiRefSeqPsl") && hTableExists(db, "ncbiRefSeqCds")) + if (fromNcbiRefSeq) query = sqlDyStringCreate("select id, cds from ncbiRefSeqCds where id in (%-s)", accs->string); -/* refGene's versionless accessions get their CDS from the genbank tables, which usually - * live in hgFixed, so these names arrive database-qualified and only sqlTableExists can - * see them; hTableExists looks inside db and would always say no. */ -else if (sameString(table, "refSeqAli") && - sqlTableExists(conn, gbCdnaInfoTable) && sqlTableExists(conn, cdsTable)) + else query = sqlDyStringCreate( "select g.acc, c.name from %s g, %s c where g.cds = c.id and g.acc in (%-s)", gbCdnaInfoTable, cdsTable, accs->string); -if (query != NULL) - { struct sqlResult *sr = sqlGetResult(conn, query->string); char **row; while ((row = sqlNextRow(sr)) != NULL) { struct genbankCds *cds; AllocVar(cds); if (genbankCdsParse(row[1], cds) && cds->start < cds->end) hashAdd(cdsHash, row[0], cds); else freez(&cds); } sqlFreeResult(&sr); dyStringFree(&query); + dyStringFree(&accs); + } } hFreeConn(&conn); -dyStringFree(&accs); hashElFreeList(&elList); } static struct txAliWindow *txAliInWindow(char *db, char *table, char *chrom) /* Return the alignments overlapping this window and their transcripts' CDS. Two queries * per table per window, the first on the bin index, and the window is never wide: the * codon coloring this feeds only happens at zoomedToCdsColorLevel. */ { static struct hash *windowHash = NULL; if (windowHash == NULL) windowHash = hashNew(0); char key[1024]; safef(key, sizeof(key), "%s:%s:%s:%d:%d", db, table, chrom, winStart, winEnd); struct txAliWindow *tw = hashFindVal(windowHash, key); if (tw == NULL) @@ -1004,32 +1015,35 @@ lo = mid + 1; else { int qOff = txAli->qStarts[mid] + (t - txAli->tStarts[mid]); if (txAli->strand[0] == '-') qOff = txAli->qSize - 1 - qOff; return qOff; } } return -1; } static int txCodonIndexForCodon(struct psl *txAli, struct genbankCds *txCds, int start, int end, boolean posStrand) /* Return the codon's 1-based number counted in the transcript's own coordinates, or 0 if - * that cannot be worked out. Any of a codon's three bases gives the same answer, so a - * codon split across an intron gets one number for both of its pieces. */ + * that cannot be worked out. Measured from the codon's 5'-most base, so a codon split + * across an intron gets one number for both of its pieces. It has to be the 5'-most base + * and not just any of the three: an indel inside the codon - the very thing this numbering + * exists to account for - moves the other two bases to transcript offsets that are not + * consecutive, and they would divide out to a different codon. */ { if (txAli == NULL || txCds == NULL || txCds->start >= txCds->end) return 0; if (!txCodonNumbersEnabled()) return 0; int txOff = pslTToTxOffset(txAli, posStrand ? start : end-1); if (txOff < 0) return 0; int cdsOff = txOff - txCds->start; if (cdsOff < 0) return 0; return cdsOff/3 + 1; } boolean baseColorCodonIsShifted(struct simpleFeature *sf)