4898794edd81be5285ea6e544acbedeaeb31bf78 max Tue Nov 23 08:10:57 2021 -0800 Fixing pointers to README file for license in all source code files. refs #27614 diff --git src/hg/mouseStuff/chainNet/chainNet.c src/hg/mouseStuff/chainNet/chainNet.c index 329149f..8b56906 100644 --- src/hg/mouseStuff/chainNet/chainNet.c +++ src/hg/mouseStuff/chainNet/chainNet.c @@ -1,845 +1,845 @@ /* chainNet - Make alignment nets out of chains. */ /* Copyright (C) 2014 The Regents of the University of California - * See README in this or parent directory for licensing information. */ + * See kent/LICENSE or http://genome.ucsc.edu/license/ for licensing information. */ #include "common.h" #include "linefile.h" #include "localmem.h" #include "hash.h" #include "options.h" #include "dnautil.h" #include "rbTree.h" #include "chainBlock.h" #include "portable.h" int minSpace = 25; /* Minimum gap size to fill. */ int minFill; /* Minimum fill to record. */ double minScore = 2000; /* Minimum chain score to look at. */ boolean inclHap = FALSE; /* include haplotype pseudochromosome queries */ /* command line option specifications */ static struct optionSpec optionSpecs[] = { {"minSpace", OPTION_INT}, {"minFill", OPTION_INT}, {"minScore", OPTION_DOUBLE}, {"inclHap", OPTION_BOOLEAN}, {NULL, 0} }; void usage() /* Explain usage and exit. */ { errAbort( "chainNet - Make alignment nets out of chains\n" "usage:\n" " chainNet in.chain target.sizes query.sizes target.net query.net\n" "where:\n" " in.chain is the chain file sorted by score\n" " target.sizes contains the size of the target sequences\n" " query.sizes contains the size of the query sequences\n" " target.net is the output over the target genome\n" " query.net is the output over the query genome\n" "options:\n" " -minSpace=N - minimum gap size to fill, default %d\n" " -minFill=N - default half of minSpace\n" " -minScore=N - minimum chain score to consider, default %.1lf\n" " -verbose=N - Alter verbosity (default 1)\n" " -inclHap - include query sequences name in the form *_hap*|*_alt*.\n" " Normally these are excluded from nets as being haplotype\n" " pseudochromosomes\n" , minSpace, minScore); } struct gap /* A gap in sequence alignments. */ { struct gap *next; /* Next gap in list. */ int start,end; /* Range covered in chromosome. */ int oStart, oEnd; /* Range covered in other chromosome. */ struct fill *fillList; /* List of gap fillers. */ }; struct fill /* Some alignments that hopefully fill some gaps. */ { struct fill *next; /* Next fill in list. */ int start,end; /* Range covered, always plus strand. */ int oStart,oEnd; /* Range covered in other coordinate, always plus. */ struct gap *gapList; /* List of internal gaps. */ struct chain *chain; /* Alignment chain (not all is necessarily used) */ }; struct space /* A part of a gap. */ { int start, end; /* Portion of gap covered, always plus strand. */ /* The start of this structure has to be identical * with struct range above. */ struct gap *gap; /* The actual gap. */ }; struct chrom /* A chromosome. */ { struct chrom *next; /* Next in list. */ char *name; /* Name - allocated in hash */ int size; /* Size of chromosome. */ struct gap *root; /* Root of the gap/chain tree */ struct rbTree *spaces; /* Store gaps here for fast lookup. Items are spaces. */ }; static boolean inclQuery(struct chain *chain) /* should this query be included? */ { return inclHap || ! haplotype(chain->qName); } int gapCmpStart(const void *va, const void *vb) /* Compare to sort based on start. */ { const struct gap *a = *((struct gap **)va); const struct gap *b = *((struct gap **)vb); return a->start - b->start; } int fillCmpStart(const void *va, const void *vb) /* Compare to sort based on start. */ { const struct fill *a = *((struct fill **)va); const struct fill *b = *((struct fill **)vb); return a->start - b->start; } int spaceCmp(void *va, void *vb) /* Return -1 if a before b, 0 if a and b overlap, * and 1 if a after b. */ { struct space *a = va; struct space *b = vb; if (a->end <= b->start) return -1; else if (b->end <= a->start) return 1; else return 0; } void dumpSpace(void *item, FILE *f) /* Print out range info. */ { struct space *space = item; fprintf(f, "%d,%d", space->start, space->end); } void addSpaceForGap(struct chrom *chrom, struct gap *gap) /* Given a gap, create corresponding space in chromosome's * space rbTree. */ { struct space *space; AllocVar(space); space->gap = gap; space->start = gap->start; space->end = gap->end; rbTreeAdd(chrom->spaces, space); } struct gap *gapNew(int start, int end, int oStart, int oEnd) /* Create a new gap. */ { struct gap *gap; AllocVar(gap); gap->start = start; gap->end = end; gap->oStart = oStart; gap->oEnd = oEnd; return gap; } boolean strictlyInside(int minStart, int maxEnd, int start, int end) /* Return TRUE if minStart < start <= end < maxEnd * and end - start >= minSpace */ { return (minStart < start && start + minSpace <= end && end < maxEnd); } void makeChroms(char *fileName, struct lm *lm, struct rbTreeNode **stack, struct hash **retHash, struct chrom **retList) /* Read size file and make chromosome structure for each element. */ { char *row[2]; struct lineFile *lf = lineFileOpen(fileName, TRUE); struct hash *hash = newHash(0); struct chrom *chrom, *chromList = NULL; while (lineFileRow(lf, row)) { char *name = row[0]; if (hashLookup(hash, name) != NULL) errAbort("Duplicate %s in %s", name, fileName); AllocVar(chrom); slAddHead(&chromList, chrom); hashAddSaveName(hash, name, chrom, &chrom->name); chrom->size = lineFileNeedNum(lf, row, 1); chrom->spaces = rbTreeNewDetailed(spaceCmp, lm, stack); chrom->root = gapNew(0, chrom->size, 0, 0); addSpaceForGap(chrom, chrom->root); } lineFileClose(&lf); slReverse(&chromList); *retHash = hash; *retList = chromList; } boolean innerBounds(struct cBlock *startBlock, boolean isQ, int inStart, int inEnd, int *retStart, int *retEnd) /* Return the portion of chain which is between inStart and * inEnd. Only considers blocks inbetween inStart and inEnd. * Return NULL if chain has no blocks between these. */ { struct cBlock *b; int start = BIGNUM, end = -BIGNUM; for (b = startBlock; b != NULL; b = b->next) { int s, e; if (isQ) { s = b->qStart; e = b->qEnd; } else { s = b->tStart; e = b->tEnd; } if (e <= inStart) continue; if (s >= inEnd) break; if (s < inStart) s = inStart; if (e > inEnd) e = inEnd; if (start > s) start = s; if (end < e) end = e; } if ((end < 0) || ((end - start) < minFill)) return FALSE; *retStart = start; *retEnd = end; return TRUE; } static void qFillOtherRange(struct fill *fill) /* Given bounds of fill in q coordinates, calculate * oStart/oEnd in t coordinates, and refine * start/end to reflect parts of chain actually used. */ { struct chain *chain = fill->chain; int clipStart = fill->start; int clipEnd = fill->end; boolean isRev = (chain->qStrand == '-'); int tMin = BIGNUM, tMax = -BIGNUM; int qMin = BIGNUM, qMax = -BIGNUM; struct cBlock *b; if (isRev) reverseIntRange(&clipStart, &clipEnd, chain->qSize); for (b = chain->blockList; b != NULL; b = b->next) { int qs, qe, ts, te; /* Clipped bounds of block */ if ((qe = b->qEnd) <= clipStart) continue; if ((qs = b->qStart) >= clipEnd) break; ts = b->tStart; te = b->tEnd; if (qs < clipStart) { ts += (clipStart - qs); qs = clipStart; } if (qe > clipEnd) { te -= (qe - clipEnd); qe = clipEnd; } if (qMin > qs) qMin = qs; if (qMax < qe) qMax = qe; if (tMin > ts) tMin = ts; if (tMax < te) tMax = te; } if (isRev) reverseIntRange(&qMin, &qMax, chain->qSize); fill->start = qMin; fill->end = qMax; fill->oStart = tMin; fill->oEnd = tMax; assert(tMin < tMax); } static void tFillOtherRange(struct fill *fill) /* Given bounds of fill in t coordinates, calculate * oStart/oEnd in q coordinates, and refine * start/end to reflect parts of chain actually used. */ { struct chain *chain = fill->chain; int clipStart = fill->start; int clipEnd = fill->end; boolean isRev = (chain->qStrand == '-'); int tMin = BIGNUM, tMax = -BIGNUM; int qMin = BIGNUM, qMax = -BIGNUM; struct cBlock *b; for (b = chain->blockList; b != NULL; b = b->next) { int qs, qe, ts, te; /* Clipped bounds of block */ if ((te = b->tEnd) <= clipStart) continue; if ((ts = b->tStart) >= clipEnd) break; qs = b->qStart; qe = b->qEnd; if (ts < clipStart) { qs += (clipStart - ts); ts = clipStart; } if (te > clipEnd) { qe -= (te - clipEnd); te = clipEnd; } if (qMin > qs) qMin = qs; if (qMax < qe) qMax = qe; if (tMin > ts) tMin = ts; if (tMax < te) tMax = te; } if (isRev) reverseIntRange(&qMin, &qMax, chain->qSize); fill->start = tMin; fill->end = tMax; fill->oStart = qMin; fill->oEnd = qMax; assert(qMin < qMax); } struct fill *fillSpace(struct chrom *chrom, struct space *space, struct chain *chain, struct cBlock *startBlock, boolean isQ) /* Fill in space with chain, remove existing space from chrom, * and add smaller spaces on either side if big enough. */ { struct fill *fill; int s, e; struct space *lSpace, *rSpace; if (!innerBounds(startBlock, isQ, space->start, space->end, &s, &e)) return NULL; assert(s < e); AllocVar(fill); fill->start = s; fill->end = e; fill->chain = chain; rbTreeRemove(chrom->spaces, space); if (s - space->start >= minSpace) { AllocVar(lSpace); lSpace->gap = space->gap; lSpace->start = space->start; lSpace->end = s; rbTreeAdd(chrom->spaces, lSpace); } if (space->end - e >= minSpace) { AllocVar(rSpace); rSpace->gap = space->gap; rSpace->start = e; rSpace->end = space->end; rbTreeAdd(chrom->spaces, rSpace); } slAddHead(&space->gap->fillList, fill); return fill; } static struct slRef *fsList; void fsAdd(void *item) /* Add item to fsList. */ { refAdd(&fsList, item); } struct slRef *findSpaces(struct rbTree *tree, int start, int end) /* Return a list of spaces that intersect interval between start * and end. */ { static struct space space; space.start = start; space.end = end; fsList = NULL; rbTreeTraverseRange(tree, &space, &space, fsAdd); slReverse(&fsList); return fsList; } void reverseBlocksQ(struct cBlock **pList, int qSize) /* Reverse qside of blocks. */ { struct cBlock *b; slReverse(pList); for (b = *pList; b != NULL; b = b->next) reverseIntRange(&b->qStart, &b->qEnd, qSize); } void addChainT(struct chrom *chrom, struct chrom *otherChrom, struct chain *chain) /* Add T side of chain to fill/gap tree of chromosome. * This is the easier case since there are no strand * issues to worry about. */ { struct slRef *spaceList; struct slRef *ref; struct cBlock *startBlock, *block, *nextBlock; struct gap *gap; spaceList = findSpaces(chrom->spaces,chain->tStart,chain->tEnd); startBlock = chain->blockList; for (ref = spaceList; ref != NULL; ref = ref->next) { struct space *space = ref->val; struct fill *fill; int gapStart, gapEnd; for (;;) { nextBlock = startBlock->next; if (nextBlock == NULL) break; gapEnd = nextBlock->tStart; if (gapEnd > space->start) break; startBlock = nextBlock; } if ((fill = fillSpace(chrom, space, chain, startBlock, FALSE)) != NULL) { for (block = startBlock; ; block = nextBlock) { nextBlock = block->next; if (nextBlock == NULL) break; gapStart = block->tEnd; gapEnd = nextBlock->tStart; if (strictlyInside(space->start, space->end, gapStart, gapEnd)) { int qs = block->qEnd; int qe = nextBlock->qStart; if (chain->qStrand == '-') reverseIntRange(&qs, &qe, chain->qSize); gap = gapNew(gapStart, gapEnd, qs, qe); addSpaceForGap(chrom, gap); slAddHead(&fill->gapList, gap); } } freez(&ref->val); /* aka space */ } } slFreeList(&spaceList); } void addChainQ(struct chrom *chrom, struct chrom *otherChrom, struct chain *chain) /* Add Q side of chain to fill/gap tree of chromosome. * For this side we have to cope with reverse strand * issues. */ { struct slRef *spaceList; struct slRef *ref; struct cBlock *startBlock, *block, *nextBlock; int gapStart, gapEnd; struct gap *gap; boolean isRev = (chain->qStrand == '-'); int qStart = chain->qStart, qEnd = chain->qEnd; if (isRev) { reverseIntRange(&qStart, &qEnd, chain->qSize); reverseBlocksQ(&chain->blockList, chain->qSize); } spaceList = findSpaces(chrom->spaces,qStart,qEnd); startBlock = chain->blockList; for (ref = spaceList; ref != NULL; ref = ref->next) { struct space *space = ref->val; struct fill *fill; for (;;) { nextBlock = startBlock->next; if (nextBlock == NULL) break; gapEnd = nextBlock->qStart; if (gapEnd > space->start) break; startBlock = nextBlock; } if ((fill = fillSpace(chrom, space, chain, startBlock, TRUE)) != NULL) { for (block = startBlock; ; block = nextBlock) { nextBlock = block->next; if (nextBlock == NULL) break; gapStart = block->qEnd; gapEnd = nextBlock->qStart; if (strictlyInside(space->start, space->end, gapStart, gapEnd)) { int ts, te; if (chain->qStrand == '+') { ts = block->tEnd; te = nextBlock->tStart; } else { ts = nextBlock->tStart; te = block->tEnd; } gap = gapNew(gapStart, gapEnd, ts, te); addSpaceForGap(chrom, gap); slAddHead(&fill->gapList, gap); } } freez(&ref->val); /* aka space */ } } slFreeList(&spaceList); if (isRev) reverseBlocksQ(&chain->blockList, chain->qSize); } void addChain(struct chrom *qChrom, struct chrom *tChrom, struct chain *chain) /* Add as much of chain as possible to chromosomes. */ { addChainQ(qChrom, tChrom, chain); addChainT(tChrom, qChrom, chain); } boolean chromHasData(struct chrom *chrom) /* Return TRUE if chrom is non-empty */ { return chrom != NULL && chrom->root != NULL && chrom->root->fillList != NULL; } void sortNet(struct gap *gap) /* Recursively sort lists. */ { struct fill *fill; slSort(&gap->fillList, fillCmpStart); for (fill = gap->fillList; fill != NULL; fill = fill->next) { slSort(&fill->gapList, gapCmpStart); for (gap = fill->gapList; gap != NULL; gap = gap->next) sortNet(gap); } } void rCalcOtherFill(struct gap *gap, boolean isQ) /* Recursively fill in oStart/oEnd fields of fill * structures. */ { struct fill *fill; for (fill = gap->fillList; fill != NULL; fill = fill->next) { struct gap *g; if (isQ) qFillOtherRange(fill); else tFillOtherRange(fill); for (g = fill->gapList; g != NULL; g = g->next) rCalcOtherFill(g, isQ); } } void finishNet(struct chrom *chromList, boolean isQ) /* Fill in oStart/oEnd fields of fill structures attatched * to chromosomes, sort, and generally spiff up net. */ { struct chrom *chrom; for (chrom = chromList; chrom != NULL; chrom = chrom->next) { if (chromHasData(chrom)) { sortNet(chrom->root); rCalcOtherFill(chrom->root, isQ); } } } int rOutDepth = 0; boolean rOutQ = FALSE; static void rOutputGap(struct fill *parent, struct gap *gap, FILE *f); static void rOutputFill(struct fill *fill, FILE *f); static void rOutputGap(struct fill *parent, struct gap *gap, FILE *f) /* Recursively output gap and it's fillers. */ { struct fill *fill; struct chain *chain = parent->chain; char *oChrom = (rOutQ ? chain->tName : chain->qName); ++rOutDepth; spaceOut(f, rOutDepth); fprintf(f, "gap %d %d %s %c %d %d\n", gap->start, gap->end - gap->start, oChrom, chain->qStrand, gap->oStart, gap->oEnd - gap->oStart); for (fill = gap->fillList; fill != NULL; fill = fill->next) rOutputFill(fill, f); --rOutDepth; } int chainBaseCount(struct chain *chain) /* Return number of bases in gap free alignments in chain. */ { struct cBlock *b; int total = 0; for (b = chain->blockList; b != NULL; b = b->next) total += b->qEnd - b->qStart; return total; } int chainBaseCountSubT(struct chain *chain, int tMin, int tMax) /* Return number of bases in gap free alignments in chain between * tMin and tMax. */ { struct cBlock *b; int total = 0; for (b = chain->blockList; b != NULL; b = b->next) total += positiveRangeIntersection(b->tStart, b->tEnd, tMin, tMax); return total; } int chainBaseCountSubQ(struct chain *chain, int qMin, int qMax) /* Return number of bases in gap free alignments in chain between * tMin and tMax. */ { struct cBlock *b; int total = 0; for (b = chain->blockList; b != NULL; b = b->next) total += positiveRangeIntersection(b->qStart, b->qEnd, qMin, qMax); return total; } void subchainInfo(struct chain *chain, int start, int end, boolean isQ, int *retSize, double *retScore) /* Return score of subchain. */ { int fullSize = chainBaseCount(chain); if (isQ) { if (chain->qStrand == '-') reverseIntRange(&start, &end, chain->qSize); if (start <= chain->qStart && end >= chain->qEnd) { *retScore = chain->score; *retSize = fullSize; } else { int subSize = chainBaseCountSubQ(chain, start, end); *retScore = chain->score * subSize / fullSize; *retSize = subSize; // uglyf("subchainInfo Q%c %d,%d %d,%d ali %d, subSize %d, score %f, subScore %f\n", chain->qStrand, start, end, chain->qStart, chain->qEnd, fullSize, subSize, chain->score, *retScore); } } else { if (start <= chain->tStart && end >= chain->tEnd) { *retScore = chain->score; *retSize = fullSize; } else { int subSize = chainBaseCountSubT(chain, start, end); *retScore = chain->score * subSize / fullSize; *retSize = subSize; } } } void fillOut(FILE *f, struct fill *fill, char *oChrom, int depth, int subSize, double subScore) /* Output fill. */ { struct chain *chain = fill->chain; spaceOut(f, depth); fprintf(f, "fill %d %d %s %c %d %d id %d score %1.0f ali %d\n", fill->start, fill->end - fill->start, oChrom, chain->qStrand, fill->oStart, fill->oEnd - fill->oStart, chain->id, subScore, subSize); } static void rOutputFill(struct fill *fill, FILE *f) /* Recursively output fill and it's gaps. */ { struct gap *gap; struct chain *chain = fill->chain; int subSize; double subScore; subchainInfo(chain, fill->start, fill->end, rOutQ, &subSize, &subScore); if (subScore >= minScore && subSize >= minFill) { ++rOutDepth; if (rOutQ) fillOut(f, fill, fill->chain->tName, rOutDepth, subSize, subScore); else fillOut(f, fill, fill->chain->qName, rOutDepth, subSize, subScore); for (gap = fill->gapList; gap != NULL; gap = gap->next) rOutputGap(fill, gap, f); --rOutDepth; } } void outputNetSide(struct chrom *chromList, FILE *f, boolean isQ) /* Output one side of net */ { struct chrom *chrom; for (chrom = chromList; chrom != NULL; chrom = chrom->next) { rOutDepth = 0; rOutQ = isQ; if (chromHasData(chrom)) { struct fill *fill; fprintf(f, "net %s %d\n", chrom->name, chrom->size); for (fill = chrom->root->fillList; fill != NULL; fill = fill->next) rOutputFill(fill, f); } } } void printMem(FILE *f) /* Print out memory used and other stuff from linux. */ { struct lineFile *lf = lineFileMayOpen("/proc/self/stat", TRUE); if (lf != NULL) { char *line, *words[50]; int wordCount; if (lineFileNext(lf, &line, NULL)) { wordCount = chopLine(line, words); if (wordCount >= 23) fprintf(f, "memory usage %s, utime %s s/100, stime %s\n", words[22], words[13], words[14]); } lineFileClose(&lf); } } void chainNet(char *chainFile, char *tSizes, char *qSizes, char *tNet, char *qNet) /* chainNet - Make alignment nets out of chains. */ { struct lineFile *lf = lineFileOpen(chainFile, TRUE); struct hash *qHash, *tHash; struct chrom *qChromList, *tChromList, *tChrom, *qChrom; struct chain *chain; double lastScore = -1; struct lm *lm = lmInit(0); struct rbTreeNode **rbStack; FILE *tNetFile = mustOpen(tNet, "w"); FILE *qNetFile = mustOpen(qNet, "w"); lmAllocArray(lm, rbStack, 256); makeChroms(qSizes, lm, rbStack, &qHash, &qChromList); makeChroms(tSizes, lm, rbStack, &tHash, &tChromList); verbose(1, "Got %d chroms in %s, %d in %s\n", slCount(tChromList), tSizes, slCount(qChromList), qSizes); lineFileSetMetaDataOutput(lf, tNetFile); lineFileSetMetaDataOutput(lf, qNetFile); /* Loop through chain file building up net. */ while ((chain = chainRead(lf)) != NULL) { /* Make sure that input is really sorted. */ if (lastScore >= 0 && chain->score > lastScore) errAbort("%s must be sorted in order of score", chainFile); lastScore = chain->score; if (chain->score < minScore) { break; } verbose(2, "chain %f (%d els) %s %d-%d %c %s %d-%d\n", chain->score, slCount(chain->blockList), chain->tName, chain->tStart, chain->tEnd, chain->qStrand, chain->qName, chain->qStart, chain->qEnd); qChrom = hashMustFindVal(qHash, chain->qName); if (qChrom->size != chain->qSize) errAbort("%s is %d in %s but %d in %s", chain->qName, chain->qSize, chainFile, qChrom->size, qSizes); tChrom = hashMustFindVal(tHash, chain->tName); if (tChrom->size != chain->tSize) errAbort("%s is %d in %s but %d in %s", chain->tName, chain->tSize, chainFile, tChrom->size, tSizes); if (!inclQuery(chain)) verbose(2, "skipping chain on query %s\n", chain->qName); else { addChain(qChrom, tChrom, chain); verbose(2, "%s has %d inserts, %s has %d\n", tChrom->name, tChrom->spaces->n, qChrom->name, qChrom->spaces->n); } } /* Build up other side of fills. It's just for historical * reasons this is not done during the main build up. * It's a little less efficient this way, but to change it * some hard reverse strand issues would have to be juggled. */ verbose(1, "Finishing nets\n"); finishNet(qChromList, TRUE); finishNet(tChromList, FALSE); /* Write out basic net files. */ verbose(1, "writing %s\n", tNet); outputNetSide(tChromList, tNetFile, FALSE); verbose(1, "writing %s\n", qNet); outputNetSide(qChromList, qNetFile, TRUE); /* prevent SIGPIPE in preceding process if input is a pipe, consume remainder * of input file since we stop before EOF. */ if (isPipe(lf->fd)) { char *line; while(lineFileNext(lf, &line, NULL)) continue; } lineFileClose(&lf); if (verboseLevel() > 1) printMem(stderr); } int main(int argc, char *argv[]) /* Process command line. */ { optionInit(&argc, argv, optionSpecs); if (argc != 6) usage(); minSpace = optionInt("minSpace", minSpace); minFill = optionInt("minFill", minSpace/2); minScore = optionInt("minScore", minScore); inclHap = optionExists("inclHap"); chainNet(argv[1], argv[2], argv[3], argv[4], argv[5]); return 0; }