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/lib/psl.c src/lib/psl.c index 9a0fe14d655..4e13415465e 100644 --- src/lib/psl.c +++ src/lib/psl.c @@ -1,2276 +1,2360 @@ /* psl.c was originally generated by the autoSql program, which also * generated as_psl.h and as_psl.sql. This module links the database and the RAM * representation of objects. * * This file is copyright 2002 Jim Kent, but license is hereby * granted for all use - public, private or commercial. */ #include "common.h" #include "sqlNum.h" #include "sqlList.h" #include "localmem.h" #include "psl.h" #include "hash.h" #include "linefile.h" #include "dnaseq.h" #include "dystring.h" #include "fuzzyFind.h" #include "aliType.h" #include "binRange.h" #include "rangeTree.h" struct psl *pslxLoad(char **row) /* Load a psl from row fetched with select * from psl * from database. Dispose of this with pslFree(). */ { struct psl *ret = pslLoad(row); int retSize; sqlStringDynamicArray(row[21],&ret->qSequence, &retSize); sqlStringDynamicArray(row[22],&ret->tSequence, &retSize); return ret; } struct psl *pslLoad(char **row) /* Load a psl from row fetched with select * from psl * from database. Dispose of this with pslFree(). */ { struct psl *ret; int sizeOne; AllocVar(ret); ret->blockCount = sqlUnsigned(row[17]); ret->match = sqlUnsigned(row[0]); ret->misMatch = sqlUnsigned(row[1]); ret->repMatch = sqlUnsigned(row[2]); ret->nCount = sqlUnsigned(row[3]); ret->qNumInsert = sqlUnsigned(row[4]); ret->qBaseInsert = sqlSigned(row[5]); ret->tNumInsert = sqlUnsigned(row[6]); ret->tBaseInsert = sqlSigned(row[7]); strcpy(ret->strand, row[8]); ret->qName = cloneString(row[9]); ret->qSize = sqlUnsigned(row[10]); ret->qStart = sqlUnsigned(row[11]); ret->qEnd = sqlUnsigned(row[12]); ret->tName = cloneString(row[13]); ret->tSize = sqlUnsigned(row[14]); ret->tStart = sqlUnsigned(row[15]); ret->tEnd = sqlUnsigned(row[16]); sqlUnsignedDynamicArray(row[18], &ret->blockSizes, &sizeOne); if (sizeOne != ret->blockCount) { printf("sizeOne bloxksizes %d bs %d block=%s\n",sizeOne, ret->blockCount,row[18]); } assert(sizeOne == ret->blockCount); sqlUnsignedDynamicArray(row[19], &ret->qStarts, &sizeOne); if (sizeOne != ret->blockCount) { printf("sizeOne qStarts %d bs %d\n",sizeOne, ret->blockCount); } assert(sizeOne == ret->blockCount); sqlUnsignedDynamicArray(row[20], &ret->tStarts, &sizeOne); if (sizeOne != ret->blockCount) { printf("sizeOne tStarts %d bs %d\n",sizeOne, ret->blockCount); } assert(sizeOne == ret->blockCount); return ret; } struct psl *pslCommaIn(char **pS, struct psl *ret) /* Create a psl out of a comma separated string. * This will fill in ret if non-null, otherwise will * return a new psl */ { char *s = *pS; int i; if (ret == NULL) AllocVar(ret); ret->match = sqlUnsignedComma(&s); ret->misMatch = sqlUnsignedComma(&s); ret->repMatch = sqlUnsignedComma(&s); ret->nCount = sqlUnsignedComma(&s); ret->qNumInsert = sqlUnsignedComma(&s); ret->qBaseInsert = sqlSignedComma(&s); ret->tNumInsert = sqlUnsignedComma(&s); ret->tBaseInsert = sqlSignedComma(&s); sqlFixedStringComma(&s, ret->strand, sizeof(ret->strand)); ret->qName = sqlStringComma(&s); ret->qSize = sqlUnsignedComma(&s); ret->qStart = sqlUnsignedComma(&s); ret->qEnd = sqlUnsignedComma(&s); ret->tName = sqlStringComma(&s); ret->tSize = sqlUnsignedComma(&s); ret->tStart = sqlUnsignedComma(&s); ret->tEnd = sqlUnsignedComma(&s); ret->blockCount = sqlUnsignedComma(&s); s = sqlEatChar(s, '{'); AllocArray(ret->blockSizes, ret->blockCount); for (i=0; i<ret->blockCount; ++i) { ret->blockSizes[i] = sqlUnsignedComma(&s); } s = sqlEatChar(s, '}'); s = sqlEatChar(s, ','); s = sqlEatChar(s, '{'); AllocArray(ret->qStarts, ret->blockCount); for (i=0; i<ret->blockCount; ++i) { ret->qStarts[i] = sqlUnsignedComma(&s); } s = sqlEatChar(s, '}'); s = sqlEatChar(s, ','); s = sqlEatChar(s, '{'); AllocArray(ret->tStarts, ret->blockCount); for (i=0; i<ret->blockCount; ++i) { ret->tStarts[i] = sqlUnsignedComma(&s); } s = sqlEatChar(s, '}'); s = sqlEatChar(s, ','); *pS = s; return ret; } void pslFree(struct psl **pEl) /* Free a single dynamically allocated psl such as created * with pslLoad(). */ { struct psl *el; if ((el = *pEl) == NULL) return; freeMem(el->qName); freeMem(el->tName); freeMem(el->blockSizes); freeMem(el->qStarts); freeMem(el->tStarts); if (el->qSequence) { freeMem(el->qSequence[0]); freeMem(el->qSequence); } if (el->tSequence) { freeMem(el->tSequence[0]); freeMem(el->tSequence); } freez(pEl); } void pslFreeList(struct psl **pList) /* Free a list of dynamically allocated psl's */ { struct psl *el, *next; for (el = *pList; el != NULL; el = next) { next = el->next; pslFree(&el); } *pList = NULL; } void pslOutputJson(struct psl *el, FILE *f) /* print out psl as a json array */ { fputs("[", f); fprintf(f, "%u,", el->match); fprintf(f, "%u,", el->misMatch); fprintf(f, "%u,", el->repMatch); fprintf(f, "%u,", el->nCount); fprintf(f, "%u,", el->qNumInsert); fprintf(f, "%d,", el->qBaseInsert); fprintf(f, "%u,", el->tNumInsert); fprintf(f, "%d,", el->tBaseInsert); fprintf(f, "\"%s\",", el->strand); fprintf(f, "\"%s\",", el->qName); fprintf(f, "%u,", el->qSize); fprintf(f, "%u,", el->qStart); fprintf(f, "%u,", el->qEnd); fprintf(f, "\"%s\",", el->tName); fprintf(f, "%u,", el->tSize); fprintf(f, "%u,", el->tStart); fprintf(f, "%u,", el->tEnd); fprintf(f, "%u,", el->blockCount); fputs("\"", f); for (int i=0; i<el->blockCount; ++i) { fprintf(f, "%u", el->blockSizes[i]); if (i < el->blockCount-1) fputc(',', f); } fputs("\",", f); fputs("\"", f); for (int i=0; i<el->blockCount; ++i) { fprintf(f, "%u", el->qStarts[i]); if (i < el->blockCount-1) fputc(',', f); // json does not allow trailing commas } fputs("\",", f); fputs("\"", f); for (int i=0; i<el->blockCount; ++i) { fprintf(f, "%u", el->tStarts[i]); if (i < el->blockCount-1) fputc(',', f); } fputs("\"", f); if (el->qSequence) { fputc(',',f); fputc('[',f); for (int i=0; i<el->blockCount; ++i) { fprintf(f, "'%s'", el->qSequence[i]); if (i < el->blockCount-1) fputc(',', f); } fputc(']',f); fputc(',',f); fputc('[',f); for (int i=0; i<el->blockCount; ++i) { fprintf(f, "\"%s\"", el->tSequence[i]); if (i < el->blockCount-1) fputc(',', f); } fputc(']',f); } if (ferror(f)) { perror("Error writing psl file\n"); errAbort("\n"); } fputs("]", f); } void pslOutput(struct psl *el, FILE *f, char sep, char lastSep) /* Print out psl. Separate fields with sep. Follow last field with lastSep. */ { int i; fprintf(f, "%u", el->match); fputc(sep,f); fprintf(f, "%u", el->misMatch); fputc(sep,f); fprintf(f, "%u", el->repMatch); fputc(sep,f); fprintf(f, "%u", el->nCount); fputc(sep,f); fprintf(f, "%u", el->qNumInsert); fputc(sep,f); fprintf(f, "%d", el->qBaseInsert); fputc(sep,f); fprintf(f, "%u", el->tNumInsert); fputc(sep,f); fprintf(f, "%d", el->tBaseInsert); fputc(sep,f); if (sep == ',') fputc('"',f); fprintf(f, "%s", el->strand); if (sep == ',') fputc('"',f); fputc(sep,f); if (sep == ',') fputc('"',f); fprintf(f, "%s", el->qName); if (sep == ',') fputc('"',f); fputc(sep,f); fprintf(f, "%u", el->qSize); fputc(sep,f); fprintf(f, "%u", el->qStart); fputc(sep,f); fprintf(f, "%u", el->qEnd); fputc(sep,f); if (sep == ',') fputc('"',f); fprintf(f, "%s", el->tName); if (sep == ',') fputc('"',f); fputc(sep,f); fprintf(f, "%u", el->tSize); fputc(sep,f); fprintf(f, "%u", el->tStart); fputc(sep,f); fprintf(f, "%u", el->tEnd); fputc(sep,f); fprintf(f, "%u", el->blockCount); fputc(sep,f); if (sep == ',') fputc('{',f); for (i=0; i<el->blockCount; ++i) { fprintf(f, "%u", el->blockSizes[i]); fputc(',', f); } if (sep == ',') fputc('}',f); fputc(sep,f); if (sep == ',') fputc('{',f); for (i=0; i<el->blockCount; ++i) { fprintf(f, "%u", el->qStarts[i]); fputc(',', f); } if (sep == ',') fputc('}',f); fputc(sep,f); if (sep == ',') fputc('{',f); for (i=0; i<el->blockCount; ++i) { fprintf(f, "%u", el->tStarts[i]); fputc(',', f); } if (sep == ',') fputc('}',f); if (el->qSequence) { fputc(sep,f); if (sep == ',') fputc('{',f); for (i=0; i<el->blockCount; ++i) { fprintf(f, "%s", el->qSequence[i]); fputc(',', f); } if (sep == ',') fputc('}',f); fputc(sep,f); if (sep == ',') fputc('{',f); for (i=0; i<el->blockCount; ++i) { fprintf(f, "%s", el->tSequence[i]); fputc(',', f); } if (sep == ',') fputc('}',f); } fputc(lastSep,f); if (ferror(f)) { perror("Error writing psl file\n"); errAbort("\n"); } } /* ----- end autoSql generated part --------------- */ void pslOutputShort(struct psl *el, FILE *f, char sep, char lastSep) /* Print out psl. Separate fields with sep. Follow last field with lastSep. */ { fprintf(f, "%u", el->match); fputc(sep,f); fprintf(f, "%u", el->misMatch); fputc(sep,f); fprintf(f, "%u", el->repMatch); fputc(sep,f); fprintf(f, "%u", el->qNumInsert); fputc(sep,f); fprintf(f, "%d", el->qBaseInsert); fputc(sep,f); fprintf(f, "%u", el->tNumInsert); fputc(sep,f); fprintf(f, "%d", el->tBaseInsert); fputc(sep,f); if (sep == ',') fputc('"',f); fprintf(f, "%s", el->strand); if (sep == ',') fputc('"',f); fputc(sep,f); if (sep == ',') fputc('"',f); fprintf(f, "%s", el->qName); if (sep == ',') fputc('"',f); fputc(sep,f); fprintf(f, "%u", el->qStart); fputc(sep,f); fprintf(f, "%u", abs(el->qEnd - el->qStart)); fputc(sep,f); if (sep == ',') fputc('"',f); fprintf(f, "%s", el->tName); if (sep == ',') fputc('"',f); fputc(sep,f); fprintf(f, "%u", el->tStart); fputc(sep,f); fprintf(f, "%u", abs(el->tEnd - el->tStart)); fputc(sep,f); fprintf(f, "%u", el->blockCount); fputc(sep,f); if (sep == ',') fputc('{',f); fputc(lastSep,f); if (ferror(f)) { perror("Error writing psl file\n"); errAbort("\n"); } } void pslOutFormat(struct psl *el, FILE *f, char sep, char lastSep) /* Print out selected psl values. Separate fields with sep. Follow last field with lastSep. */ /* Prints out a better format with bold field headings followed by value */ /* Requires further upstream work to ensure that only the field headers */ /* declared here are printed if replacing an existing psl print function*/ { const char *headers[] = {"Matches", "Mismatches", "Matches in repeats", "Number of N bases", "Query name", "Size", "Start", "End", "Chromosome", "Strand", "Start", "End"}; char *hformat = "<B>%s:</B> "; /* string for formatted print for headers */ char *uformat = "<B>%s:</B> %u%c"; /* string for formatted print for unsigned variable */ char *targName; fprintf(f, uformat, headers[0], el->match, sep); fprintf(f, uformat, headers[1], el->misMatch, sep); fprintf(f, uformat, headers[2], el->repMatch, sep); fprintf(f, uformat, headers[3], el->nCount, sep); fprintf(f, hformat, headers[4]); if (sep == ',') fputc('"',f); fprintf(f, "%s", el->qName); if (sep == ',') fputc('"',f); fputc(sep,f); fprintf(f, uformat, headers[5], el->qSize, sep); fprintf(f, uformat, headers[6], el->qStart, sep); fprintf(f, uformat, headers[7], el->qEnd, sep); fprintf(f, hformat, headers[8]); if (sep == ',') fputc('"',f); /* skip leading 'chr' in string to get only chromosome part */ targName = el->tName; if (startsWith("chr", el->tName)) targName += 3; fprintf(f, "%s", targName); if (sep == ',') fputc('"',f); fputc(sep,f); fprintf(f, hformat, headers[9]); if (sep == ',') fputc('"',f); fprintf(f, "%s", el->strand); if (sep == ',') fputc('"',f); fputc(sep,f); fprintf(f, uformat, headers[10], el->tStart, sep); fprintf(f, uformat, headers[11], el->tEnd, sep); fputc(lastSep,f); if (ferror(f)) { perror("Error writing psl file\n"); errAbort("\n"); } } struct psl *pslLoadAll(char *fileName) /* Load all psl's in file. */ { struct lineFile *lf = pslFileOpen(fileName); struct psl *pslList = NULL, *psl; while ((psl = pslNext(lf)) != NULL) { slAddHead(&pslList, psl); } slReverse(&pslList); lineFileClose(&lf); return pslList; } int pslCmpQuery(const void *va, const void *vb) /* Compare to sort based on query start. */ { const struct psl *a = *((struct psl **)va); const struct psl *b = *((struct psl **)vb); int dif; dif = strcmp(a->qName, b->qName); if (dif == 0) dif = a->qStart - b->qStart; return dif; } int pslCmpTarget(const void *va, const void *vb) /* Compare to sort based on target start. */ { const struct psl *a = *((struct psl **)va); const struct psl *b = *((struct psl **)vb); int dif; dif = strcmp(a->tName, b->tName); if (dif == 0) dif = a->tStart - b->tStart; return dif; } int pslCmpTargetAndStrand(const void *va, const void *vb) /* Compare to sort based on target, strand, tStart. */ { const struct psl *a = *((struct psl **)va); const struct psl *b = *((struct psl **)vb); int dif; dif = strcmp(a->tName, b->tName); if (dif == 0) dif = strcmp(a->strand, b->strand); if (dif == 0) dif = a->tStart - b->tStart; return dif; } int pslCmpScore(const void *va, const void *vb) /* Compare to sort based on score (descending). */ { const struct psl *a = *((struct psl **)va); const struct psl *b = *((struct psl **)vb); return pslScore(b) - pslScore(a); } int pslCmpQueryScore(const void *va, const void *vb) /* Compare to sort based on query then score (descending). */ { const struct psl *a = *((struct psl **)va); const struct psl *b = *((struct psl **)vb); int diff = strcmp(a->qName, b->qName); if (diff == 0) diff = pslScore(b) - pslScore(a); return diff; } int pslCmpMatch(const void *va, const void *vb) /* Compare to sort based on match */ { const struct psl *a = *((struct psl **)va); const struct psl *b = *((struct psl **)vb); return b->match - a->match; } static void pslLabelColumns(FILE *f) /* Write column info. */ { fputs("\n" "match\tmis- \trep. \tN's\tQ gap\tQ gap\tT gap\tT gap\tstrand\tQ \tQ \tQ \tQ \tT \tT \tT \tT \tblock\tblockSizes \tqStarts\t tStarts\n" " \tmatch\tmatch\t \tcount\tbases\tcount\tbases\t \tname \tsize\tstart\tend\tname \tsize\tstart\tend\tcount\n" "---------------------------------------------------------------------------------------------------------------------------------------------------------------\n", f); } static void pslLabelColumnsJson(FILE *f) /* Write column info as a JSON array */ { fputs("[\"matches\", \"misMatches\", \"repMatches\", \"nCount\", \"qNumInsert\", \"qBaseInsert\", " "\"tNumInsert\", \"tBaseInsert\", \"strand\", " "\"qName\", \"qSize\", \"qStart\", \"qEnd\", " "\"tName\", \"tSize\", \"tStart\", \"tEnd\", " "\"blockCount\", \"blockSizes\", \"qStarts\", \"tStarts\"]", f); } void pslxWriteHead(FILE *f, enum gfType qType, enum gfType tType) /* Write header for extended (possibly protein) psl file. */ { fprintf(f, "psLayout version 4 %s %s\n", gfTypeName(qType), gfTypeName(tType)); pslLabelColumns(f); } void pslWriteHead(FILE *f) /* Write head of psl. */ { fputs("psLayout version 3\n", f); pslLabelColumns(f); } void pslWriteAll(struct psl *pslList, char *fileName, boolean writeHeader) /* Write a psl file from list. */ { FILE *f; struct psl *psl; f = mustOpen(fileName, "w"); if (writeHeader) pslWriteHead(f); for (psl = pslList; psl != NULL; psl = psl->next) pslTabOut(psl, f); fclose(f); } void pslWriteAllJson(struct psl *pslList, FILE *f, char *db, boolean writeHeader) /* Write a psl file from list as a json array . */ { fputs("{\n", f); if (writeHeader) { fputs("\"track\": \"blat\",\n", f); fprintf(f, "\"genome\": \"%s\",\n", db); fputs("\"fields\": ", f); pslLabelColumnsJson(f); fputs(",\n", f); } fputs("\"blat\" : [\n", f); for (struct psl *psl = pslList; psl; psl = psl->next) { pslOutputJson(psl, f); if (psl->next) fputs(",\n", f); } puts("\n]\n}\n"); } void pslxFileOpen(char *fileName, enum gfType *retQueryType, enum gfType *retTargetType, struct lineFile **retLf) /* Read header part of psl and make sure it's right. Return * sequence types and file handle. */ { char *line; int lineSize; char *words[30]; char *version; int wordCount; int i; enum gfType qt = gftRna, tt = gftDna; struct lineFile *lf = lineFileOpen(fileName, TRUE); if (!lineFileNext(lf, &line, &lineSize)) warn("%s is empty", fileName); else { if (startsWith("psLayout version", line)) { wordCount = chopLine(line, words); if (wordCount < 3) errAbort("%s is not a psLayout file", fileName); version = words[2]; if (sameString(version, "3")) { } else if (sameString(version, "4")) { qt = gfTypeFromName(words[3]); tt = gfTypeFromName(words[4]); } else { errAbort("%s is version %s of psLayout, this program can only handle through version 4", fileName, version); } for (i=0; i<4; ++i) { if (!lineFileNext(lf, &line, &lineSize)) errAbort("%s severely truncated", fileName); } } else lineFileReuse(lf); } *retQueryType = qt; *retTargetType = tt; *retLf = lf; } static void pslxFileOpenWithMetaConfig(char *fileName, bool isMetaUnique, enum gfType *retQueryType, enum gfType *retTargetType, struct lineFile **retLf, FILE *f) /* Read header part of psl and make sure it's right. Return * sequence types and file handle and send meta data to output file f */ { char *line; int lineSize; char *words[30]; char *version; int wordCount; int i; enum gfType qt = gftRna, tt = gftDna; struct lineFile *lf = lineFileOpen(fileName, TRUE); lineFileSetMetaDataOutput(lf, f); if (isMetaUnique) lineFileSetUniqueMetaData(lf); if (!lineFileNext(lf, &line, &lineSize)) warn("%s is empty", fileName); else { if (startsWith("psLayout version", line)) { wordCount = chopLine(line, words); if (wordCount < 3) errAbort("%s is not a psLayout file", fileName); version = words[2]; if (sameString(version, "3")) { } else if (sameString(version, "4")) { qt = gfTypeFromName(words[3]); tt = gfTypeFromName(words[4]); } else { errAbort("%s is version %s of psLayout, this program can only handle through version 4", fileName, version); } for (i=0; i<4; ++i) { if (!lineFileNext(lf, &line, &lineSize)) errAbort("%s severely truncated", fileName); } } else { char *s = cloneString(line); boolean eof = FALSE; while ((line[0] == '#') && (!eof)) { freeMem(s); if (!lineFileNext(lf, &line, &lineSize)) eof = TRUE; s = cloneString(line); } wordCount = chopLine(s, words); if ((wordCount < 21 || wordCount > 23 || (words[8][0] != '+' && words[8][0] != '-')) && (!eof)) errAbort("%s is not a psLayout file", fileName); else lineFileReuse(lf); freeMem(s); } } *retQueryType = qt; *retTargetType = tt; *retLf = lf; } void pslxFileOpenWithMeta(char *fileName, enum gfType *retQueryType, enum gfType *retTargetType, struct lineFile **retLf, FILE *f) /* Read header part of psl and make sure it's right. Return * sequence types and file handle and send meta data to output file f */ { pslxFileOpenWithMetaConfig(fileName, FALSE, retQueryType, retTargetType, retLf, f); } void pslxFileOpenWithUniqueMeta(char *fileName, enum gfType *retQueryType, enum gfType *retTargetType, struct lineFile **retLf, FILE *f) /* Read header part of psl and make sure it's right. Return * sequence types and file handle and send only unique meta data to output f */ { pslxFileOpenWithMetaConfig(fileName, TRUE, retQueryType, retTargetType, retLf, f); } struct lineFile *pslFileOpen(char *fileName) /* Read header part of psl and make sure it's right. * Return line file handle to it. */ { enum gfType qt, tt; struct lineFile *lf; pslxFileOpen(fileName, &qt, &tt, &lf); return lf; } struct lineFile *pslFileOpenWithMeta(char *fileName, FILE *f) /* Read header part of psl and make sure it's right. * Return line file handle to it. */ { enum gfType qt, tt; struct lineFile *lf; pslxFileOpenWithMeta(fileName, &qt, &tt, &lf, f); return lf; } struct lineFile *pslFileOpenWithUniqueMeta(char *fileName, FILE *f) /* Read header part of psl and make sure it's right. * Set flag to suppress duplicate header comments. * Return line file handle to it. */ { enum gfType qt, tt; struct lineFile *lf; pslxFileOpenWithUniqueMeta(fileName, &qt, &tt, &lf, f); return lf; } struct psl *pslNext(struct lineFile *lf) /* Read next line from file and convert it to psl. Return * NULL at eof. */ { char *line; int lineSize; char *words[32]; int wordCount; static int lineAlloc = 0; static char *chopBuf = NULL; if (!lineFileNextReal(lf, &line)) { return NULL; } lineSize = strlen(line); if (lineSize >= lineAlloc) { lineAlloc = lineSize+256; chopBuf = needMoreMem(chopBuf, 0, lineAlloc); } memcpy(chopBuf, line, lineSize+1); wordCount = chopLine(chopBuf, words); if (wordCount == 21) { return pslLoad(words); } if (wordCount == 23) { return pslxLoad(words); } else { errAbort("Bad line %d of %s wordCount is %d instead of 21 or 23\n", lf->lineIx, lf->fileName, wordCount); return NULL; } } struct psl *pslxLoadLm(char **row, struct lm *lm) /* Load row into local memory pslx. */ { struct psl *ret = pslLoadLm(row, lm); ret->qSequence = lmAlloc(lm, sizeof(ret->qSequence[0]) * ret->blockCount); sqlStringArray(lmCloneString(lm,row[21]),ret->qSequence, ret->blockCount); ret->tSequence = lmAlloc(lm, sizeof(ret->tSequence[0]) * ret->blockCount); sqlStringArray(lmCloneString(lm,row[22]),ret->tSequence, ret->blockCount); return ret; } struct psl *pslLoadLm(char **row, struct lm *lm) /* Load row into local memory psl. */ { struct psl *ret; ret = lmAlloc(lm, sizeof(*ret)); ret->blockCount = sqlUnsigned(row[17]); ret->match = sqlUnsigned(row[0]); ret->misMatch = sqlUnsigned(row[1]); ret->repMatch = sqlUnsigned(row[2]); ret->nCount = sqlUnsigned(row[3]); ret->qNumInsert = sqlUnsigned(row[4]); ret->qBaseInsert = sqlSigned(row[5]); ret->tNumInsert = sqlUnsigned(row[6]); ret->tBaseInsert = sqlSigned(row[7]); strcpy(ret->strand, row[8]); ret->qName = lmCloneString(lm,row[9]); ret->qSize = sqlUnsigned(row[10]); ret->qStart = sqlUnsigned(row[11]); ret->qEnd = sqlUnsigned(row[12]); ret->tName = lmCloneString(lm, row[13]); ret->tSize = sqlUnsigned(row[14]); ret->tStart = sqlUnsigned(row[15]); ret->tEnd = sqlUnsigned(row[16]); ret->blockSizes = lmAlloc(lm, sizeof(ret->blockSizes[0]) * ret->blockCount); sqlUnsignedArray(row[18], ret->blockSizes, ret->blockCount); ret->qStarts = lmAlloc(lm, sizeof(ret->qStarts[0]) * ret->blockCount); sqlUnsignedArray(row[19], ret->qStarts, ret->blockCount); ret->tStarts = lmAlloc(lm, sizeof(ret->tStarts[0]) * ret->blockCount); sqlUnsignedArray(row[20], ret->tStarts, ret->blockCount); return ret; } boolean pslIsProtein(const struct psl *psl) /* is psl a protein psl (are it's blockSizes and scores in protein space) */ { int lastBlock = psl->blockCount - 1; return (((psl->strand[1] == '+' ) && (psl->tEnd == psl->tStarts[lastBlock] + 3*psl->blockSizes[lastBlock])) || ((psl->strand[1] == '-') && (psl->tStart == (psl->tSize-(psl->tStarts[lastBlock] + 3*psl->blockSizes[lastBlock]))))); } int pslCalcMilliBad(struct psl *psl, boolean isMrna) /* Calculate badness in parts per thousand. */ { int sizeMul = pslIsProtein(psl) ? 3 : 1; int qAliSize, tAliSize, aliSize; int milliBad = 0; int sizeDif; int insertFactor; int total; qAliSize = sizeMul * (psl->qEnd - psl->qStart); tAliSize = psl->tEnd - psl->tStart; aliSize = min(qAliSize, tAliSize); if (aliSize <= 0) return 0; sizeDif = qAliSize - tAliSize; if (sizeDif < 0) { if (isMrna) sizeDif = 0; else sizeDif = -sizeDif; } insertFactor = psl->qNumInsert; if (!isMrna) insertFactor += psl->tNumInsert; total = (sizeMul * (psl->match + psl->repMatch + psl->misMatch)); if (total != 0) milliBad = (1000 * (psl->misMatch*sizeMul + insertFactor + round(3*log(1+sizeDif)))) / total; return milliBad; } int pslScore(const struct psl *psl) /* Return score for psl. */ { int sizeMul = pslIsProtein(psl) ? 3 : 1; return sizeMul * (psl->match + ( psl->repMatch>>1)) - sizeMul * psl->misMatch - psl->qNumInsert - psl->tNumInsert; } int pslCmpScoreDesc(const void *va, const void *vb) /* Compare to sort based on score. */ { const struct psl *a = *((struct psl **)va); const struct psl *b = *((struct psl **)vb); return pslScore(b) - pslScore(a); } struct ffAli *pslToFakeFfAli(struct psl *psl, DNA *needle, DNA *haystack) /* Convert from psl to ffAli format. In some cases you can pass NULL * for needle and haystack - depending what the post-processing is going * to be. */ { struct ffAli *ffList = NULL, *ff; int blockCount = psl->blockCount; unsigned *blockSizes = psl->blockSizes; unsigned *qStarts = psl->qStarts; unsigned *tStarts = psl->tStarts; int size; int i; for (i=0; i<blockCount; ++i) { size = blockSizes[i]; AllocVar(ff); ff->left = ffList; ffList = ff; ff->nStart = ff->nEnd = needle + qStarts[i]; ff->nEnd += size; ff->hStart = ff->hEnd = haystack + tStarts[i]; ff->hEnd += size; } ffList = ffMakeRightLinks(ffList); return ffList; } struct psl *pslFromFakeFfAli(struct ffAli *ff, DNA *needle, DNA *haystack, char strand, char *qName, int qSize, char *tName, int tSize) /* This will create a basic psl structure from a sorted series of ffAli * blocks. The fields that would need actual sequence to be filled in * are left zero however - fields including match, repMatch, mismatch. */ { struct psl *psl; unsigned *blockSizes; unsigned *qStarts; unsigned *tStarts; int blockCount; int i; int nStart, hStart; int nEnd, hEnd; AllocVar(psl); psl->blockCount = blockCount = ffAliCount(ff); psl->blockSizes = AllocArray(blockSizes, blockCount); psl->qStarts = AllocArray(qStarts, blockCount); psl->tStarts = AllocArray(tStarts, blockCount); psl->qName = cloneString(qName); psl->qSize = qSize; psl->tName = cloneString(tName); psl->tSize = tSize; psl->strand[0] = strand; for (i=0; i<blockCount; ++i) { nStart = ff->nStart - needle; nEnd = ff->nEnd - needle; hStart = ff->hStart - haystack; hEnd = ff->hEnd - haystack; blockSizes[i] = nEnd - nStart; qStarts[i] = nStart; tStarts[i] = hStart; if (i == 0) { psl->qStart = nStart; psl->tStart = hStart; } if (i == blockCount-1) { psl->qEnd = nEnd; psl->tEnd = hEnd; } ff = ff->right; } if (strand == '-') { reverseIntRange(&psl->qStart, &psl->qEnd, psl->qSize); } return psl; } struct ffAli *pslToFfAli(struct psl *psl, struct dnaSeq *query, struct dnaSeq *target, int targetOffset) /* Convert from psl to ffAli format. Clip to parts that we actually * have sequence for. */ { struct ffAli *ffList = NULL, *ff; DNA *needle = query->dna; DNA *haystack = target->dna; int blockCount = psl->blockCount; unsigned *blockSizes = psl->blockSizes; unsigned *qStarts = psl->qStarts; unsigned *tStarts = psl->tStarts; int size; int i; int tMin = targetOffset; int tMax = targetOffset + target->size; int tStart, tEnd; int clipStart, clipEnd, clipOffset, clipSize; for (i=0; i<blockCount; ++i) { clipStart = tStart = tStarts[i]; size = blockSizes[i]; clipEnd = tEnd = tStart + size; if (tStart < tMax && tEnd > tMin) { if (clipStart < tMin) clipStart = tMin; if (clipEnd > tMax) clipEnd = tMax; clipOffset = clipStart - tStart; clipSize = clipEnd - clipStart; AllocVar(ff); ff->left = ffList; ffList = ff; ff->nStart = ff->nEnd = needle + qStarts[i] + clipOffset; ff->nEnd += clipSize; ff->hStart = ff->hEnd = haystack + clipStart - targetOffset; ff->hEnd += clipSize; } } ffList = ffMakeRightLinks(ffList); ffCountGoodEnds(ffList); return ffList; } int pslOrientation(struct psl *psl) /* Translate psl strand + or - to orientation +1 or -1 */ { if (psl->strand[1] != '\0') { /* translated blat */ if (psl->strand[0] != psl->strand[1]) return -1; else return 1; } else { if (psl->strand[0] == '-') return -1; else return 1; } } int pslWeightedIntronOrientation(struct psl *psl, struct dnaSeq *genoSeq, int offset) /* Return >0 if introns make it look like alignment is on + strand, * <0 if introns make it look like alignment is on - strand, * 0 if can't tell. The absolute value of the return indicates * how many splice sites we've seen supporting the orientation. * Sequence should NOT be reverse complemented. */ { int intronDir = 0; int oneDir; int i; DNA *dna = genoSeq->dna; /* code below doesn't support negative target strand (translated blat) */ if (psl->strand[1] == '-') errAbort("pslWeightedIntronOrientation doesn't support a negative target strand"); for (i=1; i<psl->blockCount; ++i) { int iStart, iEnd, blockSize = psl->blockSizes[i-1]; if (psl->qStarts[i-1] + blockSize == psl->qStarts[i]) { iStart = psl->tStarts[i-1] + psl->blockSizes[i-1] - offset; iEnd = psl->tStarts[i] - offset; oneDir = intronOrientation(dna+iStart, dna+iEnd); intronDir += oneDir; } } return intronDir; } int pslIntronOrientation(struct psl *psl, struct dnaSeq *genoSeq, int offset) /* Return 1 if introns make it look like alignment is on + strand, * -1 if introns make it look like alignment is on - strand, * 0 if can't tell. * Sequence should NOT be reverse complemented. */ { int intronDir = pslWeightedIntronOrientation(psl, genoSeq, offset); if (intronDir < 0) intronDir = -1; else if (intronDir > 0) intronDir = 1; return intronDir; } boolean pslHasIntron(struct psl *psl, struct dnaSeq *seq, int seqOffset) /* Return TRUE if there's a probable intron. Sequence should NOT be * reverse complemented.*/ { int blockCount = psl->blockCount, i; unsigned *tStarts = psl->tStarts; unsigned *blockSizes = psl->blockSizes; unsigned *qStarts = psl->qStarts; int blockSize, start, end; DNA *dna = seq->dna; for (i=1; i<blockCount; ++i) { blockSize = blockSizes[i-1]; start = qStarts[i-1]+blockSize; end = qStarts[i]; if (start == end) { start = tStarts[i-1] + blockSize; end = tStarts[i]; if (psl->strand[1] == '-') reverseIntRange(&start, &end, psl->tSize); start -= seqOffset; end -= seqOffset; if (intronOrientation(dna+start, dna+end) != 0) return TRUE; } } return FALSE; } void pslTailSizes(struct psl *psl, int *retStartTail, int *retEndTail) /* Find the length of "tails" (rather than extensions) implied by psl. */ { int orientation = pslOrientation(psl); int qFloppyStart, qFloppyEnd; int tFloppyStart, tFloppyEnd; if (orientation > 0) { qFloppyStart = psl->qStart; qFloppyEnd = psl->qSize - psl->qEnd; } else { qFloppyStart = psl->qSize - psl->qEnd; qFloppyEnd = psl->qStart; } tFloppyStart = psl->tStart; tFloppyEnd = psl->tSize - psl->tEnd; *retStartTail = min(qFloppyStart, tFloppyStart); *retEndTail = min(qFloppyEnd, tFloppyEnd); } static void rcSeqs(char **seqs, unsigned blockCount, unsigned *blockSizes) /* reverses complement sequences in list, maintain property that all strings * are in one malloc block. blockSizes should already be reversed. */ { char *buf, *next; int i, memSz = 0; /* get a new memory block for strings */ for (i = 0; i < blockCount; i++) memSz += blockSizes[i]+1; next = buf = needLargeMem(memSz); /* reverse compliment and copy to new memory block */ for (i = blockCount-1; i >= 0; i--) { int len = strlen(seqs[i]); reverseComplement(seqs[i], len); memcpy(next, seqs[i], len+1); next += len+1; } /* swap memory and update pointers */ freeMem(seqs[0]); seqs[0] = buf; next = buf; for (i = 0; i < blockCount; i++) { seqs[i] = next; next += blockSizes[i]+1; } } void pslRc(struct psl *psl) /* Reverse-complement a PSL alignment. This makes the target strand explicit. */ { unsigned tSize = psl->tSize, qSize = psl->qSize; unsigned blockCount = psl->blockCount, i; unsigned *tStarts = psl->tStarts, *qStarts = psl->qStarts, *blockSizes = psl->blockSizes; int mult = pslIsProtein(psl) ? 3 : 1; /* swap strand, forcing target to have an explict strand */ psl->strand[0] = (psl->strand[0] != '-') ? '-' : '+'; psl->strand[1] = (psl->strand[1] != '-') ? '-' : '+'; psl->strand[2] = 0; for (i=0; i<blockCount; ++i) { tStarts[i] = tSize - (tStarts[i] + mult * blockSizes[i]); qStarts[i] = qSize - (qStarts[i] + blockSizes[i]); } reverseUnsigned(tStarts, blockCount); reverseUnsigned(qStarts, blockCount); reverseUnsigned(blockSizes, blockCount); if (psl->qSequence != NULL) { rcSeqs(psl->qSequence, blockCount, blockSizes); rcSeqs(psl->tSequence, blockCount, blockSizes); } } /* macro to swap to variables */ #define swapVars(a, b, tmp) ((tmp) = (a), (a) = (b), (b) = (tmp)) static void swapBlocks(struct psl *psl) /* Swap the blocks in a psl without reverse complementing them. */ { int i; unsigned utmp; char *stmp; for (i = 0; i < psl->blockCount; i++) { swapVars(psl->qStarts[i], psl->tStarts[i], utmp); if (psl->qSequence != NULL) swapVars(psl->qSequence[i], psl->tSequence[i], stmp); } } static void swapRCBlocks(struct psl *psl) /* Swap and reverse complement blocks in a psl. Other psl fields must * be modified first */ { int i; unsigned *uatmp; char **satmp; reverseUnsigned(psl->tStarts, psl->blockCount); reverseUnsigned(psl->qStarts, psl->blockCount); reverseUnsigned(psl->blockSizes, psl->blockCount); swapVars(psl->tStarts, psl->qStarts, uatmp); /* qSize and tSize have already been swapped */ for (i = 0; i < psl->blockCount; i++) { psl->qStarts[i] = psl->qSize - (psl->qStarts[i] + psl->blockSizes[i]); psl->tStarts[i] = psl->tSize - (psl->tStarts[i] + psl->blockSizes[i]); } if (psl->qSequence != NULL) { /* note: all block sequences are stored in one malloc block, which is * entry zero */ rcSeqs(psl->qSequence, psl->blockCount, psl->blockSizes); rcSeqs(psl->tSequence, psl->blockCount, psl->blockSizes); swapVars(psl->qSequence, psl->tSequence, satmp); } } void pslSwap(struct psl *psl, boolean noRc) /* swap query and target in psl. If noRc is TRUE, don't reverse-complement * PSL if needed, instead make target strand explict. */ { int itmp; unsigned utmp; char ctmp, *stmp; swapVars(psl->qBaseInsert, psl->tBaseInsert, utmp); swapVars(psl->tNumInsert, psl->qNumInsert, utmp); swapVars(psl->qName, psl->tName, stmp); swapVars(psl->qSize, psl->tSize, utmp); swapVars(psl->qStart, psl->tStart, itmp); swapVars(psl->qEnd, psl->tEnd, itmp); /* handle strand and block copy */ if (psl->strand[1] != '\0') { /* translated */ swapVars(psl->strand[0], psl->strand[1], ctmp); swapBlocks(psl); } else if (noRc) { /* untranslated with no reverse complement */ psl->strand[1] = psl->strand[0]; psl->strand[0] = '+'; swapBlocks(psl); } else { /* untranslated */ if (psl->strand[0] == '+') swapBlocks(psl); else swapRCBlocks(psl); } } void pslTargetOffset(struct psl *psl, int offset) /* Add offset to target positions in psl. */ { int i, blockCount = psl->blockCount; unsigned *tStarts = psl->tStarts; psl->tStart += offset; psl->tEnd += offset; for (i=0; i<blockCount; ++i) tStarts[i] += offset; } void pslDump(struct psl *psl, FILE *f) /* Dump most of PSL to file - for debugging. */ { int i; fprintf(f, "<PRE>\n"); fprintf(f, "psl %s:%d-%d %s %s:%d-%d %d\n", psl->qName, psl->qStart, psl->qEnd, psl->strand, psl->tName, psl->tStart, psl->tEnd, psl->blockCount); for (i=0; i<psl->blockCount; ++i) fprintf(f, " size %d, qStart %d, tStart %d\n", psl->blockSizes[i], psl->qStarts[i], psl->tStarts[i]); fprintf(f, "</PRE>"); } void pslRecalcBounds(struct psl *psl) /* Calculate qStart/qEnd tStart/tEnd at top level to be consistent * with blocks. */ { int qStart, qEnd, tStart, tEnd, size; int last = psl->blockCount - 1; qStart = psl->qStarts[0]; tStart = psl->tStarts[0]; size = psl->blockSizes[last]; qEnd = psl->qStarts[last] + size; tEnd = psl->tStarts[last] + size; if (psl->strand[0] == '-') reverseIntRange(&qStart, &qEnd, psl->qSize); if (psl->strand[1] == '-') reverseIntRange(&tStart, &tEnd, psl->tSize); psl->qStart = qStart; psl->qEnd = qEnd; psl->tStart = tStart; psl->tEnd = tEnd; } void pslRecalcMatchCounts(struct psl *psl) /* Update the match/mismatch counts in PSL, assuming everything is a match. */ { psl->match = psl->misMatch = psl->repMatch = psl->nCount = 0; for (int iBlk = 0; iBlk < psl->blockCount; iBlk++) psl->match += psl->blockSizes[iBlk]; } struct psl *pslTrimToTargetRange(struct psl *oldPsl, int tMin, int tMax) /* Return psl trimmed to fit inside tMin/tMax. Note this does not * update the match/misMatch and related fields. */ { int newSize; int oldBlockCount = oldPsl->blockCount; boolean tIsRc = (oldPsl->strand[1] == '-'); int newBlockCount = 0, completeBlockCount = 0; int i; struct psl *newPsl = NULL; int tMn = tMin, tMx = tMax; /* tMin/tMax adjusted for strand. */ /* Deal with case where we're completely trimmed out quickly. */ newSize = rangeIntersection(oldPsl->tStart, oldPsl->tEnd, tMin, tMax); if (newSize <= 0) return NULL; if (tIsRc) reverseIntRange(&tMn, &tMx, oldPsl->tSize); /* Count how many blocks will survive trimming. */ oldBlockCount = oldPsl->blockCount; for (i=0; i<oldBlockCount; ++i) { int s = oldPsl->tStarts[i]; int e = s + oldPsl->blockSizes[i]; int sz = e - s; int overlap; if ((overlap = rangeIntersection(s, e, tMn, tMx)) > 0) ++newBlockCount; if (overlap == sz) ++completeBlockCount; } if (newBlockCount == 0) return NULL; /* Allocate new psl and fill in what we already know. */ AllocVar(newPsl); strcpy(newPsl->strand, oldPsl->strand); newPsl->qName = cloneString(oldPsl->qName); newPsl->qSize = oldPsl->qSize; newPsl->tName = cloneString(oldPsl->tName); newPsl->tSize = oldPsl->tSize; newPsl->blockCount = newBlockCount; AllocArray(newPsl->blockSizes, newBlockCount); AllocArray(newPsl->qStarts, newBlockCount); AllocArray(newPsl->tStarts, newBlockCount); /* Fill in blockSizes, qStarts, tStarts with real data. */ newBlockCount = completeBlockCount = 0; for (i=0; i<oldBlockCount; ++i) { int oldSz = oldPsl->blockSizes[i]; int sz = oldSz; int tS = oldPsl->tStarts[i]; int tE = tS + sz; int qS = oldPsl->qStarts[i]; int qE = qS + sz; if (rangeIntersection(tS, tE, tMn, tMx) > 0) { int diff; if ((diff = (tMn - tS)) > 0) { tS += diff; qS += diff; sz -= diff; } if ((diff = (tE - tMx)) > 0) { tE -= diff; qE -= diff; sz -= diff; } newPsl->qStarts[newBlockCount] = qS; newPsl->tStarts[newBlockCount] = tS; newPsl->blockSizes[newBlockCount] = sz; ++newBlockCount; if (sz == oldSz) ++completeBlockCount; } } pslRecalcBounds(newPsl); return newPsl; } struct psl *pslTrimToQueryRange(struct psl *oldPsl, int qMin, int qMax) /* Return psl trimmed to fit inside qMin/qMax. Note this does not * update the match/misMatch and related fields. */ { int newSize; int oldBlockCount = oldPsl->blockCount; boolean qIsRc = (oldPsl->strand[0] == '-'); int newBlockCount = 0, completeBlockCount = 0; int i; struct psl *newPsl = NULL; int qMn = qMin, qMx = qMax; /* qMin/qMax adjusted for strand. */ /* Deal with case where we're completely trimmed out quickly. */ newSize = rangeIntersection(oldPsl->qStart, oldPsl->qEnd, qMin, qMax); if (newSize <= 0) return NULL; if (qIsRc) reverseIntRange(&qMn, &qMx, oldPsl->qSize); /* Count how many blocks will survive trimming. */ oldBlockCount = oldPsl->blockCount; for (i=0; i<oldBlockCount; ++i) { int s = oldPsl->qStarts[i]; int e = s + oldPsl->blockSizes[i]; int sz = e - s; int overlap; if ((overlap = rangeIntersection(s, e, qMn, qMx)) > 0) ++newBlockCount; if (overlap == sz) ++completeBlockCount; } if (newBlockCount == 0) return NULL; /* Allocate new psl and fill in what we already know. */ AllocVar(newPsl); strcpy(newPsl->strand, oldPsl->strand); newPsl->qName = cloneString(oldPsl->qName); newPsl->qSize = oldPsl->qSize; newPsl->tName = cloneString(oldPsl->tName); newPsl->tSize = oldPsl->tSize; newPsl->blockCount = newBlockCount; AllocArray(newPsl->blockSizes, newBlockCount); AllocArray(newPsl->qStarts, newBlockCount); AllocArray(newPsl->tStarts, newBlockCount); /* Fill in blockSizes, qStarts, tStarts with real data. */ newBlockCount = completeBlockCount = 0; for (i=0; i<oldBlockCount; ++i) { int oldSz = oldPsl->blockSizes[i]; int sz = oldSz; int qS = oldPsl->qStarts[i]; int qE = qS + sz; int tS = oldPsl->tStarts[i]; int tE = tS + sz; if (rangeIntersection(qS, qE, qMn, qMx) > 0) { int diff; if ((diff = (qMn - qS)) > 0) { tS += diff; qS += diff; sz -= diff; } if ((diff = (qE - qMx)) > 0) { tE -= diff; qE -= diff; sz -= diff; } newPsl->qStarts[newBlockCount] = qS; newPsl->tStarts[newBlockCount] = tS; newPsl->blockSizes[newBlockCount] = sz; ++newBlockCount; if (sz == oldSz) ++completeBlockCount; } } pslRecalcBounds(newPsl); return newPsl; } static void printPslDesc(char* pslDesc, FILE* out, struct psl* psl) /* print description of a PSL on first error */ { fprintf(out, "Error: invalid PSL: %s:%u-%u %s:%u-%u %s %s\n", psl->qName, psl->qStart, psl->qEnd, psl->tName, psl->tStart, psl->tEnd, psl->strand, pslDesc); } static void chkError(char* pslDesc, FILE* out, struct psl* psl, int* errCount, char* format, ...) /* forward needed to specify printf signature for gcc checking */ #if defined(__GNUC__) __attribute__((format(printf, 5, 6))) #endif ; static void chkError(char* pslDesc, FILE* out, struct psl* psl, int* errCount, char* format, ...) /* error handling on an pslCheck error, counting error and issuing description * of PSL on the first error. */ { if (*errCount == 0) printPslDesc(pslDesc, out, psl); va_list args; va_start(args, format); vfprintf(out, format, args); va_end(args); (*errCount)++; } static void chkBlkRanges(char* pslDesc, FILE* out, struct psl* psl, char* pName, char* pLabel, char pCLabel, char pStrand, unsigned pSize, unsigned pStart, unsigned pEnd, unsigned iBlk, unsigned* pBlockStarts, int* errCount) /* check the target or query ranges in a PSL incrementing errorCnt */ { unsigned blkStart = pBlockStarts[iBlk]; unsigned blkEnd = blkStart+psl->blockSizes[iBlk]; /* translate stand to genomic coords */ unsigned gBlkStart = (pStrand == '+') ? blkStart : (pSize - blkEnd); unsigned gBlkEnd = (pStrand == '+') ? blkEnd : (pSize - blkStart); if ((pSize > 0) && (blkEnd > pSize)) chkError(pslDesc, out, psl, errCount, "\t%s %s block %u end %u > %cSize %u\n", pName, pLabel, iBlk, blkEnd, pCLabel, pSize); if (gBlkStart < pStart) chkError(pslDesc, out, psl, errCount, "\t%s %s block %u start %u < %cStart %u\n", pName, pLabel, iBlk, gBlkStart, pCLabel, pStart); if (gBlkStart >= pEnd) chkError(pslDesc, out, psl, errCount, "\t%s %s block %u start %u >= %cEnd %u\n", pName, pLabel, iBlk, gBlkStart, pCLabel, pEnd); if (gBlkEnd < pStart) chkError(pslDesc, out, psl, errCount, "\t%s %s block %u end %u < %cStart %u\n", pName, pLabel, iBlk, gBlkEnd, pCLabel, pStart); if (gBlkEnd > pEnd) chkError(pslDesc, out, psl, errCount, "\t%s %s block %u end %u > %cEnd %u\n", pName, pLabel, iBlk, gBlkEnd, pCLabel, pEnd); if (iBlk > 0) { unsigned prevBlkEnd = pBlockStarts[iBlk-1]+psl->blockSizes[iBlk-1]; if (blkStart < prevBlkEnd) chkError(pslDesc, out, psl, errCount, "\t%s %s block %u start %u < previous block end %u\n", pName, pLabel, iBlk, blkStart, prevBlkEnd); } } static void chkRanges(char* pslDesc, FILE* out, struct psl* psl, char* pName, char* pLabel, char pCLabel, char pStrand, unsigned pSize, unsigned pStart, unsigned pEnd, unsigned* pBlockStarts, int blockSizeMult, int* errCount) /* check the target or query ranges in a PSL, increment errorCnt */ { unsigned iBlk; if (pStart >= pEnd) chkError(pslDesc, out, psl, errCount, "\t%s %cStart %u >= %cEnd %u\n", pName, pCLabel, pStart, pCLabel, pEnd); if (pEnd > pSize) chkError(pslDesc, out, psl, errCount, "\t%s %cEnd %u >= %cSize %u\n", pName, pCLabel, pEnd, pCLabel, pSize); // check that block start/end matches overall start end unsigned pStartStrand = pStart, pEndStrand = pEnd; if (pStrand != '+') reverseUnsignedRange(&pStartStrand, &pEndStrand, pSize); unsigned lastBlkEnd = pBlockStarts[psl->blockCount-1] + (blockSizeMult * psl->blockSizes[psl->blockCount-1]); if ((pStartStrand != pBlockStarts[0]) || (pEndStrand != lastBlkEnd)) chkError(pslDesc, out, psl, errCount, "\t%s strand \"%c\" adjusted %cStart-%cEnd range %u-%u != block range %u-%u\n", pName, pStrand, pCLabel, pCLabel, pStartStrand, pEndStrand, pBlockStarts[0], lastBlkEnd); for (iBlk = 0; iBlk < psl->blockCount; iBlk++) chkBlkRanges(pslDesc, out, psl, pName, pLabel, pCLabel, pStrand, pSize, pStart, pEnd, iBlk, pBlockStarts, errCount); } static void chkInsertCounts(char* pslDesc, FILE* out, struct psl* psl, char* pName, char pCLabel, unsigned* pBlockStarts, unsigned pNumInsert, unsigned pBaseInsert, int* errCount) /* check the insert counts, incrementing errorCnt */ { unsigned numInsert = 0, baseInsert = 0; int iBlk; for (iBlk = 1; iBlk < psl->blockCount; iBlk++) { unsigned gapSize = pBlockStarts[iBlk] - (pBlockStarts[iBlk-1]+psl->blockSizes[iBlk-1]); if (gapSize > 0) { numInsert++; baseInsert += gapSize; } } if (numInsert != pNumInsert) chkError(pslDesc, out, psl, errCount, "\t%s %cNumInsert %u != expected %u\n", pName, pCLabel, pNumInsert, numInsert); if (baseInsert != pBaseInsert) chkError(pslDesc, out, psl, errCount, "\t%s %cBaseInsert %u != expected %u\n", pName, pCLabel, pBaseInsert, baseInsert); } int pslCheck2(unsigned opts, char *pslDesc, FILE* out, struct psl* psl) /* Validate a PSL for consistency. pslDesc is printed the error messages to * file out (open /dev/null to discard). Return count of errors. Option * PSL_CHECK_IGNORE_INSERT_CNTS doesn't validate problems insert counts fields * in each PSL. Useful because protein PSL doesn't seen to compute these in a * consistent way. */ { static char* VALID_STRANDS[] = { "+", "-", "++", "+-", "-+", "--", NULL }; int i, errCount = 0; int tBlockSizeMult = pslIsProtein(psl) ? 3 : 1; /* check strand value */ for (i = 0; VALID_STRANDS[i] != NULL; i++) { if (strcmp(psl->strand, VALID_STRANDS[i]) == 0) break; } if (VALID_STRANDS[i] == NULL) chkError(pslDesc, out, psl, &errCount, "\tinvalid PSL strand: \"%s\"\n", psl->strand); /* check query */ chkRanges(pslDesc, out, psl, psl->qName, "query", 'q', pslQStrand(psl), psl->qSize, psl->qStart, psl->qEnd, psl->qStarts, 1, &errCount); if ((opts & PSL_CHECK_IGNORE_INSERT_CNTS) == 0) chkInsertCounts(pslDesc, out, psl, psl->qName, 'q', psl->qStarts, psl->qNumInsert, psl->qBaseInsert, &errCount); /* check target */ chkRanges(pslDesc, out, psl, psl->tName, "target", 't', pslTStrand(psl), psl->tSize, psl->tStart, psl->tEnd, psl->tStarts, tBlockSizeMult, &errCount); if ((opts & PSL_CHECK_IGNORE_INSERT_CNTS) == 0) chkInsertCounts(pslDesc, out, psl, psl->tName, 't', psl->tStarts, psl->tNumInsert, psl->tBaseInsert, &errCount); return errCount; } int pslCheck(char *pslDesc, FILE* out, struct psl* psl) /* Validate a PSL for consistency. pslDesc is printed the error messages * to file out (open /dev/null to discard). Return count of errors. */ { return pslCheck2(0, pslDesc, out, psl); } struct hash *readPslToBinKeeper(char *sizeFileName, char *pslFileName) /* read a list of psls and return results in hash of binKeeper structure for fast query*/ { struct binKeeper *bk; struct psl *psl; struct lineFile *sf = lineFileOpen(sizeFileName, TRUE); struct lineFile *pf = lineFileOpen(pslFileName , TRUE); struct hash *hash = newHash(0); char *chromRow[2]; char *row[21] ; while (lineFileRow(sf, chromRow)) { char *name = chromRow[0]; int size = lineFileNeedNum(sf, chromRow, 1); if (hashLookup(hash, name) != NULL) warn("Duplicate %s, ignoring all but first\n", name); else { bk = binKeeperNew(0, size); assert(size > 1); hashAdd(hash, name, bk); } } while (lineFileNextRow(pf, row, ArraySize(row))) { psl = pslLoad(row); bk = hashMustFindVal(hash, psl->tName); binKeeperAdd(bk, psl->tStart, psl->tEnd, psl); } lineFileClose(&pf); lineFileClose(&sf); return hash; } static boolean isDelChar(char c) /* is this a indel character? */ { return (c == '-') || (c == '.') || (c == '=') || (c == '_'); } static void trimAlignment(struct psl* psl, char** qStrPtr, char** tStrPtr, int* aliSizePtr) /* remove leading or trailing indels from alignment */ { char* qStr = *qStrPtr; char* tStr = *tStrPtr; int aliSize = *aliSizePtr; // skip lending indels while ((aliSize > 0) && (isDelChar(*qStr) || isDelChar(*tStr))) { if (!isDelChar(*qStr)) psl->qStart++; else if (!isDelChar(*tStr)) psl->tStart++; qStr++; tStr++; aliSize--; } // skip trailing indels while ((aliSize > 0) && (isDelChar(qStr[aliSize-1]) || isDelChar(tStr[aliSize-1]))) { if (!isDelChar(qStr[aliSize-1])) psl->qEnd--; else if (!isDelChar(tStr[aliSize-1])) psl->tEnd--; aliSize--; } *qStrPtr = qStr; *tStrPtr = tStr; *aliSizePtr = aliSize; } static void addBlock(struct psl* psl, int qs, int qe, int ts, int te, int *blockSpace) /* add a block to the psl, growing if necessary */ { assert((qe-qs) == (te-ts)); /* same lengths? */ if (psl->blockCount == *blockSpace) pslGrow(psl, blockSpace); psl->blockSizes[psl->blockCount] = qe - qs; psl->qStarts[psl->blockCount] = qs; psl->tStarts[psl->blockCount] = ts; psl->blockCount++; } static void accumCounts(struct psl *psl, char prevQ, char prevT, char q, char t, unsigned options) /* accumulate block and base counts */ { if (!isDelChar(q) && !isDelChar(t)) { /* aligned column. */ char qu = toupper(q); char tu = toupper(t); if ((q == 'N') || (t == 'N')) psl->nCount++; else if (qu == tu) { if ((options & PSL_IS_SOFTMASK) && ((qu != q) || (tu != t))) psl->repMatch++; else psl->match++; } else psl->misMatch++; } else if (isDelChar(q) && !isDelChar(t)) { /* target insert */ psl->tBaseInsert++; if (!isDelChar(prevQ)) psl->tNumInsert++; } else if (isDelChar(t) && !isDelChar(q)) { /* query insert */ psl->qBaseInsert++; if (!isDelChar(prevT)) psl->qNumInsert++; } } struct psl* pslFromAlign(char *qName, int qSize, int qStart, int qEnd, char *qString, char *tName, int tSize, int tStart, int tEnd, char *tString, char* strand, unsigned options) /* Create a PSL from an alignment. Options PSL_IS_SOFTMASK if lower case * bases indicate repeat masking. Returns NULL if alignment is empty after * triming leading and trailing indels.*/ { /* Note, the unit tests for these programs exercise this function: * hg/embossToPsl * hg/mouseStuff/axtToPsl * hg/spideyToPsl * utils/est2genomeToPsl */ int blockSpace = 16; struct psl* psl = NULL; int aliSize = strlen(qString); boolean eitherInsert = FALSE; /* True if either in insert state. */ int i, qs,qe,ts,te; char prevQ = '\0', prevT = '\0'; AllocVar(psl); if (strlen(tString) != aliSize) errAbort("query and target alignment strings are different lengths"); psl = pslNew(qName, qSize, qStart, qEnd, tName, tSize, tStart, tEnd, strand, blockSpace, 0); trimAlignment(psl, &qString, &tString, &aliSize); /* Don't create if either query or target is zero length after trim */ if ((psl->qStart == psl->qEnd) || (psl->tStart == psl->tEnd)) { pslFree(&psl); return NULL; } /* Get range of alignment in strand-specified coordinates */ qs = psl->qStart; qe = psl->qEnd; if (strand[0] == '-') reverseIntRange(&qs, &qe, psl->qSize); ts = psl->tStart; te = psl->tEnd; if (strand[1] == '-') reverseIntRange(&ts, &te, psl->tSize); eitherInsert = FALSE; qe = qs; /* current block coords */ te = ts; for (i=0; i<aliSize; ++i) { char q = qString[i]; char t = tString[i]; if (isDelChar(q) && isDelChar(t)) { continue; /* nothing in this column, just ignore it */ } else if (isDelChar(q) || isDelChar(t)) { /* insert in one of the columns */ if (!eitherInsert) { /* end of a block */ addBlock(psl, qs, qe, ts, te, &blockSpace); eitherInsert = TRUE; } if (!isDelChar(q)) qe += 1; if (!isDelChar(t)) te += 1; } else { /* columns aligned */ if (eitherInsert) { /* start new block */ qs = qe; ts = te; eitherInsert = FALSE; } qe += 1; te += 1; } accumCounts(psl, prevQ, prevT, q, t, options); prevQ = q; /* will not include skipped empty columns */ prevT = t; } addBlock(psl, qs, qe, ts, te, &blockSpace); return psl; } struct psl* pslNew(char *qName, unsigned qSize, int qStart, int qEnd, char *tName, unsigned tSize, int tStart, int tEnd, char *strand, unsigned blockSpace, unsigned opts) /* create a new psl with space for the specified number of blocks allocated. * pslGrow maybe used to expand this space if needed. Valid options are * PSL_XA_FORMAT. */ { struct psl *psl; AllocVar(psl); assert(blockSpace > 0); psl->qName = cloneString(qName); psl->qSize = qSize; psl->qStart = qStart; psl->qEnd = qEnd; psl->tName = cloneString(tName); psl->tSize = tSize; psl->tStart = tStart; psl->tEnd = tEnd; strncpy(psl->strand, strand, 2); AllocArray(psl->blockSizes, blockSpace); AllocArray(psl->qStarts, blockSpace); AllocArray(psl->tStarts, blockSpace); if (opts & PSL_XA_FORMAT) { AllocArray(psl->qSequence, blockSpace); AllocArray(psl->tSequence, blockSpace); } return psl; } void pslGrow(struct psl *psl, int *blockSpacePtr) /* Increase memory allocated to a psl to hold more blocks. blockSpacePtr * should point the the current maximum number of blocks and will be * updated to with the new amount of space. */ { int blockSpace = *blockSpacePtr; int newSpace = 2 * blockSpace; ExpandArray(psl->blockSizes, blockSpace, newSpace); ExpandArray(psl->qStarts, blockSpace, newSpace); ExpandArray(psl->tStarts, blockSpace, newSpace); if (psl->qSequence != NULL) { ExpandArray(psl->qSequence, blockSpace, newSpace); ExpandArray(psl->tSequence, blockSpace, newSpace); } *blockSpacePtr = newSpace; } void pslComputeInsertCounts(struct psl *psl) /* compute numInsert and baseInsert fields from the blocks */ { psl->qNumInsert = psl->qBaseInsert = 0; psl->tNumInsert = psl->tBaseInsert = 0; int iBlk; for (iBlk = 1; iBlk < psl->blockCount; iBlk++) { unsigned qGapSize = psl->qStarts[iBlk] - (psl->qStarts[iBlk-1]+psl->blockSizes[iBlk-1]); if (qGapSize != 0) { psl->qNumInsert++; psl->qBaseInsert += qGapSize; } unsigned tGapSize = psl->tStarts[iBlk] - (psl->tStarts[iBlk-1]+psl->blockSizes[iBlk-1]); if (tGapSize != 0) { psl->tNumInsert++; psl->tBaseInsert += tGapSize; } } } static boolean getNextCigarOp(char *startPtr, boolean reverse, char **ptr, char *op, int *size) /* gets one cigar op out of the CIGAR string. Reverse the order if asked */ { char *str = *ptr; if (str == NULL) return FALSE; if ((!reverse && (*str == 0)) || (reverse && (str == startPtr))) return FALSE; // between each cigar op there could be nothing, or a space, or a plus if (reverse) { char *end = str - 1; for(;*end ; end--) { if (isalpha(*end)) break; } str = end; *ptr = end; } else { char *end = str + 1; for(;*end ; end++) { if (! (isdigit(*end) || (*end == ' ') || (*end == '+'))) break; } *ptr = end; } *op = *str++; *size = atoi(str); return TRUE; } struct psl* pslFromGff3Cigar(char *qName, int qSize, int qStart, int qEnd, char *tName, int tSize, int tStart, int tEnd, char* strand, char *cigar) /* create a PSL from a GFF3-style cigar formatted alignment */ { // this is currently tested by the gff3ToPsl int blocksAlloced = 4; struct psl *psl = pslNew(qName, qSize, qStart, qEnd, tName, tSize, tStart, tEnd, strand, blocksAlloced, 0); char op; int size; int totalSize = 0; int qNext = qStart, qBlkEnd = qEnd; if (strand[0] == '-') reverseIntRange(&qNext, &qBlkEnd, qSize); int tNext = tStart, tBlkEnd = tEnd; if (strand[1] == '-') reverseIntRange(&tNext, &tBlkEnd, tSize); if (cigar == NULL) { // no cigar means one block size = qEnd - qStart; totalSize += size; psl->blockSizes[psl->blockCount] = size; psl->qStarts[psl->blockCount] = qNext; psl->tStarts[psl->blockCount] = tNext; psl->blockCount++; tNext += size; qNext += size; } else { if (strand[0] == '-' && strand[1] == '-') errAbort("GFF3 spec is vague about Gap when both strands are '-'; not implemented yet."); char cigarSpec[strlen(cigar)+1]; // copy since parsing is destructive safecpy(cigarSpec, sizeof cigarSpec, cigar); char *cigarNext = cigarSpec; while(getNextCigarOp(cigarSpec, FALSE, &cigarNext, &op, &size)) { switch (op) { case 'M': // match or mismatch (gapless aligned block) if (psl->blockCount == blocksAlloced) pslGrow(psl, &blocksAlloced); totalSize += size; psl->blockSizes[psl->blockCount] = size; psl->qStarts[psl->blockCount] = qNext; psl->tStarts[psl->blockCount] = tNext; psl->blockCount++; tNext += size; qNext += size; break; case 'I': // inserted in target tNext += size; break; case 'D': // deleted from target qNext += size; break; default: errAbort("unrecognized CIGAR op %c in %s", op, cigar); } } } /* CIGARs starting/ending with indels require adjusting of query/target ranges, * as PSL starts/ends with matches */ psl->qStart = psl->qStarts[0]; psl->qEnd = pslQEnd(psl, psl->blockCount-1); if (strand[0] == '-') reverseIntRange(&psl->qStart, &psl->qEnd, qSize); psl->tStart = psl->tStarts[0]; psl->tEnd = pslTEnd(psl, psl->blockCount-1); if (strand[1] == '-') reverseIntRange(&psl->tStart, &psl->tEnd, tSize); /* sanity check */ if (qNext != qBlkEnd) errAbort("CIGAR query length does not match specified query range %s:%d-%d", qName, qStart, qEnd); if (tNext != tBlkEnd) errAbort("CIGAR target length does not match specified target range %s:%d-%d", tName, tStart, tEnd); psl->match = totalSize; pslComputeInsertCounts(psl); return psl; } int pslRangeTreeOverlap(struct psl *psl, struct rbTree *rangeTree) /* Return amount that psl overlaps (on target side) with rangeTree. */ { int i; int overlap = 0; boolean isRc = (psl->strand[1] == '-'); for (i=0; i<psl->blockCount; ++i) { int start = psl->tStarts[i]; int end = start + psl->blockSizes[i]; if (isRc) reverseIntRange(&start, &end, psl->tSize); overlap += rangeTreeOverlapSize(rangeTree, start, end); } return overlap; } float pslIdent(struct psl *psl) /* computer fraction identity */ { float aligned = psl->match + psl->misMatch + psl->repMatch; if (aligned == 0.0) return 0.0; else return ((float)(psl->match + psl->repMatch))/aligned; } float pslQueryAligned(struct psl *psl) /* compute fraction of query that was aligned */ { float aligned = psl->match + psl->misMatch + psl->repMatch; return aligned/(float)psl->qSize; } struct psl* pslClone(struct psl *psl) /* clone a psl */ { struct psl* pslCp = pslNew(psl->qName, psl->qSize, psl->qStart, psl->qEnd, psl->tName, psl->tSize, psl->tStart, psl->tEnd, psl->strand, psl->blockCount, ((psl->tSequence != NULL) ? PSL_XA_FORMAT : 0)); pslCp->match = psl->match; pslCp->misMatch = psl->misMatch; pslCp->repMatch = psl->repMatch; pslCp->nCount = psl->nCount; pslCp->qNumInsert = psl->qNumInsert; pslCp->qBaseInsert = psl->qBaseInsert; pslCp->tNumInsert = psl->tNumInsert; pslCp->tBaseInsert = psl->tBaseInsert; int iBlk; for (iBlk = 0; iBlk < psl->blockCount; iBlk++) { pslCp->blockSizes[iBlk] = psl->blockSizes[iBlk]; pslCp->qStarts[iBlk] = psl->qStarts[iBlk]; pslCp->tStarts[iBlk] = psl->tStarts[iBlk]; if (psl->qSequence != NULL) pslCp->qSequence[iBlk] = cloneString(psl->qSequence[iBlk]); if (psl->tSequence != NULL) pslCp->tSequence[iBlk] = cloneString(psl->tSequence[iBlk]); pslCp->blockCount++; } return pslCp; } int cmpChrom(char *a, char *b) /* Compare two chromosomes. */ { return cmpStringsWithEmbeddedNumbers(a, b); } int pslCmpTargetScore(const void *va, const void *vb) /* Compare to sort based on target then score. */ { const struct psl *a = *((struct psl **)va); const struct psl *b = *((struct psl **)vb); int diff = cmpChrom(a->tName, b->tName); if (diff == 0) diff = pslScore(b) - pslScore(a); return diff; } int pslCmpTargetStart(const void *va, const void *vb) /* Compare to sort based on target start. */ { const struct psl *a = *((struct psl **)va); const struct psl *b = *((struct psl **)vb); int diff = cmpChrom(a->tName, b->tName); if (diff == 0) diff = a->tStart - b->tStart; return diff; } char *pslSortList[] = {"query,score", "query,start", "chrom,score", "chrom,start", "score"}; void pslSortListByVar(struct psl **pslList, char *sort) /* Sort a list of psls using the method definied in the sort string. */ { if (sameString(sort, "query,start")) { slSort(pslList, pslCmpQuery); } else if (sameString(sort, "query,score")) { slSort(pslList, pslCmpQueryScore); } else if (sameString(sort, "score")) { slSort(pslList, pslCmpScore); } else if (sameString(sort, "chrom,start")) { slSort(pslList, pslCmpTargetStart); } else if (sameString(sort, "chrom,score")) { slSort(pslList, pslCmpTargetScore); } else { slSort(pslList, pslCmpQueryScore); } } void pslRemoveFrameShifts(struct psl *psl) /* Remove any frameshits if present. Changes in place, doesn't update statistics in first nine fields. */ { unsigned prevEnd = psl->tStarts[0]; int ii; for(ii = 0; ii < psl->blockCount; ii++) { int diff = prevEnd - psl->tStarts[ii]; if (diff > 0) { if (diff > psl->blockSizes[ii]) errAbort("frame shift (%d) larger than block size (%d)", diff, psl->blockSizes[ii]); psl->blockSizes[ii] -= diff; psl->tStarts[ii] += diff; psl->qStarts[ii] += diff; } prevEnd = psl->tStarts[ii] + psl->blockSizes[ii]; } } + +struct psl *pslProtFromNaLike(struct psl *psl, int protSize, int *retDroppedAa) +/* An "na-like" protein psl counts its query in bases, three to a residue, rather than in residues. + * That is what pslProtCnv writes, and what the UniProt alignments carry: their proteins reach the + * genome through transcripts, so they come out in transcript coordinates. Return a copy counted + * in residues, which is what the protein alignment display expects, or NULL if the psl is not that + * shape. A minus-strand query is left alone: its coordinates are reverse complemented, but qSize + * is a multiple of three, so codon boundaries survive that and the arithmetic below is the same. + * + * Splicing ends a block wherever an exon ends, and about one codon in a hundred straddles one. + * Such a residue belongs to two blocks at once, in two places in the genome, so no single column + * can hold it: blocks are trimmed to whole codons and those residues are left out, counted in + * retDroppedAa so the page can say how many are missing rather than quietly dropping them. */ +{ +if (retDroppedAa != NULL) + *retDroppedAa = 0; +if (protSize <= 0 || psl->qSize != 3 * protSize || psl->blockCount < 1) + return NULL; +/* The query being three times the protein is not enough on its own: a file that still declares + * seqType amino_acid reaches here with its block sizes already divided by three and its query + * coordinates left in bases, which looks the same by that test but is in two units at once. + * Require the blocks to measure the target the way the target is measured, which is true of a + * psl whose query side is in bases and false once pslFromBigPsl() has divided them. */ +int last = psl->blockCount - 1; +if (psl->tStarts[last] + psl->blockSizes[last] != psl->tEnd) + return NULL; + +/* Cloned rather than built with pslNew(), which keeps the names, strand and counts without + * restating them, and gives an array big enough for the blocks that survive trimming. */ +struct psl *prot = pslClone(psl); +prot->qSize = protSize; +prot->qStart = psl->qStart / 3; +prot->qEnd = (psl->qEnd + 2) / 3; +int i, j = 0, alignedBases = 0, alignedAa = 0; +for (i = 0; i < psl->blockCount; ++i) + { + int qs = psl->qStarts[i], sz = psl->blockSizes[i]; + int firstAa = (qs + 2) / 3; /* first codon that begins at or after the block start */ + int endAa = (qs + sz) / 3; /* one past the last codon that ends at or before its end */ + alignedBases += sz; + if (endAa <= firstAa) + continue; /* the block is too short to hold a whole codon */ + prot->qStarts[j] = firstAa; + prot->blockSizes[j] = endAa - firstAa; + prot->tStarts[j] = psl->tStarts[i] + (firstAa * 3 - qs); + alignedAa += prot->blockSizes[j]; + ++j; + } +prot->blockCount = j; +if (j == 0) + { + pslFree(&prot); + return NULL; + } +/* Trimming can pull the first and last block in off the ends, so restate the bounds from the + * blocks that are left. A protein psl is recognised by its last block reaching tEnd three bases + * to a residue, and pslRc() below reads that to decide whether to step the target in codons, so + * these have to agree before anything else looks at them. */ +prot->qStart = prot->qStarts[0]; +prot->qEnd = prot->qStarts[j-1] + prot->blockSizes[j-1]; +prot->tStart = prot->tStarts[0]; +prot->tEnd = prot->tStarts[j-1] + 3 * prot->blockSizes[j-1]; + +/* The counts are in bases in the psl we were handed, and the page prints them next to a query + * measured in residues, so bring them into the same units. */ +prot->match = psl->match / 3; +prot->misMatch = psl->misMatch / 3; +prot->repMatch = psl->repMatch / 3; +prot->nCount = psl->nCount / 3; + +/* A bigPsl keeps the reference strand in the field a psl reads as the query strand, so an + * alignment to the minus strand arrives claiming the protein is reversed. Proteins have no + * strand; it is the genome that is reverse complemented here, and a reader that believes + * otherwise reverse complements the protein and prints nucleotide ambiguity codes. pslRc() + * moves the alignment to the convention where the query runs forward and the target carries the + * strand, which is the shape blat gives a protein alignment and what the display expects. */ +if (prot->strand[0] == '-') + pslRc(prot); + +pslComputeInsertCounts(prot); +if (retDroppedAa != NULL) + *retDroppedAa = alignedBases / 3 - alignedAa; +return prot; +}