bd53b8dfe0a86f9d1c8ed57289fd409a6c52bd6b chmalee Fri Jun 26 10:54:41 2020 -0700 Add option to shade alleles by inconsitent transmission between child and parent haplotypes diff --git src/hg/hgTracks/vcfTrack.c src/hg/hgTracks/vcfTrack.c index baa9560..d513525 100644 --- src/hg/hgTracks/vcfTrack.c +++ src/hg/hgTracks/vcfTrack.c @@ -1936,56 +1936,103 @@ } } } // now fill out a distance matrix based on the cwaDistance between all the pairs in hapArray struct hapDistanceMatrix *hapDistMatrix = fillOutDistanceMatrix(hapArray, vcff, sample, &helper, gtCount); unsigned int hapCount = 2 * gtCount; unsigned int *hapOrder = needMem(sizeof(unsigned int) * hapCount); setHapOrderFromMatrix(hapOrder, hapCount, hapDistMatrix, hapArray, vcff, sample, sampleDrawOrder); return hapOrder; } #define VCFPHASED_UNPHASED_COLOR MG_RED +static int getChildAlleleColor(struct track *track, struct vcfRecord *rec, int childHapIx, + int alleleIx, int nameIx, int gtHapCount, unsigned int *gtHapOrder, + char *sampleOrder[]) +/* Return the color we should use for this variant, depending on whether the allele + * of the child matches what the parent transmitted or if there was some kind of + * 'mendelian inconsistency' */ +{ +// if there are no parents to draw then no concept of transmitted vs unstransmitted +if (gtHapCount <= 2) + return MG_BLACK; +int parIx = 0; +if (nameIx == 0) + { + if (alleleIx % 2 == 0) + parIx = alleleIx + 2; + else + parIx = alleleIx + 3; + } +else if (nameIx == 1) + { + if (alleleIx % 2 == 0) + parIx = alleleIx - 1; + else + parIx = alleleIx + 1; + } +else + { + if (alleleIx % 2 == 0) + parIx = alleleIx - 3; + else + parIx = alleleIx - 2; + } +struct vcfGenotype *parentGt = &(rec->genotypes[gtHapOrder[parIx] >> 1]); +if (parentGt->isPhased || (parentGt->hapIxA == 1 && parentGt->hapIxB == 1)) + { + int parentAlleleIx = gtHapOrder[parIx] & 1 ? parentGt->hapIxB : parentGt->hapIxA; + if (rec->alleles[childHapIx] != rec->alleles[parentAlleleIx]) // if they disagree mark as red + return MG_RED; + else + return MG_BLACK; + } +return MG_BLACK; +} + void vcfPhasedDrawOneRecord(struct track *track, struct hvGfx *hvg, struct vcfRecord *rec, void *item, unsigned int *gtHapOrder, int gtHapCount, int xOff, int yOffsets[], - char *sampleOrder[], double scale) + char *sampleOrder[], char *childSample, boolean highlightChildDiffs, + double scale) // Draw a record's haplotypes on the appropriate lines { int i; int x1 = round((double)(rec->chromStart-winStart)*scale) + xOff; int x2 = round((double)(rec->chromEnd-winStart)*scale) + xOff; struct pgSnpVcfStartEnd *psvs = (struct pgSnpVcfStartEnd *)item; int w = x2-x1; if (w < 1) w = 1; int tickHeight = track->itemHeight(track, track->items); for (i = 0; i < gtHapCount ; i++) { struct vcfGenotype *gt = &(rec->genotypes[gtHapOrder[i] >> 1]); int nameIx = stringArrayIx(gt->id, sampleOrder, track->customInt); struct dyString *mouseover = dyStringNew(0); if (gt->isPhased || (gt->hapIxA == 1 && gt->hapIxB == 1)) // if phased or homozygous alt { int alIx = gtHapOrder[i] & 1 ? gt->hapIxB : gt->hapIxA; int tickColor = MG_BLACK; + if (sameString(childSample, gt->id) && highlightChildDiffs) + tickColor = getChildAlleleColor(track, rec, alIx, i, nameIx, gtHapCount, gtHapOrder, sampleOrder); dyStringPrintf(mouseover, "%s ", gt->id); if (alIx != 0) // non-reference allele { - if (i > 1 && i < 4) // we're drawing child ticks + if (sameString(childSample, gt->id)) // we're drawing child ticks { dyStringPrintf(mouseover, "allele: %s", rec->alleles[alIx]); } else { dyStringPrintf(mouseover, "%s allele: %s", i == 0 || i == 5 ? "untransmitted" : "transmitted", rec->alleles[alIx]); } int y = yOffsets[i] - (tickHeight/2); hvGfxBox(hvg, x1, y, w, tickHeight, tickColor); vcfPhasedAddMapBox(hvg, rec, psvs, mouseover->string, x1, y, w, tickHeight, track); } } else if (gt->hapIxA != 0 || gt->hapIxB != 0)// draw the tick between the two haplotype lines { int yMid = ((yOffsets[2*nameIx] + yOffsets[(2*nameIx)+1]) / 2); // midpoint of two haplotype lines @@ -2076,35 +2123,36 @@ int i; for (pair = sampleNames, i = 0; pair != NULL && i < gtCount; pair = pair->next, i++) sampleOrder[i] = pair->name; // child sample is required to be in the middle/bottom of the trio char *childSample = gtCount > 2 ? sampleOrder[1] : sampleOrder[0]; // set up the "haplotype" lines and the transparent yellow box to delineate samples vcfPhasedSetupHaplotypesLines(track, hvg, xOff, yOff, width, yOffsets, sampleNames, font); // sort the variants by haplotype then draw ticks int nRecords = slCount(vcff->records); int centerIx = getCenterVariantIx(track, seqStart, seqEnd, vcff->records); int startIx = 0; int endIx = nRecords; unsigned int *hapOrder = computeHapDist(vcff, centerIx, startIx, endIx, childSample, gtCount, sampleOrder); +boolean highlightChildDiffs = cartUsualBooleanClosestToHome(cart, track->tdb, FALSE, VCF_PHASED_HIGHLIGHT_INCONSISTENT, FALSE); struct vcfRecord *rec = NULL; struct slList *item = NULL; for (rec = vcff->records, item = track->items; rec != NULL && item != NULL; rec = rec->next, item = item->next) { - vcfPhasedDrawOneRecord(track, hvg, rec, item, hapOrder, gtCount * 2, xOff, yOffsets, sampleOrder, scale); + vcfPhasedDrawOneRecord(track, hvg, rec, item, hapOrder, gtCount * 2, xOff, yOffsets, sampleOrder, childSample, highlightChildDiffs, scale); } } static int vcfPhasedItemHeight(struct track *tg, void *item) { return tg->heightPer; } int vcfPhasedTrackHeight(struct track *tg, enum trackVisibility vis) { const struct vcfFile *vcff = tg->extraUiData; // when doing composite track height, vcfPhasedLoadItems won't have been called yet! if (vis == tvDense) return pgSnpHeight(tg, vis); if (!vcff || vcff->records == NULL)