b49e61be4ad54a46e01be904fa8a8985e9850f0d angie Tue Nov 12 12:27:30 2019 -0800 dbSnp153: add a bigBed4 subtrack of coordinate ranges for mappings that we dropped due to inconsistent SPDI. refs #23283 Overall counts increased because we used to bail on an entire variant when we discovered an inconsistent SPDI, losing some valid mappings. Now we go through all mappings, and the bad ones are stored instead of dropped. diff --git src/hg/snp/dbSnpJsonToTab/dbSnpJsonToTab.c src/hg/snp/dbSnpJsonToTab/dbSnpJsonToTab.c index 9540b8e..e42d34c 100644 --- src/hg/snp/dbSnpJsonToTab/dbSnpJsonToTab.c +++ src/hg/snp/dbSnpJsonToTab/dbSnpJsonToTab.c @@ -21,34 +21,35 @@ #include "twoBit.h" void usage() /* Explain usage and exit. */ { errAbort( "dbSnpJsonToTab - Extract fields from dbSNP JSON files, write out as tab-separated track files\n" "usage:\n" " dbSnpJsonToTab assemblyList refSeqToUcsc refsnp-chrN.json[.gz] outRoot\n" "assemblyList is a comma-separated set of assembly_name values in JSON, e.g. 'GRCh37.p13,GRCh38.p12'\n" "refSeqToUcsc is tab-sep, 4 columns: RefSeq acc (NC_...), assembly_name, twoBitFile, ucscName (chr...)\n" "Track files starting with outRoot will be created:\n" " outRoot.<assembly>.bigDbSnp: BED4+ for main dbSNP track in each assembly in assemblyList\n" " (final columns need to be added later by bedJoinTabOffset)\n" " outRootDetails.tab: tab-separated allele frequency counts and other details (for bedJoinTabOffset)\n" + " outRoot.<assembly>.badCoords.bed: plain BED4 with range of inconsistent SPDI coords\n" "Additional files are written out for error reporting:\n" - " outRootFailed.json: lines of JSON input that encountered an error, for debugging/reprocessing\n" " outRootMerged.tab: two columns, obsolete dbSNP rs# ID and current rs# ID\n" - " outRootErrors.tab: lines of error info\n" + " outRootErrors.tab: lines of error info for rejected JSON objects\n" + " outRootWarnings.tab: warnings to report to dbSNP for rejected mappings or frequencies\n" "\n" "\n" "options:\n" " -freqSourceOrder=LIST LIST is an ordered, comma-sep list of project names (no spaces)\n" " for frequency data submitted to dbSNP, e.g. '1000genomes,GnomAD'\n" " If given, then the bigDbSnp minorAlleleFreq array uses the order.\n" " Otherwise minorAlleleFreq contains only one frequency, from the\n" " project with the highest sample count reported for that variant.\n" " -equivRegions=FILE Each line of FILE contains space-separated acc:start:end regions\n" " that are equivalent, for example PAR1 or PAR2 on X and Y,\n" " or an alt sequence and the corresponding portion of its chrom.\n" " These use RefSeq accesions not UCSC chr names, so that multiple\n" " assemblies can be described in the same file.\n" ); } @@ -60,66 +61,66 @@ struct hash *ncGrcToChrom = NULL; struct hash *ncToTwoBitChrom = NULL; struct hash *ncToSeqWin = NULL; /* Command line validation table. */ static struct optionSpec options[] = { {"freqSourceOrder", OPTION_STRING}, {"equivRegions", OPTION_STRING}, {NULL, 0}, }; struct outStreams /* Container for various output files that are not assembly-specific */ { FILE *details; // outRootDetails.tab - FILE *failJson; // outRootFailed.json FILE *merged; // outRootMerged.tab FILE *err; // outRootErrors.tab + FILE *warn; // outRootWarnings.tab }; static FILE *openOutFile(char *format, char *outRoot) /* Open a file for writing, named by format (which must have exactly one %s) and outRoot. */ { struct dyString *name = dyStringCreate(format, outRoot); FILE *outF = mustOpen(name->string, "w"); dyStringFree(&name); return outF; } static struct outStreams *outStreamsOpen(char *outRoot) /* Open output file handles */ { struct outStreams *outStreams; AllocVar(outStreams); outStreams->details = openOutFile("%sDetails.tab", outRoot); -outStreams->failJson = openOutFile("%sFailed.json", outRoot); outStreams->merged = openOutFile("%sMerged.tab", outRoot); outStreams->err = openOutFile("%sErrors.tab", outRoot); +outStreams->warn = openOutFile("%sWarnings.tab", outRoot); return outStreams; } static void outStreamsClose(struct outStreams **pOutStreams) /* Close file handles and free struct. */ { if (*pOutStreams != NULL) { struct outStreams *outStreams = *pOutStreams; carefulClose(&outStreams->details); - carefulClose(&outStreams->failJson); carefulClose(&outStreams->merged); carefulClose(&outStreams->err); + carefulClose(&outStreams->warn); freez(pOutStreams); } } struct spdiBed /* dbSNP 2.0 variant representation: {Sequence, Position, Deletion, Insertion}, expanded to * cover ambiguous placement region for indels when applicable. Add chromEnd to get bed3+. */ { struct spdiBed *next; char *chrom; // Sequence ID uint chromStart; // Position: 0-based start uint chromEnd; // Position + number of deleted bases char *del; // Deleted sequence (starting at pos) -- or number of deleted bases char *ins; // Inserted sequence }; @@ -178,79 +179,86 @@ static boolean seqIsNucleotide(const char *seq) /* Return TRUE unless seq contains a character that is not in the IUPAC nucleotide set. */ { if (seq == NULL) return FALSE; const char *p; for (p = seq; *p != '\0'; p++) { if (!isIupac(*p)) return FALSE; } return TRUE; } static struct spdiBed *parseSpdis(struct slRef *spdiList, char *name, struct lm *lm) -/* Make sure that all SPDI objects have the consistent seq_id, deleted_sequence and position. */ +/* Transform a list of SPDI JSON objects into a list of spdiBed. errAbort on invalid values. */ { struct spdiBed *spdiBList = NULL; -char *firstSeq = NULL; -int firstPos = -1; -char *firstDel = NULL; struct slRef *spdiRef; for (spdiRef = spdiList; spdiRef != NULL; spdiRef = spdiRef->next) { struct jsonElement *spdiEl = spdiRef->val; char *seq = jsonQueryString(spdiEl, "spdi", "seq_id", lm); int pos = jsonQueryInt(spdiEl, "spdi", "position", -1, lm); char *del = jsonQueryString(spdiEl, "spdi", "deleted_sequence", lm); char *ins = jsonQueryString(spdiEl, "spdi", "inserted_sequence", lm); if (seq == NULL) errAbort("Missing or null seq_id in spdi for %s", name); if (pos < 0) errAbort("Invalid or missing position in spdi for %s", name); if (del == NULL) errAbort("Missing or null deleted_sequence in spdi for %s", name); if (ins == NULL) errAbort("Missing or null inserted_sequence in spdi for %s", name); if (! seqIsNucleotide(del)) errAbort("Invalid deleted_sequence '%s' in spdi for %s", del, name); if (! seqIsNucleotide(ins)) errAbort("Invalid inserted_sequence '%s' in spdi for %s", ins, name); - if (firstSeq == NULL) + slAddHead(&spdiBList, spdiBedNewLm(seq, pos, del, ins, lm)); + } +slReverse(&spdiBList); +return spdiBList; +} + +static boolean spdiBedListSameRef(struct spdiBed *spdiBList, char *rsId) +/* Return TRUE if all spdiBed objects have the same chrom, chromStart and del. + * errAbort if diff chrom, return FALSE for diff start/del. */ +{ +if (spdiBList && spdiBList->next) + { + struct spdiBed *spdiB; + for (spdiB = spdiBList->next; spdiB != NULL; spdiB = spdiB->next) + { + if (differentString(spdiBList->chrom, spdiB->chrom)) + errAbort("Mismatching SPDI seq ('%s' vs '%s') for %s", + spdiBList->chrom, spdiB->chrom, rsId); + if (spdiBList->chromStart != spdiB->chromStart) { - // First spdi - firstSeq = seq; - firstDel = del; - firstPos = pos; + warn("Mismatching SPDI pos (%s: %d vs %d) for %s", + spdiB->chrom, spdiBList->chromStart, spdiB->chromStart, rsId); + return FALSE; } - else + if (differentString(spdiBList->del, spdiB->del)) { - if (differentString(firstSeq, seq)) - errAbort("Mismatching seq_id ('%s' vs '%s') in spdis for %s", - firstSeq, seq, name); - if (differentString(firstDel, del)) - errAbort("Mismatching deleted sequence ('%s' vs '%s') in %s spdis for %s", - firstDel, del, seq, name); - if (firstPos != pos) - errAbort("Mismatching position (%d vs %d) in %s spdis for %s", - firstPos, pos, seq, name); + warn("Mismatching SPDI del (%s:%d: '%s' vs '%s') for %s", + spdiB->chrom, spdiB->chromStart, spdiBList->del, spdiB->del, rsId); + return FALSE; } - slAddHead(&spdiBList, spdiBedNewLm(seq, pos, del, ins, lm)); } -slReverse(&spdiBList); -return spdiBList; + } +return TRUE; } 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; @@ -1250,30 +1258,35 @@ static struct bigDbSnp *parsePlacement(struct jsonElement *pl, struct sharedProps *props, boolean isMultiMapper, struct dyString *dy, struct lm *lm) /* Return a bigDbSnp with information about one placement (mapping). */ { struct bigDbSnp *bds; lmAllocVar(lm, bds); struct slRef *allSpdiRef = jsonQueryElement(pl, "placement", "alleles[*].allele.spdi", lm); bds->chrom = jsonQueryString(pl, "placement", "seq_id", lm); if (bds->chrom == NULL) errAbort("Missing or NULL seq_id for %s", props->name); struct spdiBed *spdiBList = parseSpdis(allSpdiRef, props->name, lm); if (!sameString(bds->chrom, spdiBList->chrom)) errAbort("seq_id mismatch: placement has '%s' but spdi has '%s for %s'", bds->chrom, spdiBList->chrom, props->name); +if (! spdiBedListSameRef(spdiBList, props->name)) + { + // Bad coords/alleles for this placement (usually due to mapping bug) -- don't proceed. + return NULL; + } bds->chromStart = spdiBList->chromStart; bds->chromEnd = spdiBList->chromEnd; bds->name = props->name; bds->ref = spdiBList->del; // spdiBList may or may not include no-change allele (del == ins) lmAllocArray(lm, bds->alts, slCount(spdiBList)); struct spdiBed *spdiB; for (spdiB = spdiBList; spdiB != NULL; spdiB = spdiB->next) { if (!sameString(bds->ref, spdiB->ins)) bds->alts[bds->altCount++] = spdiB->ins; } calcShiftBases(bds); boolean isRc = jsonQueryBoolean(pl, "placement", "placement_annot.is_aln_opposite_orientation", FALSE, lm); @@ -1602,38 +1615,30 @@ appendFrequencies(props, dy); dyStringPrintf(dy, "\t%d\t", slCount(props->soTerms)); dyStringPrintSlIntList(dy, props->soTerms, ","); dyStringPrintf(dy, "\t%d\t", slCount(props->clinVarAccs)); dyStringPrintSlNameList(dy, props->clinVarAccs, ","); dyStringAppendC(dy, '\t'); dyStringPrintSlNameList(dy, props->clinVarSigs, ","); dyStringPrintf(dy, "\t%d\t", slCount(props->submitters)); dyStringPrintSlNameList(dy, props->submitters, ","); dyStringPrintf(dy, "\t%d\t", slCount(props->pubMedIds)); dyStringPrintSlIntList(dy, props->pubMedIds, ","); dyStringAppendC(dy, '\n'); fputs(dy->string, outPropsF); } -static void dumpLineToFile(char *line, FILE *jsonF) -/* Write line to jsonF, making sure it ends in \n. */ -{ -fputs(line, jsonF); -if (line[strlen(line)-1] != '\n') - fputc('\n', jsonF); -} - static void updateSeqIsRc(struct slRef *placements, struct lm *lm) /* Note is_aln_opposite_orientation flag for each mapped sequence, regardless of whether it's * in the assembly that we're working on; later we can look up sequences from freq reports. */ { hashIntReset(seqIsRc); struct slRef *plRef; for (plRef = placements; plRef != NULL; plRef = plRef->next) { char *seqId = jsonQueryString(plRef->val, "placement", "seq_id", lm); boolean isRc = jsonQueryBoolean(plRef->val, "placement", "placement_annot.is_aln_opposite_orientation", FALSE, lm); struct hashEl *hel = hashLookup(seqIsRc, seqId); if (hel) { if (isRc != ptToInt(hel->val)) @@ -1642,30 +1647,31 @@ } else hashAddInt(seqIsRc, seqId, isRc); } } struct perAssembly // For each assembly whose mappings we're extracting from JSON, we need the name, // a jsonQuery path to extract just the mappings for that assembly, and an output file. { struct perAssembly *next; char *name; // assembly_name e.g. GRCh38.p12 char *placementPath; // jsonQuery path to placements on assembly char *plIsTopLevelPath; // jsonQuery path to find whether placement is top level in assembly FILE *outF; // bigDbSnp output file + FILE *outBad; // buggy coords BED4 output file }; static void perAssemblyFree(struct perAssembly **pPa) /* Close *pPa's outF and free mem. */ { if (pPa) { struct perAssembly *pa = *pPa; freeMem(pa->name); carefulClose(&pa->outF); freeMem(pa->placementPath); freeMem(pa->plIsTopLevelPath); freez(pPa); } } @@ -1695,30 +1701,50 @@ for (mergedId = mergedIds; mergedId != NULL; mergedId = mergedId->next) fprintf(outMergedF, "rs%s\trs%s\n", mergedId->name, rsId); } } static char *getChrom(char *refSeqAcc, char *assembly) /* Look up "refSeqAcc assembly" in ncGrcToChrom; if found, return chrom, otherwise return * refSeqAcc for translation later (e.g. NC_012920.1 liftOver to hg19 chrM). Don't free ret. */ { char ncGrc[strlen(refSeqAcc) + 1 + strlen(assembly) + 1]; safef(ncGrc, sizeof ncGrc, "%s %s", refSeqAcc, assembly); char *chrom = hashFindVal(ncGrcToChrom, ncGrc); return chrom ? chrom : refSeqAcc; } +static void writeBadCoords(struct jsonElement *pl, char *rsId, char *assembly, FILE *outBad, + struct lm *lm) +/* Write BED4 with range encompassing inconsistent SPDI coords. */ +{ +struct slRef *allSpdiRef = jsonQueryElement(pl, "placement", "alleles[*].allele.spdi", lm); +struct spdiBed *spdiBList = parseSpdis(allSpdiRef, rsId, lm); +char *chrom = getChrom(spdiBList->chrom, assembly); +uint minChromStart = BIGNUM; +uint maxChromEnd = 0; +struct spdiBed *spdiB; +for (spdiB = spdiBList; spdiB != NULL; spdiB = spdiB->next) + { + if (spdiB->chromStart < minChromStart) + minChromStart = spdiB->chromStart; + if (spdiB->chromEnd > maxChromEnd) + maxChromEnd = spdiB->chromEnd; + } +fprintf(outBad, "%s\t%u\t%u\t%s\n", chrom, minChromStart, maxChromEnd, rsId); +} + static void parseDbSnpJson(char *line, struct perAssembly *assemblyProps, struct outStreams *outStreams) /* Each line of dbSNP's file contains one JSON blob describing an rs# variant. Parse JSON, * extract the bits that we want to keep, print out BED+ and accompanying files. */ { struct lm *lm = lmInit(1<<16); struct jsonElement *top = jsonParseLm(line, lm); struct sharedProps *props = NULL; struct bigDbSnp *bds = NULL; struct slRef *allPlacements = jsonQueryElement(top, "top", "primary_snapshot_data.placements_with_allele[*]", lm); updateSeqIsRc(allPlacements, lm); struct dyString *dyScratch = dyStringNew(0); struct perAssembly *ap; @@ -1746,46 +1772,52 @@ //#*** but who knows, there could be more. Ultimately we will need to compare assembly //#*** sequences... unless/until dbSNP fixes their projection onto alts (JIRA VR-176). if (errCatchStart(errCatch)) { for (plRef = placements; plRef != NULL; plRef = plRef->next) { struct jsonElement *pl = plRef->val; if (props == NULL) { // Write props only if we have at least one placement. props = extractProps(top, lm); writeDetails(props, dyScratch, outStreams->details); writeMerged(top, outStreams->merged, lm); } bds = parsePlacement(pl, props, multiMapper, dyScratch, lm); + if (bds) + { bds->chrom = getChrom(bds->chrom, ap->name); bigDbSnpTabOut(bds, ap->outF); // Prevent stale bds in errCatch below: bds = NULL; } + else + { + writeBadCoords(pl, props->name, ap->name, ap->outBad, lm); + } + } } errCatchEnd(errCatch); if (errCatch->gotError) { - dumpLineToFile(line, outStreams->failJson); char *rsId = jsonQueryString(top, "top", "refsnp_id", lm); char *seqId = bds ? bds->chrom : "NOSEQ"; fprintf(outStreams->err, "%s\t%s\t%s", rsId, seqId, errCatch->message->string); } else if (isNotEmpty(errCatch->message->string)) - warn("%s", errCatch->message->string); + fprintf(outStreams->warn, "%s", errCatch->message->string); errCatchFree(&errCatch); } dyStringFree(&dyScratch); lmCleanup(&lm); } static struct slName *initFreqSourceOrder() /* If -freqSourceOrder option was given, extract source names from the comma-sep list * into globals freqSourceCount and freqSourceOrder[]. */ { char *freqSourceOrderStr = optionVal("freqSourceOrder", NULL); struct slName *sources = NULL; if (freqSourceOrderStr) { sources = slNameListFromComma(freqSourceOrderStr); @@ -1808,30 +1840,31 @@ { struct perAssembly *ap; AllocVar(ap); ap->name = cloneString(assembly->name); // Fancy path to get all placements that are on the desired assembly. struct dyString *dy = dyStringCreate("primary_snapshot_data.placements_with_allele" "[placement_annot.seq_id_traits_by_assembly[*].assembly_name=%s]", assembly->name); ap->placementPath = dyStringCannibalize(&dy); dy = dyStringCreate("placement_annot.seq_id_traits_by_assembly[assembly_name=%s].is_top_level", assembly->name); ap->plIsTopLevelPath = dyStringCannibalize(&dy); char prefix[2048]; safef(prefix, sizeof prefix, "%s.%s", outRoot, assembly->name); ap->outF = openOutFile("%s.bigDbSnp", prefix); + ap->outBad = openOutFile("%s.badCoords.bed", prefix); slAddHead(&assemblyProps, ap); } return assemblyProps; } static void parseRefSeqToUcsc(char *refSeqToUcsc, struct hash *ncGrcToChrom, struct hash *ncToTwoBitChrom) /* Extract mappings of {NC_*, GRC*} to chrom and NC_* to {twoBitFile, chrom}. */ { struct lineFile *lf = lineFileOpen(refSeqToUcsc, TRUE); char *line; while (lineFileNextReal(lf, &line)) { char *row[5]; int wordCount = chopTabs(line, row);