6974f14d41a28ed1198110de47946370cfc2bb73 angie Fri Nov 22 10:29:44 2019 -0800 dbSnp153: When different freq alleles can be expanded by different amounts, pad out the smaller range & alleles for consistency. refs #23283 One variant (rs782394990) was dropped because of that quirk, but now no variants are dropped. :) A few ucscNotes increased by 1 or 2. Some of the counts of various ucscNotes were a little out of date due to recent work on adding freq notes instead of dropping; up to date now. I was all proud of this but it's a dead end. One b153 variant, rs782394990, errored out because its freq alleles couldn't be expanded by the same amount. The insertions of pure A's could be expanded to a larger range than an ins that included a T in addition to A's. So I siezed upon the problem of independently trimmed and expanded SPDIs resulting in inconsistent dels, and set about trimming by list instead of trimming each SPDI separately. I would restrict the expansion to the minimum range by which any allele could be shifted. That would fix the problem of independent del's for rs782394990. However, there was a problem. A variant that includes both an ins and a del that could be expanded to the same range, e.g. rs201454468 with alleles ref=T, alt="", alt=TT, cannot be trimmed as a list because it already has "". T/"" is already minimal. However, T/TT is not already minimal and needs to be minimized before expansion. So no expansion was performed, and we ended up with inconsistent alleles between different sources because one was expanded (it had only the del) and the other wasn't. No good. I think that's why, although I was trimming listwise already, I did another trim just before trying to expand. Is there even a point in trimming listwise? So... back to individual trimming of each SPDI. And individual expansion. But what do we do about the 1 in 680M case of rs782394990? I think we should detect inconsistent dels -- and then pad out all dels and inses to the maximal range. So for the allele that could be shifted 2 bases less than the others because it included a T, just pad its del and ins with genomic bases. OK, now it works. Squash branch into master... diff --git src/hg/snp/dbSnpJsonToTab/dbSnpJsonToTab.c src/hg/snp/dbSnpJsonToTab/dbSnpJsonToTab.c index cb46c1d..06d244d 100644 --- src/hg/snp/dbSnpJsonToTab/dbSnpJsonToTab.c +++ src/hg/snp/dbSnpJsonToTab/dbSnpJsonToTab.c @@ -239,30 +239,47 @@ warn("Mismatching SPDI pos (%s: %d vs %d) for %s", spdiB->chrom, spdiBList->chromStart, spdiB->chromStart, rsId); return FALSE; } if (differentString(spdiBList->del, spdiB->del)) { warn("Mismatching SPDI del (%s:%d: '%s' vs '%s') for %s", spdiB->chrom, spdiB->chromStart, spdiBList->del, spdiB->del, rsId); return FALSE; } } } return TRUE; } +static void spdiBedListAssertSameRef(struct spdiBed *spdiBList) +/* Sanity check -- raise assertion if spdiBList elements do not all have the same seq, pos, del. */ +{ +if (spdiBList && spdiBList->next) + { + assert(spdiBList->chromEnd == spdiBList->chromStart + strlen(spdiBList->del)); + struct spdiBed *spdiB; + for (spdiB = spdiBList->next; spdiB != NULL; spdiB = spdiB->next) + { + assert(spdiB->chrom == spdiBList->chrom || sameString(spdiB->chrom, spdiBList->chrom)); + assert(spdiB->chromStart == spdiBList->chromStart); + assert(spdiB->chromEnd == spdiBList->chromEnd); + assert(spdiB->del == spdiBList->del || sameString(spdiB->del, spdiBList->del)); + } + } +} + static enum bigDbSnpClass varTypeToClass(char *varType) /* Convert JSON variant_type to bigDbSnp enum form. */ { if (sameString(varType, "snv")) return bigDbSnpSnv; else if (sameString(varType, "mnv")) return bigDbSnpMnv; else if (sameString(varType, "ins")) return bigDbSnpIns; else if (sameString(varType, "del")) return bigDbSnpDel; else if (sameString(varType, "delins")) return bigDbSnpDelins; else if (sameString(varType, "identity")) return bigDbSnpIdentity; @@ -517,38 +534,40 @@ if (twoBitChrom != NULL) { struct dnaSeq *dnaSeq = twoBitLoadAll(twoBitChrom); if (dnaSeq != NULL) { // Use memSeqWindow instead of twoBitSeqWindow to avoid keeping a bunch of open // twoBitFile handles. // Hmm, a dnaSeqWindow would be more efficient, no cloning and strlen... // or at the least, memSeqWindowNew could take a size arg. seqWin = memSeqWindowNew(twoBitChrom, dnaSeq->dna); dnaSeqFree(&dnaSeq); hashAdd(ncToSeqWin, nc, seqWin); } } } +if (seqWin == NULL) + errAbort("No twoBit file and chrom found in refSeqToUcsc file for seq '%s'", nc); return seqWin; } static boolean maybeExpandRange(struct spdiBed *spdiBTrim, struct seqWindow *seqWin, - struct spdiBed *spdiB, struct alObs *obs, struct lm *lm) + struct spdiBed *spdiB, struct lm *lm) /* If spdiBTrim is a pure insertion or deletion within a repetitive region, - * then expand spdiB (and obs if non-NULL) to the full repetitive range and return TRUE. - * Trim spdiBTrim to minimal representation (see maybeTrimSpdi) before calling this. */ + * then expand spdiB to the full repetitive range and return TRUE. + * spdiBTrim is assumed to already have minimal representation (see trimSpdi). */ { boolean changed = FALSE; uint refTrimStart = spdiBTrim->chromStart; uint refTrimEnd = spdiBTrim->chromEnd; char *altTrim = spdiBTrim->ins; int refTrimLen = refTrimEnd - refTrimStart; int altTrimLen = strlen(altTrim); if (indelShiftIsApplicable(refTrimLen, altTrimLen)) { // Now shift the minimal representation left and right and compare with spdiB coords. uint leftStart = refTrimStart, leftEnd = refTrimEnd; char leftAlt[altTrimLen+1]; safecpy(leftAlt, sizeof leftAlt, altTrim); int leftShifted = indelShift(seqWin, &leftStart, &leftEnd, leftAlt, INDEL_SHIFT_NO_MAX, isdLeft); @@ -560,263 +579,213 @@ if (leftShifted || rightShifted) { changed = TRUE; // Expand spdiB->del to cover the whole range., int newRefLen = rightEnd - leftStart; spdiB->del = lmAlloc(lm, newRefLen+1); seqWindowCopy(seqWin, leftStart, newRefLen, spdiB->del, newRefLen+1); // Concatenate ref bases to the left, altTrim, ref bases to the right -> new alt char refLeft[leftShifted+1]; seqWindowCopy(seqWin, leftStart, leftShifted, refLeft, sizeof refLeft); char refRight[rightShifted+1]; seqWindowCopy(seqWin, refTrimEnd, rightShifted, refRight, sizeof refRight); int newAltLen = leftShifted + altTrimLen + rightShifted; spdiB->ins = lmAlloc(lm, newAltLen+1); safef(spdiB->ins, newAltLen+1, "%s%s%s", refLeft, altTrim, refRight); - if (obs) - obs->allele = spdiB->ins; spdiB->chromStart = leftStart; spdiB->chromEnd = rightEnd; } } return changed; } -static void checkSpdisAfterExpansion(struct spdiBed *spdiBList, struct alObs *obsList) -/* Check spdiBList and obsList for inconsistencies; if the reference allele spdiB+obs missed out - * on the expansion then expand them. */ +static void checkObsVsSpdi(struct alObs *obsList, struct spdiBed *spdiBList) +/* assert that obsList and spdiBList have same length and each obs->allele is the same as the + * corresponding spdiB->ins. */ { -boolean alreadyExpandedRefAllele = FALSE; +assert(slCount(spdiBList) == slCount(obsList)); struct spdiBed *spdiB; struct alObs *obs; for (spdiB = spdiBList, obs = obsList; spdiB != NULL; spdiB = spdiB->next, obs = obs->next) { - if (strlen(spdiB->del) != spdiB->chromEnd - spdiB->chromStart) - errAbort("Inconsistent freq spdi coords & del: %s|%d|%d %s/%s", - spdiB->chrom, spdiB->chromStart, spdiB->chromEnd, spdiB->del, spdiB->ins); - if (spdiB != spdiBList) - { - if (differentString(spdiB->chrom, spdiBList->chrom)) - errAbort("Differing freq spdiB chrom: '%s' vs. '%s'", - spdiB->chrom, spdiBList->chrom); - if (differentString(spdiB->del, spdiBList->del)) - { - // Ref allele can't be shifted, so expand it to alt's shifted range... but only once, - // so we'll catch differences between alt alleles' expanded ranges. - struct spdiBed *spdiBRef = NULL, *spdiBAlt = NULL; - struct alObs *obsRef = NULL; - if (sameString(spdiBList->del, spdiBList->ins)) - { - spdiBRef = spdiBList; - spdiBAlt = spdiB; - obsRef = obsList; - } - else if (sameString(spdiB->del, spdiB->ins)) - { - spdiBRef = spdiB; - spdiBAlt = spdiBList; - obsRef = obs; - } - if (spdiBRef != NULL && !alreadyExpandedRefAllele) - { - // Copy spdiBAlt's expanded coords and del into spdiBRef - spdiBRef->chromStart = spdiBAlt->chromStart; - spdiBRef->chromEnd = spdiBAlt->chromEnd; - spdiBRef->del = spdiBAlt->del; - // The reference allele's ins is the same as its del, and the corresponding - // obs->allele must match it. - spdiBRef->ins = spdiBRef->del; - obsRef->allele = spdiBRef->ins; - alreadyExpandedRefAllele = TRUE; - } - else - errAbort("Differing freq spdiB: %s|%d|%d %s/%s vs %s|%d|%d %s/%s", - spdiBList->chrom, spdiBList->chromStart, spdiBList->chromEnd, - spdiBList->del, spdiBList->ins, - spdiB->chrom, spdiB->chromStart, spdiB->chromEnd, - spdiB->del, spdiB->ins); - } - } + assert(sameString(obs->allele, spdiB->ins)); } } static struct spdiBed *trimSpdi(struct spdiBed *spdiB, struct lm *lm) /* Return a new spdiBed that is the minimal representation of spdiB. */ { struct spdiBed *spdiBTrim = spdiBedNewLm(spdiB->chrom, spdiB->chromStart, lmCloneString(lm, spdiB->del), lmCloneString(lm, spdiB->ins), lm); int refTrimLen = 0, altTrimLen = 0; trimRefAlt(spdiBTrim->del, spdiBTrim->ins, &spdiBTrim->chromStart, &spdiBTrim->chromEnd, &refTrimLen, &altTrimLen); return spdiBTrim; } -static boolean spdiListTrimLeftBase(struct spdiBed *spdiList, struct lm *lm) -/* If all alleles in spdiList begin with the same base at the same pos, then modify them - * all to start one base to the right and trim the left base from all alleles. - * errAbort if the don't all have the same pos and del. */ +static void spdiBedListNormalizeRange(struct spdiBed *spdiBList, struct seqWindow *seqWin, + struct lm *lm) +/* If members of spdiBList have inconsistent chromStart/chromEnd/del then expand chromStart, + * chromEnd, del and ins to cover the full range of all elements in list. Ref allele is + * not included when finding range -- it is reset using an alt's final del. */ { -if (spdiList && spdiList->next) +if (spdiBList && spdiBList->next) { - boolean doTrim = TRUE; - char *del = spdiList->del; - char leftBase = del[0]; - uint start = spdiList->chromStart; - struct spdiBed *spdi; - for (spdi = spdiList; spdi != NULL; spdi = spdi->next) + struct spdiBed *spdiBRef = NULL, *spdiBAlt = NULL; + char *firstChrom = NULL; + int minChromStart = BIGNUM; + int maxChromEnd = 0; + boolean inconsistent = FALSE; + struct spdiBed *spdiB; + for (spdiB = spdiBList; spdiB != NULL; spdiB = spdiB->next) { - if (spdi->chromStart != start) - errAbort("spdiListTrimLeftBase: items have different starts (%s:%u, %s:%u)", - spdiList->chrom, start, spdi->chrom, spdi->chromStart); - if (differentString(spdi->del, del)) - errAbort("spdiListTrimLeftBase: items have different del (%s, %s)", del, spdi->del); - if (spdi->del[0] == 0 || spdi->ins[0] == 0 || - spdi->ins[0] != leftBase) + if (firstChrom == NULL) + firstChrom = spdiB->chrom; + else if (differentString(spdiB->chrom, firstChrom)) + errAbort("spdiBedListNormalizeRange: inconsistent chrom '%s' vs. '%s'", + firstChrom, spdiB->chrom); + if (spdiB->del == spdiB->ins || sameString(spdiB->del, spdiB->ins)) { - doTrim = FALSE; - break; - } + if (spdiBRef != NULL) + errAbort("spdiBedListNormalizeRange: multiple ref allele SPDIs not supported."); + spdiBRef = spdiB; } - if (doTrim) - { - for (spdi = spdiList; spdi != NULL; spdi = spdi->next) + else { - spdi->chromStart++; - // Allocate new strings because the spdi->del pointer can be shared all over. - spdi->del = lmCloneString(lm, spdi->del+1); - spdi->ins = lmCloneString(lm, spdi->ins+1); - } - } - return doTrim; + if (spdiBAlt == NULL) + spdiBAlt = spdiB; + else if (spdiB->chromStart != spdiBAlt->chromStart || + spdiB->chromEnd != spdiBAlt->chromEnd) + inconsistent = TRUE; + else if (! (spdiB->del == spdiBAlt->del || sameString(spdiB->del, spdiBAlt->del))) + errAbort("spdiBedListNormalizeRange: same coords %s|%d|%d, different del " + "'%s' / '%s'", + spdiB->chrom, spdiB->chromStart, spdiB->chromEnd, + spdiB->del, spdiBAlt->del); + if (spdiB->chromStart < minChromStart) + minChromStart = spdiB->chromStart; + if (spdiB->chromEnd > maxChromEnd) + maxChromEnd = spdiB->chromEnd; } -return FALSE; } - -static boolean spdiListTrimRightBase(struct spdiBed *spdiList, struct lm *lm) -/* If all alleles in spdiList end with the same base, then trim the right base from all alleles. - * errAbort if not all at same pos. */ -{ -if (spdiList && spdiList->next) - { - boolean doTrim = TRUE; - char *del = spdiList->del; - char rightBase = lastChar(del); - uint start = spdiList->chromStart; - struct spdiBed *spdi; - for (spdi = spdiList; spdi != NULL; spdi = spdi->next) - { - if (spdi->chromStart != start) - errAbort("spdiListTrimRightBase: items have different starts (%s:%u, %s:%u)", - spdiList->chrom, start, spdi->chrom, spdi->chromStart); - if (differentString(spdi->del, del)) - errAbort("spdiListTrimRightBase: items have different del (%s, %s)", del, spdi->del); - if (spdi->del[0] == 0 || spdi->ins[0] == 0 || - lastChar(spdi->ins) != rightBase) + if (inconsistent) { - doTrim = FALSE; - break; - } - } - if (doTrim) + int newDelLen = maxChromEnd - minChromStart; + char *newDel = lmAlloc(lm, newDelLen+1); + seqWindowCopy(seqWin, minChromStart, newDelLen, newDel, newDelLen+1); + for (spdiB = spdiBList; spdiB != NULL; spdiB = spdiB->next) { - for (spdi = spdiList; spdi != NULL; spdi = spdi->next) + // Skip ref allele for now, it is updated at end. + if (sameString(spdiB->del, spdiB->ins)) + continue; + // If needed, expand chromStart/chromEnd and add bases from the reference to the + // left and/or right of del and ins. + int padLeft = spdiB->chromStart - minChromStart; + int padRight = maxChromEnd - spdiB->chromEnd; + if (padLeft > 0 || padRight > 0) + { + spdiB->chromStart = minChromStart; + spdiB->chromEnd = maxChromEnd; + spdiB->del = newDel; + // Concatenate left padding, old ins value, and right padding to make new ins. + int oldInsLen = strlen(spdiB->ins); + int newInsLen = padLeft + oldInsLen + padRight; + char *newIns = lmAlloc(lm, newInsLen+1); + safencpy(newIns, newInsLen+1, newDel, padLeft); + safencpy(newIns+padLeft, newInsLen+1-padLeft, spdiB->ins, oldInsLen); + safencpy(newIns+padLeft+oldInsLen, newInsLen+1-padLeft-oldInsLen, + newDel+newDelLen-padRight, padRight); + spdiB->ins = newIns; + } + } + } + // Now that alts are settled, fix ref allele. + if (spdiBRef && spdiBAlt && + (spdiBRef->chromStart != spdiBAlt->chromStart || spdiBRef->chromEnd != spdiBAlt->chromEnd)) { - spdi->chromEnd--; - // Modify clones because the spdi->del pointer can be shared all over. - spdi->del = lmCloneString(lm, spdi->del); - spdi->ins = lmCloneString(lm, spdi->ins); - trimLastChar(spdi->del); - trimLastChar(spdi->ins); + spdiBRef->chromStart = spdiBAlt->chromStart; + spdiBRef->chromEnd = spdiBAlt->chromEnd; + spdiBRef->del = spdiBRef->ins = spdiBAlt->del; } } - return doTrim; - } -return FALSE; } -static boolean trimVcfSpdiList(struct spdiBed *spdiList, struct alObs *obsList, struct lm *lm) -/* SPDIs from VCF indel allele frequency submissions often still have the extra left base... and - * sometimes they have an extra right base. Detect and trim those cases. */ +static void setObsToSpdi(struct alObs *obsList, struct spdiBed *spdiBList) +/* Set alleles in obsList to match ins in spdiBList. */ { -boolean changed = FALSE; -while (spdiListTrimLeftBase(spdiList, lm)) - changed = TRUE; -while (spdiListTrimRightBase(spdiList, lm)) - changed = TRUE; -if (changed) - { - struct spdiBed *spdi; +assert(slCount(spdiBList) == slCount(obsList)); +struct spdiBed *spdiB; struct alObs *obs; - // Propagate changes to obsList - for (spdi = spdiList, obs = obsList; - spdi != NULL; - spdi = spdi->next, obs = obs->next) +for (spdiB = spdiBList, obs = obsList; + spdiB != NULL; + spdiB = spdiB->next, obs = obs->next) { - obs->allele = spdi->ins; - } + obs->allele = spdiB->ins; } -return changed; } static void spdiNormalizeFreq(struct sharedProps *props, struct hash *ncToTwoBitChrom, struct hash *ncToSeqWin, struct lm *lm) /* Allele frequency counts come directly from VCF, so although they are encoded as SPDI, they do * not cover the entire repetitive range. Expand variant alleles on the sequence on which they * were originally reported; if there is an indel difference between the reporting assembly and * the mapping assembly then this will resolve the actual alleles correctly (sometimes ref and alt * may be different for different assemblies). */ { int sIx; for (sIx = 0; sIx < props->freqSourceCount; sIx++) { - struct spdiBed *spdiB, *spdiBList = props->freqSourceSpdis[sIx]; - struct alObs *obs, *obsList = props->freqSourceObs[sIx]; - // Sanity check that spdiBList and obsList are consistent before we start tweaking: - assert(slCount(spdiBList) == slCount(obsList)); - for (spdiB = spdiBList, obs = obsList; - spdiB != NULL; - spdiB = spdiB->next, obs = obs->next) - { - assert(sameString(obs->allele, spdiB->ins)); - } - boolean changed = trimVcfSpdiList(spdiBList, obsList, lm); - // For each allele, see if it can be slid around. - for (spdiB = spdiBList, obs = obsList; - spdiB != NULL; - spdiB = spdiB->next, obs = obs->next) + struct spdiBed *spdiBList = props->freqSourceSpdis[sIx]; + if (spdiBList != NULL) { - assert(sameString(obs->allele, spdiB->ins)); - // If spdiB->del can be expanded, then expand both ref and alt in spdiB and obs accordingly. - if (isEmpty(spdiB->del) || isEmpty(spdiB->ins) || - stringIn(spdiB->del, spdiB->ins) || stringIn(spdiB->ins, spdiB->del)) + spdiBedListAssertSameRef(spdiBList); + struct seqWindow *seqWin = getChromSeq(ncToSeqWin, spdiBList->chrom, ncToTwoBitChrom); + struct alObs *obsList = props->freqSourceObs[sIx]; + // Sanity check that spdiBList and obsList are consistent before we start tweaking: + checkObsVsSpdi(obsList, spdiBList); + boolean changed = FALSE; + struct spdiBed *spdiB; + for (spdiB = spdiBList; spdiB != NULL; spdiB = spdiB->next) { - // First trim ref and alt to their minimal version for indelShift to work with. - // There may be nothing to trim because they're usually minimal already and - // stringIn catches false positive like A/TAAC, and that's fine. + // The reference allele (del == ins) can't be shifted in a meaningful way -- + // skip it for now and see how the alt alleles (the actual changes) might be shifted. + if (sameString(spdiB->del, spdiB->ins)) + continue; struct spdiBed *spdiBTrim = trimSpdi(spdiB, lm); - // If possible, expand the minimal representation to full SPDI range. - struct seqWindow *seqWin = getChromSeq(ncToSeqWin, spdiB->chrom, ncToTwoBitChrom); - if (seqWin == NULL) - errAbort("No twoBit file and chrom found for freq seq '%s'", spdiB->chrom); - changed |= maybeExpandRange(spdiBTrim, seqWin, spdiB, obs, lm); + if (maybeExpandRange(spdiBTrim, seqWin, spdiB, lm)) + changed = TRUE; + else if (spdiBTrim->chromStart != spdiB->chromStart || + spdiBTrim->chromEnd != spdiB->chromEnd) + { + // Use the trimmed version + changed = TRUE; + spdiB->chromStart = spdiBTrim->chromStart; + spdiB->chromEnd = spdiBTrim->chromEnd; + spdiB->del = spdiBTrim->del; + spdiB->ins = spdiBTrim->ins; } } if (changed) - checkSpdisAfterExpansion(spdiBList, obsList); + { + spdiBedListNormalizeRange(spdiBList, seqWin, lm); + setObsToSpdi(obsList, spdiBList); + } + spdiBedListAssertSameRef(spdiBList); + } } } static double mafFromSource(struct sharedProps *props, int sourceIx, char **retMajorAllele, char **retMinorAllele, struct lm *lm) /* Compare allele counts from source and return minor allele frequency (2nd-highest). * Set *retMajorAllele to allele with highest count (for comparison with mapped ref alleles). * If there are no observations from source, return NAN (divide by zero). */ { char *majorAllele = NULL, *minorAllele = NULL; int majorAlleleCount = -1, minorAlleleCount = -1; double totalCount = 0; struct alObs *obs; for (obs = props->freqSourceObs[sourceIx]; obs != NULL; obs = obs->next) { @@ -1105,31 +1074,31 @@ * bds's SPDI-based coords cover the whole repetitive region as expected. */ { int i; for (i = 0; i < bds->altCount; i++) { char *alt = bds->alts[i]; if (isEmpty(bds->ref) || isEmpty(alt) || stringIn(bds->ref, alt) || stringIn(alt, bds->ref)) { // First trim ref and alt to their minimal version for indelShift to work with. // There may be nothing to trim because they're usually minimal already and // stringIn catches false positive like A/TAAC, and that's fine. struct spdiBed *spdiB = spdiBedNewLm(bds->chrom, bds->chromStart, bds->ref, alt, lm); struct spdiBed *spdiBTrim = trimSpdi(spdiB, lm); - maybeExpandRange(spdiBTrim, seqWin, spdiB, NULL, lm); + maybeExpandRange(spdiBTrim, seqWin, spdiB, lm); if (spdiB->chromStart != bds->chromStart || spdiB->chromEnd != bds->chromEnd) errAbort("Range of %s (%s|%d|%d ref='%s', alt='%s') " "could be expanded to %s|%d|%d\n", bds->name, bds->chrom, bds->chromStart, bds->chromEnd, bds->ref, alt, bds->chrom, spdiB->chromStart, spdiB->chromEnd); } } } static void addClinVarSigs(struct dyString *dyUcscNotes, struct sharedProps *props) /* If clinVarSigs indicate benign, pathogenic, or both (conflicting), add ucscNote. */ { boolean isBenign = FALSE, isPathogenic = FALSE; struct slName *sig; for (sig = props->clinVarSigs; sig != NULL; sig = sig->next)