ef7f905f8875763947ebe09110b0c3205e567cfa
ceisenhart
  Wed May 20 11:15:40 2015 -0700
REFS #15377 fixing an issue with the usage statement

diff --git src/utils/bigWigCluster/bigWigCluster.c src/utils/bigWigCluster/bigWigCluster.c
index 0572a9c..f2d0a4b 100644
--- src/utils/bigWigCluster/bigWigCluster.c
+++ src/utils/bigWigCluster/bigWigCluster.c
@@ -1,416 +1,416 @@
 /* bigWigCluster - Cluster bigWigs using a hactree. */
 #include "sqlNum.h"
 #include "common.h"
 #include "linefile.h"
 #include "hash.h"
 #include "options.h"
 #include "hacTree.h"
 #include "rainbow.h"
 #include "portable.h"
 
 int gThreadCount = 10;
 char *gTmpDir = ".";
 
 void usage()
 /* Explain usage and exit. */
 {
 errAbort(
   "bigWigCluster - Cluster bigWigs using a hacTree\n"
   "usage:\n"
   "   bigWigCluster input.list chrom.sizes output.json output.tab\n"
   "where: input.list is a list of bigWig file names\n"
   "       chrom.sizes is tab separated <chrom><size> for assembly for bigWigs\n"
   "       output.json is json formatted output suitable for graphing with D3\n"
   "       output.tab is tab-separated file of  of items ordered by tree with the fields\n"
   "           label - label from -labels option or from file name with no dir or extention\n"
   "           pos - number from 0-1 representing position according to tree and distance\n"
   "           red - number from 0-255 representing recommended red component of color\n"
   "           green - number from 0-255 representing recommended green component of color\n"
   "           blue - number from 0-255 representing recommended blue component of color\n"
   "           path - file name from input.list including directory and extension\n" 
   "options:\n"
   "   -labels=fileName - label files from tabSeparated file with fields\n"
   "           path - path to bigWig file\n"
   "           label - a string with no tabs\n"
   "   -precalc=precalc.tab - tab separated file with <file1> <file2> <distance>\n"
   "            columns.\n"
   "   -threads=N - number of threads to use, default %d\n"
   "   -tmpDir=/tmp/path - place to put temp files, default current dir\n"
   , gThreadCount
   );
 }
 
 char tmpPrefix[PATH_LEN] = "bwc_tmp_";  /* Prefix for temp file names */
 
 /* Command line validation table. */
 static struct optionSpec options[] = {
    {"precalc", OPTION_STRING},
    {"threads", OPTION_INT},
    {"labels", OPTION_STRING},
    {"tmpDir", OPTION_STRING},
    {NULL, 0},
    {"threads", OPTION_INT},
    {"hacTree", OPTION_STRING},
 };
 
 
 double longest = 0;
 char* clHacTree = "fromItems";
 int clThreads = 10; // The number of threads to run with the multiThreads option
 int nameCount = 0;
 
 int bigWigNextId;
 
 struct bigWig
 /* Information on a bigWig */
     {
     struct bigWig *next;  //next item in series
     int id;	// Unique numerical identifier
     char *fileName; // Associated file name
     char *label;    // Associated labels if any
     struct rgbColor color; //for coloring
     double pos;	    // Position between 0.0 and 1.0 when ordered by tree and distance
     };
 
 struct bigWig *readBigWigList(char* input)
 // get the bigWig files 
 {
 struct bigWig *list = NULL;
 struct lineFile *lf = lineFileOpen(input,TRUE);
 char *line = NULL;
 while(lineFileNextReal(lf, &line))
     {
     /* Allocate new item and initialize id and fileName fields */
     struct bigWig *bw;
     AllocVar(bw);
     bw->id = ++bigWigNextId;
     bw->fileName = cloneString(trimSpaces(line));
 
     // Make up default label - same as file name without directory and extension
     char root[FILENAME_LEN];
     splitPath(line, NULL, root, NULL);
     bw->label = cloneString(root);
 
     /* Add to list */
     slAddHead(&list,bw);
     }
 
 /* Clean up and go home */
 lineFileClose(&lf);
 slReverse(&list);
 return list;
 }
 
 
 static void rAddLeaf(struct hacTree *tree, struct slRef **pList)
 /* Recursively add itemOrCluster from leaf nodes of hacTree to list */
 {
 if (tree->left == NULL && tree->right == NULL)
     refAdd(pList, tree->itemOrCluster);
 else
     {
     rAddLeaf(tree->left, pList);
     rAddLeaf(tree->right, pList);
     }
 }
 
 struct slRef *getOrderedLeafList(struct hacTree *tree)
 /* Return list of references to bigWigs from leaf nodes
  * ordered by position in hacTree */
 {
 struct slRef *leafList = NULL;
 rAddLeaf(tree, &leafList);
 slReverse(&leafList);
 return leafList;
 }
 
 static void rPrintHierarchicalJson(FILE *f, struct hacTree *tree, int level, double distance,
 		int normConstant, int cgConstant)
 /* Recursively prints out the elements of the hierarchical .json file. */
 {
 struct bigWig *bw = (struct bigWig *)tree->itemOrCluster;
 char *label = bw->label;
 struct rgbColor colors = bw->color;
 int i;
 for (i = 0;  i < level;  i++)
     fputc(' ', f); // correct spacing for .json format
 if (tree->left == NULL && tree->right == NULL)
     {
     // Prints out the leaf objects
     fprintf(f, "{\"%s\"%s \"%s\"%s\"%s\"%s %f %s\"%s\"%s \"rgb(%i,%i,%i)\"}", "name", ":", label, ", ",
     		"similarity", ":", distance, "," , "colorGroup", ":", colors.r, colors.g, colors.b);
     return;
     }
 else if (tree->left == NULL || tree->right == NULL)
     errAbort("\nHow did we get a node with one NULL kid??");
 
 // Prints out the node object and opens a new children block
 fprintf(f, "{\"%s\"%s \"%s\"%s", "name", ":", " ", ",");
 fprintf(f, "\"colorGroup\": \"rgb(%i,%i,%i)\",", colors.r, colors.g, colors.b );
 fprintf(f, "\"%s\"%s \"%f\"%s\n", "distance", ":", 100*tree->childDistance, ",");
 for (i = 0;  i < level + 1;  i++)
     fputc(' ', f);
 fprintf(f, "\"%s\"%s\n", "children", ": [");
 distance = tree->childDistance; 
 rPrintHierarchicalJson(f, tree->left, level+1, distance, normConstant, cgConstant);
 fputs(",\n", f);
 rPrintHierarchicalJson(f, tree->right, level+1, distance, normConstant, cgConstant);
 fputc('\n', f);
 // Closes the children block for node objects
 for (i=0;  i < level + 1;  i++)
     fputc(' ', f);
 fputs("]\n", f);
 for (i = 0;  i < level;  i++)
     fputc(' ', f);
 fputs("}", f);
 }
 
 void printHierarchicalJson(char *fileName, struct hacTree *tree, int normConstant, int cgConstant)
 /* Prints out the binary tree into .json format intended for d3
  * hierarchical layouts */
 {
 FILE *f = mustOpen(fileName, "w");
 if (tree == NULL)
     {
     fputs("Empty tree.\n", f);
     return;
     }
 rPrintHierarchicalJson(f, tree, 0, 0, normConstant, cgConstant);
 fputc('\n', f);
 carefulClose(&f);
 }
 
 void printOrderedTab(char *fileName, struct slRef *refList)
 /* Print list of references to bigWigs */
 {
 FILE *f = mustOpen(fileName, "w");
 struct slRef *ref;
 for (ref = refList; ref != NULL; ref = ref->next)
     {
     struct bigWig *bw = ref->val;
     fprintf(f, "%s\t%g\t%d\t%d\t%d\t%s\n", bw->label, bw->pos,
 	bw->color.r, bw->color.g, bw->color.b, bw->fileName);
     }
 carefulClose(&f);
 }
 
 double slBigWigDistance(const struct slList *item1, const struct slList *item2, void *extraData)
 /* Return the absolute difference between the two kids' values. 
  * Designed for HAC tree use*/
 {
 verbose(1,"Calculating Distance...\n");
 const struct bigWig *kid1 = (const struct bigWig *)item1;
 const struct bigWig *kid2 = (const struct bigWig *)item2;
 char tmpName[PATH_LEN];
 safef(tmpName, sizeof(tmpName), "%s%d_%d.txt", tmpPrefix, kid1->id, kid2->id);
 char cmd[1024];
 safef(cmd, 1024, "bigWigCorrelate %s %s > %s", kid1->fileName, kid2->fileName, tmpName);
 double diff = 0;
 mustSystem(cmd);
 struct lineFile *lf = lineFileOpen(tmpName,TRUE);
 char *line = NULL;
 if (!lineFileNext(lf, &line, NULL))
     errAbort("no difference output, check bigWigCorrelate");
 lineFileClose(&lf);
 diff = sqlDouble(line);
 remove(tmpName);
 return  1 - diff;
 }
 
 
 struct slList *slBigWigMerge(const struct slList *item1, const struct slList *item2,
 				void *extraData)
 /* Make a new slPair where the name is the children names concattenated and the 
  * value is the average of kids' values.
  * Designed for HAC tree use*/
 {
 char *chromSizesFile = extraData;
 verbose(1,"Merging...\n");
 struct bigWig *result;
 AllocVar(result);
 const struct bigWig *kid1 = (const struct bigWig *)item1;
 const struct bigWig *kid2 = (const struct bigWig *)item2;
 char tmpName[PATH_LEN];
 safef(tmpName, sizeof(tmpName), "%s%d_%d.bedGraph", tmpPrefix, kid1->id, kid2->id);
 char cmd1[1024];
 char cmd2[1024];
 safef(cmd1, 1024, "bigWigMerge %s %s %s -verbose=0", kid1->fileName, kid2->fileName, tmpName);
 char fileName[PATH_LEN];
 safef(fileName, sizeof(fileName), "%s%d_%d.bw", tmpPrefix, kid1->id, kid2->id);
 safef(cmd2, 1024, "bedGraphToBigWig %s %s %s", tmpName, chromSizesFile, fileName);
 mustSystem(cmd1);
 mustSystem(cmd2);
 remove(tmpName);
 result->id = ++bigWigNextId;
 result->fileName = cloneString(fileName);
 result->label = " ";
 return (struct slList *)result;
 }
 
 double findCachedDistance(struct hash *distanceHash, struct bigWig *bw1, struct bigWig *bw2)
 /* Return cached distance - from hash if possible */
 {
 double distance;
 double *cached = hacTreeDistanceHashLookup(distanceHash, bw1, bw2);
 if (cached != NULL)
    distance = *cached;
 else
    {
    distance = slBigWigDistance((struct slList *)bw1, (struct slList *)bw2, NULL);
    hacTreeDistanceHashAdd(distanceHash, bw1, bw2, distance);
    }
 return distance;
 }
 
 
 void colorLeaves(struct slRef *leafList, struct hash *distanceHash)
 /* Assign colors of rainbow to leaves. */
 {
 /* Loop through list once to figure out total, since we need to
  * normalize */
 double total = 0;
 double purplePos = 0.80;
 struct slRef *el, *nextEl;
 for (el = leafList; el != NULL; el = nextEl)
    {
    nextEl = el->next;
    if (nextEl == NULL)
        break;
    double distance = findCachedDistance(distanceHash, el->val, nextEl->val);
    total += distance;
    }
 
 /* Loop through list a second time to generate actual colors. */
 double soFar = 0;
 for (el = leafList; el != NULL; el = nextEl)
    {
    nextEl = el->next;
    if (nextEl == NULL)
        break;
    struct bigWig *bw = nextEl->val;
    double distance = findCachedDistance(distanceHash, el->val, bw);
    soFar += distance;
    double normalized = soFar/total;
    bw->pos = normalized;
    bw->color = saturatedRainbowAtPos(normalized * purplePos);
    }
 
 /* Set first color to correspond to 0, since not set in above loop */
 struct bigWig *bw = leafList->val;
 bw->color = saturatedRainbowAtPos(0);
 }
 
 void preloadDistances(char *distanceFile, struct bigWig *itemList, struct hash *distanceHash)
 /* Load tab-separated distanceFile into a hash, making sure all items
  * in distanceFile correspond to file in list */
 {
 /* Create hash of items in list keyed by file name */
 struct hash *itemHash = hashNew(0);
 struct bigWig *item;
 for (item = itemList; item != NULL; item = item->next)
     hashAdd(itemHash, item->fileName, item);
 
 /* Loop through distance file's three columns building up distance hash. */
 struct lineFile *lf = lineFileOpen(distanceFile, TRUE);
 char *row[3];
 while (lineFileRow(lf, row))
     {
     char *a = row[0], *b = row[1];
     double correlation = sqlDouble(row[2]);
     double distance = 1.0 - correlation;
     struct bigWig *aItem = hashMustFindVal(itemHash, a);
     struct bigWig *bItem = hashMustFindVal(itemHash, b);
     hacTreeDistanceHashAdd(distanceHash, aItem, bItem, distance);
     hacTreeDistanceHashAdd(distanceHash, bItem, aItem, distance);
     }
 
 /* Report, clean up and go home. */
 verbose(1, "Preloaded %d distances on %d items from %s\n", distanceHash->elCount, itemHash->elCount,
     distanceFile);
 lineFileClose(&lf);
 hashFree(&itemHash);
 }
 
 void addLabels(char *labelFile, struct bigWig *list)
 /* Update label field of bigWig list according to labels in tab-separated file of format
  *    <fileName> <labels> */
 {
 /* Build up hash of labels keyed by file name */
 struct hash *hash = hashNew(0);
 struct lineFile *lf = lineFileOpen(labelFile, TRUE);
 char *row[2];
 while (lineFileRowTab(lf, row))
      hashAdd(hash, row[0], cloneString(row[1]));
 lineFileClose(&lf);
 
 /* Loop through list looking up labels in hash */
 struct bigWig *bw;
 for (bw = list; bw != NULL; bw = bw->next)
      {
      bw->label = hashFindVal(hash, bw->fileName);
      if (bw->label == NULL)
          errAbort("No label for %s\n", bw->fileName);
      }
 
 /* Clean up and go home */
 freeHash(&hash);
 }
 
 void bigWigCluster(char *inputList, char *chromSizes, char *outputJson, char *outputTab)
 /* bigWigCluster - Cluster bigWigs using a hactree. */
 {
 /* Read in input hash */
 struct bigWig *list = readBigWigList(inputList);
 
 /* Set up distance cache, preloading it if possible from file */
 struct hash *distanceHash = hashNew(0);
 char *precalcFile = optionVal("precalc", NULL);
 if (precalcFile)
     preloadDistances(precalcFile, list, distanceHash);
 
 /* Supply better labels if they exist */
 char *labelFile = optionVal("labels", NULL);
 if (labelFile)
     addLabels(labelFile, list);
 
 struct hacTree *clusters;
 struct lm *localMem = lmInit(0);
 if (gThreadCount > 1)
     {
     clusters = hacTreeMultiThread(gThreadCount, (struct slList *)list, localMem,
 					    slBigWigDistance, slBigWigMerge, chromSizes, distanceHash);
     }
 else
     {
     /* Use older code for single threaded case, just to be able to compare results to
      * parallelized version */
     clusters = hacTreeFromItems((struct slList *)list, localMem,
 					    slBigWigDistance, slBigWigMerge, NULL, chromSizes);
     }
 
 /* Convert tree to ordered list, do coloring, and make outputs */
 struct slRef *orderedList = getOrderedLeafList(clusters);
 colorLeaves(orderedList, distanceHash);
 printHierarchicalJson(outputJson, clusters, 20, 20);
 printOrderedTab(outputTab, orderedList);
 
 // Clean up remaining temp files 
 char cmd[1024];
 safef(cmd, sizeof(cmd), "rm %s*", tmpPrefix);
 mustSystem(cmd);
 }
 
 int main(int argc, char *argv[])
 /* Process command line. */
 {
+optionInit(&argc, argv, options);
 clThreads = optionInt("threads", clThreads);
 clHacTree = optionVal("hacTree", clHacTree);
-optionInit(&argc, argv, options);
 if (argc != 5)
     usage();
 gThreadCount = optionInt("threads", gThreadCount);
 gTmpDir = optionVal("tmpDir", gTmpDir);
 safef(tmpPrefix, sizeof(tmpPrefix), "%s/bwc_tmp_", gTmpDir);
 bigWigCluster(argv[1], argv[2], argv[3], argv[4]);
 return 0;
 }