fa5b31d305066e1938953374d80b3782ef87e239 max Mon Sep 7 23:23:40 2026 -0700 Position box: accept a bare codon number, and a range of codon numbers "KAT6A p.495_533" used to land on codon 495 and silently drop the end of the range, and a bare codon number after a transcript accession was not understood at all, so "ENST00000265713.8 p.495" fell through the HGVS code and ended up on an unrelated locus. Nucleotide ranges already worked. The pseudo-HGVS layer now takes an optional _end on a bare codon number, and accepts a bare codon number or range after an NM_ or ENST accession as well as after a gene symbol, looking up the reference amino acids that HGVS wants and the user did not type. The accession forms require a literal "p", so "NM_006766.5 1483" keeps meaning what it meant. A hyphen is still not a range separator: c.1483-1599 is the HGVS intronic position and stays that way. Also fixes a read past the end of the protein sequence when the codon number was larger than the protein, and documents codon ranges in query.html. refs #38285 diff --git src/hg/lib/hgHgvs.c src/hg/lib/hgHgvs.c index b603b35a2d1..de78af5a6a6 100644 --- src/hg/lib/hgHgvs.c +++ src/hg/lib/hgHgvs.c @@ -391,35 +391,62 @@ // 1................... gene symbol // 2..... original sequence // 3...... 1-based position // 4...... replacement sequence #define pseudoHgvsGeneSymbolProtRangeExp "^" geneSymbolExp maybePDot hgvsAaRangeExp "\\)?" // 0..................................................... whole matching string // 1................... gene symbol // 2... original start AA // 3... 1-based start position // 4................ optional range sep and AA+pos // 5... original end AA // 6... 1-based end position // 7..... change description -// As above but omitting the protein change -#define pseudoHgvsGeneSymbolProtPosExp "^" geneSymbolExp maybePDot posIntExp "\\)?" +// As above but omitting the protein change, and allowing a range of codon numbers. +// Someone reading about a mutation usually has the codon number but not the amino acid, +// so "KAT6A p.495" and "KAT6A p.495_533" have to work as well as "KAT6A p.Lys495". +#define posIntRangeExp posIntExp "(_" posIntExp ")?" +#define pseudoHgvsGeneSymbolProtPosExp "^" geneSymbolExp maybePDot posIntRangeExp "\\)?" // 0.......................... whole matching string // 1................... gene symbol -// 2..... 1-based position +// 2..... 1-based start position +// 3....... optional range sep and end position +// 4..... 1-based end position + +// The same bare codon number or range, but after a transcript accession rather than a gene +// symbol. Here the "p" is required: without it "NM_006766.5 1483" would silently become a +// codon number, and a bare number after an accession is far more likely to be something else. +#define pDot "[ :]+p\\.?\\(?" +#define pseudoHgvsNMPDotPosExp "^" versionedRefSeqNMExp pDot posIntRangeExp "\\)?" +// 0.......................... whole matching string +// 1............... acc & optional dot version +// 2........ optional dot version +// 3..... optional gene sym in ()s +// 4... optional gene symbol +// 5..... 1-based start position +// 6....... optional range sep and end position +// 7..... 1-based end position + +#define pseudoHgvsENSPDotPosExp "^" ensTranscriptExp pDot posIntRangeExp "\\)?" +// 0.......................... whole matching string +// 1..................................... ENS transcript ID including optional lift suffix +// 2... optional non-human species code e.g. MUS for mouse +// 3..... 1-based start position +// 4....... optional range sep and end position +// 5..... 1-based end position // Gene symbol, maybe punctuation, and a clear "c." position (and possibly change) #define pseudoHgvsGeneSymbolCDotPosExp "^" geneSymbolExp "[: ]+" hgvsCDotPosExp // 0..................................................... whole matching string // 1................... gene symbol // 2..... optional beginning of position exp // 3..... beginning of position exp // Gene symbol, maybe punctuation, and a clear "n." position (and possibly change) #define pseudoHgvsGeneSymbolNDotPosExp "^" geneSymbolExp "[: ]+" hgvsNDotPosExp // 0..................................................... whole matching string // 1................... gene symbol // 2..... optional beginning of position exp // 3..... beginning of position exp @@ -927,37 +954,60 @@ "where name = '%s'", acc); seq = sqlQuickString(conn, query); } hFreeConn(&conn); } else { aaSeq *aaSeq = hGenBankGetPep(db, acc, NULL); if (aaSeq) seq = aaSeq->dna; } } return seq; } -static char refBaseForNp(char *db, char *npAcc, int pos) -// Get the amino acid base in NP_'s sequence at 1-based offset pos. +static struct hgvsVariant *hgvsFromBareProtPos(char *db, char *protAcc, char *geneSymbol, + int startPos, int endPos) +/* Turn a bare codon number, or range of codon numbers, into a parsed HGVS p. term. + * HGVS wants the reference amino acid with the number, so look it up. endPos may be 0 + * for a single position. Returns NULL if the protein is unknown or the numbers are off + * the end of it. */ { -char *seq = getProteinSeq(db, npAcc); -char base = seq ? seq[pos-1] : '\0'; +struct hgvsVariant *hgvs = NULL; +char *seq = getProteinSeq(db, protAcc); +if (seq == NULL) + return NULL; +int protLen = strlen(seq); +if (endPos == 0) + endPos = startPos; +if (startPos >= 1 && startPos <= protLen && endPos >= startPos && endPos <= protLen) + { + struct dyString *dy = dyStringNew(0); + dyStringAppend(dy, protAcc); + if (isNotEmpty(geneSymbol)) + dyStringPrintf(dy, "(%s)", geneSymbol); + if (endPos > startPos) + dyStringPrintf(dy, ":p.%c%d_%c%d", seq[startPos-1], startPos, seq[endPos-1], endPos); + else + // a single position with no change is spelled as a synonymous substitution + dyStringPrintf(dy, ":p.%c%d=", seq[startPos-1], startPos); + hgvs = hgvsParseTerm(dy->string); + dyStringFree(&dy); + } freeMem(seq); -return base; +return hgvs; } static struct hgvsVariant* hgvsPseudoToRealHgvs(regmatch_t substrs[], char* term, char* db, int geneSymbolIx, char *prefix) /* rewrite a pseudo-HGVS with symbol:<hgvs> to a real NM_xxx HGVS */ { int len = substrs[geneSymbolIx].rm_eo - substrs[geneSymbolIx].rm_so; char geneSymbol[len+1]; safencpy(geneSymbol, sizeof(geneSymbol), term, len); char *nmAcc = nmForGeneSymbol(db, geneSymbol); struct hgvsVariant *hgvs = NULL; if (isNotEmpty(nmAcc)) { // Make it a real HGVS term with the NM and pass that on to the usual parser. int descStartIx = regexSubstrMatched(substrs[2]) ? 2 : 3; char *description = term + substrs[descStartIx].rm_so; @@ -989,30 +1039,46 @@ char ensAcc[len+1]; safencpy(ensAcc, sizeof(ensAcc), term, len); char *enspAcc = enspForEnst(db, ensAcc); if (isNotEmpty(enspAcc)) { // Make it a real HGVS term with the ENSP and pass that on to the usual parser. int descStartIx = 3; char *description = term + substrs[descStartIx].rm_so; struct dyString *enspTerm; enspTerm = dyStringCreate("%s:p.%s", enspAcc, description); hgvs = hgvsParseTerm(enspTerm->string); dyStringFree(&enspTerm); freeMem(enspAcc); } } +else if (regexMatchSubstr(term, pseudoHgvsENSPDotPosExp, substrs, ArraySize(substrs))) + { + // User gave an ENST_ accession and a bare codon number or range of codon numbers. + int ensAccIx = 1, startPosIx = 3, endPosIx = 5; + char *ensAcc = regexSubstringClone(term, substrs[ensAccIx]); + char *enspAcc = enspForEnst(db, ensAcc); + if (isNotEmpty(enspAcc)) + { + int startPos = regexSubstringInt(term, substrs[startPosIx]); + int endPos = regexSubstrMatched(substrs[endPosIx]) ? + regexSubstringInt(term, substrs[endPosIx]) : 0; + hgvs = hgvsFromBareProtPos(db, enspAcc, NULL, startPos, endPos); + } + freeMem(enspAcc); + freeMem(ensAcc); + } else if ((isSubst = regexMatchSubstr(term, pseudoHgvsNMPDotSubstExp, substrs, ArraySize(substrs))) || regexMatchSubstr(term, pseudoHgvsNMPDotRangeExp, substrs, ArraySize(substrs))) { // User gave an NM_ accession but a protein change -- swap in the right NP_. int nmAccIx = 1; int geneSymbolIx = 4; int len = substrs[nmAccIx].rm_eo - substrs[nmAccIx].rm_so; char nmAcc[len+1]; safencpy(nmAcc, sizeof(nmAcc), term, len); char *npAcc = npForNm(db, nmAcc); if (isNotEmpty(npAcc)) { // Make it a real HGVS term with the NP and pass that on to the usual parser. int descStartIx = 5; @@ -1020,30 +1086,49 @@ struct dyString *npTerm; if (regexSubstrMatched(substrs[geneSymbolIx])) { len = substrs[geneSymbolIx].rm_eo - substrs[geneSymbolIx].rm_so; char geneSymbol[len+1]; safencpy(geneSymbol, sizeof(geneSymbol), term, len); npTerm = dyStringCreate("%s(%s):p.%s", npAcc, geneSymbol, description); } else npTerm = dyStringCreate("%s:p.%s", npAcc, description); hgvs = hgvsParseTerm(npTerm->string); dyStringFree(&npTerm); freeMem(npAcc); } } +else if (regexMatchSubstr(term, pseudoHgvsNMPDotPosExp, substrs, ArraySize(substrs))) + { + // User gave an NM_ accession and a bare codon number or range of codon numbers. + int nmAccIx = 1, nmGeneSymbolIx = 4, startPosIx = 5, endPosIx = 7; + char *nmAcc = regexSubstringClone(term, substrs[nmAccIx]); + char *npAcc = npForNm(db, nmAcc); + if (isNotEmpty(npAcc)) + { + char *geneSymbol = regexSubstrMatched(substrs[nmGeneSymbolIx]) ? + regexSubstringClone(term, substrs[nmGeneSymbolIx]) : NULL; + int startPos = regexSubstringInt(term, substrs[startPosIx]); + int endPos = regexSubstrMatched(substrs[endPosIx]) ? + regexSubstringInt(term, substrs[endPosIx]) : 0; + hgvs = hgvsFromBareProtPos(db, npAcc, geneSymbol, startPos, endPos); + freeMem(geneSymbol); + } + freeMem(npAcc); + freeMem(nmAcc); + } else if ((isSubst = regexMatchSubstr(term, pseudoHgvsGeneSymbolProtSubstExp, substrs, ArraySize(substrs))) || regexMatchSubstr(term, pseudoHgvsGeneSymbolProtRangeExp, substrs, ArraySize(substrs))) { int len = substrs[geneSymbolIx].rm_eo - substrs[geneSymbolIx].rm_so; char geneSymbol[len+1]; safencpy(geneSymbol, sizeof(geneSymbol), term, len); struct slName *npAccList = npForGeneSymbol(db, geneSymbol); if (npAccList != NULL) { struct slName *npItem = NULL; for (npItem = npAccList; npItem != NULL; npItem = npItem->next) { char *npAcc = npItem->name; // Make it a real HGVS term with the NP and pass that on to the usual parser. @@ -1058,41 +1143,39 @@ } } } else if (regexMatchSubstr(term, pseudoHgvsGeneSymbolProtPosExp, substrs, ArraySize(substrs))) { int len = substrs[geneSymbolIx].rm_eo - substrs[geneSymbolIx].rm_so; char geneSymbol[len+1]; safencpy(geneSymbol, sizeof(geneSymbol), term, len); struct slName *npAccList = npForGeneSymbol(db, geneSymbol); if (npAccList != NULL) { struct slName *npItem = NULL; for (npItem = npAccList; npItem != NULL; npItem = npItem->next) { char *npAcc = npItem->name; - // Only position was provided, no change. Look up ref base and make a synonymous subst - // so it's parseable HGVS. - int posIx = 2; - int pos = regexSubstringInt(term, substrs[posIx]); - char refBase = refBaseForNp(db, npAcc, pos); - struct dyString *npTerm = dyStringCreate("%s(%s):p.%c%d=", - npAcc, geneSymbol, refBase, pos); - struct hgvsVariant *newTerm = hgvsParseTerm(npTerm->string); + // Only the codon number(s) were provided, no amino acid and no change. + int startPosIx = 2, endPosIx = 4; + int startPos = regexSubstringInt(term, substrs[startPosIx]); + int endPos = regexSubstrMatched(substrs[endPosIx]) ? + regexSubstringInt(term, substrs[endPosIx]) : 0; + struct hgvsVariant *newTerm = hgvsFromBareProtPos(db, npAcc, geneSymbol, + startPos, endPos); if (newTerm) slAddHead(&hgvs, newTerm); - dyStringFree(&npTerm); } } } else if (regexMatchSubstr(term, pseudoHgvsGeneSymbolCDotPosExp, substrs, ArraySize(substrs))) hgvs = hgvsPseudoToRealHgvs(substrs, term, db, geneSymbolIx, "c."); else if (regexMatchSubstr(term, pseudoHgvsGeneSymbolNDotPosExp, substrs, ArraySize(substrs))) hgvs = hgvsPseudoToRealHgvs(substrs, term, db, geneSymbolIx, "n."); else if (regexMatchSubstr(term, pseudoHgvsChrGDotExp, substrs, ArraySize(substrs))) { int chrIx = 1; int startPosIx = 3; int endPosIx = 5; int changeIx = 6; AllocVar(hgvs); hgvs->type = hgvstGenomic;