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/phyloPlace.c src/hg/hgPhyloPlace/phyloPlace.c index 15a5b24d791..521b1cc42b2 100644 --- src/hg/hgPhyloPlace/phyloPlace.c +++ src/hg/hgPhyloPlace/phyloPlace.c @@ -2254,30 +2254,125 @@ { if (gotClades) printf("<td>n/a</td>"); if (gotLineages) printf("<td>n/a</td>"); } // pangolin-assigned lineage char *lineage = lineageForSample(sms, si->name); printLineageTd(lineage, "n/a", db); // Maybe also #mutations with mouseover to show mutation path? printSubtreeTd(results->subtreeInfoList, jsonTns, si->name, subtreeSize); } puts("</tbody></table><p></p>"); } +static void summarizeRecombinants(struct recombinantInfo *recombinants, struct hash *descendants, struct slName *sampleIds) +/* Show a table describing potential recombinants found by ripples search. */ +{ +puts("<p><b>Potential recombinants found: </b><br>"); +puts("<table class='seqSummary'>"); +puts("<thead><tr><th>Potential recombinant node" + TOOLTIP("Sequence or internal node for which potential recombination was found. " + "If internal node, then the uploaded sequence(s) descended from the node are listed, " + "and if there are additional sequences already at that node then their number is " + "given") + "</th>"); +puts("<th>Placement clade" + TOOLTIP("Nextstrain clade in which the potential recombinant was placed") + "</th>"); +puts("<th>Placement lineage" + TOOLTIP("Pango lineage in which the potential recombinant was placed") + "</th>"); +puts("<th>Parsimony score improvement" + TOOLTIP("The number of mutations that could be explained by recombination, " + "given along with the original number of mutations on the branch " + "to the potential recombinant node") + "</th>"); +puts("<th>Donor node" + TOOLTIP("The parent node that explains mutations before the first breakpoint range " + "(and, if applicable, after the second breakpoint range)") + "</th>"); +puts("<th>Donor descendant count" + TOOLTIP("The number of sequences at the tips/leaves of the branch that starts " + "with the donor node") + "</th>"); +puts("<th>Donor clade" + TOOLTIP("Nextstrain clade in which the donor node was placed") + "</th>"); +puts("<th>Donor lineage" + TOOLTIP("Pango lineage in which the donor node was placed") + "</th>"); +puts("<th>Acceptor node" + TOOLTIP("The parent node that explains mutations after the first breakpoint range " + "(and, if applicable, before the second breakpoint range)") + "</th>"); +puts("<th>Acceptor descendant count" + TOOLTIP("The number of sequences at the tips/leaves of the branch that starts " + "with the acceptor node") + "</th>"); +puts("<th>Acceptor clade" + TOOLTIP("Nextstrain clade in which the acceptor node was placed") + "</th>"); +puts("<th>Acceptor lineage" + TOOLTIP("Pango lineage in which the donor node was placed") + "</th>"); +puts("<th>Breakpoint 1 range" + TOOLTIP("Coordinate range, bounded by mutations explained by the donor node at lower " + "coordinates and by the acceptor node at higher coordinates") + "</th>"); +puts("<th>Breakpoint 2 range" + TOOLTIP("Coordinate range, bounded by mutations explained by the acceptor node at lower " + "coordinates and by the donor node at higher coordinates") + "</th></tr></thead><tbody>"); +struct recombinantInfo *ri; +for (ri = recombinants; ri != NULL; ri = ri->next) + { + printf("<tr><td>%s", ri->recombNodeId); + if (startsWith("node_", ri->recombNodeId)) + { + printf(" ("); + struct slName *descendantList = hashFindVal(descendants, ri->recombNodeId), *desc; + int sampleCount = 0; + for (desc = descendantList; desc != NULL; desc = desc->next) + { + if (slNameInList(sampleIds, desc->name)) + { + if (sampleCount > 0) + printf(", "); + printf("%s", desc->name); + sampleCount++; + } + } + int others = ri->recombNumDesc - sampleCount; + if (others) + printf(" + %d", others); + printf(")"); + } + printf("</td><td>%s</td><td>%s</td>", ri->recombClade, ri->recombLineage); + printf("<td>%d (from %d)</td>", ri->parsimonyImprovement, ri->originalParsimony); + printf("<td>%s</td><td>%d</td><td>%s</td><td>%s</td>", + ri->donorNodeId, ri->donorNumDesc, ri->donorClade, ri->donorLineage); + printf("<td>%s</td><td>%d</td><td>%s</td><td>%s</td>", + ri->acceptorNodeId, ri->acceptorNumDesc, ri->acceptorClade, ri->acceptorLineage); + printf("<td>(%d, %d)</td><td>(%d, %d)</td></tr>\n", + ri->bp1Min, ri->bp1Max, ri->bp2Min, ri->bp2Max); + } +puts("</tbody></table>"); +puts("</p>"); +} + static struct singleNucChange *sncListFromSampleMutsAndImputed(struct slName *sampleMuts, struct baseVal *imputedBases, struct seqWindow *gSeqWin) /* Convert a list of "<ref><pos><alt>" names to struct singleNucChange list. * However, if <alt> is ambiguous, skip it because variantProjector doesn't like it. * Add imputed base predictions. */ { struct singleNucChange *sncList = NULL; struct slName *mut; for (mut = sampleMuts; mut != NULL; mut = mut->next) { char ref = mut->name[0]; if (ref < 'A' || ref > 'Z') errAbort("sncListFromSampleMuts: expected ref base value, got '%c' in '%s'", ref, mut->name); @@ -3420,41 +3515,44 @@ reportUnmatched(unmatched, (sampleIds == NULL)); slNameFreeList(&unmatched); reportTiming(&startTime, "look up uploaded sample names"); } lineFileClose(&lf); if (sampleIds == NULL) { return; } // Kick off child thread to load metadata simultaneously with running usher or matUtils. pthread_t *metadataPthread = mayStartLoaderPthread(metadataFile, loadMetadataWorker); struct usherResults *results = NULL; char *anchorFile = phyloPlaceRefSettingPath(org, refName, "anchorSamples"); +char *ripplesEnabledSetting = phyloPlaceOrgSetting(org, "ripplesEnabled"); +boolean ripplesEnabled = (isNotEmpty(ripplesEnabledSetting) && SETTING_IS_ON(ripplesEnabledSetting)); +boolean doRipples = (ripplesEnabled && slCount(sampleIds) <= MAX_RIPPLES_SEARCH); if (vcfTn) { fflush(stdout); results = runUsher(org, usherPath, protobufPath, vcfTn->forCgi, subtreeSize, &sampleIds, - treeChoices, anchorFile, &startTime); + treeChoices, anchorFile, doRipples, refGenome->size, &startTime); } else if (subtreesOnly) { char *matUtilsPath = getMatUtilsPath(TRUE); results = runMatUtilsExtractSubtrees(org, matUtilsPath, protobufPath, subtreeSize, - sampleIds, treeChoices, anchorFile, &startTime); + sampleIds, treeChoices, anchorFile, refGenome->size, &startTime); } struct sampleMetadataStore *sampleMetadata = NULL; if (metadataPthread) { pthreadJoin(metadataPthread, (void **)(&sampleMetadata)); reportTiming(&startTime, "wait for sample metadata loading thread to finish"); } else { // We really need metadata -- load it the slow way. sampleMetadata = getSampleMetadata(metadataFile); reportTiming(&startTime, "load sample metadata without pthread"); } @@ -3523,57 +3621,73 @@ "</p>\n"); puts("<p><em>Note: " "The Nextstrain subtree views, and Download files below, are temporary files and will " "expire within two days. " "Please download the Nextstrain subtree JSON files if you will want to view them " "again in the future. The JSON files can be drag-dropped onto " "<a href='https://auspice.us/' target=_blank>https://auspice.us/</a>." "</em></p>"); struct tempName *tsvTn = NULL, *sTsvTn = NULL; struct tempName *zipTn = makeSubtreeZipFile(results, jsonTns, singleSubtreeJsonTn, &startTime); struct tempName *ctTn = NULL; if (subtreesOnly) { + if (ripplesEnabled) + puts("<p><b>Did not search for potential recombinants (currently supported only " + "for uploaded sequences, not names/IDs).</b></p>"); summarizeSubtrees(sampleIds, results, sampleMetadata, jsonTns, refName, subtreeSize); reportTiming(&startTime, "describe subtrees"); } else { findNearestNeighbors(results, sampleMetadata); reportTiming(&startTime, "find nearest neighbors"); char *singleSubtreeFile = results->singleSubtreeInfo->subtreeTn->forCgi; struct phyloTree *sampleTree = uploadedSamplesTree(singleSubtreeFile, sampleIds); if (canDoCustomTracks) { // Make custom tracks for uploaded samples and subtree(s). ctTn = writeCustomTracks(org, refName, ctDb, vcfTn, results, sampleIds, source, tl->fontHeight, sampleTree, &startTime); } // Make a sample summary TSV file and accumulate S gene changes struct hash *seqInfoHash = hashFromSeqInfoListAndIds(seqInfoList, sampleIds); addSampleMutsFromSeqInfo(results->samplePlacements, seqInfoHash); struct hash *spikeChanges = hashNew(0); tsvTn = writeTsvSummary(results, sampleTree, sampleIds, seqInfoHash, geneInfoList, gSeqWin, spikeChanges, &startTime); sTsvTn = writeSpikeChangeSummary(spikeChanges, slCount(sampleIds)); downloadsRow(results->bigTreePlusTn->forHtml, tsvTn->forHtml, sTsvTn->forHtml, zipTn->forHtml); + if (results->recombinants) + { + summarizeRecombinants(results->recombinants, results->recombinantDescendants, sampleIds); + } + else if (ripplesEnabled) + { + if (doRipples) + puts("<p><b>No potential recombinants found.</b></p>"); + else + printf("<p><b>Did not search for potential recombinants, too many input sequences " + "(maximum %d for recombinant search).</b></p>\n", MAX_RIPPLES_SEARCH); + } + int seqCount = slCount(seqInfoList); if (seqCount <= MAX_SEQ_DETAILS) { char *refAcc = cloneString(refGenome->name); if (regexMatch(refAcc, "v[0-9]+$")) { char *v = strrchr(refAcc, 'v'); assert(v != NULL); *v = '.'; } summarizeSequences(seqInfoList, isFasta, results, jsonTns, refAcc, refName, subtreeSize); reportTiming(&startTime, "write summary table (including reading in lineages)"); for (ix = 0, ti = results->subtreeInfoList; ti != NULL; ti = ti->next, ix++) { int subtreeUserSampleCount = slCount(ti->subtreeUserSampleIds);