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>&nbsp;&nbsp;&nbsp;</th><th>Tissue</th><th>Effect &nbsp;&nbsp;</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);
 }