0525b1b26752ba51653ce488a3d141ccf3eee4aa max Mon Jan 15 03:24:50 2024 -0800 adding note to VCF track if filtering or wiggling. Had to restructure filterRecords() a little to get rid of the early return, as otherwise labelTrackAsFilteredNumber() would not have been called if there are no filters but the track is in coverage mode. This seemed like the most minimal change, refs #32864 diff --git src/hg/hgTracks/vcfTrack.c src/hg/hgTracks/vcfTrack.c index b779bc9..568d4c6 100644 --- src/hg/hgTracks/vcfTrack.c +++ src/hg/hgTracks/vcfTrack.c @@ -186,55 +186,61 @@ nextRecord = rec->next; boolean allRef = TRUE; for (sample = sampleOrder; sample != NULL; sample = sample->next) { gt = vcfRecordFindGenotype(rec, sample->name); if (gt && !(gt->hapIxA == 0 && gt->hapIxB == 0) ) allRef = FALSE; } if (!allRef) slAddHead(&retList, rec); } slReverse(&retList); vcff->records = retList; } -static void filterRecords(struct vcfFile *vcff, struct trackDb *tdb) -/* If a filter is specified in the cart, remove any records that don't pass filter. */ +static void filterRecords(struct vcfFile *vcff, struct track *tg) +/* If a filter is specified in the cart, remove any records that don't pass filter. Adapt longLabel if something was filtered. */ { +struct trackDb *tdb = tg->tdb; double minQual = VCF_DEFAULT_MIN_QUAL; struct slName *filterValues = NULL; double minFreq = VCF_DEFAULT_MIN_ALLELE_FREQ; boolean gotQualFilter = getMinQual(tdb, &minQual); boolean gotFilterFilter = getFilterValues(tdb, &filterValues); boolean gotMinFreqFilter = getMinFreq(tdb, &minFreq); -if (! (gotQualFilter || gotFilterFilter || gotMinFreqFilter) ) - return; - +int filtOut = 0; +if (gotQualFilter || gotFilterFilter || gotMinFreqFilter) + { struct vcfRecord *rec, *nextRec, *newList = NULL; for (rec = vcff->records; rec != NULL; rec = nextRec) { nextRec = rec->next; if (! ((gotQualFilter && minQualFail(rec, minQual)) || (gotFilterFilter && filterColumnFail(rec, filterValues)) || (gotMinFreqFilter && minFreqFail(rec, minFreq)) )) slAddHead(&newList, rec); + else + filtOut++; } slReverse(&newList); vcff->records = newList; } +labelTrackAsFilteredNumber(tg, filtOut); +} + struct pgSnpVcfStartEnd /* This extends struct pgSnp by tacking on an original VCF chromStart and End at the end, * for use by indelTweakMapItem below. This can be cast to pgs. */ { struct pgSnp pgs; unsigned int vcfStart; unsigned int vcfEnd; }; static struct pgSnp *vcfFileToPgSnp(struct vcfFile *vcff, struct trackDb *tdb) /* Convert vcff's records to pgSnp; don't free vcff until you're done with pgSnp * because it contains pointers into vcff's records' chrom. If the trackDb setting * sampleName is present, then check whether all the records are phased or not */ { struct pgSnp *pgsList = NULL; @@ -2278,60 +2284,62 @@ hFreeConn(&conn); } if (isEmpty(fileOrUrl)) return; int vcfMaxErr = -1; struct vcfFile *vcff = NULL; /* protect against temporary network error */ struct errCatch *errCatch = errCatchNew(); if (errCatchStart(errCatch)) { vcff = vcfTabixFileAndIndexMayOpen(fileOrUrl, tbiFileOrUrl, chromName, winStart, winEnd, vcfMaxErr, -1); if (vcff != NULL) { - filterRecords(vcff, tg->tdb); + filterRecords(vcff, tg); filterRefOnlyAlleles(vcff, tg->tdb); // remove items that don't differ from reference // TODO: in multi-region mode, different windows end up with different sets of variants where // all are phased or not, which throws off track heights in each window. Similar to hapCluster // mode, just switch to pgSnp view when in multi-region for now. if (slCount(windows) > 1 || tg->visibility == tvDense) pgSnpMethods(tg); tg->items = vcfFileToPgSnp(vcff, tg->tdb); // pgSnp bases coloring/display decision on count of items: tg->extraUiData = vcff; tg->customInt = slCount(tg->items); // Don't vcfFileFree here -- we are using its string pointers! } else { if (tbiFileOrUrl) errAbort("Unable to open VCF file/URL '%s' with tabix index '%s'", fileOrUrl, tbiFileOrUrl); else errAbort("Unable to open VCF file/URL '%s'", fileOrUrl); } } errCatchEnd(errCatch); if (errCatch->gotError || vcff == NULL) { if (isNotEmpty(errCatch->message->string)) + { tg->networkErrMsg = cloneString(errCatch->message->string); tg->drawItems = bigDrawWarning; tg->totalHeight = bigWarnTotalHeight; } + } errCatchFree(&errCatch); } static void vcfPhasedAddMapBox(struct hvGfx *hvg, struct vcfRecord *rec, struct pgSnpVcfStartEnd *psvs, char *text, int x, int y, int width, int height, struct track *track) // Add mapbox for a tick in the vcfPhased track { mapBoxHgcOrHgGene(hvg, psvs->vcfStart, psvs->vcfEnd, x, y, width, height, track->track, rec->name, text, NULL, TRUE, NULL); } struct hapDistanceMatrixCell { /* The pair of alleles and how distant they are */ struct hapDistanceMatrixCell *next; char *sampleId; // name of this sample char *otherId; // name of other sample @@ -3027,31 +3035,31 @@ if (isEmpty(fileOrUrl)) return; fileOrUrl = hReplaceGbdb(fileOrUrl); int vcfMaxErr = -1; struct vcfFile *vcff = NULL; boolean hapClustEnabled = cartOrTdbBoolean(cart, tg->tdb, VCF_HAP_ENABLED_VAR, TRUE); if (slCount(windows)>1) hapClustEnabled = FALSE; // haplotype sorting display not currently available with multiple windows. /* protect against temporary network error */ struct errCatch *errCatch = errCatchNew(); if (errCatchStart(errCatch)) { vcff = vcfTabixFileAndIndexMayOpen(fileOrUrl, tbiFileOrUrl, chromName, winStart, winEnd, vcfMaxErr, -1); if (vcff != NULL) { - filterRecords(vcff, tg->tdb); + filterRecords(vcff, tg); int vis = tdbVisLimitedByAncestors(cart,tg->tdb,TRUE,TRUE); if (hapClustEnabled && vcff->genotypeCount > 1 && (vis == tvPack || vis == tvSquish)) vcfHapClusterOverloadMethods(tg, vcff); else { tg->items = vcfFileToPgSnp(vcff, tg->tdb); // pgSnp bases coloring/display decision on count of items: tg->customInt = slCount(tg->items); } // Don't vcfFileFree here -- we are using its string pointers! } else { if (tbiFileOrUrl) @@ -3101,31 +3109,31 @@ conn = hAllocConnTrack(database, tg->tdb); else { conn = hAllocConn(CUSTOM_TRASH); table = ct->dbTableName; } char *vcfFile = bbiNameFromSettingOrTable(tg->tdb, conn, table); hFreeConn(&conn); /* protect against parse error */ struct errCatch *errCatch = errCatchNew(); if (errCatchStart(errCatch)) { vcff = vcfFileMayOpen(vcfFile, chromName, winStart, winEnd, vcfMaxErr, -1, TRUE); if (vcff != NULL) { - filterRecords(vcff, tg->tdb); + filterRecords(vcff, tg); int vis = tdbVisLimitedByAncestors(cart,tg->tdb,TRUE,TRUE); if (hapClustEnabled && vcff->genotypeCount > 1 && vcff->genotypeCount < 3000 && (vis == tvPack || vis == tvSquish)) vcfHapClusterOverloadMethods(tg, vcff); else { tg->items = vcfFileToPgSnp(vcff, tg->tdb); // pgSnp bases coloring/display decision on count of items: tg->customInt = slCount(tg->items); } // Don't vcfFileFree here -- we are using its string pointers! } else errAbort("Unable to open VCF file '%s'", vcfFile); }