Index: src/hg/instinct/hgHeatmap2/drawingCode.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/instinct/hgHeatmap2/drawingCode.c,v
retrieving revision 1.75
retrieving revision 1.76
diff -b -B -U 1000000 -r1.75 -r1.76
--- src/hg/instinct/hgHeatmap2/drawingCode.c	5 Oct 2009 21:32:30 -0000	1.75
+++ src/hg/instinct/hgHeatmap2/drawingCode.c	6 Oct 2009 23:40:06 -0000	1.76
@@ -1,3361 +1,3368 @@
 /********************************************************************************/
 /* Copyright 2007-2009 -- The Regents of the University of California           */
 /********************************************************************************/
 
 /* mainPage - drs the main hgHeatmap page, including some controls
  * on the top and the graphic. */
 
 #define EXPR_DATA_SHADES 16
 #define DEFAULT_MAX_DEVIATION 0.7
 #define COLOR_SCALE 1
 #define MAX_PCA_COMPONENTS 2
 //#define DBL_MAX 1000000.0  // Changed for 32-bit computers compatibility
 #define PROBE_MARK_HEIGHT 3
 #define PROBE_MARK_BUFFER 2
 
 #include <limits.h>
 #include "common.h"
 #include "bed.h"
 #include "cart.h"
 #include "customTrack.h"
 #include "errCatch.h"
 #include "genoLay.h"
 #include "trackLayout.h"
 #include "hash.h"
 #include "hdb.h"
 #include "hgColors.h"
 #include "hPrint.h"
 #include "htmshell.h"
 #include "jsHelper.h"
 #include "psGfx.h"
 #include "trashDir.h"
 #include "vGfx.h"
 #include "hvGfx.h"
 #include "web.h"
 #include "cytoBand.h"
 #include "hCytoBand.h"
 #include "hgHeatmapLib.h"
 #include "heatmapUtility.h"
 #include "featuresLib.h"
 #include "hgStatsLib.h"
 #include "hgHeatmap2.h"
 #include "filterFeatures.h"
 #include "chromGraph.h"
 #include "hgChromGraph.h"
 
 static char const rcsid[] = "$Id$";
 
 static char *heatMapDbProfile = "localDb";  // database profile to use
 
 struct gifList {
     struct gifList *next;
 
     int expId;
     char *filename;
 };
 
 struct chromZoom {
     struct chromZoom *next;
     
     char *fullName;
     int size;
     int baseStart;
     int baseEnd;
 };
 
 struct trackLayout tl;           /* Dimensions of things, fonts, etc. */
 char *blockStatTest[3];          /* compatability with library code */
 
 
 
 
 #define CLIP(p,limit) if (p < 0) p = 0; if (p >= (limit)) p = (limit)-1;  
 void verticalTextCentered(struct vGfx *vg, int x, int y,
                                  int width, int height, int colorIx, MgFont *font, char *string)
 /* Draw a vertical line of text in middle of the box.
  *      The string will read from bottom to top
  * This is not perfect by any means, but seems to draw the labels fairly nice
  * for the pathways.
  */
 {
 /*      don't let this run wild */
 CLIP(width, vg->width);
 CLIP(height, vg->height);
 
 if ((width <= 0) || (height <= 0))
     return;
 
 struct vGfx *vgHoriz;
 int i, j;
 /*  reversed meanings of width and height here since this is going
  *  to rotate 90 degrees
  */
 int offset = 0;
 int cutoff = height / 7;  // TODO, approximately correct.
 int lineSpacing = 15;
 
 vgHoriz = vgOpenGif(height, width, "/dev/null", FALSE);
 struct slName *sl, *slList = slNameListFromString(string, '_');
 
 struct dyString *dy = newDyString(0), *dyList = NULL;
 for (sl = slList; sl; sl = sl->next)
     { /* loop through words in text, separated by underscores */
     if ( (dy->stringSize > 0) && (offset < width) &&
 	 (dy->stringSize + strlen(sl->name) > cutoff))
 	{ /* enough string for a single line is ready to be drawn and can be drawn */
 	slAddHead(&dyList, dy);
 	dy = newDyString(0);
 	offset += lineSpacing;
 	}
     dyStringAppend(dy, sl->name);
     dyStringAppend(dy, " ");
     }
 
 if ( (dy->stringSize > 0) && (offset < width) )
     /* Clean up remaining string */
     slAddHead(&dyList, dy);
 
 slReverse(&dyList);
 int count = slCount(dyList);
 
 offset = 0;
 for (dy = dyList; dy ; dy = dy->next)
     { /* -3 is needed to approximately recenter text */
     vgTextCentered(vgHoriz, 0, -3*(count - 1), height, width+offset,
 		   colorIx, font, dy->string);
     offset += lineSpacing;
     }
 dyStringFreeList(&dyList);
 
 /*  now, blit from the horizontal to the vertical, rotate -90 (CCW) */
 for (i = 0; i < height; ++i)        /* xSrc -> yDest */
     {
     int yDest = height - i;
     for (j = 0; j < width; ++j)     /* ySrc -> xDest */
 	vgDot(vg, j+x, yDest+y, vgGetDot(vgHoriz, i, j));
     }
 vgClose(&vgHoriz);
 }
 
 void verticalTextRight(struct vGfx *vg, int x, int y,
 		       int width, int height, int colorIx, MgFont *font, char *string)
 /* Draw a vertical line of text in middle of the box.
  *      The string will read from bottom to top
  * This is not perfect by any means, but seems to draw the labels fairly nice
  * for the pathways.
  */
 {
 /*      don't let this run wild */
 CLIP(width, vg->width);
 CLIP(height, vg->height);
 
 if ((width <= 0) || (height <= 0))
     return;
 
 struct vGfx *vgHoriz;
 int i, j;
 /*  reversed meanings of width and height here since this is going
  *  to rotate 90 degrees
  */
 int offset = 0;
 int cutoff = height / 7;  // TODO, approximately correct.
 int lineSpacing = 15;
 
 vgHoriz = vgOpenGif(height, width, "/dev/null", FALSE);
 string = replaceChars(string, "_", " ");
 struct slName *sl, *slList = slNameListFromString(string, ' ');
 
 struct dyString *dy = newDyString(0), *dyList = NULL;
 for (sl = slList; sl; sl = sl->next)
     { /* loop through words in text, separated by underscores */
     if ( (dy->stringSize > 0) && (offset < width) &&
 	 (dy->stringSize + strlen(sl->name) > cutoff))
 	{ /* enough string for a single line is ready to be drawn and can be drawn */
 	slAddHead(&dyList, dy);
 	dy = newDyString(0);
 	offset += lineSpacing;
 	}
     dyStringAppend(dy, sl->name);
     dyStringAppend(dy, " ");
     }
 
 if ( (dy->stringSize > 0) && (offset < width) )
     /* Clean up remaining string */
     slAddHead(&dyList, dy);
 
 slReverse(&dyList);
 int count = slCount(dyList);
 
 offset = 0;
 for (dy = dyList; dy ; dy = dy->next)
     { /* -3 is needed to approximately recenter text */
     vgTextRight(vgHoriz, 0, -3*(count - 1), height + 2, width+offset,
 		colorIx, font, dy->string);
     offset += lineSpacing;
     }
 dyStringFreeList(&dyList);
 
 /*  now, blit from the horizontal to the vertical, rotate -90 (CCW) */
 for (i = 0; i < height; ++i)        /* xSrc -> yDest */
     {
     int yDest = height - i;
     for (j = 0; j < width; ++j)     /* ySrc -> xDest */
 	vgDot(vg, j+x, yDest+y, vgGetDot(vgHoriz, i, j));
     }
 vgClose(&vgHoriz);
 }
 
 
 boolean singleChromImage(struct heatmapLay *hl)
 {
 if (!hl)
     return FALSE;
 
 if (slCount(hl->elements) == 1)
     return TRUE;
 
 return FALSE;
 }
 
 /* return an array for reordering the experiments in a chromosome */
 double maxDeviation(char* heatmap)
 {
 struct hashEl *e = hashLookup(ghHash, heatmap);
 
 if (!e)
     errAbort("Could not found heatmap %s", heatmap);
 
 struct genoHeatmap *gh = e->val;
 if (gh->expScale >0)
     return gh->expScale;
 
 return DEFAULT_MAX_DEVIATION;
 }
 
 void vgMakeColorGradient(struct vGfx *vg,
     struct rgbColor *start, struct rgbColor *end,
     int steps, Color *colorIxs)
 /* Make a color gradient that goes smoothly from start
  * to end colors in given number of steps.  Put indicesgl->chromLis
  * in color table in colorIxs */
 {
 double scale = 0, invScale;
 double invStep;
 int i; int r,g,b;
 
 steps -= 1;     /* Easier to do the calculation in an inclusive way. */
 invStep = 1.0/steps;
 for (i=0; i<=steps; ++i)
     {
     invScale = 1.0 - scale;
     r = invScale * start->r + scale * end->r;
     g = invScale * start->g + scale * end->g;
     b = invScale * start->b + scale * end->b;
     colorIxs[i] = vgFindColorIx(vg, r, g, b);
     scale += invStep;
     }
 }
 
 char *getId(struct sqlConnection *conn, char *table, char *key, char *sample, char *value)
 /* get patient ID from sample (or experiment) Id */
 {
 char query[512]; 
 safef(query, sizeof(query), "select %s from %s where %s = '%s' ", key, table, value, sample);
 return sqlQuickString(conn, query);
 }
 
 char *getPathwayDb()
 /* Get pathway database */
 {
 return "pathway";
 }
 
 struct sqlConnection *getPathwayDbConn()
 /* Get connection to pathway database */
 {
 struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, getPathwayDb());
 return conn;
 }
                     
 struct sqlConnection *getFeatureDbConn(struct genoHeatmap *gh)
 /* Get feature database connection */
 {
 char *featureDb = gh->patDb;
 char *raName = gh->raFile;
 
 if ((featureDb) && (raName))
     return hAllocConnProfile(heatMapDbProfile, featureDb);
 else
     return NULL;
 } 
 
 void heatmapLaySetColor(struct heatmapLay *hlList, struct rgbColor *low, struct rgbColor *zero, 
 			struct rgbColor *high)
 {
 struct heatmapLay *hl;
 
 for (hl = hlList; hl; hl = hl->next)
     {
     hl->lowColor = low;
     hl->zeroColor = zero;
     hl->highColor = high;
     }
 }
 
 void heatmapLaySetDefaultColor(struct heatmapLay *hl, char *platform)
 {
 struct rgbColor *low = NULL;
 struct rgbColor *zero = NULL;
 struct rgbColor *high = NULL;
 
 // default, expression data color scheme:
 if (sameString(platform, "CNV") || sameString(platform, "SNP"))
     { // CNV, SNP data
     AllocVar(zero);  // white
     zero->r = 255;
     zero->g = 255;
     zero->b = 255;
     
     AllocVar(high);  // red
     high->r = 255;
     high->g = 0;
     high->b = 0;
     
     AllocVar(low);  // blue
     low->r = 0;
     low->g = 0;
     low->b = 255;
     }
 else
     { // expression data, others...
     AllocVar(zero);  // black
     zero->r = 0;
     zero->g = 0;
     zero->b = 0;
     
     AllocVar(high);  // red
     high->r = 255;
     high->g = 0;
     high->b = 0;
     
     AllocVar(low);  // green
     low->r = 0;
     low->g = 255;
     low->b = 0;
     }
 
 heatmapLaySetColor(hl, low, zero, high);
 }
 
 void heatmapLayResetMinMax(struct heatmapLay *hlList)
 {
 if (!hlList)
     return;
 
 struct heatmapLay *hl;
 for (hl = hlList; hl; hl = hl->next)
     {
     hl->minVal = INT_MAX;
     hl->minValPos = -1;
     hl->maxVal = INT_MIN;
     hl->maxValPos = -1;
     }
 }
 
 void drawLayoutLines(struct vGfx *vg, struct heatmapLay *hlList, int width)
 {
 if (!hlList)
     return;
 
 struct heatmapLay *hl;
 for (hl = hlList; hl; hl = hl->next)
     {
     int leftX = hl->pixelStart - 1;
     if (leftX > 0)
 	vgBox(vg, leftX, 0, 1, hlList->height, MG_GRAY);
     }
 }
 
 struct heatmapLay *chromLayout(struct sqlConnection *conn, struct genoHeatmap *gh, 
 			       char *chromStart, int startBase,
 			       char *chromEnd, int endBase,
 			       int width, int sampleHeight)
 {
 /* Get list of chromosomes and lay them out. */
 struct sqlConnection *hg18conn = hAllocConnProfile(heatMapDbProfile, database);
 struct genoLayChrom *chrom, *allChromList = genoLayDbChroms(hg18conn, FALSE);
 
 if (!allChromList)
     return NULL;
 
 unsigned long totalBases = 0;
 int chromCount = 0;
 int keepChroms = 0;
 
 if (!chromStart)
     chromStart = cloneString(allChromList->fullName);
 if (!chromEnd)
     for (chrom = allChromList; chrom; chrom = chrom->next)
 	if (!chrom->next)
 	    chromEnd = cloneString(chrom->fullName);
 
 struct chromZoom *cz, *czList = NULL; 
 for (chrom = allChromList; chrom != NULL; chrom = chrom->next)
     {
     if (sameString(chrom->fullName, chromStart))
 	keepChroms = 1;
 	
     if (keepChroms)
 	{
 	AllocVar(cz);
 	cz->fullName = cloneString(chrom->fullName);
 
 	cz->baseStart = -1;
 	cz->baseEnd = -1;
 	if (sameString(chrom->fullName, chromStart) && (startBase != -1))
 	    cz->baseStart = startBase;
 
 	if (sameString(chrom->fullName, chromEnd) && (endBase != -1))
 	    {
 	    if (endBase > chrom->size)
 		endBase = chrom->size;
 	    cz->baseEnd = endBase + 1;
 	    }
 
 	if (cz->baseStart < 0)
 	    cz->baseStart = 0;
 	if (cz->baseEnd < 0)
 	    cz->baseEnd = chrom->size;
 
 	cz->size = cz->baseEnd - cz->baseStart;
 
 	slAddHead(&czList, cz);
 	totalBases += cz->size;
 	chromCount += 1;
 	}
 
     if (sameString(chrom->fullName, chromEnd))
 	keepChroms = 0;
     }
 
 slReverse(&czList);
 
 int availablePixels = width - (chromCount - 1);
 double basesPerPixel = (double) totalBases / availablePixels;
 
 struct heatmapLay *hm = AllocA(struct heatmapLay);
 struct hmElement *hEl, *heList = NULL;
 
 hm->name = cloneString("chromInfo");
 hm->height = 0;
 hm->sampleHeight = sampleHeight;
 hm->width = width;
 hm->basesPerPixel = basesPerPixel;
 
 trackLayoutInit(&tl, cart);
 hm->fontHeight = tl.fontHeight;
 hm->font = tl.font;
 
 hm->minVal = INT_MAX;
 hm->minValPos = -1;
 hm->maxVal = INT_MIN;
 hm->maxValPos = -1;
 
 int xOff = 0;
 for (cz = czList; cz != NULL; cz = cz->next)
     {
     unsigned long bStart = cz->baseStart;
     unsigned long bEnd = cz->baseEnd;
 
     int pixStart = round(0 / basesPerPixel) + xOff;
     int pixEnd = round((bEnd-bStart) / basesPerPixel) + xOff;
 
     if (pixEnd >= width)
 	pixEnd = width;
 
     AllocVar(hEl);
 
     hEl->name = cloneString(cz->fullName);
     hEl->pixelStart = pixStart;
     hEl->pixelEnd = pixEnd;
 
     hEl->baseStart = bStart;
     hEl->baseEnd = bEnd;
     
     hEl->shortLabel = NULL;
     hEl->longLabel = NULL;
     hEl->probeName = NULL;
 
     slAddHead(&heList, hEl);
     xOff = pixEnd + 1;
     }
 
 slReverse(&heList);
 hm->elements = heList;
 
 if (gh)
     heatmapLaySetDefaultColor(hm, gh->platform);
 
 return hm;
 }  
 
 struct heatmapLay *genesetLayout(struct genoHeatmap *gh,
 				 char *geneSetNames,
 				 int width, int sampleHeight)
 {
 struct geneSet *geneSets = getAllPathways(cart, getPathwayDb(), geneSetNames);
 
 geneSetsChromLayout(geneSets); 
 
 boolean isSNP = FALSE;
 struct hash *geneHash = NULL;
 if (gh)
     {
     getChromHeatmapHash(&geneHash, gh->database, gh->probeTable,
 			gh->name, NULL, geneSets);
     if (!geneHash)
 	return NULL;
 
     if (sameString(gh->platform, "SNP"))
 	isSNP = TRUE;
     }
 else
     { /* for geneset labels, gh == NULL, but still want layout) */
     geneHash = hashNew(0);
     }
 
 
 geneSetsAddPathwayLayout(geneSets, geneHash);
 trackLayoutInit(&tl, cart);
 
 struct geneSet *gs;
 struct heatmapLay *hl, *hlList = NULL;
 for (gs = geneSets; gs; gs = gs->next)
     {
     AllocVar(hl);
     hl->name = cloneString(gs->name);
     hl->height = 0;
     hl->sampleHeight = sampleHeight;    
     hl->width = gs->width;
 
     hl->fontHeight = tl.fontHeight;
     hl->font = tl.font;
 
     hl->pixelStart = gs->x;
     hl->pixelEnd = gs->x + gs->width;
 
     hl->minVal = INT_MAX;
     hl->minValPos = -1;
     hl->maxVal = INT_MIN;
     hl->maxValPos = -1;
 
     struct slName *sl;
     double count = 0.0;
     int gsStart = hl->pixelStart;
 
     struct hmElement *hEl, *heList = NULL;
 
     for (sl = gs->genes; sl; sl = sl->next)
 	{
 	int numProbes = 0;
         struct hashEl *el = hashLookup(geneHash, sl->name);
         while (el)
             {
 	    el = hashLookupNext(el);
 	    numProbes++;
 	    }
 
 	int geneStart = (int) ( gs->pixelsPerGene * count + gsStart );
 	if (numProbes == 0)
 	    {
 	    int pixStart = geneStart;
 	    int pixEnd = pixStart + ceil(gs->pixelsPerGene);
 	    
 	    AllocVar(hEl);	    
 	    hEl->name = cloneString(sl->name);
 	    hEl->pixelStart = pixStart;
 	    hEl->pixelEnd = pixEnd;
 	    
 	    hEl->baseStart = 0;
 	    hEl->baseEnd = 0;
 	    
 	    hEl->shortLabel = NULL;
 	    hEl->longLabel = NULL;
 	    hEl->probeName = cloneString("N/A");
 	    hEl->bed = NULL;
 	    
 	    slAddHead(&heList, hEl);
 	    }
 	else 
 	    {
 	    el = hashLookup(geneHash, sl->name);
 	    double probeCount = 0.0;
 	    double pixelsPerProbe = gs->pixelsPerGene / (double) numProbes;
 	    while (el)
 		{
 		struct bed *nb = el->val;
 		el = hashLookupNext(el);
 		
 		int pixStart, pixEnd;
 		if (isSNP)
 		    {
 		    pixStart = geneStart;
 		    pixEnd = pixStart + ceil(gs->pixelsPerGene);
 		    }
 		else
 		    {
 		    pixStart = (int) ( pixelsPerProbe * probeCount + geneStart );
 		    pixEnd = pixStart + ceil(pixelsPerProbe);
 		    }
 		AllocVar(hEl);
 		
 		hEl->name = cloneString(sl->name);
 		hEl->pixelStart = pixStart;
 		hEl->pixelEnd = pixEnd;
 
 		hEl->baseStart = nb->chromStart;
 		hEl->baseEnd = nb->chromEnd; 
 		
 		hEl->shortLabel = NULL;
 		hEl->longLabel = NULL;
 		hEl->probeName = cloneString(nb->name);
 		hEl->bed = nb;
 		
 		slAddHead(&heList, hEl);
 		probeCount += 1.0;
 		}
 	    }
 	count += 1.0;
 	}
 
     slReverse(&heList);
     hl->elements = heList;
     slAddHead(&hlList, hl);
     }
 
 slReverse(&hlList);
 
 if (gh)
     heatmapLaySetDefaultColor(hlList, gh->platform);
 
 return hlList;
 }  
 
 struct heatmapLay *featureLayout(struct genoHeatmap *gh, int width, int sampleHeight)
 /* Layout Feature Sorter based on previously layed out heatmaps */
 {
 struct sqlConnection *conn = getFeatureDbConn(gh);
 if (!conn)
     return NULL;
 
 char *raName = gh->raFile;
 struct column *col, *colList = getColumns(conn, raName, gh->patDb);
 hFreeConn(&conn);
 
 int numVis = 0;
 for (col = colList; col != NULL; col = col->next)
     if (col->on)
         numVis++;
 
 int availablePixels = width;
 double basesPerPixel = (double) numVis / availablePixels;
 
 struct heatmapLay *hm = AllocA(struct heatmapLay);
 struct hmElement *hEl, *heList = NULL;
 
 hm->name = cloneString("featureInfo");
 hm->height = 0;
 hm->sampleHeight = sampleHeight;
 hm->width = width;
 hm->basesPerPixel = basesPerPixel;
 
 trackLayoutInit(&tl, cart);
 hm->font = tl.font;
 hm->fontHeight = tl.fontHeight;
     
 double pixelsPerBase = 0.0;
 if (basesPerPixel > 0)
     pixelsPerBase = 1.0/basesPerPixel;
 
 int colIx = 0;
 for (col = colList; col != NULL; col = col->next)
     {
     if (!col->on)
 	continue;
 
     AllocVar(hEl);
     hEl->name = cloneString(col->name);
     hEl->shortLabel = cloneString(col->shortLabel);
     hEl->longLabel = cloneString(col->longLabel);
 
     hEl->pixelStart = round(pixelsPerBase * colIx);
     hEl->pixelEnd = round(pixelsPerBase * (colIx + 1));
 
     hEl->baseStart = -1;
     hEl->baseEnd = -1;
     hEl->probeName = NULL;
 
     colIx++;
     slAddHead(&heList, hEl);
     }
 slReverse(&heList);
 
 hm->elements = heList;
 return hm;
 } 
 
 void geneSetsAddPathwayLayout(struct geneSet *geneSets, struct hash *geneHash)
 {
 /* lay out of all pathways, add setting the gs->pixelsPerGene gs->activeNum for each gene */
 if (geneSets == NULL || geneHash == NULL)
     return;
 
 struct geneSet *gs;
 for (gs = geneSets; gs; gs = gs->next)
     {
     gs->numGenesActive = slCount(gs->genes);//numActive;
     gs->pixelsPerGene = (double) gs->width / (double) gs->numGenesActive;
     }
 }
 
 void geneSetsChromLayout(struct geneSet *geneSets)
 /* lay out of all pathways similar to chromsomes,
    without setting the gs->pixelsPerGene gs->activeNum for each gene */
 {
 if (geneSets == NULL)
     return;
 
 int totalGenes = 0;
 struct geneSet *gs;
 int numSets = 0;
 
 for (gs = geneSets; gs; gs = gs->next)
     {
     totalGenes += gs->numGenes;
     numSets++;
     }
 
 if (totalGenes == 0)
     return;
 
 int picWidth = cartUsualInt(cart, hghImageWidth, hgHeatmapDefaultPixWidth);
 double xOff = 0;
 int minWidth = 20;
 
 int totalW = (double) (picWidth - numSets);
 for (gs = geneSets; gs; gs = gs->next)
     {
     int w = gs->numGenes * (picWidth - numSets) / totalGenes;
     if (w < minWidth) 
 	{
         totalGenes -= gs->numGenes;
         totalW -= minWidth;
         }
     }
 
 /* 5 is hard-coded , it works well*/
 double pixelsPerGene = (totalW - 5 ) / (double) totalGenes;
 for (gs = geneSets; gs; gs = gs->next)
     {
     gs->x = (int) xOff;
     gs->y = 0;
     gs->width = (int) (gs->numGenes * pixelsPerGene);
 
     if (gs->width < minWidth)
         gs->width = minWidth;
     xOff += (double) gs->width + 1.0;
     }
 }
 
 
 void drawGeneSetSingleStat(struct vGfx *vg, struct hash *statHash, struct heatmapLay *hl,
                            double gMin, double gMax, float colorCutoff,
                            char *chromHeatmap)
 {
 if (!statHash | !hl)
     return;
 
 Color valCol;
 double height = hghBed5Height;
 double baseline = (float) height /2.0;
 double gScale = baseline;
 if (gMax != gMin)
     gScale /= (gMax - gMin);
 
 vgBox(vg, hl->pixelStart, baseline-1, hl->width, 1, MG_GRAY);
 
 int leftX = hl->pixelStart - 1;
 if (leftX < 0)
     leftX = 0;
 int rightX = hl->pixelStart + hl->width;
 
 vgBox(vg, leftX, 0, 1, height, MG_GRAY);
 vgBox(vg, rightX, 0, 1, height, MG_GRAY);
 
 float val, dirVal;
 int x;
 int y1,y2;
 int w, h;
 
 struct hmElement *hEl;
 for (hEl = hl->elements; hEl ; hEl = hEl->next)
     {
     struct hgStats *nb = NULL;
     struct hashEl *el;
 
     el = hashLookup(statHash, hEl->name);
 
     while (el)                                                                                
 	{  // TODO:  do this in a smarter way.
 	nb = el->val;
 	if (sameString(nb->name, hEl->probeName))
 	    break;
 	nb = NULL;
 	el = hashLookupNext(el);
 	}
     if (!nb)
 	continue;
 
     nb = el->val;
     x = hEl->pixelStart; // get the X position
     w = hEl->pixelEnd - hEl->pixelStart;
 
     val = fabs(nb->prob);
     dirVal = nb->stats;
     
     if (val>gMax)
 	val = gMax;
     if (val <gMin)
 	val = gMin;
     
     if (val < colorCutoff)
 	valCol = MG_GRAY;
     else if (dirVal > 0)
 	valCol = MG_RED;
     else
 	valCol = MG_GREEN;
 
     double tmpV;
     double tmpH;
     
     if (dirVal>0)
 	{
 	y1 = baseline  - val *gScale; //smaller value
 	y2 = baseline; //greater value
 	h = y2-y1;
 	
 	/* sometimes if the bed segment is too small, resulting in width=0, 
 	     * reset to 1 pixel length */
 	if ((w ==0) && (h != 0))
 	    w=1;
 	if (w && h)
 	    vgBox(vg,x,y1, w, h,valCol);
 
 	tmpV = val;
 	tmpH = y1;
 	}
     else
 	{
 	y1 = baseline + val *gScale; //greater value
 	y2 = baseline; //smaller value
 	h = y1-y2;
 	/* sometimes if the bed segment is too small, resulting in width=0, 
 	 * reset to 1 pixel length */
 	if ((w ==0) && (h != 0))
 	    w=1;
 	if (w && h)
 	    vgBox(vg,x,y2,w,h,valCol);
 
 	tmpV = -1.0 * val;
 	tmpH = y2 + h;
 	}
 
     /* For scaling */
     if (tmpV > hl->maxVal)
 	{
 	hl->maxVal = tmpV;
 	hl->maxValPos = tmpH;
 	}
     if (tmpV < hl->minVal)
 	{
 	hl->minVal = tmpV;
 	hl->minValPos = tmpH;
 	}
 
     el = hashLookupNext(el);
     }
 }
 
 /* Draw hgStats results along gene sets  */
 void drawGeneSetStats(struct vGfx *vg, struct hash *statHash,
 		      struct heatmapLay *hlList, double gMin, double gMax, float colorCutoff,
 		      char *tableName)
 /* Drawing code */
 {
 if (!statHash)
     errAbort("No hgStats results");
 
 struct heatmapLay *hl;
 for (hl = hlList; hl ; hl = hl->next)
     drawGeneSetSingleStat(vg, statHash, hl, gMin, gMax, colorCutoff, tableName);
 }
 
 
 void drawChromStats(struct vGfx *vg, struct hgStats *bg, struct heatmapLay *hl,
 		    double gMin, double gMax, double colorCutoff)
 {
 if (!bg)
     errAbort("No hgStats data");
 
 double height = hghBed5Height;
 double baseline = (float) height / 2.0;
 vgBox(vg, 0, baseline-1, hl->width, 1, MG_GRAY);
 
 /* Chromosome layout hashs */
 struct hmElement *hEl = NULL;
 struct hash *chromHashX = newHash(0);
 struct hash *pixelHashX = newHash(0);
 
 for(hEl = hl->elements; hEl; hEl = hEl->next)
     {
     int chromX = hEl->baseStart;
     int pixelX = hEl->pixelStart;
     char *name = hEl->name;
     hashAddInt(chromHashX, name, chromX);
     hashAddInt(pixelHashX, name, pixelX);
 
     int leftX = hEl->pixelStart - 1;
     if (leftX < 0)
         leftX = 0; 
     int rightX = hEl->pixelEnd;
     if (rightX > hl->width - 1)
         rightX = hl->width - 1;
 
     vgBox(vg, leftX, 0, 1, height, MG_GRAY);
     vgBox(vg, rightX, 0, 1, height, MG_GRAY);
     }
 
 double pixelsPerBase = 1.0/hl->basesPerPixel;
 Color valCol;
 double gScale = height / 2.0;
 if (gMax != gMin)
     gScale /= (gMax - gMin);
 
 int start,end;
 char *chromName;
 float val, dirVal;
 int x1,x2;
 int y1;
 int w, h;
 
 char pixelStr[128];
 struct hash *pixelHash = hashNew(0);
 struct hmPixel *hm = NULL;
 
 struct hgStats *ibg = NULL;
 for(ibg = bg; ibg ; ibg = ibg->next)
     {
     chromName = ibg->chrom;
     start = ibg->chromStart + 1;
     end = ibg->chromEnd + 1;
 
     if (!hashLookup(chromHashX, chromName))
         continue;
 
     int chromX = hashIntVal(chromHashX, chromName);
     int pixelX = hashIntVal(pixelHashX, chromName);
 
     x1 = pixelsPerBase*(start - chromX) + pixelX;
     x2 = pixelsPerBase*(end - chromX) + pixelX;
 
     w = x2-x1;
 
     val = fabs(ibg->prob); // probability score control the absolute height of the drawing code
     if (val > gMax)
         val = gMax;
     if (val < gMin)
         val = gMin;
 
     dirVal = ibg->stats;  // statisitcs control drawing above or below baseline and color
     if (dirVal > 0)
         y1 = baseline - val *gScale; //smaller value
     else
         y1 = baseline; //smaller value
 
     h = val*gScale;
 
     double tmpV = val;
     double tmpH = y1;
     if (dirVal < 0)
 	{
 	tmpV = -1.0 * tmpV;
 	tmpH = tmpH + h;
 	}
     if (tmpV > hl->maxVal)
 	{
 	hl->maxVal = tmpV;
 	hl->maxValPos = tmpH;
 	}
     if (tmpV < hl->minVal)
 	{
 	hl->minVal = tmpV;
 	hl->minValPos = tmpH;
 	}
 
     if (dirVal > 0)
 	safef(pixelStr,sizeof (pixelStr),"%d,%d+",x1,x2);
     else
 	safef(pixelStr,sizeof (pixelStr),"%d,%d-",x1,x2);
 
     struct hashEl *el = hashLookup(pixelHash, pixelStr);
     if (!el)
         {
         hm = AllocA(struct hmPixel);
         hm->x = x1;
         hm->y = y1;
         hm->w = w;
         hm->h = h;
 	hm->val = dirVal;
         hashAdd(pixelHash, pixelStr, hm);
         }
     else
         {
         hm = el->val;
         //TO DOS: the current implementation is taking the most significant value, 
 	// a more sensible way to do it probably should be taking the average. 
 	// Implement in the future.
 	
 	if (hm->h <h)
 	    {
 	    hm->y = y1;
 	    hm->w = w;
 	    hm->h = h;
 	    hm->val = dirVal;
 	    }
         }
     }
 
 struct hashEl *elList = hashElListHash(pixelHash);
 struct hashEl *el;
 for (el = elList; el; el=el->next)
     {
     hm = el->val;
     w = hm->w;
     h = hm->h;
     x1 = hm->x;
     y1 = hm->y;
     val = hm->val;
     
     /* sometimes if the bed segment is too small, resulting in width=0, 
      * reset to 1 pixel length
      */
     if ((w ==0) && (h != 0))
 	w = 1;
     
     if (w < 1 && h < 1)
 	continue;
     
     if ( h < colorCutoff*gScale)
 	valCol= MG_GRAY;
     else if (val > 0)
 	valCol = MG_RED;
     else
 	valCol = MG_GREEN;
     
     vgBox(vg, x1, y1, w, h, valCol);
     }
 hashElFreeList(&elList);
 }
 
 
 
 void drawSubgroups(struct vGfx *vg, char *chromHeatmap, int width, 
 		   int height, int sampleHeight)
 /* Draw features to right of chromosomes in layout at given
  * y offset and height. */
 {
 if (!chromHeatmap)
     return;
 
 struct hashEl *el = hashLookup(ghHash, chromHeatmap);
 struct genoHeatmap *gh =NULL;
 if (el)
     gh = el->val;
 else
     errAbort("No heatmap %s", chromHeatmap);
 
 char *raName = gh->raFile;
 int subsetNum = cartUsualInt(cart, hgh2SubgroupNum, hgh2SubgroupDefaultNum);
 
 /* This is the key function to call for subgrouping */
 struct slName **ptSubsets = getSubsets(gh, subsetNum, raName );
 
 if (!ptSubsets)
     return;
 
 int subset = 0;
 Color valCol;
 char *sample;
 
 double pixelsPerSubgroup = (double) width / (double) subsetNum;
 
 for (subset = 0; subset < subsetNum; subset++)
     {
     if (!ptSubsets[subset])
 	continue;
 
     struct slName *subSampleList  = ptSubsets[subset];
 
     int x = round(pixelsPerSubgroup * subset);
     int w = ceil(pixelsPerSubgroup);
 
     vgSetClip(vg, x, 0, w, height);
     vgBox(vg, x, 0, w, height, MG_WHITE);
 
     struct slName *sl = NULL;
     for (sl = gh->sampleList; sl ; sl = sl->next)
         {
         sample = sl->name;
         int orderId = hashIntValDefault(gh->sampleOrder, sample, -1);
 
         if (orderId == -1)
             continue;
 
         if (slNameInList( subSampleList, sample))
             valCol = advFilterColor (subset);
         else
             continue;
 
         int h = sampleHeight;
         int y = orderId * h;
 
         vgBox(vg, x, y, w, h, valCol);
         }
     }
 
 for (subset=0; subset < subsetNum; subset++)
     if (ptSubsets[subset])
 	slNameFreeList(ptSubsets[subset]);
 free(ptSubsets);
 
 vgUnclip(vg);
 }
 
 
 int ifDrawFeatureLabel(struct heatmapLay *hl)
 {
 if (!hl)
     return 0;
 
 struct hmElement *hEl = hl->elements;
 if (!hEl)
     return 0;
 
 int width = hEl->pixelEnd - hEl->pixelStart;
 if (width < (hl->fontHeight * 2))
     return 0;
 
 // Feature width is wide enough to accomodate two lines of text.
 return 1;
 }
 
   
 void drawFeatures(struct vGfx *vg, char *chromHeatmap, struct heatmapLay *hl, int yOff, 
 		 boolean reverse)
 /* Draw features to right of chromosomes in layout at given
  * y offset and height. */
 {
 if (!hl)
     return;
 
 struct hashEl *el = hashLookup(ghHash, chromHeatmap);
 struct genoHeatmap *gh =NULL;
 if (el)
     gh = el->val;
 else
     errAbort("No heatmap %s", chromHeatmap);
 
 /* set up connection to table with both patient and sample information */
 char *labTable = gh->patTable;
 char *value = gh->sampleField;
 char *key = gh->patField;
 char *db = gh->patDb;
 
 if ((labTable == NULL) || (key == NULL) || (value == NULL) || (db==NULL))
     return;
 
 struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, db);
 
 if (!conn)
     return;
 
 double colorScale = 0.0;
 double val;
 double absVal,  minVal, maxVal;
 char *minCutVal, *maxCutVal;
 double offset;
 int reverseColor;
 char *colorScheme;
 
 Color valCol;
 Color upShades[EXPR_DATA_SHADES];
 Color downShades[EXPR_DATA_SHADES];
 
 static struct rgbColor black = {0, 0, 0};
 static struct rgbColor red = {255, 0, 0};
 static struct rgbColor green = {0, 255, 0};
 static struct rgbColor yellow = {255, 255, 0};
 
 struct slName *id = NULL;
 char *raName = gh->raFile;
 struct column *col, *colList = getColumns(conn, raName, gh->patDb);
 
 int sampleHeight = hl->sampleHeight;
 int height = heatmapHeight(gh) * sampleHeight;
 
 int nCols=0; //number of features 
 int nSamples=0; //number of samples
 struct hmElement *hEl=NULL;
 for (hEl = hl->elements; hEl; hEl = hEl->next)
     nCols++;
 struct slName *sl = NULL;
 for (sl = gh->sampleList; sl ; sl = sl->next)
     nSamples++;
 char ***featureArray= NULL; //an array with dim(nSamples , nCols+2)
 
 boolean output = FALSE; //TRUE; /* variable used to control generate numerical output of the features */
 if (output)
     {
     featureArray=  (char***)(malloc(sizeof (char**)*nSamples));
     int i,j;
     for (i=0; i<nSamples; i++)
 	{
 	featureArray[i]=(char**)(malloc(sizeof (char*)*(nCols+2)));
 	for (j=0; j<nCols+2; j++)
 	    featureArray[i][j]=NULL;
 	}
     }
 nCols=1;
 
 for (hEl = hl->elements; hEl; hEl = hEl->next)
     {
     int chromX = hEl->pixelStart, chromY = 0;
     int width = hEl->pixelEnd - hEl->pixelStart;
     int h = sampleHeight;
     nCols++;
 
     for (col = colList; col != NULL; col = col->next)
 	if (sameString(col->name, hEl->name))
 	    break;
     if (!col)
 	continue;
 
     colorScheme = col->cellColorSchemeVal(col);
     if sameWord(colorScheme,"redGreen")
       {
 	vgMakeColorGradient(vg, &black, &red, EXPR_DATA_SHADES, upShades);
 	vgMakeColorGradient(vg, &black, &green, EXPR_DATA_SHADES, downShades);
       }
     if sameWord(colorScheme,"default")
       {
 
 	vgMakeColorGradient(vg, &black, &yellow, EXPR_DATA_SHADES, upShades);
 	vgMakeColorGradient(vg, &black, &green, EXPR_DATA_SHADES, downShades);
       }
 
     if (ifDrawFeatureLabel(hl))
 	verticalTextRight(vg, chromX, height, width, hghFeatureLabel, 
 			     MG_BLACK, hl->font, col->shortLabel);
 
     vgSetClip(vg, chromX, chromY, width, height);
     vgBox(vg, chromX, chromY, width, height, MG_GRAY);
     
     offset = atof(col->cellOffsetVal(col));
     reverseColor = atoi(col->cellColorReverseVal(col));
     
     char *minValStr = col->cellMinVal(col, conn);
     if (!minValStr)
 	minVal = 0.0;
     else
 	minVal = atof(minValStr) + offset;
     
     char *maxValStr = col->cellMaxVal(col, conn);
 
     if (!maxValStr)
 	maxVal = 0.0;
     else 
 	maxVal = atof(maxValStr) + offset;
     
     if (reverseColor == -1)
 	{      
 	maxVal = reverseColor*maxVal;
 	minVal = reverseColor*maxVal;
 	}
 
     minCutVal = col->cellMinCutVal(col, conn);
     maxCutVal = col->cellMaxCutVal(col, conn);
     
     /* double color: consistent with the scale of microarray data display */
     if ((minVal<0) && (maxVal>0))
 	{
 	double md = max(fabs(minVal), maxVal);
 	colorScale = COLOR_SCALE  / md;
 	}
     else    /* single color */
 	colorScale = COLOR_SCALE / (maxVal - minVal);
 
     struct slName *sl = NULL;
     nSamples=0;
     for (sl = gh->sampleList; sl ; sl = sl->next)
 	{
 	int orderId = hashIntValDefault(gh->sampleOrder, sl->name, -1);
 
 	if (orderId == -1)
 	    continue;
 	id = slNameNew(getId(conn, labTable, key, sl->name, value));
       
 	char *cellVal = col->cellVal(col, id, conn);
 	if (output)
 	    {
 	    if (cellVal)
 		{
 		featureArray[nSamples][nCols]=(char*)(malloc(sizeof(int)*(strlen(cellVal)+1)));
 		strcpy(featureArray[nSamples][nCols],cellVal);
 		featureArray[nSamples][nCols][strlen(cellVal)]='\0';
 		}
 	    if (nCols==2)
 		{
 		featureArray[nSamples][0]=(char*)(malloc(sizeof(int)*8));
 		sprintf(featureArray[nSamples][0],"%d", orderId);
 
 		featureArray[nSamples][1]=(char*)(malloc(sizeof(int)*(strlen(sl->name)+1)));
 		strcpy(featureArray[nSamples][1],sl->name);
 		featureArray[nSamples][1][strlen(sl->name)]='\0';
 		}
 	    nSamples++;
 	    }
 
 	if (!cellVal)
 	    continue;
 
 	valCol = MG_GRAY;
 	val = atof(cellVal);
 	if (minCutVal)
 	    if (val < atof(minCutVal))
 		continue;
 	if (maxCutVal)
 	    if (val > atof(maxCutVal))
 		continue;
 	absVal = fabs(val);
 	absVal = reverseColor *( absVal +offset);
 	
 	/* double color: consistant with the scale of microarray data display */
 	int colorIndex ;
 	if ((minVal<0) && (maxVal>0))
 	    colorIndex = (int)((absVal) * (EXPR_DATA_SHADES-1.0) * colorScale);
 	else
 	    colorIndex = (int)((absVal-minVal) * (EXPR_DATA_SHADES-1.0) * colorScale);
 	
 	/* Clip color index to fit inside of array, since we may have brightened it. */
 	if (colorIndex < 0) colorIndex = 0;
 	if (colorIndex >= EXPR_DATA_SHADES)
 	    colorIndex = EXPR_DATA_SHADES-1;
 	
 	if (reverse)
 	    val = val * -1.0;
 	
 	if (val >= 0.0)
 	    valCol = upShades[colorIndex];
 	else
 	    valCol = downShades[colorIndex];
 	    
 	int y = chromY + orderId * h;
 	
 	vgBox(vg, chromX, y, width, h, valCol);
 	freez(&id);
 	}
     vgUnclip(vg);
     }
 
 /* write out feartureArray to output file, can be used for feature download function */
 if (output)
     {
     int i,j;
     // the feature value output file name in tmp directory , consider change the file name
     FILE* fout=fopen("/tmp/featureOutput.txt","w");
     for (i=0; i<nSamples; i++)
 	{
 	for (j=0; j<nCols+1; j++)
 	    {
 	    if (featureArray[i][j])
 		fprintf(fout,"%s",featureArray[i][j]);
 	    fprintf(fout,"\t");
 	    }
 	fprintf(fout,"\n");
 	}
     fclose(fout);
     }
 
 /* clean up memory of featureArray */
 if (output) 
     {
     int i;
     for (i=0; i<nSamples; i++)
 	free(featureArray[i]);
     free(featureArray);
     }
 
 hFreeConn(&conn);
 }
 
 int hmPixelCmpVal(const void *va, const void *vb)
 /* Compare to sort columns based on priority. */
 {
 const struct hmPixel *a = *((struct hmPixel **)va);
 const struct hmPixel *b = *((struct hmPixel **)vb);
 float dif = a->val/a->count - b->val/b->count;
 if (dif < 0)
     return -1;
 else if (dif > 0)
     return 1;
 else
     return 0;
 } 
 
 void drawFeatureSubgroupBar(struct vGfx *vg, struct slName *sampleList, int sampleHeight,
-			    boolean isBottomSubgroup)
+			    boolean subset)
 {
 if (!sampleList)
     return;
 
 int height = hgHeatmapDefaultSummaryHeight; 
 int width = hghSubgroupDefaultPixWidth / 2;
 vgSetClip(vg, 0, 0, hghSubgroupDefaultPixWidth, height);
-if (isBottomSubgroup) // only draw labels on bottom/green subgroup
+if (subset ==0) // only draw labels on bottom/green subgroup
     vgBox(vg, width, 0, width, height, MG_GREEN);
-else
+else if (subset ==1)
     vgBox(vg, 0, 0, width, height, MG_RED);
 vgUnclip(vg);
 }
 
 void drawFeatureSummary(struct vGfx *vg, char *chromHeatmap, struct heatmapLay *hl, int yOff, 
 			boolean reverse, struct slName *sampleList, boolean isBottomSubgroup)
 /* Draw feature summary to right of chromosomes in layout at given
  * y offset and height... actually a feature heatmap where every column is sorted ascending
  * breaking all relationships between horizontal rows
  */
 {
 if (!hl)
     return;
 
 struct hashEl *el = hashLookup(ghHash, chromHeatmap);
 struct genoHeatmap *gh =NULL;
 if (el)
     gh = el->val;
 else
     errAbort("No heatmap %s", chromHeatmap);
 
 /* set up connection to table with both patient and sample information */
 char *labTable = gh->patTable;
 char *value = gh->sampleField;
 char *key = gh->patField;
 char *db = gh->patDb;
 
 struct slName *sl;
 struct slInt *si, *siList = NULL;
 for (sl = sampleList; sl; sl = sl->next)
     {
     int i = slNameFindIx(gh->sampleList, sl->name);
     if (i < 0)
 	continue;
     si = slIntNew(i);
     slAddHead(&siList, si);
     }
 slReverse(&siList);
 
 if ((labTable == NULL) || (key == NULL) || (value == NULL) || (db==NULL))
     return;
 
 struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, db);
 
 if (!conn)
     return;
 
 double colorScale = 0.0;
 double val;
 double absVal,  minVal, maxVal;
 char *minCutVal, *maxCutVal;
 double offset;
 int reverseColor;
 char *colorScheme;
 
 Color valCol;
 Color upShades[EXPR_DATA_SHADES];
 Color downShades[EXPR_DATA_SHADES];
 
 static struct rgbColor black = {0, 0, 0};
 static struct rgbColor red = {255, 0, 0};
 static struct rgbColor green = {0, 255, 0};
 static struct rgbColor yellow = {255, 255, 0};
 
 struct slName *id = NULL;
 char *raName = gh->raFile;
 struct column *col, *colList = getColumns(conn, raName, gh->patDb);
 
 int height = hl->height;
 int buffer = 1;
 
 double h = 1.0;
 if (slCount(sampleList))
     h = (double) hl->height / (double) slCount(sampleList);
 
 struct hmElement *hEl;
 for (hEl = hl->elements; hEl; hEl = hEl->next)
     {
     int chromX = hEl->pixelStart;
     int width = hEl->pixelEnd - hEl->pixelStart;
     
     for (col = colList; col != NULL; col = col->next)
 	if (sameString(col->name, hEl->name))
 	    break;
     if (!col)
 	continue; 
 
     colorScheme = col->cellColorSchemeVal(col);
     if sameWord(colorScheme,"redGreen")
       {
 	vgMakeColorGradient(vg, &black, &red, EXPR_DATA_SHADES, upShades);
 	vgMakeColorGradient(vg, &black, &green, EXPR_DATA_SHADES, downShades);
       }
     if sameWord(colorScheme,"default")
       {
 
 	vgMakeColorGradient(vg, &black, &yellow, EXPR_DATA_SHADES, upShades);
 	vgMakeColorGradient(vg, &black, &green, EXPR_DATA_SHADES, downShades);
       }
 
     if (ifDrawFeatureLabel(hl) && isBottomSubgroup)
 	verticalTextRight(vg, chromX, height, width, hghFeatureLabel, 
 			     MG_BLACK, hl->font, col->shortLabel);
 
     /* To differentiate between normal feature heatmaps and summary heatmaps 
      * keep a one-pixel white buffer between each feature column, done by 
      * reducing width of clipped area */
     vgSetClip(vg, chromX-buffer, 0, width-buffer, height);
     vgBox(vg, chromX, 0, width, height, MG_GRAY);
 
     offset = atof(col->cellOffsetVal(col));
     reverseColor = atoi(col->cellColorReverseVal(col));
     
     char *minValStr = col->cellMinVal(col, conn);
     if (!minValStr)
 	minVal = 0.0;
     else
 	minVal = atof(minValStr) + offset;
     
     char *maxValStr = col->cellMaxVal(col, conn);
     if (!maxValStr)
 	maxVal = 0.0;
     else 
 	maxVal = atof(maxValStr) + offset;
     
     if (reverseColor == -1)
 	{      
 	maxVal = reverseColor*maxVal;
 	minVal = reverseColor*maxVal;
 	}
 
     minCutVal = col->cellMinCutVal(col, conn);
     maxCutVal = col->cellMaxCutVal(col, conn);
     
     /* double color: consistent with the scale of microarray data display */
     if ((minVal<0) && (maxVal>0))
 	{
 	double md = max(fabs(minVal), maxVal);
 	colorScale = COLOR_SCALE  / md;
 	}
     else    /* single color */
 	colorScale = COLOR_SCALE / (maxVal - minVal);
     
     struct slName *sl;
     struct hmPixel *hm, *hmList = NULL;
     for (sl = sampleList; sl ; sl = sl->next)
 	{
 	int orderId = hashIntValDefault(gh->sampleOrder, sl->name, -1);
 	
 	if (orderId == -1)
 	    continue;
 
 	int y = round((double) orderId * h);
 	id = slNameNew(getId(conn, labTable, key, sl->name, value));
 	
 	char *cellVal = col->cellVal(col, id, conn);
 		
 	if (!cellVal)
 	    continue;
 	val = atof(cellVal);
 
 	if (minCutVal)
 	    if (val < atof(minCutVal))
 		continue;
 	if (maxCutVal)
 	    if (val > atof(maxCutVal))
 		continue;
 
 	absVal = fabs(val);
 	absVal = reverseColor *( absVal +offset);
 	
 	/* double color: consistant with the scale of microarray data display */
 	int colorIndex ;
 	if ((minVal<0) && (maxVal>0))
 	    colorIndex = (int)((absVal) * (EXPR_DATA_SHADES-1.0) * colorScale);
 	else
 	    colorIndex = (int)((absVal-minVal) * (EXPR_DATA_SHADES-1.0) * colorScale);
 	
 	/* Clip color index to fit inside of array, since we may have brightened it. */
 	if (colorIndex < 0) colorIndex = 0;
 	if (colorIndex >= EXPR_DATA_SHADES)
 	    colorIndex = EXPR_DATA_SHADES-1;
 	
 	if (reverse)
 	    val = val * -1.0;
 	
 	if (val >= 0.0)
 	    valCol = upShades[colorIndex];
 	else
 	    valCol = downShades[colorIndex];
 
 	AllocVar(hm);
 	hm->val = val;
 	hm->count = 1;
 	hm->x = chromX;
 	hm->y = y;
 	hm->w = width;
 	hm->h = ceil(h);
 	hm->color = valCol;
 	slAddHead(&hmList, hm);
 	freez(&id);
 	}
 
     /* Sort column by value */
     slSort(&hmList, hmPixelCmpVal);
     float yOff = 0;    
 
     /* Draw sorted data */
     for (hm = hmList; hm; hm = hm->next)
 	{
 	vgBox(vg, hm->x, ceil(yOff), hm->w, hm->h, hm->color);
 	yOff +=h;
 	}
     vgUnclip(vg);
     slFreeList(&hmList);
     }
 hFreeConn(&conn);
 }
 
 
 void drawGeneSetHeatmapsByPixel(struct vGfx *vg, char* database, 
 				char *chromHeatmap, struct heatmapLay *hlList) 
 /* Draw chromosome graph on all chromosomes in layout at given
  * y offset and height. */
 {
 if (!hlList)
     return;
 
 struct hashEl *el = hashLookup(ghHash, chromHeatmap);
 struct genoHeatmap *gh =NULL;
 if (el)
     gh = el->val;
 else
     errAbort("No heatmap %s", chromHeatmap);
 
 int *chromOrder = getBedOrder(gh);
 struct bed *nb = NULL;
 
 float gain = gh->gainSet;
 double md = maxDeviation(chromHeatmap);
 double colorScale = COLOR_SCALE / md;
 double val;
 double absVal;
 int valId;
 Color valCol;
 
 Color upShades[EXPR_DATA_SHADES];
 Color downShades[EXPR_DATA_SHADES];
 
 vgMakeColorGradient(vg, hlList->zeroColor, hlList->highColor, EXPR_DATA_SHADES, upShades);
 vgMakeColorGradient(vg, hlList->zeroColor, hlList->lowColor, EXPR_DATA_SHADES, downShades);
 
 int sampleHeight = hlList->sampleHeight;
 int height = heatmapHeight(gh) * sampleHeight;
 int maxWidth = 0;
 struct heatmapLay *hl;
 for (hl = hlList; hl ; hl = hl->next)
     maxWidth = hl->pixelEnd;
 
 Color zeroColor = vgFindColorIx(vg, hlList->zeroColor->r, 
 			       hlList->zeroColor->g, hlList->zeroColor->b);
 
 boolean isSNP = FALSE;
 if (sameString(gh->platform, "SNP"))
     isSNP = TRUE;
 
 char pixelStr[128];
 struct hash *pixelHash = hashNew(0);
 struct hmPixel *hm, *hmList = NULL;
 
 for (hl = hlList; hl ; hl = hl->next)
     {
     int width = hl->pixelEnd - hl->pixelStart;
 
     vgSetClip(vg, hl->pixelStart, 0, width, height);
     vgBox(vg, hl->pixelStart, 0, maxWidth, height, zeroColor);
 
     struct hmElement *hEl;
 
     for (hEl = hl->elements; hEl ; hEl = hEl->next)
         {
 	nb = hEl->bed;
 	
 	if (!nb)
 	    continue;
 
 	int x = hEl->pixelStart;
 	int w = hEl->pixelEnd - hEl->pixelStart;
 	int h = sampleHeight;
 
 	int i;
 	for(i = 0; i < nb->expCount; ++i)
 	    {
 	    val = nb->expScores[i];
 	    valId = nb->expIds[i];
 	    int orderId = chromOrder[valId];
 	    if (orderId == -1)
 		continue;
 
 	    int y = orderId * h;
 	    safef(pixelStr, sizeof(pixelStr), "%d,%d", x, y);
 
 	    struct hashEl *pixelEl = hashLookup(pixelHash, pixelStr);
 	    if (!pixelEl)
 		{
 		hm = AllocA(struct hmPixel);
 		hm->x = x;
 		hm->y = y;
 		hm->w = w;
 		hm->h = h;
 		hm->val = 0.0;
 		hm->count = 0;
 		slAddHead(&hmList, hm);
 		hashAdd(pixelHash, pixelStr, hm);
 		}
 	    else
 		hm = pixelEl->val;
 	    
 	    if (isSNP)
 		{
 		if (val > 0.0)
 		    {
 		    hm->val += 1.0;
 		    hm->count += 1;
 		    }
 		}
 	    else 
 		{
 		hm->val += val;
 		hm->count += 1;
 		}
 	    }
 	}
 
 
     for (hm = hmList; hm ; hm = hm->next)
 	{
 	val = hm->val / (double) hm->count;
 	
         if(val > 0)
             absVal = val;
         else
             absVal = -val;
 
         absVal = absVal *gain;
         int colorIndex = (int)(absVal * (EXPR_DATA_SHADES-1.0) * colorScale);
         /* Clip color index to fit inside of array, since we may have brightened it. */
         if (colorIndex < 0) 
 	    colorIndex = 0;
         if (colorIndex >= EXPR_DATA_SHADES)
             colorIndex = EXPR_DATA_SHADES-1;
 
         if(val >= 0)
             valCol = upShades[colorIndex];
         else
             valCol = downShades[colorIndex];
 
         vgBox(vg, hm->x, hm->y, hm->w, hm->h, valCol);
         }
     vgUnclip(vg);
     
     freeHash(&pixelHash);
     slFreeList(&hmList);
     hmList = NULL;
     pixelHash = hashNew(0);
     }
 freeHash(&pixelHash);
 }
 
 void drawTufteTukeyBox(struct vGfx *vg, struct heatmapLay *hl, 
 		       struct hmPixel *hm, double colorScale, 
 		       Color *upShades, Color *downShades, double absMax, int height)
 /* Draws a dense, high data-to-ink ratio Tukey Boxplot, conceived by Edward Tufte 
  * Dot for the median, and then lines connecting the quartile to the whisker */
 {
 double val, absVal;
 Color valCol;
 /* Line connecting low quartile to low whisker, colored according to their average */
 val = (hm->lowW + hm->lowQ)/2.0;
 absVal = fabs(val);
 
 int colorIndex = (int)(absVal * (EXPR_DATA_SHADES-1.0) * colorScale);    
 /* Clip color index to fit inside of array, since we may have brightened it. */
 if (colorIndex < 0) colorIndex = 0;
 if (colorIndex >= EXPR_DATA_SHADES)
     colorIndex = EXPR_DATA_SHADES-1;
 if(val > 0)
     valCol = upShades[colorIndex];
 else
     valCol = downShades[colorIndex];
 
 val = hm->lowQ/absMax * height;
 int lQ = height - round(val);  // Low Quartile
 
 val = hm->lowW/absMax * height;
 int lW = height - round(val);  // Low Whisker
 
 int h = lW - lQ;
 vgBox(vg, hm->x + hm->w/2, lQ, 1, h, valCol); 
 
 /* Line connecting upper quartile to upper whisker, colored according to their average */
 val = (hm->highW + hm->highQ)/2.0;
 absVal = fabs(val);
 
 colorIndex = (int)(absVal * (EXPR_DATA_SHADES-1.0) * colorScale);    
 /* Clip color index to fit inside of array, since we may have brightened it. */
 if (colorIndex < 0) colorIndex = 0;
 if (colorIndex >= EXPR_DATA_SHADES)
     colorIndex = EXPR_DATA_SHADES-1;
 if(val > 0)
     valCol = upShades[colorIndex];
 else
     valCol = downShades[colorIndex];
 
 val = hm->highQ/absMax * height;
 int hQ = height - round(val);   // High Quartile
 val = hm->highW/absMax * height;
 int hW = height - round(val);   // High Whisker 
 h = hQ - hW;
 vgBox(vg, hm->x + hm->w/2, hW, 1, h, valCol); 
 
 /* Draw a dot for the median */
 val = hm->median/absMax * height;
 int median = height - round(val);
 vgBox(vg, hm->x + hm->w/2, median, 1, 1, MG_BLACK); 
 
 /* For scaling, check high and low whisker points */
 if (hm->highW > hl->maxVal)
     { 
     hl->maxVal = hm->highW;
     hl->maxValPos = hW;
     }
 if (hm->lowW < hl->minVal)
     {
     hl->minVal = hm->lowW;
     hl->minValPos = lW;
     }
 }
 
 void drawStandardTukeyBox(struct vGfx *vg, struct heatmapLay *hl, 
 			  struct hmPixel *hm, double colorScale, 
 			  Color *upShades, Color *downShades, double absMax, int height)
 /* Standard Tukey Boxplot, a box from upper->lower quartile, lines to whiskers, and a line
  * for the median
  */
 {
 Color valCol;
 double val = hm->median/absMax * height;
 double absVal = fabs(hm->median);
 
 /* Median a line, with 1 pixel buffer*/
 int median = height - round(val);
 vgBox(vg, hm->x+1, median, hm->w-2, 1, MG_BLACK); //valCol); 
 
 /* Lower whisker line */
 val = (hm->lowW + hm->lowQ)/2.0;
 absVal = fabs(val);
 
 int colorIndex = (int)(absVal * (EXPR_DATA_SHADES-1.0) * colorScale);    
 /* Clip color index to fit inside of array, since we may have brightened it. */
 if (colorIndex < 0) colorIndex = 0;
 if (colorIndex >= EXPR_DATA_SHADES)
     colorIndex = EXPR_DATA_SHADES-1;
 if(val > 0)
     valCol = upShades[colorIndex];
 else
     valCol = downShades[colorIndex];
 
 val = hm->lowQ/absMax * height;
 int lQ = height - round(val);  // Low Quartile
 
 val = hm->lowW/absMax * height;
 int lW = height - round(val);  // Low Whisker
 
 int h = lW - lQ;
 vgBox(vg, hm->x + hm->w/2, lQ, 1, h, valCol); 
 vgBox(vg, hm->x+1, lW, hm->w-2, 1, valCol); 
 /* Upper whisker line*/
 val = (hm->highW + hm->highQ)/2.0;
 absVal = fabs(val);
 
 colorIndex = (int)(absVal * (EXPR_DATA_SHADES-1.0) * colorScale);    
 /* Clip color index to fit inside of array, since we may have brightened it. */
 if (colorIndex < 0) colorIndex = 0;
 if (colorIndex >= EXPR_DATA_SHADES)
     colorIndex = EXPR_DATA_SHADES-1;
 if(val > 0)
     valCol = upShades[colorIndex];
 else
     valCol = downShades[colorIndex];
 
 val = hm->highQ/absMax * height;
 int hQ = height - round(val);   // High Quartile
 val = hm->highW/absMax * height;
 int hW = height - round(val);   // High Whisker 
 
 h = hQ - hW;    
 vgBox(vg, hm->x + hm->w/2, hW, 1, h, valCol);
 vgBox(vg, hm->x+1, hW, hm->w-2, 1, valCol); 
 
 /* Upper/Lower Quartile box */
 vgBox(vg, hm->x+1, lQ, hm->w-2, 1, MG_BLACK); 
 vgBox(vg, hm->x+1, hQ, hm->w-2, 1, MG_BLACK);  
 vgBox(vg, hm->x+1, hQ, 1, lQ - hQ, MG_BLACK);
 vgBox(vg, hm->x + hm->w-2, hQ, 1, lQ - hQ, MG_BLACK);
 
 /* For scaling, check high and low whisker points */
 if (hm->highW > hl->maxVal)
     { 
     hl->maxVal = hm->highW;
     hl->maxValPos = hW;
     }
 if (hm->lowW < hl->minVal)
     {
     hl->minVal = hm->lowW;
     hl->minValPos = lW;
     }
 }
 
 double hmPixelComputeTukeySummary(struct genoHeatmap *gh, struct hmPixel *hmList)
 /* Summary of data -- for drawing Tukey Box plot, returns absolute max for plotting purposes */
 {
 float gain = gh->gainFull;
  
 struct hmPixel *hm;
 double percentile = 0.25;
 double absMax = INT_MIN;
 for (hm = hmList; hm; hm = hm->next)
     {
     if (!hm->valList)
 	continue;
     double count = (double) slCount(hm->valList);
    
     slSort(&hm->valList, slDoubleCmp);
     struct slDouble *sld = slElementFromIx(hm->valList, floor(percentile * count));
     hm->lowQ = sld->val * gain;
     sld = slElementFromIx(hm->valList, ceil((1.0-percentile) * (count - 1.0)));
     hm->highQ = sld->val * gain;
 
     hm->median = gain * slDoubleMedian(hm->valList); 
     
     double iqr = hm->highQ - hm->lowQ;
     double lowBound = hm->lowQ - 1.5*iqr;
     double highBound = hm->highQ + 1.5*iqr;
     double tmp;
     hm->lowW  = INT_MAX;
     hm->highW = INT_MIN;
     for (sld = hm->valList; sld; sld = sld->next)
 	{
 	tmp = sld->val*gain;
 	if (tmp > lowBound && tmp < hm->lowW)
 	    hm->lowW = tmp;
 
 	if (tmp < highBound && tmp > hm->highW)
 	    hm->highW = tmp;
 	}
 
     tmp = -1;
     if (hm->lowW != INT_MAX && hm->highW != INT_MIN)
 	tmp = max(fabs(hm->lowW), fabs(hm->highW));
     else if (hm->lowW != INT_MAX)
 	tmp = fabs(hm->lowW);
     else if (hm->highW != INT_MIN)
 	tmp = fabs(hm->highW);
     
     if (tmp > absMax)
 	absMax = tmp;
     }
 
 return absMax;
 }
 
 void drawTukeyBoxPlot(struct vGfx *vg, struct heatmapLay *hl, struct genoHeatmap *gh, 
 		      struct hmPixel *hmList, int width, int height)
 {
 double md = maxDeviation(gh->name);
 double colorScale = COLOR_SCALE / md;
 
 Color upShades[EXPR_DATA_SHADES];
 Color downShades[EXPR_DATA_SHADES];
 
 // default, expression data color scheme:
 vgMakeColorGradient(vg, hl->zeroColor, hl->highColor, EXPR_DATA_SHADES, upShades);
 vgMakeColorGradient(vg, hl->zeroColor, hl->lowColor, EXPR_DATA_SHADES, downShades);
 
 double absMax = hmPixelComputeTukeySummary(gh, hmList);
 
 struct hmPixel *hm;
 int midHeight = height/2;
 vgSetClip(vg, 0, 0, width, height);
 //vgBox(vg, 0, 0, width, height, MG_WHITE);
 for (hm = hmList; hm ; hm = hm->next)
     {
     if (!hm->valList)
 	continue;
     if (hm->w < 5)  // dense version of the box, courtesy of Edward Tufte.
 	drawTufteTukeyBox(vg, hl, hm, colorScale, upShades, downShades, absMax, midHeight);
     else // standard box plot
 	drawStandardTukeyBox(vg, hl, hm, colorScale, upShades, downShades, absMax, midHeight);
     }
 vgUnclip(vg);
 }
 
 
 void drawGenesetAverageLine(struct vGfx *vg, struct hmPixel *hmList, 
 			    double absMax, int height)
 /* Draws a bargraph, specific for SNP data (for now) */
 {
 double absVal = fabs(absMax);
 
 struct hmPixel *hm;
 for (hm = hmList; hm; hm = hm->next)
     {
     int h = round(hm->avg/absVal * (double) height );
     if (h == 0)
 	h = 1;
     vgBox(vg, hm->x, height - h, hm->w, 2, MG_GRAY); 
     }
 
 }
 
 void drawBarGraph(struct vGfx *vg, struct heatmapLay *hl, 
 		  struct hmPixel *hm, double colorScale, 
 		  Color *upShades, Color *downShades, double absMax, int height)
 /* Draws a bargraph, specific for SNP data (for now) */
 {
 Color valCol;
 double absVal = fabs(absMax);
 
 int colorIndex = EXPR_DATA_SHADES-1;
 if(hm->avg > 0)
     valCol = upShades[colorIndex];
 else
     valCol = downShades[colorIndex];
 
 int h = round(hm->avg/absVal * (double) height );
 if (h == 0)
     h = 1;
 vgBox(vg, hm->x, height - h, hm->w, h, valCol); 
 
 /* For scaling, check high and low whisker points */
 if (hm->avg > hl->maxVal)
     { 
     hl->maxVal = hm->avg;
     hl->maxValPos = height - h;
     }
 if (hm->avg < hl->minVal)
     {
     hl->minVal = 0.0;  // Always set minimum val to zero for bar graphs.
     hl->minValPos = height;
     }
 }
 
 double hmPixelComputeBarGraphSummary(struct genoHeatmap *gh, struct hmPixel *hmList)
 /* Summary of data -- Bar graph, currently specific for SNP bed15 files */
 {
 struct hmPixel *hm;
 double absMax = INT_MIN;
 for (hm = hmList; hm; hm = hm->next)
     {
     if (!hm->valList)
 	continue;
    
     double sum = 0.0;
     struct slDouble *sld;
     for (sld = hm->valList; sld; sld = sld->next)
 	sum += sld->val;
     hm->sum = sum;
     hm->avg = sum / (double) hm->count;
 
     if (hm->avg > absMax)
 	absMax = hm->avg;
     }
 if (absMax == 0.0)
     absMax = 0.000000001;
 return absMax;
 }
 
 void drawBarGraphPlot(struct vGfx *vg, struct heatmapLay *hl, struct genoHeatmap *gh, 
 		      struct hmPixel *hmList, int width, int height)
 { /* Specifically for SNP mutation data */
 double md = maxDeviation(gh->name);
 double colorScale = COLOR_SCALE / md;
 
 Color upShades[EXPR_DATA_SHADES];
 Color downShades[EXPR_DATA_SHADES];
 
 // default, expression data color scheme:
 vgMakeColorGradient(vg, hl->zeroColor, hl->highColor, EXPR_DATA_SHADES, upShades);
 vgMakeColorGradient(vg, hl->zeroColor, hl->lowColor, EXPR_DATA_SHADES, downShades);
 
 double absMax = hmPixelComputeBarGraphSummary(gh, hmList);
 
 struct hmPixel *hm;
 vgSetClip(vg, 0, 0, width, height);
 for (hm = hmList; hm ; hm = hm->next)
     {
     if (!hm->valList)
 	continue;
     drawBarGraph(vg, hl, hm, colorScale, upShades, downShades, absMax, height);
     }
 vgUnclip(vg);
 }
 
 
 void drawGeneSetSNPSummary(struct vGfx *vg, struct genoHeatmap *gh,
 			   struct heatmapLay *hlList, struct slName *sampleList, int height) 
 /* Draw chromosome graph on all chromosomes in layout at given
  * y offset and height. */
 {
 if (!hlList)
     return;
 
 struct bed *nb = NULL;
 
 double val;
 char pixelStr[128];
 struct hash *pixelHash = hashNew(0);
 struct hmPixel *hm, *hmList = NULL;
 struct hmPixel *hmGs, *hmGsList = NULL;
 
 struct slName *sl;
 struct slInt *si, *siList = NULL;
 for (sl = sampleList; sl; sl = sl->next)
     {
     int i = slNameFindIx(gh->sampleList, sl->name);
 
     struct slName *sample = slElementFromIx(gh->sampleList, i);
     if (i < 0)
 	continue;
 
     if (!slNameInList(sampleList, sample->name))
 	continue;
   
     si = slIntNew(i);
     slAddHead(&siList, si);
     }
 slReverse(&siList);
 
 int sampleCount = slCount(siList);
 
 int maxX = 0;
 struct heatmapLay *hl;
 
 for (hl = hlList; hl ; hl = hl->next)
     {
     struct hmElement *hEl;
 
     hmGs = AllocA(struct hmPixel);
     hmGs->x = hl->pixelStart;
     hmGs->y = 0;
     hmGs->w = hl->pixelEnd - hl->pixelStart;
     hmGs->h = 0;
     hmGs->val = 0.0;
     hmGs->valList = NULL;
     hmGs->count = 0;
     hmGs->avg = -1.0;
     slAddHead(&hmGsList, hmGs);
 
     struct hash *txLength = hashNew(0);
     for (hEl = hl->elements; hEl ; hEl = hEl->next)
         {
 	nb = hEl->bed;
 
 	if (!nb)
 	    continue;
 
 	int x = hEl->pixelStart;
 	int w = hEl->pixelEnd - hEl->pixelStart;
 
 	if (hEl->pixelEnd > maxX)
 	    maxX = hEl->pixelEnd;
 
 	safef(pixelStr, sizeof(pixelStr), "%d", x);
 	
 	struct hashEl *pixelEl = hashLookup(pixelHash, pixelStr);
 	if (!pixelEl)
 	    {
 	    hm = AllocA(struct hmPixel);
 	    hm->x = x;
 	    hm->y = 0;
 	    hm->w = w;
 	    hm->h = 0;
 	    hm->val = 0.0;
 	    hm->valList = NULL;
 	    hm->count = 0;
 	    slAddHead(&hmList, hm);
 	    hashAdd(pixelHash, pixelStr, hm);
 	    }
 	else
 	    hm = pixelEl->val;
 
 	struct hashEl *el = hashLookup(txLength, hEl->name);
 	if (!el)
 	    {
 	    /* we're using score to store transcript length for SNP data */
 	    int score = nb->score;
 	    if (score == 0)
 		score = 1;
 
 	    hashAddInt(txLength, hEl->name, score);	
 	    hmGs->count += score * sampleCount;
 	    hm->count += score * sampleCount;
 	    }
 
 	for(si = siList; si; si = si->next)
 	    {
 	    int i = si->val;
 	    val = nb->expScores[i];
 
 	    /* Count observed mutations */
 	    hmGs->val += val;
 
 	    struct slDouble *sld = slDoubleNew(val);
 	    slAddHead(&hm->valList, sld);
 	    }
 
 	}
     freeHash(&txLength);
 
     if (hmGs->count)
 	hmGs->avg = hmGs->val / (double) hmGs->count;
     }
 slReverse(&hmGsList);
 
 drawBarGraphPlot(vg, hlList, gh, hmList, maxX, height);
 drawGenesetAverageLine(vg, hmGsList, hlList->maxVal, height);
 }
 
 void drawGeneSetGenomicSummary(struct vGfx *vg, struct genoHeatmap *gh,
 			       struct heatmapLay *hlList, struct slName *sampleList, int height) 
 /* Draw chromosome graph on all chromosomes in layout at given
  * y offset and height. */
 {
 if (!hlList)
     return;
 
 struct bed *nb = NULL;
 
 double val;
 char pixelStr[128];
 struct hash *pixelHash = hashNew(0);
 struct hmPixel *hm, *hmList = NULL;
 
 struct slName *sl;
 struct slInt *si, *siList = NULL;
 for (sl = sampleList; sl; sl = sl->next)
     {
     int i = slNameFindIx(gh->sampleList, sl->name);
 
     struct slName *sample = slElementFromIx(gh->sampleList, i);
     if (i < 0)
 	continue;
 
     if (!slNameInList(sampleList, sample->name))
 	continue;
   
     si = slIntNew(i);
     slAddHead(&siList, si);
     }
 slReverse(&siList);
 
 int maxX = 0;
 struct heatmapLay *hl;
 
 for (hl = hlList; hl ; hl = hl->next)
     {
     struct hmElement *hEl;
 
     for (hEl = hl->elements; hEl ; hEl = hEl->next)
         {
 	nb = hEl->bed;
 	
 	if (!nb)
 	    continue;
 
 	int x = hEl->pixelStart;
 	int w = hEl->pixelEnd - hEl->pixelStart;
 
 	if (hEl->pixelEnd > maxX)
 	    maxX = hEl->pixelEnd;
 
 	safef(pixelStr, sizeof(pixelStr), "%d", x);
 	
 	struct hashEl *pixelEl = hashLookup(pixelHash, pixelStr);
 	if (!pixelEl)
 	    {
 	    hm = AllocA(struct hmPixel);
 	    hm->x = x;
 	    hm->y = 0;
 	    hm->w = w;
 	    hm->h = 0;
 	    hm->val = 0.0;
 	    hm->valList = NULL;
 	    hm->count = 0;
 	    slAddHead(&hmList, hm);
 	    hashAdd(pixelHash, pixelStr, hm);
 	    }
 	else
 	    hm = pixelEl->val;
 	    
 	for(si = siList; si; si = si->next)
 	    {
 	    int i = si->val;
 	    val = nb->expScores[i];
 
 	    struct slDouble *sld = slDoubleNew(val);
 	    slAddHead(&hm->valList, sld);
 	    }
 	}
     }
 
 drawTukeyBoxPlot(vg, hlList, gh, hmList, maxX, height);
 }
 
 void drawGeneSetSummary(struct vGfx *vg, char* database, char *chromHeatmap, 
 			struct heatmapLay *hlList, struct slName *sampleList, int height) 
 /* Draw chromosome graph on all chromosomes in layout at given
  * y offset and height. */
 {
 if (!hlList)
     return;
 
 struct hashEl *el = hashLookup(ghHash, chromHeatmap);
 struct genoHeatmap *gh =NULL;
 if (el)
     gh = el->val;
 else
     errAbort("No heatmap %s", chromHeatmap);
 
 if (sameString(gh->platform, "SNP"))
     drawGeneSetSNPSummary(vg, gh, hlList, sampleList, height);
 else
     drawGeneSetGenomicSummary(vg, gh, hlList, sampleList, height);
 }
 
 
 void drawBedGraph(struct vGfx *vg, char* database,
 		  struct heatmapLay *hl, char *tableName)
 {
 if (!hl)
     return;
 
 struct hashEl *el = hashLookup(ghHash, tableName);
 struct genoHeatmap *gh =NULL;
 if (el)
     gh = el->val;
 else
     errAbort("No heatmap %s", tableName);
 
 struct bed *bg=NULL, *ibg=NULL;
 float val;
 
 double pixelsPerBase = 1.0/hl->basesPerPixel;
 
 int x1,x2;
 int y1,y2;
 int w, h;
 int start,end;
 
 double gMin = chromGraphMin(tableName);
 double gMax = chromGraphMax(tableName);
 double height = hl->height;
 double gScale = height/(gMax-gMin);
 
 struct rgbColor color = chromGraphColor(tableName);
 struct rgbColor black = {0,0,0};
 Color valCol;
 
 double md = 10.0;
 double colorScale = COLOR_SCALE / md;
 int nField = 4;
 
 Color shadesOfColor[EXPR_DATA_SHADES];
 if (sameWord(gh->dataType, "bed 4"))
     vgMakeColorGradient(vg, &black, &color, 1, shadesOfColor);
 else if (sameWord(gh->dataType, "bed 5"))
     {
     vgMakeColorGradient(vg, &black, &color, EXPR_DATA_SHADES, shadesOfColor);
     nField = 5;
     }
 struct hmElement *hEl;  
 for(hEl = hl->elements; hEl; hEl = hEl->next)
     {
     bg = getBedGraph(gh, hEl->name, nField);
     if (bg== NULL)
         continue;
 
     int chromX = hEl->pixelStart, chromY = 0;
     int width = hEl->pixelEnd - hEl->pixelStart;
     vgSetClip(vg, chromX, chromY, width, height);
 
     /* Draw rest of points, connecting with line to previous point
      * if not too far off. */
     for(ibg = bg; ibg ; ibg = ibg->next)
         {
         start = ibg->chromStart + 1;
         end = ibg->chromEnd + 1;
 	val = sqlFloat(ibg->name);
 	
 	valCol = shadesOfColor[0];
 	if (sameWord(gh->dataType, "bed 5"))
 	    {
 	    int prob = ibg->score/md;
 	    int colorIndex = abs(prob) * (EXPR_DATA_SHADES-1.0) * colorScale;
 	    if (colorIndex < 0) 
 		colorIndex = 0;
 	    if (colorIndex >= EXPR_DATA_SHADES)
 		colorIndex = EXPR_DATA_SHADES-1;
 	    valCol = shadesOfColor[colorIndex];
 	    }
  
         /* within a boundary */
         if( val > gMax)
             val = gMax;
         else if (val < gMin)
             val = gMin;
 
         x1 = pixelsPerBase*start + chromX;
         x2 = pixelsPerBase*end + chromX;
         w = x2-x1;
 
         if (val>0)
             {
             y1 = (height - (val - gMin) * gScale) + chromY; //smaller value
             y2 = (height - (gMin -gMin) * gScale) + chromY;//greater value
             h = y2-y1;
 
             /* sometimes if the bed segment is too small, resulting in width=0, 
 	     * reset to 1 pixel length */
             if ((w ==0) && (h != 0))
                 w=1;
             vgBox(vg,x1,y1, w, h,valCol);
             }
         if (val<0)
             {
             y1 = (height - (val - gMax - gMin) *gScale) + chromY; //greater value
             y2 = (height - (gMax - gMax - gMin)*gScale) + chromY; //smaller value
             h = y1-y2;
             /* sometimes if the bed segment is too small, resulting in width=0, 
 	     *reset to 1 pixel length */
             if ((w ==0) && (h != 0))
                 w=1;
             vgBox(vg,x1,y2,w,h,valCol);
             }
         }
     vgUnclip(vg);
     bedFree(&bg);
     }
 }  
 
 
 void drawChromSNPSummary(struct vGfx *vg, struct genoHeatmap *gh, struct heatmapLay *hl, 
 			 struct slName *sampleList, int height)
 /* Draw chromosome graph on all chromosomes in layout at given
  * y offset and height. */
 {
 struct bed *ghBed=NULL, *nb=NULL;
 
 double pixelsPerBase = 1.0/hl->basesPerPixel;
 double val;
 int start, end;
 
 char pixelStr[128];
 struct hash *pixelHash = hashNew(0);
 struct hmPixel *hm, *hmList = NULL;
 
 struct slName *sl;
 struct slInt *si, *siList = NULL;
 for (sl = sampleList; sl; sl = sl->next)
     {
     int i = slNameFindIx(gh->sampleList, sl->name);
     if (i < 0)
 	continue;
     si = slIntNew(i);
     slAddHead(&siList, si);
     }
 slReverse(&siList);
 int sampleCount = slCount(siList);
 
 struct hmElement *hEl;
 for(hEl = hl->elements; hEl; hEl = hEl->next)
     {
     ghBed = getChromHeatmapRange(gh, hEl->name, hEl->baseStart, hEl->baseEnd, FALSE);
 
     int chromX = hEl->pixelStart; 
     for(nb = ghBed; nb; nb = nb->next)
 	{
 	start = nb->chromStart + 1;  // to fix off-by-one error
 	end = nb->chromEnd + 1;      // to fix off-by-one error
 
 	/* Window overlaps a probe */
 	if (start < hEl->baseStart) 
 	    start = hEl->baseStart;
 	if ( end > hEl->baseEnd) 
 	    end = hEl->baseEnd;
 	
 	int x = pixelsPerBase * (start - hEl->baseStart) + chromX;
 	int w = pixelsPerBase * (end - start);
 	if (w == 0)
 	    w = 1;
 
 	safef(pixelStr, sizeof(pixelStr), "%d", x);
 	struct hashEl *el = hashLookup(pixelHash, pixelStr);
 	if (!el)
 	    {
 	    hm = AllocA(struct hmPixel);
 	    hm->x = x;
 	    hm->y = 0;
 	    hm->w = w;
 	    hm->h = 0;
 	    hm->val = 0.0;
 	    hm->count = 0;
 	    hm->valList = NULL;
 	    slAddHead(&hmList, hm);
 	    hashAdd(pixelHash, pixelStr, hm);
 	    }
 	else
 	    hm = el->val;
 
 	/* the score field is the number of sites used to derive mutation per site
 	 * if you want to have mutation per site per sample, use the score field
 	 * if you do not want to normalized against the size of the gene, set score =1 or 0 */
 	int score = nb->score;
 	if (score == 0)
 	    score = 1;
 
 	hm->count += score * sampleCount;
 
 	for(si = siList; si; si = si->next)
 	    {
 	    int i = si->val;
 	    val = nb->expScores[i];
 	    
 	    struct slDouble *sld = slDoubleNew(val);
 	    slAddHead(&hm->valList, sld);
 	    }
 	}
     bedFreeList(&ghBed);
     }
 
 drawBarGraphPlot(vg, hl, gh, hmList, hl->width, height);
 }
 
 void drawChromGenomicSummary(struct vGfx *vg, struct genoHeatmap *gh, struct heatmapLay *hl, 
 			     struct slName *sampleList, int height)
 /* Draw chromosome graph on all chromosomes in layout at given
  * y offset and height. */
 {
 struct bed *ghBed=NULL, *nb=NULL;
 
 double pixelsPerBase = 1.0/hl->basesPerPixel;
 double val;
 int start, end;
 
 char pixelStr[128];
 struct hash *pixelHash = hashNew(0);
 struct hmPixel *hm, *hmList = NULL;
 
 struct slName *sl;
 struct slInt *si, *siList = NULL;
 for (sl = sampleList; sl; sl = sl->next)
     {
     int i = slNameFindIx(gh->sampleList, sl->name);
     if (i < 0)
 	continue;
     si = slIntNew(i);
     slAddHead(&siList, si);
     }
 slReverse(&siList);
 
 boolean useAccessTable = TRUE;
 if (singleChromImage(hl))
     useAccessTable = FALSE;
 
 struct hmElement *hEl;
 for(hEl = hl->elements; hEl; hEl = hEl->next)
     {
     ghBed = getChromHeatmapRange(gh, hEl->name, hEl->baseStart, hEl->baseEnd, useAccessTable);
 
     int chromX = hEl->pixelStart; 
     for(nb = ghBed; nb; nb = nb->next)
 	{
 	start = nb->chromStart + 1;  // to fix off-by-one error
 	end = nb->chromEnd + 1;      // to fix off-by-one error
 
 	/* Window overlaps a probe */
 	if (start < hEl->baseStart) 
 	    start = hEl->baseStart;
 	if ( end > hEl->baseEnd) 
 	    end = hEl->baseEnd;
 
 	int x = pixelsPerBase * (start - hEl->baseStart) + chromX;
 	int w = pixelsPerBase * (end - start);
 	if (w == 0)
 	    w = 1;
 
 	safef(pixelStr, sizeof(pixelStr), "%d", x);
 	struct hashEl *el = hashLookup(pixelHash, pixelStr);
 	if (!el)
 	    {
 	    hm = AllocA(struct hmPixel);
 	    hm->x = x;
 	    hm->y = 0;
 	    hm->w = w;
 	    hm->h = 0;
 	    hm->val = 0.0;
 	    hm->count = 0;
 	    hm->valList = NULL;
 	    slAddHead(&hmList, hm);
 	    hashAdd(pixelHash, pixelStr, hm);
 	    }
 	else
 	    hm = el->val;
 
 	for(si = siList; si; si = si->next)
 	    {
 	    int i = si->val;
 	    val = nb->expScores[i];
 
 	    struct slDouble *sld = slDoubleNew(val);
 	    slAddHead(&hm->valList, sld);
 	    }
 	}
     bedFreeList(&ghBed);
     }
 
 drawTukeyBoxPlot(vg, hl, gh, hmList, hl->width, height);
 }
 
 void drawChromSummary(struct vGfx *vg, char* database, struct heatmapLay *hl, 
 		      char *chromHeatmap, struct slName *sampleList, int height)
 /* Draw chromosome graph on all chromosomes in layout at given
  * y offset and height. */
 {
 if (!hl)
     return;
 
 struct hashEl *el = hashLookup(ghHash, chromHeatmap);
 struct genoHeatmap *gh =NULL;
 if (el)
     gh = el->val;
 else
     errAbort("No heatmap %s", chromHeatmap);
 
 if (sameString(gh->platform, "SNP"))
     drawChromSNPSummary(vg, gh, hl, sampleList, height);
 else    
     drawChromGenomicSummary(vg, gh, hl, sampleList, height);
 }
 
 void drawBackgroundLines(struct vGfx *vg, int startX, int startY, int width, int height)
 {
 int x, spacing = 20;
 
 Color lightBlue = vgFindColorIx(vg, 220, 220, 255);
 vgSetClip(vg, startX, startY, width, height);
 for (x = startX; x < (startX + width); x += spacing)
     vgBox(vg, x, startY, 1, height, lightBlue);
 vgUnclip(vg);
 }
 
 void drawChromHeatmapsByPixel(struct vGfx *vg, char* database, 
 			      struct heatmapLay *hl, char *chromHeatmap)
 /* Draw chromosome graph on all chromosomes in layout at given
  * y offset and height. */
 {
 if (!hl)
     return;
 
 struct hashEl *el = hashLookup(ghHash, chromHeatmap);
 struct genoHeatmap *gh =NULL;
 if (el)
     gh = el->val;
 else
     errAbort("No heatmap %s", chromHeatmap);
 
 int *chromOrder = getBedOrder(gh);
 
 struct bed *ghBed=NULL, *nb=NULL;
 
 
 double md = maxDeviation(gh->name);
 double colorScale = COLOR_SCALE / md;
 double val;
 double absVal;
 int valId;
 Color valCol;
 int start, end;
 
 float gain = gh->gainFull;
 double pixelsPerBase = 1.0/hl->basesPerPixel; 
 Color upShades[EXPR_DATA_SHADES];
 Color downShades[EXPR_DATA_SHADES];
 
 vgMakeColorGradient(vg, hl->zeroColor, hl->highColor, EXPR_DATA_SHADES, upShades);
 vgMakeColorGradient(vg, hl->zeroColor, hl->lowColor, EXPR_DATA_SHADES, downShades);
 
 char pixelStr[128];
 struct hash *pixelHash = hashNew(0);
 struct hmPixel *hm, *hmList = NULL;
 
 boolean isSNP = FALSE;
 if (sameString(gh->platform, "SNP"))
     isSNP = TRUE;
 
 boolean useAccessTable = TRUE;
 if (singleChromImage(hl))
     useAccessTable = FALSE;
 
 int sampleHeight = hl->sampleHeight; 
 int totalHeight = sampleHeight * heatmapHeight(gh);
 
 int buffer = PROBE_MARK_HEIGHT + PROBE_MARK_BUFFER;
 struct hmElement *hEl;
 for(hEl = hl->elements; hEl; hEl = hEl->next)
     {
     ghBed = getChromHeatmapRange(gh, hEl->name, hEl->baseStart, hEl->baseEnd, useAccessTable);
  
     int chromX = hEl->pixelStart, chromY = 0;
     int width = hEl->pixelEnd - hEl->pixelStart;
     vgSetClip(vg, chromX, chromY, width, totalHeight + buffer);
     Color color = vgFindColorIx(vg, hl->zeroColor->r, hl->zeroColor->g, hl->zeroColor->b);
     vgBox(vg, chromX, chromY, width, totalHeight, color);
 
     for(nb = ghBed; nb; nb = nb->next)
 	{
 	start = nb->chromStart + 1;   // to fix off-by-one error
 	end = nb->chromEnd + 1;       // to fix off-by-one error
 
 	/* Window overlaps a probe */
 	if (start < hEl->baseStart) 
 	    start = hEl->baseStart;
 	if ( end > hEl->baseEnd) 
 	    end = hEl->baseEnd;
 
 	int x = pixelsPerBase * (start - hEl->baseStart) + chromX;
 	int w = pixelsPerBase * (end - start);
 	if (w == 0)
 	    w = 1;
 	int i;
 	for(i = 0; i < nb->expCount; ++i)
 	    {
 	    val = nb->expScores[i];
 	    valId = nb->expIds[i];
 	    int orderId = chromOrder[valId];
 	    if (orderId == -1)
 		continue;
 
 	    int y = chromY + orderId * sampleHeight;
 
 	    safef(pixelStr, sizeof(pixelStr), "%d,%d", x, y);
 	    struct hashEl *el = hashLookup(pixelHash, pixelStr);
 	    if (!el)
 		{
 		hm = AllocA(struct hmPixel);
 		hm->x = x;
 		hm->y = y;
 		hm->w = w;
 		hm->h = sampleHeight;
 		hm->val = 0.0;
 		hm->count = 0;
 		slAddHead(&hmList, hm);
 		hashAdd(pixelHash, pixelStr, hm);
 		}
 	    else
 		hm = el->val;
 
 	    if (isSNP)
 		{
 		if (val > 0)
 		    {
 		    hm->val += 1.0;
 		    hm->count += 1;
 		    }
 		}
 	    else
 		{
 		hm->val += val;
 		hm->count += 1;
 		}
 	    }
       	}
 
     for (hm = hmList; hm ; hm = hm->next)
 	{
 	val = hm->val / hm->count;
 
 	absVal = fabs(val) * gain;
 	
 	int colorIndex = (int)(absVal * (EXPR_DATA_SHADES-1.0) * colorScale);
 
 	/* Clip color index to fit inside of array, since we may have brightened it. */
 	if (colorIndex < 0) colorIndex = 0;
 	if (colorIndex >= EXPR_DATA_SHADES)
 	    colorIndex = EXPR_DATA_SHADES-1;
 
 	if(val > 0)
 	    valCol = upShades[colorIndex];
 	else
 	    valCol = downShades[colorIndex];
 
 	/* Draw probe data */
 	vgBox(vg, hm->x, hm->y, hm->w, hm->h, valCol); 
 
 	/* Draw probe mark at bottom */
 	vgBox(vg, hm->x, hl->height+PROBE_MARK_BUFFER, 
 	      hm->w, PROBE_MARK_HEIGHT, MG_BLACK); 
 	}
 
     freeHash(&pixelHash);
     pixelHash = hashNew(0);
     slFreeList(&hmList);
     hmList = NULL;
 
     bedFreeList(&ghBed);
     vgUnclip(vg);
     }
 }
 
 char *genomeGif(struct sqlConnection *conn, struct heatmapLay *hl, 
 		struct genoHeatmap *gh)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh || !hl)
     return NULL;
 
 struct hvGfx *vg;
 
 char *db = gh->database;
 
 int totalW = hgHeatmapDefaultPixWidth;
 hl->height = heatmapHeight(gh) * hl->sampleHeight;
 /* Include room for probe marks */
 int totalH = hl->height + PROBE_MARK_BUFFER + PROBE_MARK_HEIGHT;
 
 if (totalW * totalH == 0)
     return NULL;
 
 struct tempName md5Tn;
 char *strToHash = cartSettingsString(hgh2Prefix, "genomeGif");
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG)
     {
     vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
 
     if (sameWord(gh->dataType,"bed 15"))
 	drawChromHeatmapsByPixel(vg->vg, db, hl, gh->name);
     else if (sameWord(gh->dataType, "bed 4") || sameWord(gh->dataType, "bed 5"))
 	drawBedGraph(vg->vg, db, hl, gh->name);
 
     hvGfxClose(&vg);
     }
 
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 return filename;
 }
 
 char *genomeSummaryGif(struct sqlConnection *conn, struct heatmapLay *hl, 
 		       struct genoHeatmap *gh, struct slName *sampleList, 
 		       int totalH, int subset, int subsetNum)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh || !hl)
     return NULL;
 
 struct hvGfx *vg;
 
 char *db = gh->database;
 
 int totalW = hgHeatmapDefaultPixWidth;
 totalH = hgHeatmapDefaultSummaryHeight; //totalH;
 hl->height = totalH;
 
 if (totalW * totalH == 0)
     return NULL;
 
 struct tempName md5Tn;
 char modifier[128];
 safef(modifier, sizeof(modifier), "genomeSummaryGif,heatmap%dof%d", subset, subsetNum);
 
 char *strToHash = cartSettingsString(hgh2Prefix, modifier);
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG)
     {
     vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
     
     if (sameWord(gh->dataType,"bed 15"))
 	drawChromSummary(vg->vg, db, hl, gh->name, sampleList, totalH);
     else if (sameWord(gh->dataType, "bed 4") || sameWord(gh->dataType, "bed 5"))
 	drawBedGraph(vg->vg, db, hl, gh->name);
 
     hvGfxClose(&vg);
     }
 
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 return filename;
 }
 
 char *genomeStatsGif(struct sqlConnection *conn, struct heatmapLay *hl, 
 		     struct genoHeatmap *gh, float colorCutoff)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh || !hl)
     return NULL;
 
 if (!gh->anaResult)
     return NULL;
 
 struct hvGfx *vg;
 
 int totalW = hgHeatmapDefaultPixWidth;
 int totalH = hghBed5Height;
 hl->height = totalH;
 
 if (totalW * totalH == 0)
     return NULL;
 
 struct tempName md5Tn;
 char *strToHash = cartSettingsString(hgh2Prefix, "genomeStatsGif");
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG)
     {
     vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
 
 //    drawBackgroundLines(vg->vg, 0, 0, totalW, totalH);
 
     if (sameWord(gh->dataType,"bed 15"))
 	drawChromStats(vg->vg, gh->anaResult->stats, hl,                                  
 		       0,  gh->anaResult->max, colorCutoff);
     hvGfxClose(&vg);
     }
 
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 return filename;
 }
 
 char *genesetStatsGif(struct sqlConnection *conn, struct heatmapLay *hl,
 		      struct genoHeatmap *gh, float colorCutoff)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh || !hl)
     return NULL;
 
 if (!gh->anaResultHash)
     return NULL;
 
 struct hvGfx *vg;
 
 int totalW = hgHeatmapDefaultPixWidth;
 int totalH = hghBed5Height;
 hl->height = totalH;
 
 if (totalW * totalH == 0)
     return NULL;
 
 struct tempName md5Tn;
 char *strToHash = cartSettingsString(hgh2Prefix, "genesetStatsGif");
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG)
     {
     vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
     
     if (sameWord(gh->dataType,"bed 15"))
 	drawGeneSetStats(vg->vg, gh->anaResultHash->hash, hl,
 			 0, 3, colorCutoff, gh->name);
     hvGfxClose(&vg);
     }
 
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 
 return filename;
 }
 
 
 char *genesetGif(struct sqlConnection *conn, struct heatmapLay *hl,
 		 struct genoHeatmap *gh)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh || !hl)
     return NULL;
 
 struct hvGfx *vg;
 
 int totalW = hgHeatmapDefaultPixWidth;
 hl->height = heatmapHeight(gh) * hl->sampleHeight;
 int totalH = hl->height;
 
 if (totalW * totalH == 0)
     return NULL;
 
 struct tempName md5Tn;
 char *strToHash = cartSettingsString(hgh2Prefix, "genesetGif");
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG)
     {
     vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
     
     drawLayoutLines(vg->vg, hl, totalW);
     
     if (sameWord(gh->dataType,"bed 15"))
 	drawGeneSetHeatmapsByPixel(vg->vg, gh->database, gh->name, hl); //, geneHash);   
     
     hvGfxClose(&vg);
     }
 
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 return filename;
 }
 
 char *genesetSummaryGif(struct sqlConnection *conn, struct heatmapLay *hl,
 			struct genoHeatmap *gh, struct slName *sampleList, 
 			int totalH, int subset, int subsetNum)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh || !hl)
     return NULL;
 
 struct hvGfx *vg;
 
 int totalW = hgHeatmapDefaultPixWidth;
 totalH = hgHeatmapDefaultSummaryHeight;
 hl->height = totalH;
 
 if (totalW * totalH == 0)
     return NULL;
 
 struct tempName md5Tn;
 char modifier[128];
 safef(modifier, sizeof(modifier), "genomeSummaryGif,heatmap%dof%d", subset, subsetNum);
 char *strToHash = cartSettingsString(hgh2Prefix, modifier);
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG)
     {
     vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
 
     drawLayoutLines(vg->vg, hl, totalW);
     
     if (sameWord(gh->dataType,"bed 15"))
 	drawGeneSetSummary(vg->vg, gh->database, gh->name, hl, sampleList, totalH);
 
     hvGfxClose(&vg);
     }
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 
 return filename;
 }
 
 
 char *featureGif(struct sqlConnection *conn, struct heatmapLay *hl, 
 		 struct genoHeatmap *gh, char *tableName)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh || !hl)
     return NULL;
 
 struct hvGfx *vg;
 
 int totalW = hl->width;
 int totalH = heatmapHeight(gh) * hl->sampleHeight;
 hl->height = totalH;
 
 if (ifDrawFeatureLabel(hl))
     totalH += hghFeatureLabel;
 
 if (totalW * totalH == 0)
     return NULL;
 
 struct tempName md5Tn;
 char *strToHash = cartSettingsString(hgh2Prefix, "featureGif");
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG) 
     {
     vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
     
     /* draw feature sorter */
     if (sameWord(gh->dataType,"bed 15"))
 	drawFeatures(vg->vg, tableName, hl, 0, FALSE);
     
     hvGfxClose(&vg);
     }
 
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 return filename;
 }
 
 char *featureSummaryGif(struct sqlConnection *conn, struct heatmapLay *hl, 
 			struct genoHeatmap *gh, char *tableName, struct slName *sampleList, 
 			int buffer, int subset, int subsetNum)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh || !hl)
     return NULL;
 
 struct hvGfx *vg;
 
 int totalW = hl->width;
 int totalH = hgHeatmapDefaultSummaryHeight; 
 hl->height = totalH;
 
 if (totalW * totalH == 0)
     return NULL;
 
 boolean isBottomSubgroup = (subset + 1 == subsetNum);
