4decf5fbe82051b2d5aff9cabce3feae9dfafae1 max Sat Sep 26 18:48:55 2026 -0700 uniprot: stop declaring the alignments as amino acid coordinates, and show them properly The bigPsl seqType field describes the coordinates, not the letters stored beside them, and the UniProt query side is in bases: these proteins reach the genome through transcripts, so a query runs three bases to a residue. Declaring amino acids made pslFromBigPsl divide the block sizes by three and leave the query coordinates alone, so every reader got an alignment measured in two units at once. That fed a heap overflow in the alignment page, drew blocks short in hgTracks, and loaded sub-codon blocks as size 0, which aborted pslTransMap and took down the lifted SwissProt track (#38249). pslProtFromNaLike() converts such a psl to one counted in residues. Blocks are trimmed to whole codons, and a residue whose codon straddles an exon junction sits in two places in the genome at once, so it gets no column and the page says how many are missing rather than dropping them silently: 0.9% of residues, though 87% of alignments have at least one. A bigPsl also keeps the reference strand where a psl reads the query strand, so minus-strand items arrived claiming the protein was reversed and were rendered as reverse complemented nucleotide ambiguity codes; pslRc moves them to the convention blat uses, query forward and the strand on the target. Checked on hg38 against the translated genome: 210,416 residues over both strands, 99.86% identical, the remainder real protein-vs-reference variation. All 29 items of a test region render the stored protein at the right residues. Files already published still say amino acid and must keep working until they are rebuilt, so the conversion also requires the blocks to measure the target the way the target is measured; verified that separates the two shapes on 3000 records each way, and that all 29 render without crashing in the old format, where they now say plainly that the coordinates and the sequence do not match. refs #38300 diff --git src/hg/hgc/hgc.c src/hg/hgc/hgc.c index 6bfb7c2ea42..0e2126efc6c 100644 --- src/hg/hgc/hgc.c +++ src/hg/hgc/hgc.c @@ -8560,56 +8560,83 @@ "select chromStart from %s where frag = \"%s\"", goldTable, bactig->endContig); ctgStart = sqlQuickNum(conn, query); snprintf(ctgStartStr, sizeof(ctgStartStr), "%d", ctgStart); hgcAnchor("gold", bactig->endContig, ctgStartStr); } printf("%s</A><BR>\n", bactig->endContig); printPos(bactig->chrom, bactig->chromStart, bactig->chromEnd, NULL, FALSE,NULL); printTrackHtml(tdb); hFreeConn(&conn); } +/* Residues left out of the alignment because their codon spans an exon junction, set by the + * bigPsl alignment pages when they convert an na-like protein psl and reported by + * showGfAlignment(). It travels this way because showGfAlignment() is reached through + * showSomeAlignment(), whose signature is shared with every other alignment page. */ +static int gNaLikeDroppedAa = 0; + int showGfAlignment(struct psl *psl, bioSeq *qSeq, FILE *f, enum gfType qType, int qStart, int qEnd, char *qName) /* Show protein/DNA alignment or translated DNA alignment. */ { int blockCount; int tStart = psl->tStart; int tEnd = psl->tEnd; char tName[256]; struct dnaSeq *tSeq; +/* A protein psl counts its query in amino acids, so the sequence the block coordinates are about + * to be laid onto has to be qSize long. Anything else means the psl and the sequence stored with + * it do not describe each other, and drawing from coordinates that are not this sequence's reads + * off the end of it. Say so instead, the way showPartialDnaAlignment() does for the DNA case. */ +if (qType == gftProt && qSeq->size != psl->qSize) + { + fprintf(f, "<p><b>Cannot display alignment.</b> The query sequence stored for %s is %d amino " + "acids long, but the alignment was made against a query of %d. The query " + "coordinates in this file do not match the sequence it carries.\n", + psl->qName, qSeq->size, psl->qSize); + return 0; + } + /* protein psl's have a tEnd that isn't quite right */ if ((psl->strand[1] == '+') && (qType == gftProt)) tEnd = psl->tStarts[psl->blockCount - 1] + psl->blockSizes[psl->blockCount - 1] * 3; tSeq = hDnaFromSeq(database, seqName, tStart, tEnd, dnaLower); freez(&tSeq->name); tSeq->name = cloneString(psl->tName); safef(tName, sizeof(tName), "%s.%s", organism, psl->tName); if (qName == NULL) fprintf(f, "<H2>Alignment of %s and %s:%d-%d</H2>\n", psl->qName, psl->tName, psl->tStart+1, psl->tEnd); else fprintf(f, "<H2>Alignment of %s and %s:%d-%d</H2>\n", qName, psl->tName, psl->tStart+1, psl->tEnd); +if (gNaLikeDroppedAa == 1) + fprintf(f, "<p>One amino acid of %s is not shown below: its codon is split across an exon " + "junction, so it has no single position in the genome.</p>\n", psl->qName); +else if (gNaLikeDroppedAa > 1) + fprintf(f, "<p>%d amino acids of %s are not shown below: their codons are split across exon " + "junctions, so they have no single position in the genome.</p>\n", + gNaLikeDroppedAa, psl->qName); + if (!cartUsualBoolean(cart, "blatNewPage", FALSE)) /* no "frame" in the new single-page view */ fputs("Click on links in the frame to the left to navigate through " "the alignment.\n", f); blockCount = pslShowAlignment(psl, qType == gftProt, qName, qSeq, qStart, qEnd, tName, tSeq, tStart, tEnd, f); freeDnaSeq(&tSeq); return blockCount; } static struct ffAli *pslToFfAliAndSequence(struct psl *psl, struct dnaSeq *qSeq, boolean *retIsRc, struct dnaSeq **retSeq, int *retTStart) /* Given psl, dig up target sequence and convert to ffAli. * Note: if strand is -, this does a pslRc to psl! */ @@ -9043,31 +9070,41 @@ } if (psl == NULL) errAbort("item %s not found in range %s:%d-%d in bigBed %s (%s)", acc, chrom, start, end, tdb->table, fileName); if (cdsString) genbankParseCds(cdsString, &cdsStart, &cdsEnd); if (seq == NULL) { printf("Sequence for %s not available.\n", psl->qName); return; } struct dnaSeq *rnaSeq = newDnaSeq(seq, strlen(seq), acc); enum gfType type = gftRna; -if (pslIsProtein(psl)) +/* A protein whose query side is counted in bases rather than residues - the UniProt alignments, + * which reach the genome through transcripts - has to be converted before it can be shown against + * the protein it stores, or its coordinates address three times the sequence there is. */ +struct psl *protPsl = pslProtFromNaLike(psl, rnaSeq->size, &gNaLikeDroppedAa); +if (protPsl != NULL) + { + pslFree(&psl); + psl = protPsl; + type = gftProt; + } +else if (pslIsProtein(psl)) type = gftProt; showSomeAlignment(psl, rnaSeq, type, 0, rnaSeq->size, NULL, cdsStart, cdsEnd); } void htcBigPslAliInWindow(char *acc) /* Show alignment in window for accession in bigPsl file. */ { struct psl *partPsl, *wholePsl; char *aliTable; int start; unsigned int cdsStart = 0, cdsEnd = 0; struct trackDb *tdb = NULL; aliTable = cartString(cart, "aliTable"); struct quickLiftAli ali; @@ -9127,35 +9164,50 @@ } pslFree(&bbPsl); } if (wholePsl == NULL) errAbort("item %s not found in range %s:%d-%d in bigBed %s (%s)", acc, chrom, start, end, tdb->table, fileName); if (seq == NULL) { printf("Sequence for %s not available.\n", wholePsl->qName); return; } if (cdsString) genbankParseCds(cdsString, &cdsStart, &cdsEnd); +struct dnaSeq *rnaSeq = newDnaSeq(seq, strlen(seq), acc); +/* The same conversion htcBigPslAli() makes: a query counted in bases has to become one counted + * in residues before it can be shown against the protein stored with it. Convert before + * trimming, which only touches the target side, so the window still decides what is shown. */ +struct psl *protPsl = pslProtFromNaLike(wholePsl, rnaSeq->size, &gNaLikeDroppedAa); +if (protPsl != NULL) + { + pslFree(&wholePsl); + wholePsl = protPsl; + } if (wholePsl->tStart >= winStart && wholePsl->tEnd <= winEnd) partPsl = wholePsl; else partPsl = pslTrimToTargetRange(wholePsl, winStart, winEnd); -struct dnaSeq *rnaSeq = newDnaSeq(seq, strlen(seq), acc); +if (protPsl != NULL) + /* showSomePartialDnaAlignment() renders a nucleotide query. A protein one belongs on the + * translated path, which has no partial form, so show the whole alignment rather than + * laying residue coordinates out as though they were bases. */ + showSomeAlignment(wholePsl, rnaSeq, gftProt, 0, rnaSeq->size, NULL, cdsStart, cdsEnd); +else showSomePartialDnaAlignment(partPsl, wholePsl, rnaSeq, NULL, cdsStart, cdsEnd); } static struct dnaSeq *getBaseColorSequence(char *db, char *itemName, char *table) /* Grab sequence using the sequence and extFile table names out of BASE_COLOR_USE_SEQUENCE. * db is the assembly the sequence lives in, which is not the one on screen when the track * is quickLifted. */ { struct trackDb *tdb = hashMustFindVal(trackHash, table); char *spec = trackDbRequiredSetting(tdb, BASE_COLOR_USE_SEQUENCE); char *specCopy = cloneString(spec); // value is: extFile seqTbl extFileTbl // or: db [dddBbb1]