src/hg/instinct/extractData/extractData.c 1.2

Index: src/hg/instinct/extractData/extractData.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/instinct/extractData/extractData.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -b -B -U 1000000 -r1.1 -r1.2
--- src/hg/instinct/extractData/extractData.c	14 Jan 2010 05:18:36 -0000	1.1
+++ src/hg/instinct/extractData/extractData.c	14 Jan 2010 05:37:23 -0000	1.2
@@ -1,216 +1,247 @@
 /* mapProbesToGenes - Will maps probes in BED format to overlapping gene(s). */
 #include "common.h"
 #include "linefile.h"
 #include "hash.h"
 #include "options.h"
 #include "jksql.h"
 #include "bed.h"
 #include "genePred.h"
 #include "hPrint.h"
 #include "hdb.h"  
 #include "microarray.h"
 #include "ra.h"
 #include "featuresLib.h"
 #include "hgHeatmapLib.h"
 #include "cprob.h"
 #include "hgStatsLib.h" 
 
 char *hgDb = "hg18";
 char *genome = "Human";
 
 void usage()
 /* Explain usage and exit. */
 {
 errAbort(
   "extractData \n"
   "   extractData [OPTIONS] db table samplesgenes(s)\n"
   "options:\n"
   "   -median   Output median value if multiple probes\n"
   );
 }
 
 boolean median    = FALSE;
 
 static struct optionSpec options[] = {
     {"median", OPTION_BOOLEAN},
     {NULL, 0},
 };
 
 struct maGrouping *getMaGrouping(struct sqlConnection *hgConn, char *tableName)
 {
 /*microarray specific settings*/
 struct trackDb *tdb = hMaybeTrackInfo(hgConn, tableName);  
 struct microarrayGroups *maGs = maGroupings("hg18", tableName);
 trackDbFreeList(&tdb);
 if (!maGs)
     return NULL;
 return maGs->allArrays;
 }
 
 struct hash *getSettings(char *tableName)
 {
 struct column *raList = getColumns(NULL, "datasets.ra", NULL);  
 
 struct column *col;
 struct hash *settings = NULL;
 for (col = raList; col; col = col->next)
     {
     if (!sameString(col->name, tableName))
 	continue;
 
     settings = col->settings;
     break;
     }
 
 if (!settings)
     errAbort("Couldn't find datasets.ra listing for %s", tableName);
 
 return settings;
 }
 
 struct geneAlias {
     struct geneAlias *next;
 
     char *gene;
     struct slName *probes;
 }; 
  
 struct hash *getAliases(struct sqlConnection *hgConn, char *tableName, 
 			struct slName *genes)
 {
 if (!hgConn || !tableName || !genes)
     return NULL;
 
 struct slName *sl;
 struct dyString *dy = newDyString(100);
 dyStringPrintf(dy,
 	       "select * from %s where alias in (", 
 	       tableName);
 for (sl = genes; sl; sl = sl->next)
     {
     dyStringPrintf(dy, "\"%s\"", sl->name);
     if (sl->next)
 	dyStringPrintf(dy, ",");
     }
 dyStringPrintf(dy, ");");
 char *query = dyStringCannibalize(&dy);
 
 struct sqlResult *sr = sqlGetResult(hgConn, query);
 
 struct geneAlias *ga, *gaList = NULL;
 struct hash *gaHash = hashNew(0);
 
 char **row;
 while ((row = sqlNextRow(sr)) != NULL)
     {
     char *probe = cloneString(row[0]);
     char *gene = cloneString(row[1]);
 
     struct hashEl *el = hashLookup(gaHash, gene);
 
     if (!el)
 	{
 	ga = AllocA(struct geneAlias);
 	ga->gene = cloneString(gene);
 	ga->probes = NULL;
 
 	slAddHead(&gaList, ga);
 	hashAdd(gaHash, gene, ga);
 	}
     else
 	ga = el->val;
 
     slNameAddHead(&ga->probes, probe);	
     }
 
 sqlFreeResult(&sr);              
 
 return gaHash;                   
 }
 
-void getProbeData(struct sqlConnection *hgConn, struct maGrouping *allA, 
+struct slDouble *getProbeData(struct sqlConnection *hgConn, struct maGrouping *allA, 
 		  char *dataTable, struct slName *probes)
 {
 struct slName *sl;
 struct dyString *dy = newDyString(100);
 dyStringPrintf(dy,
 	       "select * from %s where name in(", 
 	       dataTable);
 for (sl = probes; sl; sl = sl->next)
     {
     dyStringPrintf(dy, "\"%s\"", sl->name);
     if (sl->next)
 	dyStringPrintf(dy, ",");
     }
 dyStringPrintf(dy, ");");
 char *query = dyStringCannibalize(&dy);
 
-fprintf(stderr, "query: %s\n", query);
 /* Get bed15 data from hg18 database */
 struct sqlResult *sr = sqlGetResult(hgConn, query);
 
 int count = 0;
 char **row = NULL;
+struct slDouble *sdList = NULL;
+
+int i = 1;
 while ((row = sqlNextRow(sr)) != NULL)
     {
     struct bed *nb = bedLoadN(row+1, 15);
   
 // nb->expCount
 // nb->expScores
     count += 1; 
 
+    if (i >= nb->expCount)
+	continue;
+
+    struct slDouble *sd = slDoubleNew(nb->expScores[i]);
+    slAddHead(&sdList, sd);
+
     bedFree(&nb);
     }
 sqlFreeResult(&sr);
 
-fprintf(stderr, "found %d probes\n", count);
+return sdList;
 }
 
 void extractData(char *tableName, char *sampleName, char *geneList)
 {
 if (!geneList)
     return;
 
 struct slName *genes = slNameListFromComma(geneList);
 
 struct hashEl *el;
 struct sqlConnection *hgConn = hAllocConnProfile("localDb", hgDb);
 
 /* Set up datasets entry */
 struct maGrouping *allA = getMaGrouping(hgConn, tableName);
 if (!allA)
     errAbort("Could not find maGrouping for %s!", tableName);
 
 struct hash *settings = getSettings(tableName);
 
 el = hashLookup(settings, "aliasTable");
 if (!el)
     errAbort("No aliasTable.\n");
 char *aliasTable = cloneString(el->val);
 
 struct hash *gaHash = getAliases(hgConn, aliasTable, genes);
 
 struct hashCookie cookie = hashFirst(gaHash);
 while ((el = hashNext(&cookie)) != NULL)
     {
+    char *name = el->name;
     struct geneAlias *ga = el->val;    
-    getProbeData(hgConn, allA, tableName, ga->probes);
+    struct slDouble *sd, *sdList = getProbeData(hgConn, allA, tableName, ga->probes);
+
+    if (!sdList)
+	continue;
+
+    printf("%s\t", name);
+    if (median)
+	{
+	double val = slDoubleMedian(sdList);
+	printf("%f\n", val);
+	}
+    else
+	{
+	struct dyString *dy = newDyString(100);
+	for (sd = sdList; sd; sd = sd->next)
+	    {
+	    dyStringPrintf(dy, "%f", sd->val);
+	    if (sd->next)
+		dyStringPrintf(dy, ",");
+	    }
+	char *val = dyStringCannibalize(&dy);
+	printf("%s\n", val);
+	}
     }
 
 hFreeConn(&hgConn);
 }
 
 
 int main(int argc, char *argv[])
 /* Process command line. */
 {
 optionInit(&argc, argv, options);
 if (argc != 4)
     usage();
 
 if (optionExists("median"))
     median = TRUE;
 
 extractData(argv[1], argv[2], argv[3]);
 return 0;
 }