src/hg/instinct/hgPathways/hgCircleMaps.c 1.15
1.15 2009/09/15 00:08:14 sbenz
Fixed gene sorting in cgi, need to add to interface
Index: src/hg/instinct/hgPathways/hgCircleMaps.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/instinct/hgPathways/hgCircleMaps.c,v
retrieving revision 1.14
retrieving revision 1.15
diff -b -B -U 1000000 -r1.14 -r1.15
--- src/hg/instinct/hgPathways/hgCircleMaps.c 18 Aug 2009 20:57:31 -0000 1.14
+++ src/hg/instinct/hgPathways/hgCircleMaps.c 15 Sep 2009 00:08:14 -0000 1.15
@@ -1,1295 +1,1352 @@
/* hgCircleMaps.c - circular heatmaps to visualize genomic high-throughput data for
* multiple samples in a concise image suitable for placement in a pathway diagram
* This file is copyright 2009 J. Zachary Sanborn. All rights reserved. */
#define EXPR_DATA_SHADES 16
#define COLOR_SCALE 1
#define DEFAULT_MAX_DEVIATION 2.5
#define PI 3.14159265358979323846
#include "common.h"
#include "bed.h"
#include "cart.h"
#include "errCatch.h"
#include "hash.h"
#include "hdb.h"
#include "hgColors.h"
#include "hPrint.h"
#include "vGfx.h"
#include "hvGfx.h"
#include "hgHeatmapLib.h"
#include "featuresLib.h"
//#include "hgHeatmap2.h"
#include "hgPathways.h"
#include "sortFeatures.h"
#include "gfxPoly.h"
#include "json.h"
#include "base64.h"
static char const rcsid[] = "$Id$";
struct ring {
struct ring *next; /* next ring */
char *name; /* name of ring (optional) */
char *dataset; /* name of dataset */
char *type; /* type of data (feature, SNP, CNV, expression, etc.) */
int numSamples; /* num samples in ring */
int numElements; /* number of elements in ring */
double minR; /* min radius for current ring */
double maxR; /* max radius for current ring */
double gain; /* gain setting for ring */
double maxDev; /* maximum deviation for ring */
double median; /* median value for ring */
boolean drawSummary; /* draw summary */
boolean visible; /* invisible rings exist, can be sorted, not displayed */
boolean sort; /* Sort this ring or not */
int sortDir; /* +1 = ASC, -1 = DESC */
double maxVal; /* maximum value, used for sorting */
struct rgbColor *zeroColor; /* colors for ring */
struct rgbColor *highColor;
struct rgbColor *lowColor;
struct hash *data; /* data for ring, keyed by sample */
};
struct circleMap {
char *name; /* name of circle (e.g. gene) */
struct slName *samples; /* global list of samples, order of samples defined
* position in ring */
struct ring *rings; /* data defining all rings for circle */
double minR; /* min radius for all rings */
double maxR; /* max radius for all rings */
int width; /* pix width */
char *position; /* position, if available */
boolean drawSummary; /* draw summary (pie chart) */
boolean subgrouped; /* if successful subgrouping applied */
Color bgColor; /* Background color (default = MG_WHITE) */
struct hash *data; /* data for center (pie chart) */
};
char *heatMapDbProfile = "localDb";
double maxDev(struct genoHeatmap *gh)
{
if (!gh)
return DEFAULT_MAX_DEVIATION;
if (gh->expScale > 0)
return gh->expScale;
return DEFAULT_MAX_DEVIATION;
}
void vgRadiusBox(struct vGfx *vg, int x0, int y0, double rStart, double rStop,
double tStart, double tStop, Color col)
{
if ((rStart == rStop) || (tStart == tStop)) // zero height/thickness
return;
struct gfxPoly *gp = gfxPolyNew();
int x, y;
double r, t, tStep;
/* rStart -> rStop @ tStart -- one point at each edge */
for (r = rStart; r <= rStop; r += (rStop - rStart))
{
x = floor(r * sin(tStart));
y = floor(r * cos(tStart));
gfxPolyAddPoint(gp, x + x0, y + y0);
}
/* tStart -> tStop @ rStop */
tStep = 0.5 / (2.0*PI);
for (t = tStart; t <= tStop; t += tStep)
{
x = floor(rStop * sin(t));
y = floor(rStop * cos(t));
gfxPolyAddPoint(gp, x + x0, y + y0);
}
/* rStop -> rStart @ tStop -- one point at each edge */
for (r = rStop; r >= rStart; r -= (rStop - rStart))
{
x = floor(r * sin(tStop));
y = floor(r * cos(tStop));
gfxPolyAddPoint(gp, x + x0, y + y0);
}
/* tStop -> tStop @ rStart */
for (t = tStop; t <= tStart; t -= tStep)
{
x = floor(rStart * sin(t));
y = floor(rStart * cos(t));
gfxPolyAddPoint(gp, x + x0, y + y0);
}
vgDrawPoly(vg, gp, col, TRUE);
gfxPolyFree(&gp);
}
struct circleMap *newCircleMap(char *name)
{ /* initialize a new circle layout */
struct circleMap *cm;
AllocVar(cm);
cm->name = cloneString(name);
cm->samples = NULL;
cm->rings = NULL;
cm->position = NULL;
cm->minR = 0.0;
cm->maxR = 0.0;
cm->width = 0;
cm->data = NULL;
cm->drawSummary = FALSE;
cm->subgrouped = FALSE;
cm->bgColor = MG_WHITE;
return cm;
}
struct ring *cloneRing(struct ring *rg)
{
if (!rg)
return NULL;
struct ring *rg2;
AllocVar(rg2);
rg2->name = cloneString(rg->name);
rg2->dataset = cloneString(rg->dataset);
rg2->type = cloneString(rg->type);
rg2->maxDev = rg->maxDev;
rg2->gain = rg->gain;
rg2->numSamples = rg->numSamples;
rg2->numElements = rg->numElements;
rg2->minR = rg->minR;
rg2->maxR = rg->maxR;
rg2->data = NULL;
rg2->drawSummary = rg->drawSummary;
rg2->median = 0.0;
rg2->visible = rg->visible;
rg2->sort = rg->sort;
rg2->sortDir = rg->sortDir;
rg2->maxVal = rg->maxVal;
return rg2;
}
//struct ring *addRingToCircleMap(struct circleMap *cm, char *name, struct genoHeatmap *gh,
struct ring *addRingToCircleMap(struct circleMap *cm, char *name, char *dataset,
boolean isVisible)
{ /* Initialize a ring, add it circleLay, and return reference to it */
struct ring *rg;
AllocVar(rg);
rg->name = cloneString(name);
rg->dataset = cloneString(dataset);
rg->type = "CNV"; // for red/blue/white color scheme
rg->maxDev = DEFAULT_MAX_DEVIATION;
rg->gain = 1.0;
//rg->dataset = cloneString(gh->name);
//rg->type = cloneString(gh->platform);
//rg->maxDev = maxDev(gh);
//rg->gain = gh->gainSet;
rg->numElements = 0;
rg->numSamples = 0;
rg->minR = 0.0;
rg->maxR = 0.0;
rg->data = NULL;
rg->drawSummary = FALSE;
rg->median = 0.0;
rg->visible = isVisible;
rg->sort = FALSE;
rg->sortDir = 0;
rg->maxVal = 0.0;
slAddHead(&cm->rings, rg);
return rg;
}
void addGeneDataToRing(struct ring *rg, struct hash *geneHash, char *gene,
struct slName *sampleList)
{
if (!rg || !geneHash)
return;
struct hashEl *el = hashLookup(geneHash, gene);
if (!el)
return;
rg->numElements = slCount(el); // number of probes matching 'gene'
if (!rg->data)
rg->data = hashNew(0);
struct bed *nb, *nbList = NULL;
for ( ; el; el = el->next) // loop through probes
{
struct bed *nb = el->val;
if (!nb)
continue;
slAddHead(&nbList, nb);
}
int i = 0;
struct slName *sl;
struct slDouble *ad, *allData = NULL;
for (sl = sampleList; sl; sl = sl->next)
{
struct slDouble *sd, *sdList = NULL;
for (nb = nbList; nb; nb = nb->next)
{
double val = nb->expScores[i];
if (val > rg->maxVal)
rg->maxVal = val;
/* Store copy of all data in list to determine median value */
ad = slDoubleNew(val);
slAddHead(&allData, ad);
sd = slDoubleNew(val);
slAddHead(&sdList, sd);
}
hashAdd(rg->data, sl->name, sdList);
rg->numSamples += 1;
i++; // increment expId pointer.
}
if(rg->numSamples > 0)
rg->median = slDoubleMedian(allData);
else
rg->median = 0;
slFreeList(&allData);
}
void addFeatureDataToRing(struct sqlConnection *conn, struct ring *rg,
struct column *col, struct genoHeatmap *gh)
{
if (!rg)
return;
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;
if (!rg->data)
rg->data = hashNew(0);
rg->numElements = 1;
char *colorScheme = col->cellColorSchemeVal(col);
rg->type = cloneString(colorScheme);
double absVal, minVal, maxVal;
char *minCutVal, *maxCutVal;
double offset;
int reverseColor;
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);
if ((minVal<0) && (maxVal>0)) /* double color: like microarray */
rg->maxDev = max(fabs(minVal), maxVal);
else /* single color */
rg->maxDev = maxVal - minVal;
int i = 0;
struct slName *sl;
double val;
for (sl = gh->sampleList; sl; sl = sl->next)
{
struct slName * 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);
if (!(minVal < 0 && maxVal > 0))
absVal = absVal - minVal;
if (val < 0)
absVal = -1.0 * absVal;
/* store final data for drawing in ring */
struct slDouble *sd = slDoubleNew(absVal);
hashAdd(rg->data, sl->name, sd);
i++; // increment expId pointer.
rg->numSamples++;
freez(&id);
}
}
void ringSetDefaultColor(struct ring *rg)
{
if (!rg)
return;
struct rgbColor *low = NULL;
struct rgbColor *zero = NULL;
struct rgbColor *high = NULL;
// default, expression data color scheme:
if (sameString(rg->type, "CNV") || sameString(rg->type, "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 if (sameString(rg->type, "default"))
{ // default feature color
AllocVar(zero); // black
zero->r = 0;
zero->g = 0;
zero->b = 0;
AllocVar(high); // yellow
high->r = 255;
high->g = 255;
high->b = 0;
AllocVar(low); // green
low->r = 0;
low->g = 255;
low->b = 0;
}
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;
}
rg->zeroColor = zero;
rg->highColor = high;
rg->lowColor = low;
}
Color getColor(double val, double gain, double maxDev,
Color *upShades, Color *downShades)
{
double absVal = fabs(val);
double colorScale = COLOR_SCALE / maxDev;
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)
return upShades[colorIndex];
else
return downShades[colorIndex];
}
void drawCenter(struct vGfx *vg, struct circleMap *cm)
{
int x0 = cm->width / 2;
int y0 = cm->width / 2;
Color valCol;
double tStart, tStop;
Color upShades[EXPR_DATA_SHADES];
Color downShades[EXPR_DATA_SHADES];
double numSamples = (double) slCount(cm->samples);
int rStart = 0;
int rStop = cm->minR - 2;
cm->data = hashNew(0);
int maxSamples = -1;
struct ring *rg;
for (rg = cm->rings; rg; rg = rg->next)
{
if (rg->numSamples == 0)
continue;
ringSetDefaultColor(rg);
vgMakeColorGradient(vg, rg->zeroColor, rg->highColor, EXPR_DATA_SHADES, upShades);
vgMakeColorGradient(vg, rg->zeroColor, rg->lowColor, EXPR_DATA_SHADES, downShades);
valCol = getColor(rg->median, rg->gain, rg->maxDev, upShades, downShades);
if (rg->numSamples > maxSamples)
{ /* Only set background color according to median color of largest subset */
maxSamples = rg->numSamples;
cm->bgColor = valCol;
}
double i = 0.0;
struct slName *sl;
for (sl = cm->samples; sl; sl = sl->next)
{
tStart = 2.0 * PI * i / numSamples;
tStop = 2.0 * PI * (i + 1.0) / numSamples;
i += 1.0;
struct hashEl *el = hashLookup(rg->data, sl->name);
if (!el)
continue; // data exists for sample in ring
el = hashLookup(cm->data, sl->name);
if (el)
continue; // only save the first
vgRadiusBox(vg, x0, y0, rStart, rStop, tStart, tStop, valCol);
hashAddInt(cm->data, sl->name, 1);
}
}
}
void drawRing(struct vGfx *vg, struct circleMap *cm, struct ring *rg)
{
if (!rg->data || rg->numSamples == 0) // nothing to draw
return;
Color valCol;
Color upShades[EXPR_DATA_SHADES];
Color downShades[EXPR_DATA_SHADES];
ringSetDefaultColor(rg);
vgMakeColorGradient(vg, rg->zeroColor, rg->highColor, EXPR_DATA_SHADES, upShades);
vgMakeColorGradient(vg, rg->zeroColor, rg->lowColor, EXPR_DATA_SHADES, downShades);
int x0 = cm->width / 2;
int y0 = cm->width / 2;
double rStart, rStop, rStep = (rg->maxR - rg->minR) / (double) rg->numElements;
double tStart, tStop; // thetas
double numSamples = (double) slCount(cm->samples);
double i = 0.0;
struct slName *sl;
struct slDouble *sd, *sdList;
Color missingColor;
if (cm->subgrouped)
missingColor = MG_WHITE; // when subgrouped, gray is distracting
else
missingColor = vgFindColorIx(vg, 225, 225, 225); // light gray
for (sl = cm->samples; sl; sl = sl->next)
{
tStart = 2.0 * PI * i / numSamples;
tStop = 2.0 * PI * (i + 1.0) / numSamples;
i += 1.0;
rStart = rg->minR;
rStop = rg->minR + rStep;
struct hashEl *el = hashLookup(rg->data, sl->name);
if (!el)
{
vgRadiusBox(vg, x0, y0, rg->minR, rg->maxR, tStart, tStop, missingColor);
continue;
}
sdList = el->val; // slDouble list of expScores
if (rg->drawSummary)
{
valCol = getColor(rg->median, rg->gain, rg->maxDev, upShades, downShades);
vgRadiusBox(vg, x0, y0, rg->minR, rg->maxR, tStart, tStop, valCol);
}
else
{
for (sd = sdList; sd; sd = sd->next)
{
valCol = getColor(sd->val, rg->gain, rg->maxDev, upShades, downShades);
vgRadiusBox(vg, x0, y0, rStart, rStop, tStart, tStop, valCol);
rStart += rStep;
rStop += rStep;
}
}
}
}
void drawRings(struct vGfx *vg, struct circleMap *cm)
{
struct ring *rg;
for (rg = cm->rings; rg; rg = rg->next)
drawRing(vg, cm, rg);
}
void drawName(struct vGfx *vg, struct circleMap *cm)
{
Color textCol = vgContrastingColor(vg, cm->bgColor);
/* Draw gene label in center */
MgFont *font = mgSmallFont();
if (mgFontStringWidth(font, cm->name) < cm->minR * 2.0)
vgTextCentered(vg, 0, 0, cm->width, cm->width, textCol, font, cm->name);
else
{
char *substr;
strncpy(substr,cm->name,5);
vgTextCentered(vg, 0, 0, cm->width, cm->width, textCol, font, substr);
}
}
void drawCircleMap(struct vGfx *vg, struct circleMap *cm)
{
if (!cm)
return;
if (cm->drawSummary)
drawCenter(vg, cm);
drawRings(vg, cm);
drawName(vg, cm);
}
struct tempName *circleMapGif(struct circleMap *cm)
/* Create genome GIF file and HT that includes it. */
{
if (!cm)
return NULL;
if (cm->width == 0)
return NULL;
struct tempName *md5Tn;
AllocVar(md5Tn);
char modifier[128];
safef(modifier, sizeof(modifier), "circleMapGif_%s", cm->name);
char *strToHash = cartSettingsString(hgh2Prefix, modifier);
trashDirMD5File(md5Tn, "hgh", ".png", strToHash);
off_t size = fileSize(md5Tn->forCgi);
if (!fileExists(md5Tn->forCgi) || (size == 0) || DEBUG_IMG)
{
struct hvGfx *vg = hvGfxOpenPng(cm->width, cm->width, md5Tn->forCgi, TRUE);
drawCircleMap(vg->vg, cm);
hvGfxClose(&vg);
}
//char *filename = replaceChars(md5Tn.forHtml, "..", "");
return md5Tn;
}
void layoutCircleMap(struct circleMap *cm, int width)
{
if (!cm)
return;
cm->width = width;
cm->maxR = (double) cm->width / 2.0;
cm->minR = 0.75 * cm->maxR;
struct ring *rg;
int numRings = 0;
for (rg = cm->rings; rg; rg = rg->next)
if (rg->visible && rg->numSamples > 0)
numRings++;
double step = (cm->maxR - cm->minR) / (double) numRings;
double buffer = 0.0;
if (cm->subgrouped)
buffer = step / 5.0; // put buffer between subgrouped rings
buffer = 0.0;
double rStart = cm->minR;
double rStop = cm->minR + (step - buffer);
for (rg = cm->rings; rg; rg = rg->next)
{
if (!rg->visible) // keep minR/maxR at zero for no display
continue;
if (rg->numSamples == 0)
{
rg->minR = 0;
rg->maxR = 0;
continue;
}
rg->minR = rStart;
rg->maxR = rStop;
rStart += step;
rStop = rStart + (step - buffer);
}
}
void setupSummary(struct circleMap *cm)
{
char *summarize = cartOptionalString(cart, hgh2Summarize);
if (!summarize)
return;
char *summaryDisplay = cartUsualString(cart, hgh2SummaryDisplay, "ring");
struct slName *names = slNameListFromComma(summarize);
if (sameString(summaryDisplay, "center"))
cm->drawSummary = TRUE;
else
{
struct ring *rg;
for (rg = cm->rings; rg; rg = rg->next)
{
if (!slNameInList(names, rg->name))
continue;
rg->drawSummary = TRUE;
}
}
}
int getSortDirection(char *sortDir)
{
if (!sortDir)
return 1;
if (sameString(sortDir, "ASC"))
return 1;
if (sameString(sortDir, "DESC"))
return -1;
return 0;
}
void setupGeneSort(struct circleMap *cm)
{
char *sortGene = cartOptionalString(cart, hgh2SortGene);
char *sortDir = cartOptionalString(cart, hgh2SortDir);
if (!sortGene)
return;
-
struct ring *rg;
for (rg = cm->rings; rg; rg = rg->next)
{
if (!sameString(rg->name, sortGene))
continue;
rg->sort = TRUE;
rg->sortDir = getSortDirection(sortDir);
}
}
void setupFeatureSort(struct circleMap *cm)
{
char *featureSort = cartOptionalString(cart, hgh2FeatureSort);
char *featureSortDir = cartOptionalString(cart, hgh2FeatureSortDir);
struct slName *s, *sorts = NULL;
struct slName *d, *sortDirs = NULL;
if (!featureSort || !featureSortDir)
return;
sorts = slNameListFromComma(featureSort);
sortDirs = slNameListFromComma(featureSortDir);
if (slCount(sorts) != slCount(sortDirs))
{
sorts = NULL;
sortDirs = NULL;
}
if (!sorts || !sortDirs)
return;
struct ring *rg;
for (rg = cm->rings; rg; rg = rg->next)
{
for (s = sorts, d = sortDirs; s && d; s = s->next, d = d->next)
{
if (!sameString(rg->name, s->name))
continue;
rg->sort = TRUE;
rg->sortDir = getSortDirection(d->name);
}
}
}
void sortCircleMap(struct circleMap *cm)
{ /* Sorting goes from inner ring out, averaging all probes within a ring for value */
if (!cm->samples)
return;
/* Setup sorting by gene value */
setupGeneSort(cm);
/* Setup sorting by features */
setupFeatureSort(cm);
struct slName *sl;
struct sortNode *sn, *snList = NULL;
struct sortList *snl;
struct ring *rg;
double val = 0.0;
for (sl = cm->samples; sl; sl = sl->next)
{
AllocVar(sn);
sn->name = cloneString(sl->name);
sn->list = NULL;
for (rg = cm->rings; rg; rg = rg->next)
{
if (!rg->sort || !rg->data)
continue;
struct hashEl *el = hashLookup(rg->data, sn->name);
if (!el)
val = rg->maxVal + 0.1; // if no val, set all to above inner ring's max
else
{
struct slDouble *sd, *sdList = el->val;
val = 0.0;
for (sd = sdList; sd; sd = sd->next)
val += sd->val;
val = val / (double) slCount(sdList);
}
AllocVar(snl);
snl->val = val;
snl->sortDirection = rg->sortDir;
slAddHead(&sn->list, snl);
}
slReverse(&sn->list);
slAddHead(&snList, sn);
}
slSort(&snList, sortNodeCmp);
if (slCount(snList) != slCount(cm->samples))
errAbort("Sorting failed");
/* Reset sample order for circleMap according to sorted list of nodes */
slNameFreeList(&cm->samples);
cm->samples = NULL;
for (sn = snList; sn; sn = sn->next)
slNameAddHead(&cm->samples, sn->name);
slReverse(&cm->samples);
}
void addFeaturesToCircleMap(struct circleMap *cm, char *datasets)
{
char *featureList = cartOptionalString(cart, hgh2FeatureList);
char *featureSort = cartOptionalString(cart, hgh2FeatureSort);
if (!featureList && !featureSort)
return;
boolean isVisible = TRUE;
if (!featureList)
{ /* in order for sorting to work, need features in rings,
* but should be invisible unless featureList was set */
isVisible = FALSE;
featureList = featureSort;
}
struct genoHeatmap *gh = NULL;
struct slName *sl, *slList = slNameListFromComma(datasets);
for (sl = slList; sl; sl = sl->next)
{
struct hashEl *el = hashLookup(ghHash, sl->name);
if (!el)
continue;
gh = el->val;
break;
}
if (!gh)
return;
struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, gh->patDb);
if (!conn)
return;
struct slName *f, *features = slNameListFromComma(featureList);
/* necessary to reverse so that features show up in order requested
* can't do a reverse of the rings in case there's already gene rings */
slReverse(&features);
char *raName = gh->raFile;
struct column *col, *colList = getColumns(conn, raName, gh->patDb);
for (f = features; f; f = f->next)
{
for (col = colList; col != NULL; col = col->next)
if (sameString(col->name, f->name))
break;
if (!col)
continue;
struct ring *rg = addRingToCircleMap(cm, col->name, gh->patDb, isVisible);
addFeatureDataToRing(conn, rg, col, gh);
}
hFreeConn(&conn);
}
boolean needSplit(struct ring *rg, struct slName **ptSubsets, int subsetNum)
{
if (!rg->data) // no data, no split!
return FALSE;
int subset;
struct slName *sl;
int needSplit = 0;
for (subset = 0; subset < subsetNum; subset++)
{
for (sl = ptSubsets[subset]; sl; sl = sl->next)
{
struct hashEl *el = hashLookup(rg->data, sl->name);
if (!el)
continue;
needSplit++;
break;
}
}
return (needSplit > 1);
}
void removeSamplesFromRing(struct circleMap *cm, struct ring *rg)
{
if (!rg->data) // nothing to remove
return;
struct slName *sl, *slList = NULL;
for (sl = cm->samples; sl; sl = sl->next)
{
struct hashEl *el = hashLookup(rg->data, sl->name);
if (!el) // not in ring, so keep sample name in circle's list
slNameAddHead(&slList, sl->name);
}
slReverse(&slList); // maintain original order
slNameFreeList(&cm->samples); // remove old sample list
cm->samples = slList; // set to new list.
}
boolean removeRing(struct circleMap *cm, struct ring *toRemove)
{
/* Delete any circleMap samples on the to-be-removed ring */
removeSamplesFromRing(cm, toRemove);
if (cm->rings == toRemove) /* ring to remove is first, set to next */
{
cm->rings = toRemove->next;
return TRUE;
}
struct ring *rg, *nextRg;
for (rg = cm->rings; rg; rg = rg->next)
{
nextRg = rg->next;
if (nextRg == toRemove)
{
rg->next = nextRg->next;
nextRg->next = NULL;
return TRUE;
}
}
return FALSE;
}
void splitRing(struct circleMap *cm, struct ring *rg,
struct slName **ptSubsets, int subsetNum)
{
int subset;
struct ring *rg2, *ptRings[subsetNum];
for (subset = 0; subset < subsetNum; subset++)
ptRings[subset] = cloneRing(rg);
struct slName *sl;
for (subset = 0; subset < subsetNum; subset++)
{
rg2 = ptRings[subset];
rg2->data = hashNew(0);
rg2->numSamples = 0;
struct slDouble *ad, *allData = NULL;
for (sl = ptSubsets[subset]; sl; sl = sl->next)
{
struct hashEl *el = hashLookup(rg->data, sl->name);
if (!el)
continue;
struct slDouble *sd, *sdList = el->val;
for (sd = sdList; sd; sd = sd->next)
{
ad = slDoubleNew(sd->val);
slAddHead(&allData, ad);
}
hashAdd(rg2->data, sl->name, sdList);
rg2->numSamples += 1;
hashRemove(rg->data, sl->name);
rg->numSamples -= 1;
}
if (!allData)
rg2->median = 0.0;
else
rg2->median = slDoubleMedian(allData);
slFreeList(&allData);
}
/* Remove empty ring that was split */
boolean removed = removeRing(cm, rg);
if (!removed)
errAbort("Problem with removing ring");
/* Add new rings */
for (subset = 0; subset < subsetNum; subset++)
{
rg = ptRings[subset];
if (rg->numSamples == 0)
continue;
slAddHead(&cm->rings, ptRings[subset]);
}
}
void splitCircleMap(struct circleMap *cm, struct slName **ptSubsets, int subsetNum)
{
struct ring *rg, *nextRg;
for (rg = cm->rings; rg; rg = nextRg)
{
nextRg = rg->next;
if (!needSplit(rg, ptSubsets, subsetNum))
continue;
splitRing(cm, rg, ptSubsets, subsetNum);
}
}
void setupSubgroups(struct circleMap *cm, char *datasets)
{
if (!ghHash) // global, make sure it exists
return;
int subsetNum = cartUsualInt(cart, hgh2SubgroupNum, hgh2SubgroupDefaultNum);
struct slName *sl, *slList = slNameListFromComma(datasets);
for (sl = slList; sl; sl = sl->next)
{
struct hashEl *el = hashLookup(ghHash, sl->name);
if (!el)
continue;
struct genoHeatmap *gh = el->val;
setSubgroups(gh);
char *raName = gh->raFile;
/* This is the key function to call for subgrouping */
struct slName **ptSubsets = NULL;
int ifSubsets = getSubsetsIfAny(gh, subsetNum, raName, &ptSubsets);
if (!ifSubsets)
continue;
splitCircleMap(cm, ptSubsets, subsetNum);
cm->subgrouped = TRUE;
}
}
void getPositionForCircleMap(struct circleMap *cm, char *gene,char *pathwayID)
{
char query[1024];
struct sqlResult *sr = NULL;
char **row = NULL;
struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, database);
if (!conn)
errAbort("Couldn't connect to database");
safef(query, sizeof(query), "select dot_pos from entities_ncipid where pathway_id = \"%s\" AND entity_name LIKE(\"%s\")",
pathwayID,gene);
sr = sqlGetResult(conn, query);
if((row = sqlNextRow(sr)) != NULL)
{
cm->position = cloneString(row[0]);
}
sqlFreeResult(&sr);
hFreeConn(&conn);
}
void getSamplesAndDataForCircleMap(struct circleMap *cm, char *datasets, char *gene,char *pathwayID,struct slName *validSamples)
{
if(!datasets || !gene || !pathwayID)
return;
struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, database);
if (!conn)
errAbort("Couldn't connect to database");
+char *sortGene = cartOptionalString(cart, hgh2SortGene);
+
struct slName *sl, *slList = slNameListFromComma(datasets);
for (sl = slList; sl; sl = sl->next)
{
char query[1024];
struct sqlResult *sr = NULL;
char **row = NULL;
+ if(sl == slList && !sameString(sortGene,"")) // first time
+ {
+ struct ring *sortRg = addRingToCircleMap(cm, sortGene, sl->name, FALSE);
+ if (!sortRg->data)
+ sortRg->data = hashNew(0);
+ if(sameStringN(sl->name,"factorGraph_",12) || sameStringN(sl->name,"fg_",3) )
+ {
+ safef(query, sizeof(query), "select DISTINCT samples.name,val from %s data inner join entities_ncipid on entities_ncipid.entity_id = data.feature_id left join samples on samples.id = data.sample_id where data.pathway_id = \"%s\" AND entities_ncipid.entity_name LIKE(\"%s\") ORDER BY samples.id",
+ sl->name,pathwayID,sortGene);
+ }
+ else
+ {
+ safef(query, sizeof(query), "select DISTINCT samples.name,val from %s data inner join analysisFeatures on analysisFeatures.id = data.feature_id left join samples on samples.id = data.sample_id where analysisFeatures.feature_name LIKE(\"%s\") ORDER BY samples.id",
+ sl->name,sortGene);
+ }
+ sr = sqlGetResult(conn, query);
+ struct slDouble *ad, *allData = NULL;
+ while ((row = sqlNextRow(sr)) != NULL)
+ {
+ if(!slNameInList(validSamples,row[0]))
+ continue;
+
+ double val = strtod(row[1],NULL);
+ struct slDouble *sd, *sdList = NULL;
+
+ if (val > sortRg->maxVal)
+ sortRg->maxVal = val;
+
+ /* Store copy of all data in list to determine median value */
+ ad = slDoubleNew(val);
+ slAddHead(&allData, ad);
+
+ sd = slDoubleNew(val);
+ slAddHead(&sdList, sd);
+
+ if (slNameInList(cm->samples, row[0]))
+ continue;
+ slNameAddHead(&cm->samples, row[0]);
+ hashAdd(sortRg->data, row[0], sdList);
+ sortRg->numSamples += 1;
+ //printf("%s %s",row[0],row[1]);
+ }
+ if(sortRg->numSamples > 0)
+ sortRg->median = slDoubleMedian(allData);
+ else
+ sortRg->median = 0;
+ //printf("%d\n",rg->numSamples);
+ sortRg->numElements = 1; // number of probes matching 'gene'
+ slFreeList(&allData);
+ slReverse(&cm->samples);
+ sqlFreeResult(&sr);
+
+ }
+
+
struct ring *rg = addRingToCircleMap(cm, gene, sl->name, TRUE);
if (!rg->data)
rg->data = hashNew(0);
if(sameStringN(sl->name,"factorGraph_",12) || sameStringN(sl->name,"fg_",3) )
{
safef(query, sizeof(query), "select DISTINCT samples.name,val from %s data inner join entities_ncipid on entities_ncipid.entity_id = data.feature_id left join samples on samples.id = data.sample_id where data.pathway_id = \"%s\" AND entities_ncipid.entity_name LIKE(\"%s\") ORDER BY samples.id",
sl->name,pathwayID,gene);
}
else
{
safef(query, sizeof(query), "select DISTINCT samples.name,val from %s data inner join analysisFeatures on analysisFeatures.id = data.feature_id left join samples on samples.id = data.sample_id where analysisFeatures.feature_name LIKE(\"%s\") ORDER BY samples.id",
sl->name,gene);
}
//errAbort(query);
//printf("%s",query);
+ slReverse(&cm->samples);
sr = sqlGetResult(conn, query);
struct slDouble *ad, *allData = NULL;
while ((row = sqlNextRow(sr)) != NULL)
{
if(!slNameInList(validSamples,row[0]))
continue;
double val = strtod(row[1],NULL);
struct slDouble *sd, *sdList = NULL;
if (val > rg->maxVal)
rg->maxVal = val;
/* Store copy of all data in list to determine median value */
ad = slDoubleNew(val);
slAddHead(&allData, ad);
sd = slDoubleNew(val);
slAddHead(&sdList, sd);
- if (slNameInList(cm->samples, row[0]))
- continue;
+ if (!slNameInList(cm->samples, row[0]))
slNameAddHead(&cm->samples, row[0]);
- hashAdd(rg->data, row[0], sdList);
+
rg->numSamples += 1;
+ hashAdd(rg->data, row[0], sdList);
//printf("%s %s",row[0],row[1]);
}
if(rg->numSamples > 0)
rg->median = slDoubleMedian(allData);
else
rg->median = 0;
//printf("%d\n",rg->numSamples);
rg->numElements = 1; // number of probes matching 'gene'
slFreeList(&allData);
slReverse(&cm->samples);
sqlFreeResult(&sr);
}
slReverse(&cm->rings);
hFreeConn(&conn);
}
void setupCircleMap(struct circleMap *cm, char *datasets, char *gene, int width, char *pathwayID, struct slName *validSamples)
{
getSamplesAndDataForCircleMap(cm,datasets,gene,pathwayID,validSamples);
getPositionForCircleMap(cm,gene,pathwayID);
/* Add feature data (if featureList is set) */
addFeaturesToCircleMap(cm, datasets);
/* Set up any subgrouping info */
setupSubgroups(cm, datasets);
/* Setup any summary (displaying only median value) */
setupSummary(cm);
/* Sort circle map data */
sortCircleMap(cm);
/* Layout circle map, setting radial widths */
layoutCircleMap(cm, width);
}
struct slName *getValidSamples(analysisID)
{
if(!analysisID)
return;
struct slName *sampleList = NULL;
struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, database);
if (!conn)
errAbort("Couldn't connect to database");
char query[1024];
struct sqlResult *sr = NULL;
char **row = NULL;
safef(query, sizeof(query), "select COUNT(*) as num from analyses,datasetCohort where analyses.id = %s and datasetCohort.cohort_id = analyses.cohort_id",
analysisID);
sr = sqlGetResult(conn, query);
char *datasetCount;
if ((row = sqlNextRow(sr)) != NULL)
{
datasetCount = cloneString(row[0]);
}
else
errAbort("Invalid AnalysisID");
sqlFreeResult(&sr);
safef(query, sizeof(query), "select COUNT(*) as num,samples.name from analyses,samples,datasetCohort where analyses.id = '%s' and datasetCohort.cohort_id = analyses.cohort_id AND datasetCohort.dataset_id = samples.dataset_id group by samples.id HAVING num = '%s'",
analysisID,datasetCount);
sr = sqlGetResult(conn, query);
while ((row = sqlNextRow(sr)) != NULL)
{
slNameAddHead(&sampleList, cloneString(row[1]));
}
sqlFreeResult(&sr);
slReverse(&sampleList);
hFreeConn(&conn);
return sampleList;
}
void modeGetCircleMap()
{
char *datasets = cartOptionalString(cart, hgh2Dataset);
int width = cartUsualInt(cart, hgh2Width, hgHeatmapDefaultCircleWidth);
char *geneList = cartOptionalString(cart, hgh2Gene);
//char *featureList = cartOptionalString(cart, hgh2FeatureList);
char *pathway = cartOptionalString(cart, hgh2PathwayID);
char *analysis = cartOptionalString(cart, hgh2AnalysisID);
struct slName *genes, *gene;
if (!datasets || !geneList)
return;
genes = slNameListFromComma(geneList);
struct json *js = newJson();
struct json *set = NULL, *group = jsonAddContainerList(js, "heatmaps");
struct slName *validSamples = getValidSamples(analysis);
for(gene = genes; gene != NULL; gene = gene->next)
{
if(set == NULL)
set = group;
else
set = jsonAddContainerToList(&group);
struct circleMap *cm = newCircleMap(gene->name);
setupCircleMap(cm, datasets, gene->name, width, pathway, validSamples);
//char *gGif = circleMapGif(cm);
struct tempName *cmTmp = circleMapGif(cm);
char *gGif = replaceChars(cmTmp->forHtml, "..", "");
// OPEN the file we just made and stream it to the output
// this saves us the multiple image openings on the client side
// cause we can just create the image using the data: URI
FILE * pFile;
long lSize;
char * buffer;
size_t result;
pFile = fopen ( cmTmp->forCgi , "rb" );
if (pFile==NULL) {fputs ("File error",stderr); exit (1);}
// obtain file size:
fseek (pFile , 0 , SEEK_END);
lSize = ftell (pFile);
rewind (pFile);
// allocate memory to contain the whole file:
buffer = (char*) malloc (sizeof(char)*lSize);
if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);}
// copy the file into the buffer:
result = fread (buffer,1,lSize,pFile);
if (result != lSize) {fputs ("Reading error",stderr); exit (3);}
char *pngFileString64 = base64Encode(buffer,lSize);
jsonAddString(set,"gene", gene->name);
jsonAddString(set, "heatmap", gGif);
jsonAddString(set, "heatmapData", pngFileString64);
if(cm->position != NULL)
jsonAddString(set, "position", cm->position);
}
hPrintf("%s", js->print(js));
}