d09b392333a2f5d2f943f2f36dad4be0a477381e kate Thu Sep 14 13:55:37 2017 -0700 Add pvalue and TSS distance for each eQTL. Input from Casey Brown, U Penn. refs #15646 diff --git src/hg/hgc/gtexEqtlClusterClick.c src/hg/hgc/gtexEqtlClusterClick.c index f25be1f..7152e66 100644 --- src/hg/hgc/gtexEqtlClusterClick.c +++ src/hg/hgc/gtexEqtlClusterClick.c @@ -1,186 +1,188 @@ /* Details page for GTEx eQTL Clusters */ /* Copyright (C) 2017 The Regents of the University of California * See README in this or parent directory for licensing information. */ #include "common.h" #include "hCommon.h" #include "web.h" #include "gtexTissue.h" #include "gtexInfo.h" #include "gtexEqtlCluster.h" #include "hgc.h" static struct gtexEqtlCluster *getGtexEqtl(char *item, char *chrom, int start, int end, char *table) /* Retrieve this item from the track table */ { char *gene = firstWordInLine(cloneString(item)); struct sqlConnection *conn = hAllocConn(database); struct gtexEqtlCluster *eqtls = NULL, *eqtl; char **row; int offset; char where[512]; sqlSafefFrag(where, sizeof(where), "target='%s'", gene); struct sqlResult *sr = hRangeQuery(conn, table, chrom, start, end, where, &offset); while ((row = sqlNextRow(sr)) != NULL) { eqtl = gtexEqtlClusterLoad(row+offset); slAddHead(&eqtls, eqtl); } 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); } static void printMinorAlleleFreq(char *rsId, struct sqlConnection *conn) /* Print minor allele frequency for a SNP (from UCSC dbSNP table) */ { #define SNP_COMMON_SUFFIX "Common" #define MAX_ALLELE_COUNT 10 char *snpTable = hFindLatestSnpTableConn(conn, SNP_COMMON_SUFFIX); if (!snpTable) return; char query[256]; sqlSafef(query, sizeof query, "SELECT alleleFreqs FROM %s WHERE name='%s'", snpTable, rsId); double freqs[MAX_ALLELE_COUNT]; int count = sqlDoubleArray(sqlQuickString(conn, query), freqs, MAX_ALLELE_COUNT); doubleSort(count, freqs); printf("
Minor allele frequency (1000 Genomes): %.0f%%\n", 100.0 * freqs[count-2]); } static void printGwasCatalogTrait(char *rsId, struct sqlConnection *conn) /* Print trait/disease for a SNP (from UCSC gwasCatalog table) */ { 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", rsId); char *trait = sqlQuickString(conn, query); printf("
GWAS disease or trait"); if (count > 1) printf(" (1 of %d)", count); printf(": %s " "GWAS Catalog\n", trait, rsId); } } 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'", table, eqtl->target); int start = sqlQuickNum(conn, query) - FLANK; sqlSafef(query, sizeof query, "SELECT MAX(chromEnd) from %s WHERE target='%s'", table, eqtl->target); int end = sqlQuickNum(conn, query) + FLANK; char posLink[1024]; safef(posLink, sizeof posLink,"%s:%d-%d", hgTracksPathAndSettings(), database, eqtl->chrom, start+1, end, eqtl->chrom, start+1, end); printf("
Region containing eQTLs for this gene: %s (%d bp, including +-%dbp flank)\n", posLink, end-start, FLANK); } 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("\n"); printf(""); -printf("\n"); +printf("\n"); int i; for (i=0; iexpCount; i++) { double effect = eqtl->expScores[i]; + double pval = eqtl->expPvals[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("\n", - color, name, effect < 0 ? "" : "+", effect, prob); + printf("\n", + color, name, effect < 0 ? "" : "+", effect, pval, prob); } printf("
   TissueEffect   Probability
   TissueEffect   P-Value (-log10)Probability
%s%s%0.2f%0.2f
%s%s%0.2f%0.2f%0.2f
"); 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("Gene: "); struct sqlConnection *conn = hAllocConn(database); char *desc = getGeneDescription(conn, geneName); if (desc == NULL) printf("%s\n", geneName); else { printf("%s
\n", hgGeneName(), database, geneName, geneName); printf("Description: %s\n", desc); } // TODO: Consider adding Ensembl gene ID, GENCODE biotype and class (as in gtexGene track) printf("
Variant: %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); +printf("
Distance from TSS: %d\n", eqtl->distance); char posLink[1024]; safef(posLink, sizeof posLink,"%s:%d-%d", hgTracksPathAndSettings(), database, eqtl->chrom, eqtl->chromStart+1, eqtl->chromEnd, eqtl->chrom, eqtl->chromStart+1, eqtl->chromEnd); printf("
Position: %s\n", posLink); printf("
Score: %d\n", eqtl->score); printEqtlRegion(eqtl, tdb->table, conn); printf("
Number of tissues with this eQTL: %d\n", eqtl->expCount); // print link to GTEx portal printf("
" "View eQTL Visualizer for this gene at the GTEx Portal\n", geneName); hFreeConn(&conn); printClusterDetails(eqtl, tdb->table); webNewEmptySection(); printTrackHtml(tdb); }