+
+struct slName **ptSubsets = NULL;
+char *raName = gh->raFile;
+int ifSubsets = getSubsetsIfAny(gh, subsetNum, raName, &ptSubsets);
+if (ifSubsets ==1)
+    isBottomSubgroup=1;
+
 if (ifDrawFeatureLabel(hl) && isBottomSubgroup)
     totalH += hghFeatureLabel;
 
 struct tempName md5Tn;
 char modifier[128];
 safef(modifier, sizeof(modifier), "featureSummaryGif,%dof%d", subset, subsetNum);
 char *strToHash = cartSettingsString(hgh2Prefix, modifier);
 				     
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG) 
     {
     vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
     
     /* draw feature sorter */
     if (sameWord(gh->dataType,"bed 15"))
 	{
 	if (subsetNum > 1)
-	    drawFeatureSubgroupBar(vg->vg, sampleList, hl->sampleHeight, isBottomSubgroup);
+	    drawFeatureSubgroupBar(vg->vg, sampleList, hl->sampleHeight, subset);
 	drawFeatureSummary(vg->vg, tableName, hl, 0, FALSE, sampleList, isBottomSubgroup);
 	}
     hvGfxClose(&vg);
     }
 
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 return filename;
 }
 
 
 char *subgroupGif(struct genoHeatmap *gh, char *tableName, int width, int sampleHeight)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh)
     return NULL;
 
 struct hvGfx *vg;
 
 int totalW = width;
 int totalH = heatmapHeight(gh) * sampleHeight;
 
 if (totalW * totalH == 0)
     return NULL;
 
 struct tempName md5Tn;
 char *strToHash = cartSettingsString(hgh2Prefix, "subgroupGif");
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG)
     {
     vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
     
     /* draw subgroups */
     if (sameWord(gh->dataType,"bed 15"))
 	drawSubgroups(vg->vg, tableName, totalW, totalH, sampleHeight);
     
     hvGfxClose(&vg);
     }
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 
 return filename;
 }
 
 
 
 void drawSimpleScale(struct vGfx *vg, struct heatmapLay *hlList, int width, int height)
 {
 if (!hlList)
     return;
 
 int minValPos = 0, maxValPos = 0;
 double minVal = INT_MAX;
 double maxVal = INT_MIN;
 
 struct heatmapLay *hl;
 for (hl = hlList; hl; hl = hl->next)
     {
     if (hl->minVal < minVal)
 	{
 	minVal = hl->minVal;
 	minValPos = hl->minValPos;
 	}
     if (hl->maxVal > maxVal)
 	{
 	maxVal = hl->maxVal;
 	maxValPos = hl->maxValPos;
 	}
     }
 
 if (minVal == INT_MAX || maxVal == INT_MIN)  // not scaled.
     return;
 
 int dashWidth = 5;
 int buffer = 2;
 
 vgSetClip(vg, 0, 0, width, height);
 vgBox(vg, 0, 0, width, height, MG_WHITE); 
 
 /* draw line on right side */
 vgBox(vg, width-1, maxValPos, 1, minValPos - maxValPos, MG_BLACK);
 
 
 int yZero = height / 2;
 if (minVal != maxVal)
     {
     double slope = ((double) (maxValPos - minValPos))/(maxVal - minVal);
     double tmp = slope * (0.0 - maxVal) + maxValPos; 
     yZero = ceil(tmp) - 1;
 
      /* Put dash at zero */
     vgBox(vg, width - dashWidth, yZero, dashWidth, 1, MG_BLACK);
     }
 
 /* Put dash at max value */
 vgBox(vg, width - dashWidth, maxValPos, dashWidth, 1, MG_BLACK);
 
 /* Label max point */
 int textHeight = maxValPos - hlList->fontHeight/2;
 if (textHeight > yZero - hlList->fontHeight)
     textHeight = yZero - hlList->fontHeight;
 if (textHeight < 0)
     textHeight = 0;
 
 char maxValStr[128];
 if (fabs(maxVal) > 1000.0)
     safef(maxValStr, sizeof(maxValStr), "%1.1e", maxVal);
 else if (fabs(maxVal) > 100.0)
     safef(maxValStr, sizeof(maxValStr), "%3.1f", maxVal);
 else if (fabs(maxVal) > 10.0)
     safef(maxValStr, sizeof(maxValStr), "%2.2f", maxVal);
 else if (fabs(maxVal) > 0 && fabs(maxVal) < 0.001)
     safef(maxValStr, sizeof(maxValStr), "%1.1e", maxVal);
 else if (fabs(maxVal) == 0)
     safef(maxValStr, sizeof(maxValStr), "0");
 else 
     safef(maxValStr, sizeof(maxValStr), "%1.3f", maxVal);
 
 vgTextRight(vg, 0, textHeight, width - (dashWidth+buffer), hlList->fontHeight,
 	    MG_BLACK, hlList->font, maxValStr);
  
 /* Put dash at min value */
 vgBox(vg, width-dashWidth, minValPos-1, dashWidth, 1, MG_BLACK);
 
 /* Label min point */
 textHeight = (minValPos-1) - hlList->fontHeight/2;
 if (textHeight < yZero + hlList->fontHeight)
     textHeight = yZero + hlList->fontHeight;
 if (textHeight > height - hlList->fontHeight)
     textHeight = height - hlList->fontHeight;
 
 char minValStr[128];
 if (fabs(minVal) > 1000.0)
     safef(minValStr, sizeof(minValStr), "%1.1e", minVal);
 else if (fabs(minVal) > 100.0)
     safef(minValStr, sizeof(minValStr), "%3.1f", minVal);
 else if (fabs(minVal) > 10.0)
     safef(minValStr, sizeof(minValStr), "%2.2f", minVal);
 else if (fabs(minVal) > 0 && fabs(minVal) < 0.001)
     safef(minValStr, sizeof(minValStr), "%1.1e", minVal);
 else if (fabs(minVal) == 0.0)
     safef(minValStr, sizeof(minValStr), "0");
 else
     safef(minValStr, sizeof(minValStr), "%1.3f", minVal);
 
 vgTextRight(vg, 0, textHeight, width - (dashWidth+buffer), hlList->fontHeight,
 	    MG_BLACK, hlList->font, minValStr); 
 
 vgUnclip(vg);
 }
 
 char *simpleScaleGif(struct heatmapLay *hl, struct genoHeatmap *gh, int subset, int subsetNum)
 /* Create genome GIF file and HT that includes it. */
 {
 if (!gh)
     return NULL;
 
 int totalW = hgDefaultScalePixWidth;
 int totalH = hl->height; 
 
 if (totalW * totalH == 0)
     return NULL;
 
 if (hl->maxVal == INT_MIN || hl->minVal == INT_MAX) // scaling not set.
     return NULL;
 
 struct tempName md5Tn;
 char modifier[128];
 safef(modifier, sizeof(modifier), "simpleScaleGif,%dof%d", subset, subsetNum);
 char *strToHash = cartSettingsString(hgh2Prefix, modifier);
 
 trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
 
 off_t size = fileSize(md5Tn.forCgi);
 if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG_IMG)
     {
     struct hvGfx *vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
     
     drawSimpleScale(vg->vg, hl, totalW, totalH);
 
     hvGfxClose(&vg);
     }
 
 char *filename = replaceChars(md5Tn.forHtml, "..", "");
 return filename;
 }
 

