4e04099939e8c5351d07cd086a3a5140cb678532 kate Mon Jul 17 16:36:35 2017 -0700 1. Print MAF (instead of allele freqs). Code cleanup. refs #15646 diff --git src/hg/hgc/gtexEqtlClusterClick.c src/hg/hgc/gtexEqtlClusterClick.c index 6933274..5d2fff2 100644 --- src/hg/hgc/gtexEqtlClusterClick.c +++ src/hg/hgc/gtexEqtlClusterClick.c @@ -30,120 +30,153 @@ slReverse(&eqtls); sqlFreeResult(&sr); hFreeConn(&conn); return eqtls; } static char *getGeneDescription(struct sqlConnection *conn, char *geneName) /* Return description from KnownGenes track */ { char query[256]; sqlSafef(query, sizeof query, "SELECT kgXref.description FROM kgXref WHERE geneSymbol='%s'", geneName); return sqlQuickString(conn, query); } -void doGtexEqtlDetails(struct trackDb *tdb, char *item) -/* Details of GTEx eQTL item */ -{ -char *chrom = cartString(cart, "c"); -int start = cartInt(cart, "o"); -int end = cartInt(cart, "t"); -struct gtexEqtlCluster *eqtl = getGtexEqtl(item, chrom, start, end, tdb->table); -genericHeader(tdb, item); -char *version = gtexVersion(tdb->table); -struct gtexTissue *tissues = gtexGetTissues(version); -struct hash *tissueHash = hashNew(0); -struct gtexTissue *tis = NULL; -for (tis = tissues; tis != NULL; tis = tis->next) - hashAdd(tissueHash, tis->name, tis); -int i; - -struct sqlConnection *conn = hAllocConn(database); -char *geneName = eqtl->target; -char *desc = getGeneDescription(conn, geneName); -printf("<b>Gene: </b>"); -if (desc == NULL) - printf("%s\n", geneName); -else +static void printMinorAlleleFreq(char *rsId, struct sqlConnection *conn) +/* Print minor allele frequency for a SNP (from UCSC dbSNP table) */ { - printf("<a target='_blank' href='%s?db=%s&hgg_gene=%s'>%s</a><br>\n", - hgGeneName(), database, geneName, geneName); - printf("<b>Description:</b> %s\n", desc); +#define ALLELE_COUNT 10 +char query[256]; +sqlSafef(query, sizeof query, "SELECT alleleFreqs FROM snp147 WHERE name='%s'", rsId); +double freqs[ALLELE_COUNT]; +int count = sqlDoubleArray(sqlQuickString(conn, query), freqs, ALLELE_COUNT); +doubleSort(count, freqs); +printf("<br><b>Minor allele frequency (1000 Genomes):</b> %.0f%%\n", 100.0 * freqs[count-2]); } -char posLink[1024]; -safef(posLink, sizeof posLink,"<a href='%s&db=%s&position=%s%%3A%d-%d'>%s:%d-%d</a>", - hgTracksPathAndSettings(), database, - eqtl->chrom, eqtl->chromStart+1, eqtl->chromEnd, - eqtl->chrom, eqtl->chromStart+1, eqtl->chromEnd); -// TODO: Consider adding Ensembl gene ID, GENCODE biotype and class (as in gtexGene track) -char query[256]; -printf("<br><b>Variant: </b>%s ", eqtl->name); -if (startsWith("rs", eqtl->name)) +static void printGwasCatalogTrait(char *rsId, struct sqlConnection *conn) +/* Print trait/disease for a SNP (from UCSC gwasCatalog table) */ { - printDbSnpRsUrl(eqtl->name, "dbSNP"); - sqlSafef(query, sizeof query, "SELECT alleleFreqs FROM snp147 WHERE name='%s'", eqtl->name); - char *freqs = sqlQuickString(conn, query); - printf("<br><b>Allele frequencies:</b> %s\n", freqs); - sqlSafef(query, sizeof query, "SELECT count(*) FROM gwasCatalog WHERE name='%s'", eqtl->name); +char query[256]; +sqlSafef(query, sizeof query, "SELECT count(*) FROM gwasCatalog WHERE name='%s'", rsId); int count = sqlQuickNum(conn, query); if (count) { - sqlSafef(query, sizeof query, "SELECT trait FROM gwasCatalog WHERE name='%s' LIMIT 1", - eqtl->name); + sqlSafef(query, sizeof query, "SELECT trait FROM gwasCatalog WHERE name='%s' LIMIT 1", rsId); char *trait = sqlQuickString(conn, query); printf("<br><b>GWAS disease or trait"); if (count > 1) printf(" (1 of %d)", count); printf(": </b>%s <a target='_blank' href='https://www.ebi.ac.uk/gwas/search?query=%s'>" - "GWAS Catalog</a>\n", trait, eqtl->name); + "GWAS Catalog</a>\n", trait, rsId); } } -else - printf("%s\n", eqtl->name); - -printf("<br><b>Position:</b> %s\n", posLink); -printf("<br><b>Score:</b> %d\n", eqtl->score); +static void printEqtlRegion(struct gtexEqtlCluster *eqtl, char *table, struct sqlConnection *conn) +/* Print position of region encompassing all identified eQTL's for this gene */ +{ #define FLANK 1000 +char query[256]; sqlSafef(query, sizeof query, "SELECT MIN(chromStart) from %s WHERE target='%s'", - tdb->table, eqtl->target); -start = sqlQuickNum(conn, query) - FLANK; + table, eqtl->target); +int start = sqlQuickNum(conn, query) - FLANK; sqlSafef(query, sizeof query, "SELECT MAX(chromEnd) from %s WHERE target='%s'", - tdb->table, eqtl->target); -end = sqlQuickNum(conn, query) + FLANK; + table, eqtl->target); +int end = sqlQuickNum(conn, query) + FLANK; +char posLink[1024]; safef(posLink, sizeof posLink,"<a href='%s&db=%s&position=%s%%3A%d-%d'>%s:%d-%d</a>", hgTracksPathAndSettings(), database, eqtl->chrom, start+1, end, eqtl->chrom, start+1, end); printf("<br><b>Region containing eQTLs for this gene:</b> %s (%d bp, including +-%dbp flank)\n", posLink, end-start, FLANK); -printf("<br><a target='_blank' href='https://www.gtexportal.org/home/bubbleHeatmapPage/%s'>" - "View eQTLs for this gene at the GTEx portal<a>\n", - geneName); - -printf("<br><b>Number of tissues with this eQTL:</b> %d\n", eqtl->expCount); -hFreeConn(&conn); +} +static void printClusterDetails(struct gtexEqtlCluster *eqtl, char *table) +/* Print details of an eQTL cluster */ +{ webNewSection("eQTL Cluster Details"); +char *version = gtexVersion(table); +struct gtexTissue *tissues = gtexGetTissues(version); +struct hash *tissueHash = hashNew(0); +struct gtexTissue *tis = NULL; +for (tis = tissues; tis != NULL; tis = tis->next) + hashAdd(tissueHash, tis->name, tis); printf("<table id='eqtls' cellspacing=1 cellpadding=3>\n"); printf("<style>#eqtls th {text-align: left; background-color: #F3E0BE;}</style>"); printf("<tr><th> </th><th>Tissue</th><th>Effect </th><th>Probability </th></tr>\n"); +int i; for (i=0; i<eqtl->expCount; i++) { double effect= eqtl->expScores[i]; double prob = eqtl->expProbs[i]; struct gtexTissue *tis = (struct gtexTissue *)hashFindVal(tissueHash, eqtl->expNames[i]); unsigned color = tis ? tis->color : 0; // BLACK char *name = tis ? tis->description : "Unknown"; printf("<tr><td bgcolor=#%06X></td><td>%s</td><td>%s%0.2f</td><td>%0.2f</td></tr>\n", color, name, effect < 0 ? "" : "+", effect, prob); } printf("</table>"); webEndSection(); +} + +void doGtexEqtlDetails(struct trackDb *tdb, char *item) +/* Details of GTEx eQTL item */ +{ +char *chrom = cartString(cart, "c"); +int start = cartInt(cart, "o"); +int end = cartInt(cart, "t"); +struct gtexEqtlCluster *eqtl = getGtexEqtl(item, chrom, start, end, tdb->table); +char *geneName = eqtl->target; + + +genericHeader(tdb, item); +printf("<b>Gene: </b>"); +struct sqlConnection *conn = hAllocConn(database); +char *desc = getGeneDescription(conn, geneName); +if (desc == NULL) + printf("%s\n", geneName); +else + { + printf("<a target='_blank' href='%s?db=%s&hgg_gene=%s'>%s</a><br>\n", + hgGeneName(), database, geneName, geneName); + printf("<b>Description:</b> %s\n", desc); + } + +// TODO: Consider adding Ensembl gene ID, GENCODE biotype and class (as in gtexGene track) +printf("<br><b>Variant: </b>%s ", eqtl->name); +if (startsWith("rs", eqtl->name)) + { + printDbSnpRsUrl(eqtl->name, "dbSNP"); + printMinorAlleleFreq(eqtl->name, conn); + printGwasCatalogTrait(eqtl->name, conn); + } +else + printf("%s\n", eqtl->name); + +char posLink[1024]; +safef(posLink, sizeof posLink,"<a href='%s&db=%s&position=%s%%3A%d-%d'>%s:%d-%d</a>", + hgTracksPathAndSettings(), database, + eqtl->chrom, eqtl->chromStart+1, eqtl->chromEnd, + eqtl->chrom, eqtl->chromStart+1, eqtl->chromEnd); +printf("<br><b>Position:</b> %s\n", posLink); + +printf("<br><b>Score:</b> %d\n", eqtl->score); + +printEqtlRegion(eqtl, tdb->table, conn); + +// print link to GTEx portal +printf("<br><a target='_blank' href='https://www.gtexportal.org/home/bubbleHeatmapPage/%s'>" + "View eQTLs for this gene at the GTEx portal<a>\n", + geneName); +printf("<br><b>Number of tissues with this eQTL:</b> %d\n", eqtl->expCount); +hFreeConn(&conn); + +printClusterDetails(eqtl, tdb->table); + webNewEmptySection(); printTrackHtml(tdb); }