2db6bab8db05291ff90bbc867784be7f7c9177f6 max Sat Jul 25 21:48:14 2026 -0700 hgLoadMafSummary: robustly split MAF sequence names into assembly and chrom Handle a pipe separator, no-dot names, GenArk accession dbs (GCA_/GCF_, which carry an accession.version dot), and a dotted chrom, without aborting on an ordinary assembly.chrom name that sits next to dotted query names. Fixes a regression from ae63ce5 (refs #36592). refs #37928 diff --git src/hg/makeDb/hgLoadMaf/hgLoadMafSummary.c src/hg/makeDb/hgLoadMaf/hgLoadMafSummary.c index 9e0355fa181..dc3c3ba685d 100644 --- src/hg/makeDb/hgLoadMaf/hgLoadMafSummary.c +++ src/hg/makeDb/hgLoadMaf/hgLoadMafSummary.c @@ -1,340 +1,362 @@ /* hgLoadMafSummary - Load a summary table of pairs in a maf into a database. */ /* Copyright (C) 2011 The Regents of the University of California * See kent/LICENSE or http://genome.ucsc.edu/license/ for licensing information. */ #include "common.h" #include "cheapcgi.h" #include "linefile.h" #include "hash.h" #include "options.h" #include "jksql.h" #include "hdb.h" #include "hgRelate.h" #include "portable.h" #include "maf.h" #include "dystring.h" #include "mafSummary.h" /* command line option specifications */ static struct optionSpec optionSpecs[] = { {"mergeGap", OPTION_INT}, {"minSize", OPTION_INT}, {"maxSize", OPTION_INT}, {"minSeqSize", OPTION_INT}, {"test", OPTION_BOOLEAN}, {NULL, 0} }; boolean test = FALSE; int mergeGap = 500; int minSize = 10000; int maxSize = 50000; int minSeqSize = 1000000; char *database = NULL; long summaryCount = 0; void usage() /* Explain usage and exit. */ { errAbort( "hgLoadMafSummary - Load a summary table of pairs in a maf into a database\n" "usage:\n" " hgLoadMafSummary database table file.maf\n" "options:\n" " -mergeGap=N max size of gap to merge regions (default %d)\n" " -minSize=N merge blocks smaller than N (default %d)\n" " -maxSize=N break up blocks larger than N (default %d)\n" " -minSeqSize=N skip alignments when reference sequence is less than N\n" " (default %d -- match with hgTracks min window size for\n" " using summary table)\n" " -test suppress loading the database. Just create .tab file(s)\n" " in current dir.\n", mergeGap, minSize, maxSize, minSeqSize ); } double scorePairwise(struct mafAli *maf) /* generate score from 0.0 to 1.0 for an alignment pair */ /* Adapted from multiz scoring in hgTracks/mafTrack.c */ { int endB; /* end in the reference (master) genome */ int deltaB; int endT; /* end in the master genome maf text (includes gaps) */ double score; double minScore = -100.0, maxScore = 100.0; double scoreScale = 1.0 / (maxScore - minScore); struct mafComp *mcMaster = maf->components; endB = mcMaster->size; deltaB = endB; for (endT = 0; endT < maf->textSize; endT++) { if (deltaB <= 0) break; if (mcMaster->text[endT] != '-') deltaB -= 1; } /* Take the score over the relevant range of text symbols in the maf, * and divide it by the bases we cover in the master genome to * get a normalized by base score. */ score = mafScoreRangeMultiz(maf, 0, endT)/endB; /* Scale the score so that it is between 0 and 1 */ score = (score - minScore) * scoreScale; if (score < 0.0) score = 0.0; if (score > 1.0) score = 1.0; return score; } void outputSummary(FILE *f, struct mafSummary *ms) /* output to .tab file */ { fprintf(f, "%u\t", hFindBin(ms->chromStart, ms->chromEnd)); mafSummaryTabOut(ms, f); summaryCount++; } double mergeScores(struct mafSummary *ms1, struct mafSummary *ms2) /* determine score for merged maf summary blocks. * Compute weighted average of block scores vs. bases * ms1 is first summary block and ms2 is second positionally */ { double total = ms1->score * (ms1->chromEnd - ms1->chromStart) + ms2->score * (ms2->chromEnd - ms2->chromStart); return total / (ms2->chromEnd - ms1->chromStart); } struct mafComp *mafMaster(struct mafAli *maf, struct mafFile *mf, char *fileName) /* Get master component from maf. Error abort if no master component */ { struct mafComp *mcMaster = mafMayFindCompPrefix(maf, database, "."); if (mcMaster == NULL) { errAbort("Couldn't find %s. sequence line %d of %s\n", database, mf->lf->lineIx, fileName); } return mcMaster; } char *mafSplitSrcGetChrom(char *src, char* database) -/* src can be in format chrom, db|chrom or db.chrom: split string on separator and return pointer to chrom. - * the db part of src can also have a dot in it, but only if the 'database' argument is not null. - * Changes 'src': The side effect of this function is that src contains only the db, not the chrom anymore. +/* src is one of: chrom, db|chrom, or db.chrom. Return a pointer to the chrom part and + * truncate src in place so it holds only the db. + * + * The hard part is that BOTH db and chrom may themselves contain dots, so the db/chrom + * split is not simply "the first dot" or "the middle dot". We resolve it like this, + * in order: + * 1. A pipe '|' is always an explicit, unambiguous separator -- use it if present. + * 2. No dot at all -> the whole string is the chrom (db is empty). + * 3. If the caller passed the reference 'database' and src begins with "<database>.", + * split right there. This nails the master component whether or not database itself + * contains a dot (e.g. a GenArk db like GCF_000001405.40). + * 4. A GenArk accession db (GCA_/GCF_) is the only db that legitimately carries an + * internal dot -- its accession.version adds exactly one -- so its db ends at the + * SECOND dot and the chrom follows. + * 5. Every other db (hg38, mm10, a species name, ...) has no dot, so split at the FIRST + * dot; anything after it (including further dots, e.g. an accession.version chrom) is + * the chrom. + * + * There is no universal way to resolve an ambiguous multi-dot name from the dots alone: if + * the db part is itself a dotted accession that is not a GenArk GCA_/GCF_ id (e.g. a bare + * INSDC accession used as the assembly), step 5 may split it in the wrong place. For those, + * use the pipe form db|chrom, which is always unambiguous. + * Changes 'src': afterwards src contains only the db, not the chrom. * */ { +/* 1. A pipe is always an explicit db|chrom separator. */ char *pipe = strchr(src, '|'); -// pipe found? It's the new format, db|chrom -if (pipe) { +if (pipe != NULL) + { *pipe = '\0'; return pipe + 1; } char *dot1 = strchr(src, '.'); -if (!dot1) - return src; // if there are no dots, assume the name is the chrom +if (dot1 == NULL) + return src; // no separator: the whole thing is the chrom -if (database) - { - // if 'database' is not NULL we can resolve a situation like GCF_1234.3.CJS12323.4 because we know that - // GCF_1234.3 is the db part - if (differentString(src, database)) +/* 3. When we know the reference database and src starts with "<database>.", split there. */ +if (database != NULL) { - // the database name isn't matching the first part of the component source, - // look to see if maybe the database has a dot in it - *dot1 = '.'; // replace the dot - char *dot2 = strchr(dot1 + 1, '.'); // look for the next dot - if (dot2 != NULL) + int dbLen = strlen(database); + if (startsWith(database, src) && src[dbLen] == '.') { - *dot2 = 0; - char *chrom = dot2 + 1; - return chrom; + src[dbLen] = '\0'; + return src + dbLen + 1; + } } - if ((dot2 == NULL) || differentString(src, database)) - errAbort("expecting first component to have assembly name with no more than one dot"); +/* 4. A GenArk accession db (GCA_/GCF_) carries one internal dot (accession.version), + * so the db ends at the second dot. */ +if (startsWith("GCA_", src) || startsWith("GCF_", src)) + { + char *dot2 = strchr(dot1 + 1, '.'); + if (dot2 != NULL) + { + *dot2 = '\0'; + return dot2 + 1; } + // accession with no chrom after the version: fall through to the first-dot split } -// if database is NULL and there is no pipe character, just split on the first dot and that's it -char* chrom = dot1 + 1; +/* 5. Ordinary db with no dot: split at the first dot. */ *dot1 = '\0'; -return chrom; +return dot1 + 1; } long processMaf(struct mafAli *maf, struct hash *componentHash, FILE *f, struct mafFile *mf, char *fileName) /* Compute scores for each pairwise component in the maf and output to .tab file */ { struct mafComp *mc = NULL, *nextMc = NULL; struct mafSummary *ms, *msPending; struct mafAli pairMaf; long componentCount = 0; struct mafComp *mcMaster = mafMaster(maf, mf, fileName); struct mafComp *oldMasterNext = mcMaster->next; char *chrom; char src[256]; strcpy(src, mcMaster->src); chrom = mafSplitSrcGetChrom(src, database); for (mc = maf->components; mc != NULL; mc = nextMc) { nextMc = mc->next; if (sameString(mcMaster->src, mc->src) || mc->size == 0) continue; /* create maf summary for this alignment component */ AllocVar(ms); ms->chrom = cloneString(chrom); /* both MAF and BED format define chromStart as 0-based */ ms->chromStart = mcMaster->start; /* BED chromEnd is start+size */ ms->chromEnd = mcMaster->start + mcMaster->size; ms->src = cloneString(mc->src); mafSplitSrcGetChrom(ms->src, database); /* construct pairwise maf for scoring */ ZeroVar(&pairMaf); pairMaf.textSize = maf->textSize; pairMaf.components = mcMaster; mcMaster->next = mc; mc->next = NULL; ms->score = scorePairwise(&pairMaf); ms->leftStatus[0] = mc->leftStatus; ms->rightStatus[0] = mc->rightStatus; /* restore component links to allow memory recovery */ mcMaster->next = oldMasterNext; mc->next = nextMc; /* output to .tab file, or save for merging with another block * if this one is too small */ /* handle pending alignment block for this species, if any */ if ((msPending = (struct mafSummary *) hashFindVal(componentHash, ms->src)) != NULL) { /* there is a pending alignment block */ /* either merge it with the current block, or output it */ if (sameString(ms->chrom, msPending->chrom) && (ms->chromStart+1 - msPending->chromEnd < mergeGap)) { /* merge pending block with current */ ms->score = mergeScores(msPending, ms); ms->chromStart = msPending->chromStart; ms->leftStatus[0] = msPending->leftStatus[0]; ms->rightStatus[0] = ms->rightStatus[0]; } else outputSummary(f, msPending); hashRemove(componentHash, msPending->src); mafSummaryFree(&msPending); } /* handle current alignment block (possibly merged) */ if (ms->chromEnd - ms->chromStart > minSize) { /* current block is big enough to output */ outputSummary(f, ms); mafSummaryFree(&ms); } else hashAdd(componentHash, ms->src, ms); componentCount++; } return componentCount; } void flushSummaryBlocks(struct hash *componentHash, FILE *f) /* flush any pending summary blocks */ { struct mafSummary *ms; struct hashCookie hc = hashFirst(componentHash); while ((ms = (struct mafSummary *)hashNextVal(&hc)) != NULL) { outputSummary(f, ms); } } void hgLoadMafSummary(char *db, char *table, char *fileName) /* hgLoadMafSummary - Load a summary table of pairs in a maf into a database. */ { long mafCount = 0, allMafCount = 0; struct mafComp *mcMaster = NULL; struct mafAli *maf; struct mafFile *mf = mafOpen(fileName); struct sqlConnection *conn; FILE *f = hgCreateTabFile(".", table); long componentCount = 0; struct hash *componentHash = newHash(0); if (!test) { conn = sqlConnect(database); mafSummaryTableCreate(conn, table, hGetMinIndexLength(db)); } verbose(1, "Indexing and tabulating %s\n", fileName); /* process mafs */ while ((maf = mafNext(mf)) != NULL) { mcMaster = mafMaster(maf, mf, fileName); allMafCount++; if (mcMaster->srcSize < minSeqSize) continue; while (mcMaster->size > maxSize) { /* break maf into maxSize pieces */ int end = mcMaster->start + maxSize; struct mafAli *subMaf = mafSubset(maf, mcMaster->src, mcMaster->start, end); verbose(3, "Splitting maf %s:%d len %d\n", mcMaster->src, mcMaster->start, mcMaster->size); componentCount += processMaf(subMaf, componentHash, f, mf, fileName); mafAliFree(&subMaf); subMaf = mafSubset(maf, mcMaster->src, end, end + (mcMaster->size - maxSize)); mafAliFree(&maf); maf = subMaf; mcMaster = mafMaster(maf, mf, fileName); } if (mcMaster->size != 0) { /* remainder of maf after splitting off maxSize submafs */ componentCount += processMaf(maf, componentHash, f, mf, fileName); } mafAliFree(&maf); mafCount++; } mafFileFree(&mf); flushSummaryBlocks(componentHash, f); verbose(1, "Created %ld summary blocks from %ld components and %ld mafs from %s\n", summaryCount, componentCount, allMafCount, fileName); if (test) return; verbose(1, "Loading into %s table %s...\n", database, table); hgLoadTabFile(conn, ".", table, &f); verbose(1, "Loading complete"); hgEndUpdate(&conn, 0, 0, "Add %ld maf summary blocks from %s\n", summaryCount, fileName); } int main(int argc, char *argv[]) /* Process command line. */ { optionInit(&argc, argv, optionSpecs); test = optionExists("test"); mergeGap = optionInt("mergeGap", mergeGap); minSize = optionInt("minSize", minSize); maxSize = optionInt("maxSize", maxSize); minSeqSize = optionInt("minSeqSize", minSeqSize); if (argc != 4) usage(); database = argv[1]; hgLoadMafSummary(database, argv[2], argv[3]); return 0; }