0750648bee41e57d7f1cbc3de8e3d9177bfd4c70 angie Fri May 1 10:02:25 2026 -0700 Initial support for ripples seach in usher-sampled-server. Enabled by config param ripplesEnabled. If enabled, there is a new checkbox for the user to click if they want the extra search. Currently hardcoded to max of 10 sequences, otherwise extra search is not done (too slow). Ripples results, if any, are displayed in a table before the usual summary table. diff --git src/hg/hgPhyloPlace/runUsher.c src/hg/hgPhyloPlace/runUsher.c index 62fe16ace9e..febad28b378 100644 --- src/hg/hgPhyloPlace/runUsher.c +++ src/hg/hgPhyloPlace/runUsher.c @@ -538,52 +538,337 @@ *ptr = '\0'; info->nextClade = hashStoreName(wordStore, words[1]); } if (wordCount > 2) { // Chop extra chars in usher-sampled clades.txt output char *ptr = strstr(words[2], "*|"); if (ptr) *ptr = '\0'; info->pangoLineage = hashStoreName(wordStore, words[2]); } } lineFileClose(&lf); } +static char *finalRecombHeaderExpected = "recomb_node_id\tdonor_node_id\tacceptor_node_id\trecombinant_num_desc\tdonor_num_desc\tacceptor_num_desc\t" + "breakpoint interval 1\tbreakpoint interval 2\trecombinant clade\trecombinant lineage\tdonor clade\tdonor lineage\t" + "acceptor clade\tacceptor lineage\trepresentative descendant\toriginal parsimony score\tparsimony score improvement"; + +static uint parseUint(char *word, struct lineFile *lf) +/* Parse non-negative int from word. Use lf for error reporting if necessary. */ +{ +if (!isAllDigits(word)) + lineFileAbort(lf, "Expected a non-negative number but got '%s'", word); +int val = atol(word); +if (val < 0) + lineFileAbort(lf, "Expected a non-negative number but got '%s'", word); +return val; +} + +#define RIPPLES_GENOME_SIZE "GENOME_SIZE" + +static void parseBpRange(char *word, uint genomeSize, uint *retMin, uint *retMax, struct lineFile *lf) +/* Parse out min and max from a recombinant breakpoint range like "(253,1432)" or "(29409,GENOME_SIZE)". + * Use lf for error reporting if necessary. */ +{ +int min = 0, max = 0; +char *p = word; +if (*p++ != '(') + lineFileAbort(lf, "Expected breakpoint range beginning with '(' but got '%s'", word); +if (startsWith(RIPPLES_GENOME_SIZE",", p)) + { + min = genomeSize; + p += strlen(RIPPLES_GENOME_SIZE); + } +else + { + char *start = p; + while (*p != '\0' && isdigit(*p)) + p++; + min = atol(start); + } +if (*p++ != ',') + lineFileAbort(lf, "Expected breakpoint range like '(1,2)' but didn't find comma, got '%s'", word); +if (sameString(p, RIPPLES_GENOME_SIZE")")) + { + max = genomeSize; + p += strlen(RIPPLES_GENOME_SIZE); + } +else + { + char *start = p; + while (*p != '\0' && isdigit(*p)) + p++; + max = atol(start); + } +if (differentString(p, ")")) + lineFileAbort(lf, "Expected breakpoint range like '(1,2)' but didn't find final ')', got '%s'", word); +if (max < min) + lineFileAbort(lf, "Expected breakpoint range like '(1,2)' but got max < min: '%s'", word); +if (min < 0 || max < 0) + lineFileAbort(lf, "Expected breakpoint range like '(1,2)' but got negative number: '%s'", word); +*retMin = min; +*retMax = max; +} + +static struct recombinantInfo *parseOneRecombinant(char *line, uint genomeSize, struct lineFile *lf) +/* Allocate and return one struct recombinantInfo with values extracted from line. Use lf for error reporting if necessary. */ +{ +struct recombinantInfo *ri = NULL; +char *words[18]; +int wordCount = chopTabs(line, words); +if (wordCount != 17) + lineFileAbort(lf, "Expected 17 tab-separated words but got %d", wordCount); +AllocVar(ri); +ri->recombNodeId = cloneString(words[0]); +ri->donorNodeId = cloneString(words[1]); +ri->acceptorNodeId = cloneString(words[2]); +ri->recombNumDesc = parseUint(words[3], lf); +ri->donorNumDesc = parseUint(words[4], lf); +ri->acceptorNumDesc = parseUint(words[5], lf); +parseBpRange(words[6], genomeSize, &(ri->bp1Min), &(ri->bp1Max), lf); +parseBpRange(words[7], genomeSize, &(ri->bp2Min), &(ri->bp2Max), lf); +ri->recombClade = cloneString(words[8]); +ri->recombLineage = cloneString(words[9]); +ri->donorClade = cloneString(words[10]); +ri->donorLineage = cloneString(words[11]); +ri->acceptorClade = cloneString(words[12]); +ri->acceptorLineage = cloneString(words[13]); +ri->representative = cloneString(words[14]); +ri->originalParsimony = parseUint(words[15], lf); +ri->parsimonyImprovement = parseUint(words[16], lf); +return ri; +} + +static void recombinantInfoFree(struct recombinantInfo **pRi) +/* Free a struct recombinantInfo and strings that it points to. */ +{ +if (pRi && *pRi) + { + struct recombinantInfo *ri = *pRi; + freeMem(ri->recombNodeId); + freeMem(ri->donorNodeId); + freeMem(ri->acceptorNodeId); + freeMem(ri->recombClade); + freeMem(ri->recombLineage); + freeMem(ri->donorClade); + freeMem(ri->donorLineage); + freeMem(ri->acceptorClade); + freeMem(ri->acceptorLineage); + freeMem(ri->representative); + freez(pRi); + } +} + +static int recombinantInfoCmp(const void *el1, const void *el2) +/* For sorting by parsimonyImprovement, descending. */ +{ +struct recombinantInfo *ri1 = *(struct recombinantInfo **)el1; +struct recombinantInfo *ri2 = *(struct recombinantInfo **)el2; +return ri2->parsimonyImprovement - ri1->parsimonyImprovement; +} + +static boolean mergeRecombinants(struct recombinantInfo **pRiOldList, struct recombinantInfo **pRiNew) +/* If riNew has a greater parsimony improvement than riOldList then replace riOldList with riNew. + * If riNew has a smaller parsimony improvement than riOldList then free up riNew. + * If riNew has the same parsimony improvement then attempt to merge its breakpoint ranges with + * each member of riOldList. + * Return TRUE if any of those were successful, FALSE only if there were breakpoint ranges that + * could not be merged. */ +{ +boolean success = FALSE; +struct recombinantInfo *riNew = *pRiNew, *riOld = *pRiOldList; +if (riNew->parsimonyImprovement > riOld->parsimonyImprovement) + { + slFreeListWithFunc(pRiOldList, recombinantInfoFree); + *pRiOldList = *pRiNew; + *pRiNew = NULL; + success = TRUE; + } +else if (riNew->parsimonyImprovement < riOld->parsimonyImprovement) + { + recombinantInfoFree(pRiNew); + success = TRUE; + } +else + { + // Attempt to merge breakpoint ranges for each member of riOldList. + for (riOld = *pRiOldList; riOld != NULL; riOld = riOld->next) + { + if (riNew->bp1Min == riOld->bp1Min && riNew->bp1Max == riOld->bp1Max) + { + // Only bp2 differs; merge riNew bp2 into riOld bp2 + riOld->bp2Min = min(riOld->bp2Min, riNew->bp2Min); + riOld->bp2Max = max(riOld->bp2Max, riNew->bp2Max); + recombinantInfoFree(pRiNew); + success = TRUE; + } + else if (riNew->bp2Min == riOld->bp2Min && riNew->bp2Max == riOld->bp2Max) + { + // Only bp1 differs; merge r1New bp1 into riOld bp1 + riOld->bp1Min = min(riOld->bp1Min, riNew->bp1Min); + riOld->bp2Max = max(riOld->bp1Max, riNew->bp1Max); + recombinantInfoFree(pRiNew); + success = TRUE; + } + if (success) + break; + } + } +return success; +} + +static struct recombinantInfo *filterRecombinants(struct recombinantInfo *riList) +/* Filter riList, which must be sorted by parsimonyImprovement, to keep at most the top 3 results + * for each potential recombinant node (plus any subsequent results that are tied for 3rd place + * with the same parsimony improvement). */ +{ +struct recombinantInfo *riListNew = NULL; +struct hash *resultCounts = hashNew(0); +struct hash *minScores = hashNew(0); +struct recombinantInfo *ri, *riNext; +for (ri = riList; ri != NULL; ri = riNext) + { + boolean keepThis = TRUE; + riNext = ri->next; + int count = hashIncInt(resultCounts, ri->recombNodeId); + if (count == 3) + hashAddInt(minScores, ri->recombNodeId, ri->parsimonyImprovement); + else if (count > 3) + { + int minScore = hashIntVal(minScores, ri->recombNodeId); + if (ri->parsimonyImprovement < minScore) + { + // Not tied for third place -- discard. + keepThis = FALSE; + } + } + ri->next = NULL; + if (keepThis) + slAddHead(&riListNew, ri); + else + recombinantInfoFree(&ri); + } +slReverse(&riListNew); +hashFree(&resultCounts); +hashFree(&minScores); +return riListNew; +} + +static char *getRecombinantKey(struct recombinantInfo *ri) +/* Make a hash key by concatenating node IDs. Not thread safe, do not free return value. */ +{ +static struct dyString *dy = NULL; +if (dy == NULL) + dy = dyStringNew(0); +else + dyStringClear(dy); +dyStringPrintf(dy, "%s %s %s", ri->recombNodeId, ri->donorNodeId, ri->acceptorNodeId); +return dy->string; +} + +static struct recombinantInfo *parseRecombinants(char *filename, uint genomeSize) +/* Parse final_recombination.tsv from ripples search and merge rows where possible. */ +{ +struct recombinantInfo *recombList = NULL; +struct lineFile *lf = lineFileOpen(filename, TRUE); +char *line; +lineFileNext(lf, &line, NULL); +if (isNotEmpty(line) && differentString(line, finalRecombHeaderExpected)) + lineFileAbort(lf, "Header fields do not match expected. Header:\n%s\nExpected:\n%s", + line, finalRecombHeaderExpected); +// Hash reported recombinants by concatenated node IDs because there tend to be many lines per +// triplet that have all the same info except for parsimony score improvement (where we want to +// keep only the highest improvement found per triplet) and breakpoints (where we want to merge +// overlapping breakpoint regions). +struct hash *mergedRecombinants = hashNew(0); +while (lineFileNext(lf, &line, NULL)) + { + struct recombinantInfo *riNew = parseOneRecombinant(line, genomeSize, lf); + char *key = getRecombinantKey(riNew); + struct hashEl *hel = hashLookup(mergedRecombinants, key); + if (hel) + { + if (! mergeRecombinants((struct recombinantInfo **)&(hel->val), &riNew)) + slAddHead(&(hel->val), riNew); + } + else + { + hashAdd(mergedRecombinants, key, riNew); + } + } +lineFileClose(&lf); +// Add all mergedRecombinants to list, and sort by parsimony improvement (highest first) +struct hashEl *hel, *helList = hashElListHash(mergedRecombinants); +for (hel = helList; hel != NULL; hel = hel->next) + { + struct recombinantInfo *riList = hel->val; + recombList = slCat(recombList, riList); + } +hashElFreeList(&helList); +slSort(&recombList, recombinantInfoCmp); +// Filter the list to keep only the top 3 findings for each potential recombinant +recombList = filterRecombinants(recombList); +return recombList; +} + +static char *descendantsHeaderExpected = "#node_id\tdescendants"; + +static struct hash *parseDescendants(char *filename) +/* Parse descendants.tsv from ripples search into hash of node ID to slName list of leaves. */ +{ +struct hash *hash = hashNew(0); +struct lineFile *lf = lineFileOpen(filename, TRUE); +char *line; +lineFileNext(lf, &line, NULL); +if (isNotEmpty(line) && differentString(line, descendantsHeaderExpected)) + lineFileAbort(lf, "Header fields do not match expected. Header:\n%s\nExpected:\n%s", + line, descendantsHeaderExpected); +char *words[3]; +int wordCount; +while ((wordCount = lineFileChopTab(lf, words)) > 0) + { + lineFileExpectWords(lf, 2, wordCount); + hashAdd(hash, words[0], slNameListFromComma(words[1])); + } +lineFileClose(&lf); +return hash; +} + static char *dirPlusFile(struct dyString *dy, char *dir, char *file) /* Write dir/file into dy and return pointer to dy->string. Do not free result! */ { dyStringClear(dy); dyStringPrintf(dy, "%s/%s", dir, file); return dy->string; } static boolean isDevHost() /* Return TRUE if it looks like we're running on hgwdev (or command line). */ { char *httpHost = getenv("HTTP_HOST"); return (httpHost == NULL || startsWith("hgwdev", httpHost) || startsWith("genome-test", httpHost)); } #define USHER_SAMPLED_TREE_PREFIX "tree-0-" static int processOutDirFiles(struct usherResults *results, char *outDir, struct tempName **retSingleSubtreeTn, struct variantPathNode **retSingleSubtreeMuts, struct tempName *subtreeTns[], struct variantPathNode *subtreeMuts[], - int maxSubtrees) + int maxSubtrees, uint genomeSize) /* Get paths to files in outDir; parse them and move files that we'll keep up to trash/ct/, * leaving behind files that we will remove when done. */ { int subtreeCount = 0; memset(subtreeTns, 0, maxSubtrees * sizeof(*subtreeTns)); memset(subtreeMuts, 0, maxSubtrees * sizeof(*subtreeMuts)); struct dyString *dyScratch = dyStringNew(0); struct slName *outDirFiles = listDir(outDir, "*"), *file; boolean isDev = isDevHost(); for (file = outDirFiles; file != NULL; file = file->next) { char *path = dirPlusFile(dyScratch, outDir, file->name); if (sameString(file->name, "uncondensed-final-tree.nh") || sameString(file->name, "uncondensed-final-tree1.nh")) { mustRename(path, results->bigTreePlusTn->forCgi); @@ -687,33 +972,42 @@ } else if (sameString(parts[2], "expanded.txt")) { // Don't need this, just remove it } else if (isDev) warn("Unexpected filename '%s' from usher, ignoring", file->name); } else if (isDev) warn("Unexpected filename '%s' from usher, ignoring", file->name); } else if (sameString(file->name, "clades.txt")) { parseClades(path, results->samplePlacements); } + else if (sameString(file->name, "final_recombination.tsv")) + { + results->recombinants = parseRecombinants(path, genomeSize); + } + else if (sameString(file->name, "descendants.tsv")) + { + results->recombinantDescendants = parseDescendants(path); + } else if (sameString(file->name, "final-tree.nh") || sameString(file->name, "current-tree.nh") || - sameString(file->name, "placement_stats.tsv")) + sameString(file->name, "placement_stats.tsv") || + sameString(file->name, "recombination.tsv")) { // Don't need this (or already parsed it elsewhere not here), just remove it. } else if (isDev) warn("Unexpected filename '%s' from usher, ignoring", file->name); } dyStringFree(&dyScratch); // Make sure we got a complete range of subtrees [0..subtreeCount-1] int i; for (i = 0; i < subtreeCount; i++) { if (subtreeTns[i] == NULL) errAbort("Missing file subtree-%d.nh in usher results", i+1); if (subtreeMuts[i] == NULL) errAbort("Missing file subtree-%d-mutations.txt in usher results", i+1); @@ -1149,48 +1443,40 @@ // Give up - fall back on regular usher command fprintf(errFile, "Second attempt to connect socket %d to path '%s' failed: %s\n", socketFd, usherServerSocketPath, strerror(errno)); socketFd = -1; } } } } return socketFd; } // Server sends ASCII EOT character (4) when done. Sadly I can't find a header file that defines EOT. #define EOT 4 static boolean sendQuery(int socketFd, char *cmd[], char *org, struct treeChoices *treeChoices, - FILE *errFile, boolean addNoIgnorePrefix, char *anchorFile) + FILE *errFile) /* Send command to socket, read response on socket, return TRUE if we get a successful response. */ { boolean success = FALSE; struct dyString *dyMessage = dyStringNew(0); int ix; for (ix = 0; cmd[ix] != NULL; ix++) { - // Don't include args from -T onward; server rejects requests with -T or --optimization_radius - if (sameString("-T", cmd[ix])) - break; dyStringPrintf(dyMessage, "%s\n", cmd[ix]); } -if (addNoIgnorePrefix) - // Needed when placing uploaded sequences, but not when finding uploaded names - dyStringPrintf(dyMessage, "--no-ignore-prefix\n"USHER_DEDUP_PREFIX"\n"); -if (isNotEmpty(anchorFile)) - dyStringPrintf(dyMessage, "--anchor-samples\n%s\n", anchorFile); dyStringAppendC(dyMessage, '\n'); boolean serverError = FALSE; int bytesWritten = write(socketFd, dyMessage->string, dyMessage->stringSize); if (bytesWritten == dyMessage->stringSize) { struct lineFile *lf = lineFileAttach("server socket", TRUE, socketFd); if (lf) { char *line; while (lineFileNext(lf, &line, NULL)) { if (startsWith("Tree", line) && endsWith(line, "not found")) { // Tell the server to reload the latest protobufs serverReloadProtobufs(org, treeChoices); @@ -1206,122 +1492,142 @@ else if (isNotEmpty(line)) fprintf(errFile, "%s\n", line); } } else fprintf(errFile, "Failed to attach linefile to socket %d.\n", socketFd); } else fprintf(errFile, "Failed to send query to socket %d: attempted to write %ld bytes, ret=%d\n", socketFd, dyMessage->stringSize, bytesWritten); dyStringFree(&dyMessage); return success; } static boolean runUsherServer(char *org, char *cmd[], char *stderrFile, - struct treeChoices *treeChoices, char *anchorFile, int *pStartTime) + struct treeChoices *treeChoices, int *pStartTime) /* Start the server if necessary, connect to it, send a query, get response and return TRUE if. * all goes well. If unsuccessful, write reasons to errFile and return FALSE. */ { boolean success = FALSE; if (serverIsConfigured(org)) { FILE *errFile = mustOpen(stderrFile, "w"); int serverSocket = getServerSocket(org, treeChoices, errFile); reportTiming(pStartTime, "get socket"); if (serverSocket > 0) { - success = sendQuery(serverSocket, cmd, org, treeChoices, errFile, TRUE, anchorFile); + success = sendQuery(serverSocket, cmd, org, treeChoices, errFile); close(serverSocket); if (success) reportTiming(pStartTime, "send query and get response (successful)"); else reportTiming(pStartTime, "send query and get response (failed)"); } carefulClose(&errFile); } return success; } static int indexOfNull(char *stringArray[]) /* Return the index of the first NULL element of stringArray. * Do not call this unless stringArray has at least one NULL! */ { int ix = 0; while (stringArray[ix] != NULL) ix++; return ix; } #define MAX_SUBTREES 1000 struct usherResults *runUsher(char *org, char *usherPath, char *usherAssignmentsPath, char *vcfFile, int subtreeSize, struct slName **pUserSampleIds, - struct treeChoices *treeChoices, char *anchorFile, int *pStartTime) + struct treeChoices *treeChoices, char *anchorFile, + boolean ripplesEnabled, uint genomeSize, int *pStartTime) /* Open a pipe from Yatish Turakhia's usher program, save resulting big trees and * subtrees to trash files, return list of slRef to struct tempName for the trash files * and parse other results out of stderr output. The usher-sampled version of usher might * modify userSampleIds, adding a prefix if a sample with the same name is already in the tree. */ { struct usherResults *results = usherResultsNew(); char subtreeSizeStr[16]; safef(subtreeSizeStr, sizeof subtreeSizeStr, "%d", subtreeSize); struct tempName tnOutDir; trashDirFile(&tnOutDir, "ct", "usher_outdir", ".dir"); char *cmd[] = { usherPath, "-v", vcfFile, "-i", usherAssignmentsPath, "-d", tnOutDir.forCgi, "-k", subtreeSizeStr, "-K", SINGLE_SUBTREE_SIZE, "-u", - "-T", USHER_NUM_THREADS, // Don't pass args from -T onward to server - // Room for extra arguments if using usher-sampled - NULL, NULL, NULL, NULL, NULL, NULL, + // Room for extra arguments if using usher-sampled or server ripples search + NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL }; +int cmdBaseEnd = indexOfNull(cmd); +// Extra options for usher-sampled-server, which we'll try first: +int ix = cmdBaseEnd; +cmd[ix++] = "--no-ignore-prefix"; +cmd[ix++] = USHER_DEDUP_PREFIX; +if (isNotEmpty(anchorFile)) + { + cmd[ix++] = "--anchor-samples"; + cmd[ix++] = anchorFile; + } +if (ripplesEnabled) + { + cmd[ix++] = "--run-ripples"; + cmd[ix++] = "--ripples-ancestor-radius"; + cmd[ix++] = RIPPLES_ANCESTOR_RADIUS; + } struct tempName tnStderr; trashDirFile(&tnStderr, "ct", "usher_stderr", ".txt"); struct tempName tnServerStderr; trashDirFile(&tnServerStderr, "ct", "usher_server_stderr", ".txt"); char *stderrFile = tnServerStderr.forCgi; -if (! runUsherServer(org, cmd, tnServerStderr.forCgi, treeChoices, anchorFile, pStartTime)) +if (! runUsherServer(org, cmd, tnServerStderr.forCgi, treeChoices, pStartTime)) { + // Querying the server didn't work out; use command line, either usher or usher-sampled. + ix = cmdBaseEnd; + cmd[ix++] = "-T"; + cmd[ix++] = USHER_NUM_THREADS; if (endsWith(usherPath, "-sampled")) { // Add --no-ignore-prefix - int ix = indexOfNull(cmd); cmd[ix++] = "--no-ignore-prefix"; cmd[ix++] = USHER_DEDUP_PREFIX; // Add --anchor-samples if configured if (isNotEmpty(anchorFile)) { cmd[ix++] = "--anchor-samples"; cmd[ix++] = anchorFile; } // Add --optimization-radius 0 unless optimization is explicitly enabled char *enableOptimization = phyloPlaceOrgSetting(org, "enableOptimization"); if (SETTING_NOT_ON(enableOptimization)) { cmd[ix++] = "--optimization_radius"; cmd[ix++] = "0"; } } + cmd[ix++] = NULL; runUsherCommand(cmd, tnStderr.forCgi, anchorFile, pStartTime); stderrFile = tnStderr.forCgi; } struct tempName *singleSubtreeTn = NULL, *subtreeTns[MAX_SUBTREES]; struct variantPathNode *singleSubtreeMuts = NULL, *subtreeMuts[MAX_SUBTREES]; parsePlacements(tnOutDir.forCgi, stderrFile, results->samplePlacements, pUserSampleIds); int subtreeCount = processOutDirFiles(results, tnOutDir.forCgi, &singleSubtreeTn, - &singleSubtreeMuts, subtreeTns, subtreeMuts, MAX_SUBTREES); + &singleSubtreeMuts, subtreeTns, subtreeMuts, MAX_SUBTREES, genomeSize); if (singleSubtreeTn) { results->subtreeInfoList = parseSubtrees(subtreeCount, singleSubtreeTn, singleSubtreeMuts, subtreeTns, subtreeMuts, *pUserSampleIds); results->singleSubtreeInfo = results->subtreeInfoList; results->subtreeInfoList = results->subtreeInfoList->next; removeOutDir(tnOutDir.forCgi); } else { results = NULL; warn("Sorry, there was a problem running usher. " "Please ask genome-www@soe.ucsc.edu to take a look at %s.", stderrFile); } reportTiming(pStartTime, "parse results from usher"); @@ -1367,74 +1673,80 @@ struct pipeline *pl = pipelineOpen(cmds, pipelineRead, NULL, tnStderr.forCgi, 0); pipelineClose(&pl); reportTiming(pStartTime, "run matUtils command"); } static boolean runMatUtilsServer(char *org, char *protobufPath, char *subtreeSizeStr, struct tempName *tnSamples, struct tempName *tnOutDir, struct treeChoices *treeChoices, char *anchorFile, int *pStartTime) /* Cheng Ye added a 'matUtils mode' to usher-sampled-server so we can get subtrees super-fast * for uploaded sample names too. */ { boolean success = FALSE; char *cmd[] = { "usher-sampled-server", "-i", protobufPath, "-d", tnOutDir->forCgi, "-k", subtreeSizeStr, "-K", SINGLE_SUBTREE_SIZE, "--existing_samples", tnSamples->forCgi, "-D", - NULL }; + NULL, NULL, NULL }; +if (isNotEmpty(anchorFile)) + { + int ix = indexOfNull(cmd); + cmd[ix++] = "--anchor-samples"; + cmd[ix++] = anchorFile; + } struct tempName tnErrFile; trashDirFile(&tnErrFile, "ct", "matUtils_server_stderr", ".txt"); if (serverIsConfigured(org)) { FILE *errFile = mustOpen(tnErrFile.forCgi, "w"); int serverSocket = getServerSocket(org, treeChoices, errFile); reportTiming(pStartTime, "get socket"); if (serverSocket > 0) { - success = sendQuery(serverSocket, cmd, org, treeChoices, errFile, FALSE, anchorFile); + success = sendQuery(serverSocket, cmd, org, treeChoices, errFile); close(serverSocket); if (success) reportTiming(pStartTime, "send query and get response (successful)"); else reportTiming(pStartTime, "send query and get response (failed)"); } carefulClose(&errFile); } return success; } struct usherResults *runMatUtilsExtractSubtrees(char *org, char *matUtilsPath, char *protobufPath, int subtreeSize, struct slName *sampleIds, struct treeChoices *treeChoices, char *anchorFile, - int *pStartTime) + uint genomeSize, int *pStartTime) /* Open a pipe from Yatish Turakhia and Jakob McBroome's matUtils extract to extract subtrees * containing sampleIds, save resulting subtrees to trash files, return subtree results. * Caller must ensure that sampleIds are names of leaves in the protobuf tree. */ { struct usherResults *results = usherResultsNew(); char subtreeSizeStr[16]; safef(subtreeSizeStr, sizeof subtreeSizeStr, "%d", subtreeSize); struct tempName tnSamples; trashDirFile(&tnSamples, "ct", "matUtilsExtractSamples", ".txt"); FILE *f = mustOpen(tnSamples.forCgi, "w"); struct slName *sample; for (sample = sampleIds; sample != NULL; sample = sample->next) fprintf(f, "%s\n", sample->name); carefulClose(&f); struct tempName tnOutDir; trashDirFile(&tnOutDir, "ct", "matUtils_outdir", ".dir"); if (! runMatUtilsServer(org, protobufPath, subtreeSizeStr, &tnSamples, &tnOutDir, treeChoices, anchorFile, pStartTime)) runMatUtilsExtractCommand(matUtilsPath, protobufPath, subtreeSizeStr, &tnSamples, &tnOutDir, anchorFile, pStartTime); addEmptyPlacements(sampleIds, results->samplePlacements); struct tempName *singleSubtreeTn = NULL, *subtreeTns[MAX_SUBTREES]; struct variantPathNode *singleSubtreeMuts = NULL, *subtreeMuts[MAX_SUBTREES]; int subtreeCount = processOutDirFiles(results, tnOutDir.forCgi, &singleSubtreeTn, &singleSubtreeMuts, - subtreeTns, subtreeMuts, MAX_SUBTREES); + subtreeTns, subtreeMuts, MAX_SUBTREES, genomeSize); results->subtreeInfoList = parseSubtrees(subtreeCount, singleSubtreeTn, singleSubtreeMuts, subtreeTns, subtreeMuts, sampleIds); results->singleSubtreeInfo = results->subtreeInfoList; results->subtreeInfoList = results->subtreeInfoList->next; reportTiming(pStartTime, "process results from matUtils"); return results; }