into buf. Return buf for convenience. */
{
safef(buf, bufSize, GENCODE_PREFIX"%s%s", suffix, version);
return buf;
}
boolean refGeneHasGencodeTags(struct slName *gencodeVersionList)
/* Return TRUE if this database has a wgEncodeGencodeRefSeq table as well as a ...Tag table. */
{
char *version = getLatestGencodeVersion(gencodeVersionList);
char table[PATH_LEN];
return hTableExists(database, gencodeTableName("RefSeq", version, table, sizeof(table)));
}
boolean startsWithGencodeGene(char *geneTrack)
/* Return TRUE if geneTrack starts with wgEncodeGencode{Basic,Comp,PseudoGene}.
* (There are other GENCODE genepred tracks that don't have tags associated with them.) */
{
return (startsWith(GENCODE_PREFIX"Basic", geneTrack) ||
startsWith(GENCODE_PREFIX"Comp", geneTrack) ||
startsWith(GENCODE_PREFIX"PseudoGene", geneTrack));
}
boolean isGencodeWithVersion(char *geneTrack, struct slName *versionList)
/* Return TRUE if geneTrack looks like a Gencode gene track for a supported version. */
{
if (! startsWithGencodeGene(geneTrack))
return FALSE;
struct slName *v;
for (v = versionList; v != NULL; v = v->next)
{
if (endsWith(geneTrack, v->name))
return TRUE;
}
return FALSE;
}
void selectTxStatus(boolean hasTxStatus, char *geneTrack)
/* Offer to include transcript status, e.g. whether it is in knownCanonical or has GENCODE tags.
* This makes one div per category of txStatus info; each div is visible only if its info is
* applicable to the selected gene track. If no divs are visible, display a message that
* there's nothing for the currently selected gene track. */
{
if (! hasTxStatus)
return;
startCollapsibleSection("txStatus", "Transcript status", FALSE);
boolean somethingIsVisible = FALSE;
if (hasGencodeTags())
{
struct slName *versionList = getGencodeTagVersions();
char *maybeKnownGene = knownGeneHasGencodeTags() ? "knownGene" : "";
char *maybeRefGene = refGeneHasGencodeTags(versionList) ? "refGene" : "";
char *maybeEnsGene = "";
char *versions = "";
if (versionList != NULL)
{
if (hTableExists(database, "ensGene"))
maybeEnsGene = "ensGene";
versions = slNameListToString(versionList, ' ');
}
boolean isVisible = (sameString(geneTrack, maybeKnownGene) ||
sameString(geneTrack, maybeEnsGene) ||
sameString(geneTrack, maybeRefGene) ||
isGencodeWithVersion(geneTrack, versionList));
somethingIsVisible |= isVisible;
printf("",
maybeKnownGene, maybeRefGene, maybeEnsGene, versions,
isVisible ? "block" : "none");
cartMakeCheckBox(cart, "hgva_txStatus_gencode", FALSE);
puts("Include the
GENCODE tags (if any) "
"associated with each transcript.
");
puts("
",
isVisible ? "block" : "none");
cartMakeCheckBox(cart, "hgva_txStatus_knownCanonical", FALSE);
char *desc = hTableExists(database, "knownToTag") ?
"based on
"
"APPRIS status or inclusion in "
"
GENCODE Basic subset: "
"principal > alternative > basic > longest isoform" :
"generally the longest isoform of a gene";
printf("Indicate whether each transcript is 'canonical' (%s).
\n", desc);
puts("
",
isVisible ? "block" : "none");
cartMakeCheckBox(cart, "hgva_txStatus_refSeqStatus", FALSE);
puts("Include the "
"
RefSeq status of each transcript.
");
puts("
",
somethingIsVisible ? "none" : "block");
puts("No transcript status data are available for the selected gene track.");
puts("
");
puts("regulatory_region_variant will be assigned. "
"Follow the links to description pages that explain how each dataset was "
"constructed. "
"Some datasets cover a significant portion of the genome and it may be "
"desirable to filter these annotations by cell type and/or score in order "
"to avoid an overabundance of hits. ");
// Use a table with checkboxes in one column and label/other stuff that depends on
// checkbox in other column.
puts("");
struct slRef *ref;
for (ref = trackRefList; ref != NULL; ref = ref->next)
{
struct trackDb *tdb = ref->val;
char cartVar[512];
safef(cartVar, sizeof(cartVar), "hgva_track_%s_%s", database, tdb->track);
puts("");
cartMakeCheckBox(cart, cartVar, FALSE);
puts(" ");
struct trackDb *topTdb = trackDbTopLevelSelfOrParent(tdb);
printf("%s
");
}
}
boolean isConsElTrack(struct trackDb *tdb, void *filterData)
/* This is a TdbFilterFunction to get "phastConsNwayElements" tracks. */
{
return (startsWith("phastCons", tdb->table) && stringIn("Elements", tdb->table));
}
boolean isConsScoreTrack(struct trackDb *tdb, void *filterData)
/* This is a TdbFilterFunction to get tracks that start with "phastCons" (no Elements)
* or "phyloP". */
{
return ((startsWith("phastCons", tdb->table) && !stringIn("Elements", tdb->table))
|| startsWith("phyloP", tdb->table));
}
void findCons(struct slRef **retElTrackRefList, struct slRef **retScoreTrackRefList)
/* See if this database has Conserved Elements and/or Conservation Scores */
{
tdbFilterGroupTrack(fullTrackList, fullGroupList, isConsElTrack,
NULL, NULL, retElTrackRefList);
tdbFilterGroupTrack(fullTrackList, fullGroupList, isConsScoreTrack,
NULL, NULL, retScoreTrackRefList);
}
boolean trackNameMatches(struct trackDb *tdb, void *filterData)
/* This is a TdbFilterFunction to get track(s) whose tdb->track matches name (filterData). */
{
char *name = filterData;
return sameString(tdb->track, name);
}
struct slRef *findTrackRefByName(char *name)
/* Return a reference to the named track, if found. */
{
struct slRef *trackRefList = NULL;
tdbFilterGroupTrack(fullTrackList, fullGroupList, trackNameMatches, name, NULL, &trackRefList);
return trackRefList;
}
void trackCheckBoxSection(char *sectionSuffix, char *title, struct slRef *trackRefList)
/* If trackRefList is non-NULL, make a collapsible section with a checkbox for each track,
* labelled with longLabel. */
{
if (trackRefList != NULL)
{
startCollapsibleSection(sectionSuffix, title, FALSE);
struct slRef *ref;
for (ref = trackRefList; ref != NULL; ref = ref->next)
{
struct trackDb *tdb = ref->val;
char cartVar[512];
safef(cartVar, sizeof(cartVar), "hgva_track_%s_%s", database, tdb->track);
cartMakeCheckBox(cart, cartVar, FALSE);
struct trackDb *topTdb = trackDbTopLevelSelfOrParent(tdb);
printf("%s ");
selectDbNsfp(dbNsfpTables);
selectTxStatus(hasTxStat, geneTrack);
selectDbSnp(gotSnp);
trackCheckBoxSection("Cosmic", "COSMIC", cosmicTrackRefList);
trackCheckBoxSection("ConsEl", "Conserved elements", elTrackRefList);
trackCheckBoxSection("ConsScore", "Conservation scores", scoreTrackRefList);
puts("
");
}
void selectFiltersFunc()
/* Options to restrict variants based on gene region/soTerm from gpFx */
{
startCollapsibleSection("filtersFunc", "Functional role", FALSE);
printf("Include variants annotated as");
cartMakeCheckBox(cart, "hgva_require_consEl", FALSE);
printf(" Include only the variants that overlap:");
if (slCount(elTrackRefList) == 1)
{
struct trackDb *tdb = elTrackRefList->val;
printf(" %s
");
}
else
{
puts("");
char *selected = cartUsualString(cart, "hgva_require_consEl_track", "");
struct slRef *ref;
for (ref = elTrackRefList; ref != NULL; ref = ref->next)
{
printf("");
struct trackDb *tdb = ref->val;
- cgiMakeOnClickRadioButton("hgva_require_consEl_track", tdb->track,
+ cgiMakeOnEventRadioButtonWithClass("hgva_require_consEl_track", tdb->track,
sameString(tdb->track, selected),
- "onclick=\"setCheckboxList('hgva_require_consEl', true);\"");
+ NULL, "click", "setCheckboxList('hgva_require_consEl', true);");
printf("%s ");
selectFiltersFunc();
selectFiltersKnownVar();
selectFiltersCons();
puts("
");
}
void selectOutput()
/* Just VEP text for now... should also have VEP HTML with limited # lines too.
* Next: custom track with same subset of info as we would stuff in VEP */
{
puts("output format: ");
char *selected = cartUsualString(cart, "hgva_outFormat", "vepTab");
printf("\n");
printOption("vepTab", selected, "Variant Effect Predictor (tab-separated text)");
printOption("vepHtml", selected, "Variant Effect Predictor (HTML)");
printf(" output file: ");
cgiMakeTextVar("hgva_outFile", fileName, 29);
printf(" (leave blank to keep output in browser)file type returned: ");
cgiMakeRadioButton("hgva_compressType", textOutCompressNone,
sameWord(textOutCompressNone, compressType));
printf(" plain text ");
cgiMakeRadioButton("hgva_compressType", textOutCompressGzip,
sameWord(textOutCompressGzip, compressType));
printf(" gzip compressed (ignored if output file is blank)");
puts("");
puts("This tool is for research use only. While this tool is open to the "
"public, users seeking information about a personal medical or genetic "
"condition are urged to consult with a qualified physician for "
"diagnosis and for answers to personal questions.");
puts("
\n");
printf("
, with spaces replaced by '_'s. */
{
int stringInLen = strlen(stringIn);
int stringOutLen = stringInLen + strlen("") + 1;
char *stringOut = needMem(stringOutLen);
safef(stringOut, stringOutLen, "", stringIn);
subChar(stringOut, ' ', '_');
return stringOut;
}
// dbSNP named alleles have many ways to describe a deletion from the reference,
// for example "LARGEDELETION", "LARGE DELETION", "... DELETED", "... DEL":
static const char *dbSnpDelRegex = "^\\(.*(DELET.*| DEL)\\)$";
static char **parseDbSnpAltAlleles(char *refAl, char *obsAls, boolean minusStrand,
int *retAltAlCount, boolean *retNeedLeftBase)
/* Given a non-symbolic reference allele and slash-sep observed alleles from dbSNP,
* return an array of +-strand alleles that are not the same as the reference.
* If any allele is "-" (deleted, zero-length), then set retNeedLeftBase to TRUE
* because in this case VCF requires that the reference base to the left of the indel
* must be added to all alleles, and the start coord also moves one base to the left.
* Also, if any alt allele is symbolic, padding is required.
* Note: this trashes obsAls. Resulting array can be freed but not its contents. */
{
int obsCount = countChars(obsAls, '/') + 1;
char *obsWords[obsCount];
chopByChar(obsAls, '/', obsWords, obsCount);
boolean obsHasDeletion = FALSE;
int i;
for (i = 0; i < obsCount; i++)
if (sameString(obsWords[i], "-"))
{
obsHasDeletion = TRUE;
break;
}
char **altAls;
AllocArray(altAls, obsCount);
int altCount = 0;
boolean needLeftBase = isEmpty(refAl) || sameString(refAl, "-");
for (i = 0; i < obsCount; i++)
{
char *altAl = obsWords[i];
int altAlLen = strlen(altAl);
if (minusStrand && isAllNt(altAl, altAlLen))
reverseComplement(altAl, altAlLen);
if (differentString(altAl, refAl))
{
if (sameString(altAl, "-"))
{
altAls[altCount] = "";
needLeftBase = TRUE;
}
else
{
// It would be nice to expand the "(CA)11/12/14/15/16/17/18/19/20" syntax of
// some dbSNP observed's. What are these?: "(D1S243)", "(D1S2870)"
// Unfortunately for observed="lengthTooLong" we just can't get the correct allele
// sequence. (76,130 of those in snp138)
// Hmmm, I guess we could at least stick in the right number of N's if we can
// parse "(245 BP INSERTION)". (2403 rows rlike "[0-9]+ BP ?INSERTION" in snp138)
if (!isAllNt(altAl, altAlLen))
{
// Symbolic allele: left base required, and enclose it in 's.
// But if it's one of dbSNP's LARGEDELETION kind of alleles, that is redundant
// with the reference allele, so if we know there is already a "-" allele,
// skip it.
if (obsHasDeletion && regexMatch(altAl, dbSnpDelRegex))
continue;
needLeftBase = TRUE;
altAl = encloseInAngleBracketsDbSnp(altAl);
}
altAls[altCount] = altAl;
}
altCount++;
}
}
*retAltAlCount = altCount;
*retNeedLeftBase = needLeftBase;
return altAls;
}
char *firstNCommaSep(struct slName *nameList, int n)
/* Return a comma-separated string with the first n names in nameList. */
{
struct dyString *dy = dyStringNew(0);
int i;
struct slName *el;
for (i=0, el=nameList; i < 5 && el != NULL; i++, el = el->next)
{
if (i > 0)
dyStringAppend(dy, ", ");
dyStringPrintf(dy, "'%s'", el->name);
}
return dyStringCannibalize(&dy);
}
//#*** Variant ID-matching should be metadata-driven too. termRegex -> data source.
static const char *rsIdRegex = "^rs[0-9]+$";
static void rsIdsToVcfRecords(struct annoAssembly *assembly, struct slName *rsIds,
struct vcfFile *vcff, struct vcfRecord **pRecList,
struct slName **pCommentList)
/* If possible, look up coordinates and alleles of dbSNP rs IDs. */
{
if (rsIds == NULL)
return;
char *table = hFindLatestSnpTable(database, NULL);
if (table == NULL)
return;
struct sqlConnection *conn = hAllocConn(assembly->name);
// Build a 'name in (...)' query, and build a hash of IDs so we can test whether all were found
struct dyString *dq = sqlDyStringCreate("select chrom, chromStart, chromEnd, name, strand, "
"refUCSC, observed from %s where name in (",
table);
struct hash *idHash = hashNew(0);
struct slName *id;
for (id = rsIds; id != NULL; id = id->next)
{
tolowers(id->name);
dyStringPrintf(dq, "%s'%s'", (id != rsIds ? "," : ""), id->name);
hashStoreName(idHash, id->name);
}
dyStringAppend(dq, ");");
struct sqlResult *sr = sqlGetResult(conn, dq->string);
// Construct a minimal VCF row to make a vcfRecord for each variant.
char *vcfRow[9];
vcfRow[5] = vcfRow[6] = vcfRow[7] = "."; // placeholder for qual, filter, info
// It would be cool to someday add snpNNN's exceptions column to the filter or info.
struct dyString *dyAltAlStr = dyStringNew(0);
int i;
char **row;
while ((row = sqlNextRow(sr)) != NULL)
{
char *chrom = row[0];
uint chromStart = atoll(row[1]);
uint chromEnd = atoll(row[2]);
char *name = row[3];
char *strand = row[4];
char refAl[max(strlen(row[5]), chromEnd-chromStart) + 2];
safecpy(refAl, sizeof(refAl), row[5]);
if (sameString(refAl, "-"))
refAl[0] = '\0';
else if (! isAllNt(refAl, strlen(refAl)))
{
// refUCSC is symbolic like "( 2548bp insertion )" -- just use the actual reference
// sequence for now, to keep things simple downstream.
//#*** Need something more efficient, like a sequence cache inside assembly!
struct dnaSeq *seq = twoBitReadSeqFragLower(assembly->tbf, chrom, chromStart, chromEnd);
safecpy(refAl, sizeof(refAl), seq->dna);
toUpperN(refAl, seq->size);
dnaSeqFree(&seq);
}
char *obsAls = row[6];
int altAlCount = 0;
boolean needLeftBase = FALSE;
char **altAls = parseDbSnpAltAlleles(refAl, obsAls, sameString(strand, "-"),
&altAlCount, &needLeftBase);
needLeftBase |= (chromStart == chromEnd); // should be redundant, but just in case.
hashRemove(idHash, name);
uint vcfStart = chromStart + 1;
dyStringClear(dyAltAlStr);
// If this is a non-symbolic indel, we'll have to move start one base to the left
// and add the reference base to the left of the actual alleles
if (needLeftBase)
{
vcfStart--;
//#*** Need something more efficient, like a sequence cache inside assembly!
struct dnaSeq *seq = twoBitReadSeqFragLower(assembly->tbf, chrom,
chromStart-1, chromStart);
toUpperN(seq->dna, 1);
char leftBase = seq->dna[0];
dnaSeqFree(&seq);
char refAlPlus[strlen(refAl)+2];
safef(refAlPlus, sizeof(refAlPlus), "%c%s", leftBase, refAl);
safecpy(refAl, sizeof(refAl), refAlPlus);
for (i = 0; i < altAlCount; i++)
{
if (i > 0)
dyStringAppendC(dyAltAlStr, ',');
if (isAllNt(altAls[i], strlen(altAls[i])))
dyStringPrintf(dyAltAlStr, "%c%s", leftBase, altAls[i]);
else
dyStringAppend(dyAltAlStr, altAls[i]);
}
}
else
{
// Not an indel, just make comma-sep string of alt alleles.
for (i = 0; i < altAlCount; i++)
{
if (i > 0)
dyStringAppendC(dyAltAlStr, ',');
dyStringAppend(dyAltAlStr, altAls[i]);
}
}
char vcfStartStr[64];
safef(vcfStartStr, sizeof(vcfStartStr), "%d", vcfStart);
vcfRow[0] = chrom;
vcfRow[1] = vcfStartStr;
vcfRow[2] = name;
vcfRow[3] = refAl;
vcfRow[4] = dyAltAlStr->string;
struct vcfRecord *rec = vcfRecordFromRow(vcff, vcfRow);
slAddHead(pRecList, rec);
}
dyStringFree(&dyAltAlStr);
// Check for IDs not found
struct slName *notFoundIds = hashListNames(idHash);
if (notFoundIds != NULL)
{
char *namesNotFound = firstNCommaSep(notFoundIds, 5);
struct dyString *dy = dyStringCreate("%d rs# IDs not found, e.g. %s",
slCount(notFoundIds), namesNotFound);
slAddTail(pCommentList, slNameNew(dy->string));
freeMem(namesNotFound);
dyStringFree(&dy);
}
slNameFreeList(¬FoundIds);
}
static struct vcfRecord *parseVariantIds(struct annoAssembly *assembly, char *variantIdText,
struct slName **pCommentList)
/* Return a sorted list of minimal vcfRecords (coords, name, ref and alt alleles)
* corresponding to each pasted variant. */
{
struct vcfRecord *recList = NULL;
char *p = cloneString(variantIdText), *id;
struct slName *rsIds = NULL, *unknownIds = NULL;
subChar(p, ',', ' ');
while ((id = nextWord(&p)) != NULL)
{
if (regexMatchNoCase(id, rsIdRegex))
slNameAddHead(&rsIds, id);
else
slNameAddHead(&unknownIds, id);
}
if (unknownIds != NULL)
{
slReverse(&unknownIds);
char *firstUnknownIds = firstNCommaSep(unknownIds, 5);
struct dyString *dy = dyStringCreate("%d variant identifiers are unrecognized, e.g. %s",
slCount(unknownIds), firstUnknownIds);
slAddTail(pCommentList, slNameNew(dy->string));
freeMem(firstUnknownIds);
dyStringFree(&dy);
}
struct vcfFile *vcff = vcfFileNew();
rsIdsToVcfRecords(assembly, rsIds, vcff, &recList, pCommentList);
slSort(&recList, vcfRecordCmp);
slNameFreeList(&rsIds);
slNameFreeList(&unknownIds);
return recList;
}
static void adjustRangeForVariants(struct vcfRecord *recList, char **pChrom,
uint *pStart, uint *pEnd)
/* Given a *sorted* list of VCF records, look at the first and last records and if
* they fall outside of the range, expand the range to include them. */
{
if (pChrom != NULL && *pChrom != NULL && recList != NULL)
{
// We have a non-empty sorted list of records, and search is not genome-wide;
// look at first and last variants to see if they're not already included in search range.
struct vcfRecord *first = recList, *last = recList;
while (last->next)
last = last->next;
if (differentString(*pChrom, first->chrom) || differentString(*pChrom, last->chrom))
{
// Variants span multiple chroms; force genome-wide search.
*pChrom = NULL;
*pStart = *pEnd = 0;
}
else
{
// Variant(s) are on the same chrom as search range; check starts and ends.
if (*pStart > first->chromEnd)
*pStart = first->chromEnd;
if (*pEnd < last->chromStart)
*pEnd = last->chromStart;
}
}
}
static struct annoStreamer *makeVariantIdStreamer(struct annoAssembly *assembly, int maxOutRows,
char **pChrom, uint *pStart, uint *pEnd,
struct slName **pCommentList)
/* Translate user's pasted/uploaded variant IDs into minimal VCF if possible.
* Return a VCF streamer for those, and if the current search position is too narrow
* to include all of the variants, widen it as necessary. */
{
// Hash variant text to get trash filename. Use if exists, otherwise build it.
char *variantIds = cartString(cart, hgvaVariantIds);
char *varFile = variantIdPath(assembly, variantIds);
boolean forceRebuild = cartUsualBoolean(cart, "hgva_rebuildVariantIds", FALSE);
if (! fileExists(varFile) || forceRebuild)
{
struct vcfRecord *varList = parseVariantIds(assembly, variantIds, pCommentList);
//#*** If no variants were recognized, we should probably show main page with a warning.
adjustRangeForVariants(varList, pChrom, pStart, pEnd);
FILE *f = mustOpen(varFile, "w");
writeMinimalVcfHeader(f, assembly->name);
struct vcfRecord *var;
for (var = varList; var != NULL; var = var->next)
writeMinimalVcfRow(f, var);
carefulClose(&f);
}
return annoStreamVcfNew(varFile, NULL, FALSE, assembly, maxOutRows);
}
static struct trackDb *getVariantTrackDb(char *variantTrack)
/* Get a tdb for the user's selected variant track, or warn and return NULL
* (main page will be displayed) */
{
struct trackDb *varTdb = tdbForTrack(database, variantTrack, &fullTrackList);
if (varTdb == NULL)
{
if (isHubTrack(variantTrack))
warn("Can't find hub track '%s'", variantTrack);
else
warn("Can't find tdb for variant track '%s'", variantTrack);
}
else
checkVariantTrack(varTdb);
return varTdb;
}
static char *gencodeVersionFromTrack(char *track)
/* If track is a GENCODE table, find and return a pointer to the version at the end;
* otherwise return NULL. */
{
if (startsWithGencodeGene(track))
{
char *v = strrchr(track, 'V');
return v;
}
return NULL;
}
static char *gencodeTagTableForTrack(char *db, char *track)
/* If there is a wgEncodeGencodeTag table that can be associated with track,
* return it; otherwise return NULL. */
{
struct slName *versionList = getGencodeTagVersions();
if (startsWithGencodeGene(track))
{
char *version = gencodeVersionFromTrack(track);
if (version != NULL)
{
char table[PATH_LEN];
return cloneString(gencodeTableName("Tag", version, table, sizeof(table)));
}
}
else if (sameString(track, "refGene") && refGeneHasGencodeTags(versionList))
{
char *version = getLatestGencodeVersion(versionList);
char table[PATH_LEN];
if (hTableExists(db, gencodeTableName("RefSeq", version, table, sizeof(table))))
return cloneString(gencodeTableName("Tag", version, table, sizeof(table)));
}
else if (sameString(track, "knownGene") && knownGeneHasGencodeTags())
{
if (hTableExists(db, "knownToTag"))
return cloneString("knownToTag");
}
return NULL;
}
static struct joinerDtf *getTxStatusExtras(char *db, char *track)
// Several kinds of transcript status may be enabled in the cart; if any are enabled,
// and if they apply to track, return the tables & fields to be joined with the track.
{
struct joinerDtf *txStatusExtras = NULL;
if (cartUsualBoolean(cart, "hgva_txStatus_gencode", FALSE))
{
char *gencodeTagTable = gencodeTagTableForTrack(db, track);
if (gencodeTagTable != NULL)
{
char *field = "tag";
if (sameString("knownToTag", gencodeTagTable))
field = "value";
slAddHead(&txStatusExtras, joinerDtfNew(db, gencodeTagTable, field));
}
}
if (cartUsualBoolean(cart, "hgva_txStatus_knownCanonical", FALSE) &&
sameString(track, "knownGene") &&
hTableExists(db, "knownCanonical"))
{
slAddHead(&txStatusExtras, joinerDtfNew(db, "knownCanonical", "transcript"));
}
if (cartUsualBoolean(cart, "hgva_txStatus_refSeqStatus", FALSE) &&
sameString(track, "refGene") &&
genbankTableExists(db, refSeqStatusTable))
{
char *genbankDb=NULL, *tableName=NULL;
genbankGetDbTable(db, refSeqStatusTable, &genbankDb, &tableName);
slAddHead(&txStatusExtras, joinerDtfNew(genbankDb, tableName, "status"));
}
return txStatusExtras;
}
static void configAddDtf(struct dyString *dy, struct joinerDtf *dtf, boolean *pIsFirst)
/* Add a JSON object with [db].table and (list of one) field. */
{
if (! *pIsFirst)
dyStringAppend(dy, ", ");
// The convention for related tables in same db as track is to prepend "." instead of ".":
char *db = (sameString(dtf->database, database) ? "" : dtf->database);
dyStringPrintf(dy, "{ \"table\": \"%s.%s\", \"fields\": [\"%s\"] }", db, dtf->table, dtf->field);
*pIsFirst = FALSE;
}
static void configAddTableField(struct dyString *dy, char *table, char *field, boolean *pIsFirst)
/* Add a JSON object with current database, table and (list of one) field. */
{
struct joinerDtf dtf;
dtf.database = database;
dtf.table = table;
dtf.field = field;
configAddDtf(dy, &dtf, pIsFirst);
}
static struct jsonElement *configForStreamer(char *db, struct trackDb *tdb,
struct joinerDtf *txStatusExtras)
/* Add VAI-specific config options, if applicable. */
{
struct jsonElement *config = NULL;
char *track = tdb->track;
struct dyString *dyConfig = dyStringCreate("{ \"naForMissing\": false,"
" \"relatedTables\": [ ");
boolean isFirst = TRUE;
// If track is sql-based knownGene and we have kgXref, then add kgXref.geneSymbol after
// the columns of knownGene.
if (sameString(track, "knownGene") &&
!isCustomTrack(track) && !isHubTrack(track) &&
!trackDbSetting(tdb, "bigDataUrl") &&
hTableExists(db, "kgXref"))
{
configAddTableField(dyConfig, "kgXref", "geneSymbol", &isFirst);
}
struct joinerDtf *txStatDtf;
for (txStatDtf = txStatusExtras; txStatDtf != NULL; txStatDtf = txStatDtf->next)
configAddDtf(dyConfig, txStatDtf, &isFirst);
// If any of the above apply, close the relatedTables list and config object
// and parse into jsonElements.
if (! isFirst)
{
dyStringAppend(dyConfig, " ] }");
config = jsonParse(dyConfig->string);
dyStringFree(&dyConfig);
}
return config;
}
static void adjustGpVarOverlapRule(struct annoGrator *gpVarGrator, boolean haveRegulatory)
/* If we're able to detect regulatory elements, and want to keep those annotations, loosen up
* gpVarGrator's overlap rule from the default (must overlap). */
{
if (haveRegulatory && cartUsualBoolean(cart, "hgva_include_regulatory", TRUE))
gpVarGrator->setOverlapRule(gpVarGrator, agoNoConstraint);
}
static void addTxStatusExtras(struct annoFormatter *vepOut, char *geneTrack,
struct annoGrator *gpVarGrator,
struct joinerDtf *txStatusExtras)
/* Given a list of tables and fields that will be joined with geneTrack to provide transcript
* status info, configure vepOut to put them in the EXTRAs column. */
{
struct joinerDtf *txStatDtf;
for (txStatDtf = txStatusExtras; txStatDtf != NULL; txStatDtf = txStatDtf->next)
{
char *tag = NULL, *description = NULL;
boolean isBoolean = FALSE;
// Double-check that we're not joining in some table from a different assembly database:
if (differentString(txStatDtf->database, database) &&
differentString(txStatDtf->database, "hgFixed"))
errAbort("addTxStatusExtras: Expected db=%s or hgFixed in txStatDtf but got %s",
database, txStatDtf->database);
if ((startsWith(GENCODE_PREFIX"Tag", txStatDtf->table) &&
sameString(txStatDtf->field, "tag")) ||
(sameString(txStatDtf->table, "knownToTag") &&
sameString(txStatDtf->field, "value")))
{
tag = "GENCODE_TAG";
description = "GENCODE tags for the transcript";
}
else if (sameString(txStatDtf->table, "knownCanonical") &&
sameString(txStatDtf->field, "transcript"))
{
tag = "CANONICAL";
description = "If present, the transcript is the 'canonical' transcript of the gene "
"(generally the longest isoform of the gene)";
isBoolean = TRUE;
}
else if (sameString(txStatDtf->table, "refSeqStatus") &&
sameString(txStatDtf->field, "status"))
{
tag = "REFSEQ_STATUS";
description = "RefSeq status of the transcript";
}
else
{
errAbort("addTxStatusExtras: Unrecognized {table,field}: {%s,%s}",
txStatDtf->table, txStatDtf->field);
}
char *column = annoStreamDbColumnNameFromDtf(database, geneTrack, txStatDtf);
annoFormatVepAddExtraItem(vepOut, (struct annoStreamer *)gpVarGrator,
tag, description, column, isBoolean);
}
}
void doQuery()
/* Translate simple form inputs into anno* components and execute query. */
{
dyInfo = dyStringNew(0);
char *chrom = NULL;
uint start = 0, end = 0;
if (sameString(regionType, hgvaRegionTypeRange))
getCartPosOrDie(&chrom, &start, &end);
struct annoAssembly *assembly = hAnnoGetAssembly(database);
boolean isCommandLine = (cgiOptionalString("cgiSpoof") != NULL);
char *geneTrack = cartString(cart, "hgva_geneTrack");
struct trackDb *geneTdb = tdbForTrack(database, geneTrack, &fullTrackList);
if (geneTdb == NULL)
{
warn("Can't find tdb for gene track %s", geneTrack);
if (! isCommandLine)
doUi();
return;
}
int maxVarRows = cartUsualInt(cart, "hgva_variantLimit", 10);
struct annoStreamer *primary = NULL;
char *primaryLongLabel = NULL;
char *variantTrack = cartString(cart, "hgva_variantTrack");
struct slName *commentList = NULL;
if (sameString(variantTrack, hgvaSampleVariants))
{
primary = makeSampleVariantsStreamer(assembly, geneTdb, maxVarRows);
primaryLongLabel = hgvaSampleVariantsLabel;
// Sample variants can't always be made within the currently selected position range,
// so just for these, force search to be genome-wide.
chrom = NULL;
start = 0;
end = 0;
}
else if (sameString(variantTrack, hgvaUseVariantIds))
{
// Override search position if cart position doesn't include all variants:
primary = makeVariantIdStreamer(assembly, maxVarRows, &chrom, &start, &end, &commentList);
primaryLongLabel = hgvaVariantIdsLabel;
}
else if (sameString(variantTrack, hgvaUseVariantFileOrUrl))
{
char *fileOrUrl = cartString(cart, hgvaVariantFileOrUrl);
char *type = cartOptionalString(cart, hgvaVariantFileOrUrlType);
primary = hAnnoStreamerFromBigFileUrl(fileOrUrl, NULL, assembly, maxVarRows, type);
primaryLongLabel = hgvaVariantFileOrUrlLabel;
}
else
{
struct trackDb *varTdb = getVariantTrackDb(variantTrack);
if (varTdb == NULL)
{
if (! isCommandLine)
doUi();
return;
}
primary = hAnnoStreamerFromTrackDb(assembly, varTdb->table, varTdb, chrom, maxVarRows, NULL);
primaryLongLabel = varTdb->longLabel;
}
enum annoGratorOverlap geneOverlapRule = agoMustOverlap;
struct joinerDtf *txStatusExtras = getTxStatusExtras(database, geneTrack);
struct jsonElement *gpConfig = configForStreamer(database, geneTdb, txStatusExtras);
struct annoGrator *gpVarGrator = hAnnoGratorFromTrackDb(assembly, geneTdb->table, geneTdb, chrom,
ANNO_NO_LIMIT, primary->asObj,
geneOverlapRule, gpConfig);
setGpVarFuncFilter(gpVarGrator);
// Some grators may be used as both filters and output values. To avoid making
// multiple grators for the same source, hash them by trackName:
struct hash *gratorsByName = hashNew(8);
struct annoGrator *snpGrator = NULL;
char *snpDesc = NULL;
if (cartUsualBoolean(cart, "hgva_rsId", TRUE))
snpGrator = gratorForSnpBed4(gratorsByName, "", assembly, chrom, agoNoConstraint, &snpDesc);
// Now construct gratorList in the order in which annoFormatVep wants to see them,
// i.e. first the gpVar, then the snpNNN, then whatever else:
struct annoGrator *gratorList = NULL;
slAddHead(&gratorList, gpVarGrator);
if (snpGrator != NULL)
slAddHead(&gratorList, snpGrator);
// Text or HTML output?
char *outFormat = cartUsualString(cart, "hgva_outFormat", "vepTab");
boolean doHtml = sameString(outFormat, "vepHtml");
// Initialize VEP formatter:
struct annoFormatter *vepOut = annoFormatVepNew("stdout", doHtml,
primary, primaryLongLabel,
(struct annoStreamer *)gpVarGrator,
geneTdb->longLabel,
(struct annoStreamer *)snpGrator,
snpDesc, assembly);
addTxStatusExtras(vepOut, geneTrack, gpVarGrator, txStatusExtras);
boolean haveRegulatory = FALSE;
addOutputTracks(&gratorList, gratorsByName, vepOut, assembly, chrom, doHtml, &haveRegulatory);
adjustGpVarOverlapRule(gpVarGrator, haveRegulatory);
addFilterTracks(&gratorList, gratorsByName, assembly, chrom);
slReverse(&gratorList);
if (doHtml)
{
webStart(cart, database, "Annotated Variants in VEP/HTML format");
}
else if (! isCommandLine)
{
// Undo the htmlPushEarlyHandlers() because after this point they make ugly text:
popWarnHandler();
popAbortHandler();
textOpen();
webStartText();
}
if (isNotEmpty(dyInfo->string))
puts(dyInfo->string);
struct annoGratorQuery *query = annoGratorQueryNew(assembly, primary, gratorList, vepOut);
struct slName *comment;
for (comment = commentList; comment != NULL; comment = comment->next)
vepOut->comment(vepOut, comment->name);
if (chrom != NULL)
annoGratorQuerySetRegion(query, chrom, start, end);
annoGratorQueryExecute(query);
annoGratorQueryFree(&query);
if (doHtml)
webEnd();
else if (! isCommandLine)
textOutClose(&compressPipeline, NULL);
}
int main(int argc, char *argv[])
/* Process command line. */
{
long enteredMainTime = clock1000();
if (hIsPrivateHost())
pushCarefulMemHandler(LIMIT_2or6GB);
cgiSpoof(&argc, argv);
boolean isCommandLine = (cgiOptionalString("cgiSpoof") != NULL);
if (!isCommandLine)
htmlPushEarlyHandlers(); /* Make errors legible during initialization. */
oldVars = hashNew(10);
boolean startQuery = (cgiUsualString("hgva_startQuery", NULL) != NULL);
if (startQuery)
{
if (isCommandLine)
// No HTTP header for command-line use.
cart = cartForSession(hUserCookie(), excludeVars, oldVars);
else
cart = cartAndCookieNoContent(hUserCookie(), excludeVars, oldVars);
}
else
cart = cartAndCookie(hUserCookie(), excludeVars, oldVars);
// Try to deal with virt chrom position used by hgTracks.
if (startsWith("virt:", cartUsualString(cart, "position", "")))
cartSetString(cart, "position", cartUsualString(cart, "nonVirtPosition", ""));
/* Set up global variables. */
getDbAndGenome(cart, &database, &genome, oldVars);
initGenbankTableNames(database);
regionType = cartUsualString(cart, hgvaRegionType, hgvaRegionTypeGenome);
if (isEmpty(cartOptionalString(cart, hgvaRange)))
cartSetString(cart, hgvaRange, hDefaultPos(database));
int timeout = cartUsualInt(cart, "udcTimeout", 300);
if (udcCacheTimeout() < timeout)
udcSetCacheTimeout(timeout);
knetUdcInstall();
cartTrackDbInit(cart, &fullTrackList, &fullGroupList, TRUE);
if (lookupPosition(cart, hgvaRange))
{
if (startQuery)
doQuery();
else if (! isCommandLine)
doUi();
}
else
{
// Revert to lastPosition if we have multiple matches or warnings,
// especially in case user manually edits browser location as in #13009:
char *position = cartUsualString(cart, "lastPosition", hDefaultPos(database));
cartSetString(cart, hgvaRange, position);
if (webGotWarnings())
{
// We land here when lookupPosition pops up a warning box.
// Reset the problematic position and show the main page.
doMainPage();
}
// If lookupPosition returned FALSE and didn't report warnings,
// then it wrote HTML showing multiple position matches & links.
}
cartCheckout(&cart);
if (! isCommandLine)
cgiExitTime("hgVai", enteredMainTime);
return 0;
}