src/hg/instinct/hgBamBam/hgAnnotations.c 1.4
1.4 2010/05/31 00:55:19 jsanborn
minor tweakpdate
Index: src/hg/instinct/hgBamBam/hgAnnotations.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/instinct/hgBamBam/hgAnnotations.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -b -B -U 1000000 -r1.3 -r1.4
--- src/hg/instinct/hgBamBam/hgAnnotations.c 31 May 2010 00:02:34 -0000 1.3
+++ src/hg/instinct/hgBamBam/hgAnnotations.c 31 May 2010 00:55:19 -0000 1.4
@@ -1,1413 +1,1414 @@
/********************************************************************************/
/* Copyright 2007-2009 -- The Regents of the University of California */
/********************************************************************************/
/* hgAnnotations.c
* Work in progress -- to draw small region to full genome annotation tracks (i.e. RefSeq)
* in concise space.
*/
#define EXPR_DATA_SHADES 16
#define DEFAULT_MAX_DEVIATION 0.7
#define COLOR_SCALE 1
#define GENE_HEIGHT 10
#include <limits.h>
#include "common.h"
#include "cart.h"
#include "vGfx.h"
#include "hvGfx.h"
#include "hCommon.h"
#include "hdb.h"
#include "cytoBand.h"
#include "hCytoBand.h"
#include "spaceSaver.h"
#include "hgBamBam.h"
#include "chromColors.h"
static char const rcsid[] = "$Id$";
static char *heatMapDbProfile = "localDb"; // database profile to use
struct bed *getTable(struct sqlConnection *conn, char *tableName, char *chromName)
{
char **row = NULL;
char query[256];
if (chromName)
safef(query, sizeof(query), "select * from %s where chrom = \"%s\"",
tableName, chromName);
else
safef(query, sizeof(query), "select * from %s", tableName);
if (!sqlExists(conn, query))
return NULL;
struct sqlResult *sr = sqlGetResult(conn, query);
struct bed *tuple = NULL;
struct bed *tupleList = NULL;
while ((row = sqlNextRow(sr)) != NULL)
{
struct genePred *gp = genePredLoad(row+1);
tuple = bedFromGenePred(gp);
if (row[12])
tuple->name = cloneString(row[12]); // switch to name2;
genePredFree(&gp);
slAddHead(&tupleList, tuple);
}
sqlFreeResult(&sr);
return tupleList;
}
struct bed *getSNPs(struct sqlConnection *conn, char *tableName, char *chromName,
int start, int stop)
{
char **row = NULL;
char query[256];
if (chromName)
safef(query, sizeof(query),
"select * from %s where chrom = \"%s\" and chromStart > %d and chromEnd <= %d",
tableName, chromName, start, stop);
else
safef(query, sizeof(query), "select * from %s", tableName);
if (!sqlExists(conn, query))
return NULL;
struct sqlResult *sr = sqlGetResult(conn, query);
struct bed *bed = NULL;
struct bed *bedList = NULL;
while ((row = sqlNextRow(sr)) != NULL)
{
AllocVar(bed);
bed->chrom = cloneString(row[0]);
bed->chromStart = atoi(row[1]);
bed->chromEnd = atoi(row[2]);
bed->thickStart = bed->chromStart;
bed->thickEnd = bed->chromEnd;
bed->name = cloneString(row[4]);
slAddHead(&bedList, bed);
}
sqlFreeResult(&sr);
return bedList;
}
int hmPixelCmpCount(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->count - b->count;
if (dif < 0)
return -1;
else if (dif > 0)
return 1;
else
return 0;
}
void drawBox(struct vGfx *vg, double x1, double x2, int pStart, int y, int h, Color valCol)
{
int x = round(x1) + pStart;
int w = round(x2) - round(x1);
if (w == 0)
w = 1;
vgBox(vg, x, y, w, h, valCol);
}
void addHmPixel(struct hmPixel **hmList, struct hash *pixelHash,
char mod, double x1, double x2, int pStart, int y, int h)
{
/* calculate starting position and width of box to draw */
int x = round(x1) + pStart;
int w = round(x2) - round(x1);
if (w == 0)
w = 1;
char pixelStr[32];
safef(pixelStr, sizeof(pixelStr), "%d,%c", x, mod);
struct hmPixel *hm;
struct hashEl *el = hashLookup(pixelHash, pixelStr);
if (!el)
{
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 = el->val;
if (w > hm->w)
hm->w = w;
hm->count += 1;
}
void drawTable(struct vGfx *vg, struct settings *settings,
char *db, char *tableName, char *profile, int height)
/* Draw chromosome graph on all chromosomes in layout at given
* y offset and height. */
{
if (!settings)
return;
/* get remote database connection */
struct sqlConnection *conn;
if (!profile)
conn = hAllocConn(db);
else
conn = hAllocConnProfile(profile, db);
struct bed *nb, *ghBed = NULL;
double pixPerBase = (double) settings->dataWidth / (double) settings->totalBases;
Color upShades[EXPR_DATA_SHADES];
static struct rgbColor zeroColor = {0, 0, 0}; // black
static struct rgbColor highColor = {255, 0, 0}; // red
vgMakeColorGradient(vg, &zeroColor, &highColor, EXPR_DATA_SHADES, upShades);
struct hash *pixelHash = hashNew(0);
struct hmPixel *hm, *hmList = NULL;
struct chromLay *cl;
for(cl = settings->layoutList; cl; cl = cl->next)
{
ghBed = getTable(conn, tableName, cl->chrom);
if (!ghBed)
continue;
struct bed *filterList =
bedFilterListInRange(ghBed, NULL, cl->chrom, cl->baseStart, cl->baseEnd);
int pStart = cl->pxStart;
int baseStart = cl->baseStart;
for(nb = filterList; nb; nb = nb->next)
{
int cStart = nb->chromStart;
int cEnd = nb->chromEnd;
double x1, x2;
int baseline = cStart - baseStart;
x1 = pixPerBase * (double) (cStart - baseStart);
x2 = pixPerBase * (double) (cEnd - baseStart);
if ( (x2 - x1) > 2.0)
{ /* width > two pixels, attempt to draw exon/utr structure
* similar to genome browser */
addHmPixel(&hmList, pixelHash, 'i', x1, x2, pStart, 4, 1);
int tStart = nb->thickStart - nb->chromStart;
int tEnd = nb->thickEnd - nb->chromStart;
int i, bStart, bEnd;
for (i = 0; i < nb->blockCount; i++)
{
bStart = nb->chromStarts[i];
bEnd = bStart + nb->blockSizes[i];
if (bStart < tStart)
{ /* 5' utr region */
x1 = pixPerBase * (double) (bStart + baseline);
if (bEnd < tStart)
x2 = pixPerBase * (double) (bEnd + baseline);
else
x2 = pixPerBase * (double) (tStart + baseline);
addHmPixel(&hmList, pixelHash, 'u', x1, x2, pStart, 2, 5);
if (bEnd > tEnd)
{ /* coding region inside single block */
x1 = pixPerBase * (double) (tStart + baseline);
x2 = pixPerBase * (double) (tEnd + baseline);
addHmPixel(&hmList, pixelHash, 'c', x1, x2, pStart, 0, 9);
x1 = pixPerBase * (double) (tEnd + baseline);
x2 = pixPerBase * (double) (bEnd + baseline);
addHmPixel(&hmList, pixelHash, 'u', x1, x2, pStart, 2, 5);
}
else if (bEnd > tStart)
{ /* coding region downstream of 5' utr */
x1 = pixPerBase * (double) (tStart + baseline);
x2 = pixPerBase * (double) (bEnd + baseline);
addHmPixel(&hmList, pixelHash, 'c', x1, x2, pStart, 0, 9);
}
}
else if (bEnd > tEnd)
{ /* 3' utr region */
if (bStart > tEnd)
x1 = pixPerBase * (double) (bStart + baseline);
else
x1 = pixPerBase * (double) (tEnd + baseline);
x2 = pixPerBase * (double) (bEnd + baseline);
addHmPixel(&hmList, pixelHash, 'u', x1, x2, pStart, 2, 5);
if (bStart < tEnd)
{ /* coding region upstream of 3' utr */
x1 = pixPerBase * (double) (bStart + baseline);
x2 = pixPerBase * (double) (tEnd + baseline);
addHmPixel(&hmList, pixelHash, 'c', x1, x2, pStart, 0, 9);
}
}
else
{ /* block is all coding region, no utr */
x1 = pixPerBase * (double) (bStart + baseline);
x2 = pixPerBase * (double) (bEnd + baseline);
addHmPixel(&hmList, pixelHash, 'c', x1, x2, pStart, 0, 9);
}
}
}
else
{ /* width < 2 pixels -- gene too dense, so just draw one large block */
x1 = pixPerBase * (double) (cStart - baseStart);
x2 = pixPerBase * (double) (cEnd - baseStart);
addHmPixel(&hmList, pixelHash, 'f', x1, x2, pStart, 0, 9);
}
}
bedFreeList(&ghBed);
}
hFreeConn(&conn);
if (!hmList) // nothing to draw.
return;
slSort(&hmList, hmPixelCmpCount);
struct hmPixel *last = slLastEl(hmList);
int maxVal = last->count;
vgSetClip(vg, settings->dataOffsetX, 0, settings->dataWidth, settings->height);
Color valCol;
for (hm = hmList; hm ; hm = hm->next)
{
double val = (double) hm->count / (double) maxVal;
int colorIndex = (int)(val * (EXPR_DATA_SHADES-1.0) );
/* 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;
valCol = upShades[colorIndex];
vgBox(vg, hm->x, hm->y + settings->dataOffsetY, hm->w, hm->h, valCol);
}
vgUnclip(vg);
}
struct spaceSaver *bedToSS(struct settings *settings, char *db,
char *tableName, boolean drawLabels)
{
if (!settings)
return NULL;
/* get remote database connection */
struct sqlConnection *conn = hAllocConn(db);
struct bed *nb, *ghBed = NULL;
struct chromLay *cl = settings->layoutList;
ghBed = getTable(conn, tableName, cl->chrom);
if (!ghBed)
return NULL;
struct bed *filterList =
bedFilterListInRange(ghBed, NULL, cl->chrom, cl->baseStart, cl->baseEnd);
struct spaceSaver *ss = spaceSaverNew(cl->pxStart, cl->pxStart + cl->pxWidth, 100);
double pixPerBase = (double) cl->pxWidth / (double) settings->totalBases;
int baseStart = cl->baseStart;
for(nb = filterList; nb; nb = nb->next)
{
int cStart = nb->chromStart;
int cEnd = nb->chromEnd;
int x1 = round(pixPerBase * (double) (cStart - baseStart)) + cl->pxStart;
int x2 = round(pixPerBase * (double) (cEnd - baseStart)) + cl->pxStart;
if (drawLabels)
x1 -= mgFontStringWidth(settings->font, nb->name) + 4;
spaceSaverAdd(ss, x1, x2, nb);
}
spaceSaverFinish(ss);
return ss;
}
struct spaceSaver *snpsToSS(struct settings *settings, char *db,
char *tableName, boolean drawLabels)
{
if (!settings)
return NULL;
/* get remote database connection */
struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, db);
struct bed *nb, *ghBed = NULL;
struct chromLay *cl = settings->layoutList;
ghBed = getSNPs(conn, tableName, cl->chrom, cl->baseStart, cl->baseEnd);
if (!ghBed)
return NULL;
struct spaceSaver *ss = spaceSaverNew(cl->pxStart, cl->pxStart + cl->pxWidth, 100);
double pixPerBase = (double) settings->dataWidth / (double) settings->totalBases;
int baseStart = cl->baseStart;
for(nb = ghBed; nb; nb = nb->next)
{
int cStart = nb->chromStart;
int cEnd = nb->chromEnd;
int x1 = round(pixPerBase * (double) (cStart - baseStart)) + cl->pxStart;
int x2 = round(pixPerBase * (double) (cEnd - baseStart)) + cl->pxStart;
if (drawLabels)
x1 -= mgFontStringWidth(settings->font, nb->name) + 4;
spaceSaverAdd(ss, x1, x2, nb);
}
spaceSaverFinish(ss);
return ss;
}
void drawBedTable(struct vGfx *vg, struct settings *settings,
struct spaceSaver *ss, boolean drawLabels)
/* Draw chromosome graph on all chromosomes in layout at given
* y offset and height. */
{
double pixPerBase = (double) settings->dataWidth / (double) settings->totalBases;
Color upShades[EXPR_DATA_SHADES];
static struct rgbColor zeroColor = {0, 0, 0}; // black
static struct rgbColor highColor = {255, 0, 0}; // red
vgMakeColorGradient(vg, &zeroColor, &highColor, EXPR_DATA_SHADES, upShades);
struct chromLay *cl = settings->layoutList;
int pStart = cl->pxStart;
int baseStart = cl->baseStart;
vgSetClip(vg, settings->dataOffsetX, 0, settings->dataWidth, settings->height);
struct spaceNode *sn;
for (sn = ss->nodeList; sn; sn = sn->next)
{
int row = sn->row * GENE_HEIGHT + settings->dataOffsetY;
struct bed *nb = (struct bed *) sn->val;
int cStart = nb->chromStart;
int cEnd = nb->chromEnd;
double x1, x2;
int baseline = cStart - baseStart;
x1 = pixPerBase * (double) (cStart - baseStart);
x2 = pixPerBase * (double) (cEnd - baseStart);
if (drawLabels)
vgTextRight(vg, 0, row, x1-2+pStart, settings->fontHeight,
MG_BLACK, settings->font, nb->name);
if ( (x2 - x1) > 2.0)
{ /* width > two pixels, attempt to draw exon/utr structure
* similar to genome browser */
drawBox(vg, x1, x2, pStart, 4 + row, 1, MG_BLACK);
int tStart = nb->thickStart - nb->chromStart;
int tEnd = nb->thickEnd - nb->chromStart;
int i, bStart, bEnd;
for (i = 0; i < nb->blockCount; i++)
{
bStart = nb->chromStarts[i];
bEnd = bStart + nb->blockSizes[i];
if (bStart < tStart)
{ /* 5' utr region */
x1 = pixPerBase * (double) (bStart + baseline);
if (bEnd < tStart)
x2 = pixPerBase * (double) (bEnd + baseline);
else
x2 = pixPerBase * (double) (tStart + baseline);
drawBox(vg, x1, x2, pStart, 2 + row, 5, MG_BLACK);
if (bEnd > tEnd)
{ /* coding region inside single block */
x1 = pixPerBase * (double) (tStart + baseline);
x2 = pixPerBase * (double) (tEnd + baseline);
drawBox(vg, x1, x2, pStart, 0 + row, 9, MG_BLACK);
x1 = pixPerBase * (double) (tEnd + baseline);
x2 = pixPerBase * (double) (bEnd + baseline);
drawBox(vg, x1, x2, pStart, 0 + row, 9, MG_BLACK);
}
else if (bEnd > tStart)
{ /* coding region downstream of 5' utr */
x1 = pixPerBase * (double) (tStart + baseline);
x2 = pixPerBase * (double) (bEnd + baseline);
drawBox(vg, x1, x2, pStart, 0 + row, 9, MG_BLACK);
}
}
else if (bEnd > tEnd)
{ /* 3' utr region */
if (bStart > tEnd)
x1 = pixPerBase * (double) (bStart + baseline);
else
x1 = pixPerBase * (double) (tEnd + baseline);
x2 = pixPerBase * (double) (bEnd + baseline);
drawBox(vg, x1, x2, pStart, 2 + row, 5, MG_BLACK);
if (bStart < tEnd)
{ /* coding region upstream of 3' utr */
x1 = pixPerBase * (double) (bStart + baseline);
x2 = pixPerBase * (double) (tEnd + baseline);
drawBox(vg, x1, x2, pStart, 0 + row, 9, MG_BLACK);
}
}
else
{ /* block is all coding region, no utr */
x1 = pixPerBase * (double) (bStart + baseline);
x2 = pixPerBase * (double) (bEnd + baseline);
drawBox(vg, x1, x2, pStart, 0 + row, 9, MG_BLACK);
}
}
}
else
{ /* width < 2 pixels -- gene too dense, so just draw one large block */
x1 = pixPerBase * (double) (cStart - baseStart);
x2 = pixPerBase * (double) (cEnd - baseStart);
drawBox(vg, x1, x2, pStart, 0 + row, 9, MG_BLACK);
}
}
vgUnclip(vg);
spaceSaverFree(&ss);
}
void drawSnpTable(struct vGfx *vg, struct settings *settings,
struct spaceSaver *ss, boolean drawLabels)
/* Draw chromosome graph on all chromosomes in layout at given
* y offset and height. */
{
if (!settings | !ss)
return;
double pixPerBase = (double) settings->dataWidth / (double) settings->totalBases;
Color upShades[EXPR_DATA_SHADES];
static struct rgbColor zeroColor = {0, 0, 0}; // black
static struct rgbColor highColor = {255, 0, 0}; // red
vgMakeColorGradient(vg, &zeroColor, &highColor, EXPR_DATA_SHADES, upShades);
struct chromLay *cl = settings->layoutList;
int pStart = cl->pxStart;
int baseStart = cl->baseStart;
struct spaceNode *sn;
for (sn = ss->nodeList; sn; sn = sn->next)
{
int row = sn->row * GENE_HEIGHT + settings->dataOffsetY;
struct bed *nb = (struct bed *) sn->val;
int cStart = nb->chromStart;
int cEnd = nb->chromEnd;
double x1, x2;
x1 = pixPerBase * (double) (cStart - baseStart);
x2 = pixPerBase * (double) (cEnd - baseStart);
if (drawLabels)
vgTextRight(vg, 0, row, x1-2, settings->fontHeight,
MG_BLACK, settings->font, nb->name);
x1 = pixPerBase * (double) (cStart - baseStart);
x2 = pixPerBase * (double) (cEnd - baseStart);
drawBox(vg, x1, x2, pStart, 0 + row, 9, MG_BLACK);
}
spaceSaverFree(&ss);
}
char *annotationWholeGenome(struct settings *settings, char *db, char *tableName)
{
int width = settings->width;
int height = settings->height; // was set to 10
if (width * height == 0)
return NULL;
struct tempName md5Tn;
char *strToHash = cartSettingsString(bbPrefix, "annotationGif");
trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
off_t size = fileSize(md5Tn.forCgi);
if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG)
{
struct hvGfx *vg = hvGfxOpenGif(width, height, md5Tn.forCgi, FALSE);
prepareImage(vg->vg, settings);
drawTable(vg->vg, settings, db, tableName, NULL, height);
hvGfxClose(&vg);
}
char *filename = replaceChars(md5Tn.forHtml, "..", "");
return filename;
}
char *annotationSingleGenome(struct settings *settings, char *db, char *tableName)
{
struct chromLay *cl = settings->layoutList;
boolean drawLabels = TRUE;
if (cl->baseEnd - cl->baseStart > 50000000)
drawLabels = FALSE;
struct spaceSaver *ss = bedToSS(settings, db, tableName, drawLabels);
int maxRow = -1;
struct spaceNode *sn;
for (sn = ss->nodeList; sn; sn = sn->next)
{
if (sn->row > maxRow)
maxRow = sn->row;
}
if (maxRow == -1)
return NULL;
settings->dataHeight = (maxRow + 1) * GENE_HEIGHT;
settings->height = settings->dataHeight + TITLE_HEIGHT;
struct tempName md5Tn;
char *strToHash = cartSettingsString(bbPrefix, "annotationGif");
trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
off_t size = fileSize(md5Tn.forCgi);
if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG)
{
struct hvGfx *vg = hvGfxOpenGif(settings->width, settings->height, md5Tn.forCgi, FALSE);
prepareImage(vg->vg, settings);
drawBedTable(vg->vg, settings, ss, drawLabels);
hvGfxClose(&vg);
}
char *filename = replaceChars(md5Tn.forHtml, "..", "");
return filename;
}
char *annotationGif(struct settings *settings, char *tableName)
/* Create genome GIF file and HT that includes it. */
{
if (!settings || !tableName)
return NULL;
char *db = "hg18";
if (slCount(settings->layoutList) > 1)
return annotationWholeGenome(settings, db, tableName);
else
return annotationSingleGenome(settings, db, tableName);
}
char *snpsWholeGenome(struct settings *settings, char *db, char *tableName)
{
int width = settings->width;
int height = 10; // was set to 10
if (width * height == 0)
return NULL;
struct tempName md5Tn;
char *strToHash = cartSettingsString(bbPrefix, "snpsGif");
trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
off_t size = fileSize(md5Tn.forCgi);
if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG)
{
struct hvGfx *vg = hvGfxOpenGif(width, height, md5Tn.forCgi, FALSE);
drawTable(vg->vg, settings, db, tableName, heatMapDbProfile, height);
hvGfxClose(&vg);
}
char *filename = replaceChars(md5Tn.forHtml, "..", "");
return filename;
}
char *snpsSingleGenome(struct settings *settings, char *db, char *tableName)
{
if (!settings)
return NULL;
struct chromLay *cl = settings->layoutList;
boolean drawLabels = TRUE;
if (cl->baseEnd - cl->baseStart > 500000)
drawLabels = FALSE;
struct spaceSaver *ss = snpsToSS(settings, db, tableName, drawLabels);
if (!ss)
return NULL;
int maxRow = -1;
struct spaceNode *sn;
for (sn = ss->nodeList; sn; sn = sn->next)
{
if (sn->row > maxRow)
maxRow = sn->row;
}
if (maxRow == -1)
return NULL;
settings->dataHeight = (maxRow + 1) * GENE_HEIGHT;
settings->height = settings->dataHeight + TITLE_HEIGHT;
struct tempName md5Tn;
char *strToHash = cartSettingsString(bbPrefix, "snpsGif");
trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
off_t size = fileSize(md5Tn.forCgi);
if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG)
{
struct hvGfx *vg = hvGfxOpenGif(settings->width, settings->height, md5Tn.forCgi, FALSE);
prepareImage(vg->vg, settings);
drawSnpTable(vg->vg, settings, ss, drawLabels);
hvGfxClose(&vg);
}
char *filename = replaceChars(md5Tn.forHtml, "..", "");
return filename;
}
char *snpsGif(struct settings *settings, char *tableName)
/* Create genome GIF file and HT that includes it. */
{
if (!settings || !tableName)
return NULL;
char *db = "pdata";
if (slCount(settings->layoutList) > 1)
return NULL; // don't do whole genome yet. snpsWholeGenome(settings, db, tableName);
else
{
struct chromLay *cl = settings->layoutList;
if (!cl)
return NULL;
if ((cl->baseEnd - cl->baseStart) > 5000000) // large region, compress.
return NULL; // snpsWholeGenome(settings, db, tableName);
else
return snpsSingleGenome(settings, db, tableName);
}
return NULL;
}
double basesInScale(double numBases)
{
double exp = floor(log(numBases)/log(10.0));
if (exp < 0.0)
exp = 0.0;
double divisor = pow(10.0, exp);
double rescaledBases = numBases / divisor;
double maxBases = rescaledBases / 2.0;
double bases = 0.0;
if (maxBases >= 5.0)
bases = 5.0 * divisor;
else if (maxBases >= 2.0)
bases = 2.0 * divisor;
else if (maxBases >= 1.0)
bases = 1.0 * divisor;
else if (maxBases >= 0.5)
bases = 0.5 * divisor;
else if (maxBases >= 0.2)
bases = 0.2 * divisor;
else if (maxBases >= 0.1)
bases = 0.1 * divisor;
bases = floor(bases);
return bases;
}
struct hash *getCytoBandHash(struct sqlConnection *conn)
{
if (!sqlTableExists(conn, "cytoBand"))
return NULL;
char **row;
char query[128];
safef(query, sizeof(query), "select * from cytoBand");
struct hash *cbHash = hashNew(0);
struct sqlResult *sr = sqlGetResult(conn, query);
while ((row = sqlNextRow(sr)) != NULL)
{
struct cytoBand *cb = cytoBandLoad(row);
struct hashEl *el = hashLookup(cbHash, cb->chrom);
struct cytoBand *cbList = NULL;
if (el)
cbList = el->val;
slAddTail(&cbList, cb);
if (!el) // add first member of list to hash
hashAdd(cbHash, cb->chrom, cbList);
}
sqlFreeResult(&sr);
return cbHash;
}
void drawIdeogram(struct hvGfx *hvg, struct settings *settings, struct chromLay *clList, int yOff)
{
if (!settings)
return;
struct vGfx *vg = hvg->vg;
struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, "hg18");
struct hash *cbHash = getCytoBandHash(conn);
hFreeConn(&conn);
if (!cbHash)
return;
/* Make grayscale for bands */
int maxShade = 9;
Color shadesOfGray[10+1];
static struct rgbColor black = {0, 0, 0};
static struct rgbColor white = {255, 255, 255};
vgMakeColorGradient(vg, &white, &black, 10, shadesOfGray);
shadesOfGray[10] = MG_RED; // to check for overflow, should never see red.
struct chromLay *cl;
vgSetClip(vg, settings->dataOffsetX, 0, settings->dataWidth, settings->height);
for (cl = clList; cl; cl = cl->next)
{
struct hashEl *el = hashLookup(cbHash, cl->chrom);
if (!el)
continue;
struct cytoBand *cb, *cbList = el->val;
int pixStart = cl->pxStart;
int pixEnd = cl->pxStart + cl->pxWidth;
int pixSize = cl->pxWidth;
unsigned bStart = cl->baseStart;
unsigned bEnd = cl->baseEnd;
unsigned bSize = bEnd - bStart;
double basesPerPixel = (double) bSize / (double) pixSize;
int xBuffer = 0; //Start = pixStart;
if (pixStart == 0)
xBuffer = 0;
/* Draw box around cytoband */
vgBox(vg, pixStart+1, yOff, pixSize-1, 1, MG_BLACK);
vgBox(vg, pixStart+1, yOff + CYTO_HEIGHT-1, pixSize-1, 1, MG_BLACK);
unsigned minBand = INT_MAX;
unsigned maxBand = 0;
for (cb = cbList; cb; cb = cb->next)
{
unsigned start = cb->chromStart; // our bases start at 1 TODO: FIX!
unsigned stop = cb->chromEnd;
if ((stop < bStart) || (start > bEnd)) // out of scope
continue;
int pStart = round(((double) start - (double) bStart) / basesPerPixel) + pixStart;
if (pStart < pixStart)
pStart = pixStart;
if (pStart <= 0)
pStart = 1;
int pEnd = round(((double) stop - (double) bStart) / basesPerPixel) + pixStart;
if (pEnd >= pixEnd)
pEnd = pixEnd;
Color col = hCytoBandColor(cb, hvg, FALSE, MG_BLACK, MG_WHITE,
shadesOfGray, maxShade);
vgBox(vg, pStart, yOff+1, (pEnd - pStart), CYTO_HEIGHT-2, col);
Color textCol = hvGfxContrastingColor(hvg, col);
char *pt = cb->name;
int fullWidth = mgFontStringWidth(settings->font, pt);
int shortWidth = mgFontStringWidth(settings->font, pt+1);
if (fullWidth < (pEnd - pStart))
vgTextCentered(vg, pStart, yOff+1, (pEnd - pStart), CYTO_HEIGHT-2,
textCol, settings->font, pt);
else if (shortWidth < (pEnd - pStart))
vgTextCentered(vg, pStart, yOff+1, (pEnd - pStart), CYTO_HEIGHT-2,
textCol, settings->font, pt+1);
if (start < minBand)
minBand = start;
if (stop > maxBand)
maxBand = stop;
}
/* If very near left/right edges, if we see the chromosome edge */
if (abs(minBand - bStart) < 100)
vgBox(vg, pixStart, yOff+1, 1, CYTO_HEIGHT-2, MG_BLACK);
if (abs(maxBand- bEnd) < 100)
vgBox(vg, pixEnd-1, yOff+1, 1, CYTO_HEIGHT-2, MG_BLACK);
for (cb = cbList; cb; cb = cb->next)
{
/* If centromere do some drawing. */
if(!sameString(cb->gieStain, "acen"))
continue;
int cenLeft, cenRight, cenTop, cenBottom;
/* Get the coordinates of the edges of the centromere. */
cenLeft = round(((double) cb->chromStart - (double) bStart) / basesPerPixel) + pixStart;
cenRight = round(((double) cb->next->chromEnd - (double) bStart)/basesPerPixel) + pixStart;
cenTop = yOff + CYTO_HEIGHT;
cenBottom = yOff;
if (cenLeft < 0 && cenRight < 0)
break;
if (cenLeft > pixEnd || cenRight < pixStart)
break;
/* Draw centromere itself. */
hCytoBandDrawCentromere(hvg, cenLeft, yOff, cenRight - cenLeft,
CYTO_HEIGHT, MG_WHITE, hCytoBandCentromereColor(hvg));
break;
}
}
}
Color getChromColor(struct vGfx *vg, char *chrom)
{
if (!chrom)
return MG_BLACK;
int r, g, b;
switch (chromToInt(chrom))
{
case 1:
r = CHROM_1_R;
g = CHROM_1_G;
b = CHROM_1_B;
break;
case 2:
r = CHROM_2_R;
g = CHROM_2_G;
b = CHROM_2_B;
break;
case 3:
r = CHROM_3_R;
g = CHROM_3_G;
b = CHROM_3_B;
break;
case 4:
r = CHROM_4_R;
g = CHROM_4_G;
b = CHROM_4_B;
break;
case 5:
r = CHROM_5_R;
g = CHROM_5_G;
b = CHROM_5_B;
break;
case 6:
r = CHROM_6_R;
g = CHROM_6_G;
b = CHROM_6_B;
break;
case 7:
r = CHROM_7_R;
g = CHROM_7_G;
b = CHROM_7_B;
break;
case 8:
r = CHROM_8_R;
g = CHROM_8_G;
b = CHROM_8_B;
break;
case 9:
r = CHROM_9_R;
g = CHROM_9_G;
b = CHROM_9_B;
break;
case 10:
r = CHROM_10_R;
g = CHROM_10_G;
b = CHROM_10_B;
break;
case 11:
r = CHROM_11_R;
g = CHROM_11_G;
b = CHROM_11_B;
break;
case 12:
r = CHROM_12_R;
g = CHROM_12_G;
b = CHROM_12_B;
break;
case 13:
r = CHROM_13_R;
g = CHROM_13_G;
b = CHROM_13_B;
break;
case 14:
r = CHROM_14_R;
g = CHROM_14_G;
b = CHROM_14_B;
break;
case 15:
r = CHROM_15_R;
g = CHROM_15_G;
b = CHROM_15_B;
break;
case 16:
r = CHROM_16_R;
g = CHROM_16_G;
b = CHROM_16_B;
break;
case 17:
r = CHROM_17_R;
g = CHROM_17_G;
b = CHROM_17_B;
break;
case 18:
r = CHROM_18_R;
g = CHROM_18_G;
b = CHROM_18_B;
break;
case 19:
r = CHROM_19_R;
g = CHROM_19_G;
b = CHROM_19_B;
break;
case 20:
r = CHROM_20_R;
g = CHROM_20_G;
b = CHROM_20_B;
break;
case 21:
r = CHROM_21_R;
g = CHROM_21_G;
b = CHROM_21_B;
break;
case 22:
r = CHROM_22_R;
g = CHROM_22_G;
b = CHROM_22_B;
break;
case 23:
r = CHROM_X_R;
g = CHROM_X_G;
b = CHROM_X_B;
break;
case 24:
r = CHROM_Y_R;
g = CHROM_Y_G;
b = CHROM_Y_B;
break;
default:
r = 0;
g = 0;
b = 0;
}
return vgFindColorIx(vg, r, g, b);
}
void drawColorIdeogram(struct hvGfx *hvg, struct settings *settings,
struct chromLay *clList, int yOff)
{
if (!settings)
return;
struct vGfx *vg = hvg->vg;
struct sqlConnection *conn = hAllocConnProfile(heatMapDbProfile, "hg18");
struct hash *cbHash = getCytoBandHash(conn);
hFreeConn(&conn);
if (!cbHash)
return;
/* Make grayscale for bands */
struct chromLay *cl;
vgSetClip(vg, 0, 0, settings->width, settings->height);
for (cl = clList; cl; cl = cl->next)
{
struct hashEl *el = hashLookup(cbHash, cl->chrom);
if (!el)
continue;
struct cytoBand *cb, *cbList = el->val;
int pixStart = cl->pxStart;
int pixEnd = cl->pxStart + cl->pxWidth;
int pixSize = cl->pxWidth;
unsigned bStart = cl->baseStart;
unsigned bEnd = cl->baseEnd;
unsigned bSize = bEnd - bStart;
double basesPerPixel = (double) bSize / (double) pixSize;
int xBuffer = 0; //Start = pixStart;
if (pixStart == 0)
xBuffer = 0;
Color col = getChromColor(vg, cl->chrom);
/* Draw box around cytoband */
vgBox(vg, pixStart+1, yOff, pixSize-1, 1, MG_BLACK);
vgBox(vg, pixStart+1, yOff + CYTO_HEIGHT-1, pixSize-1, 1, MG_BLACK);
unsigned minBand = INT_MAX;
unsigned maxBand = 0;
for (cb = cbList; cb; cb = cb->next)
{
unsigned start = cb->chromStart; // our bases start at 1 TODO: FIX!
unsigned stop = cb->chromEnd;
if ((stop < bStart) || (start > bEnd)) // out of scope
continue;
int pStart = round(((double) start - (double) bStart) / basesPerPixel) + pixStart;
if (pStart < pixStart)
pStart = pixStart;
if (pStart <= 0)
pStart = 1;
int pEnd = round(((double) stop - (double) bStart) / basesPerPixel) + pixStart;
if (pEnd >= pixEnd)
pEnd = pixEnd;
vgBox(vg, pStart, yOff+1, (pEnd - pStart), CYTO_HEIGHT-2, col);
if (start < minBand)
minBand = start;
if (stop > maxBand)
maxBand = stop;
}
/* If very near left/right edges, if we see the chromosome edge */
if (abs(minBand - bStart) < 100)
vgBox(vg, pixStart, yOff+1, 1, CYTO_HEIGHT-2, MG_BLACK);
if (abs(maxBand- bEnd) < 100)
vgBox(vg, pixEnd-1, yOff+1, 1, CYTO_HEIGHT-2, MG_BLACK);
for (cb = cbList; cb; cb = cb->next)
{
/* If centromere do some drawing. */
if(!sameString(cb->gieStain, "acen"))
continue;
int cenLeft, cenRight, cenTop, cenBottom;
/* Get the coordinates of the edges of the centromere. */
cenLeft = round(((double) cb->chromStart - (double) bStart) / basesPerPixel) + pixStart;
cenRight = round(((double) cb->next->chromEnd - (double) bStart)/basesPerPixel) + pixStart;
cenTop = yOff + CYTO_HEIGHT;
cenBottom = yOff;
if (cenLeft < 0 && cenRight < 0)
break;
if (cenLeft > pixEnd || cenRight < pixStart)
break;
/* Draw centromere itself. */
hCytoBandDrawCentromere(hvg, cenLeft, yOff, cenRight - cenLeft,
CYTO_HEIGHT, MG_WHITE, hCytoBandCentromereColor(hvg));
break;
}
char *name = cl->chrom;
Color textCol = MG_BLACK; //hvGfxContrastingColor(hvg, col);
vgTextCentered(vg, pixStart,
settings->dataOffsetY + DEFAULT_LABEL_HEIGHT - settings->fontHeight,
(pixEnd - pixStart), DEFAULT_LABEL_HEIGHT, textCol, settings->font, name);
}
}
void drawChromLabels(struct vGfx *vg, struct settings *settings, struct chromLay *clList)
{
if (!settings)
return;
int width = settings->width;
int labelHeight = DEFAULT_LABEL_HEIGHT;
struct chromLay *cl;
vgSetClip(vg, settings->dataOffsetX, 0, settings->dataWidth, settings->height);
+int yOff = CYTO_HEIGHT;
boolean singleChrom = FALSE;
if (slCount(settings->layoutList) == 1)
singleChrom = TRUE;
for(cl = clList; cl; cl = cl->next)
{
int chromWidth = cl->pxWidth;
int leftX = cl->pxStart;
if (leftX < 0)
leftX = 0;
int rightX = leftX + chromWidth;
if (rightX > width - 1)
rightX = width - 1;
- vgBox(vg, leftX, 0, 1, labelHeight, MG_GRAY);
- vgBox(vg, rightX, 0, 1, labelHeight, MG_GRAY);
+ vgBox(vg, leftX, yOff, 1, labelHeight, MG_GRAY);
+ vgBox(vg, rightX, yOff, 1, labelHeight, MG_GRAY);
int strWidth;
if (!singleChrom)
{
char *pt = cl->chrom;
strWidth = mgFontStringWidth(settings->font, pt);
if (strWidth < chromWidth)
- vgTextCentered(vg, cl->pxStart, 0, chromWidth, labelHeight,
+ vgTextCentered(vg, cl->pxStart, yOff, chromWidth, labelHeight,
MG_BLACK, settings->font, pt);
else
{
pt = cl->chrom + 3;
strWidth = mgFontStringWidth(settings->font, pt);
if (strWidth < chromWidth)
- vgTextCentered(vg, cl->pxStart, 0, chromWidth, labelHeight,
+ vgTextCentered(vg, cl->pxStart, yOff, chromWidth, labelHeight,
MG_BLACK, settings->font, pt);
}
}
else
{
char str[128];
unsigned long x;
- int y = labelHeight/2 - settings->fontHeight/4;
+ int y = labelHeight/2 - settings->fontHeight/4 + yOff;
vgText(vg, cl->pxStart + 3, y, MG_BLACK, settings->font, cl->chrom);
int minX = cl->pxStart + mgFontStringWidth(settings->font, cl->chrom);
double numBases = cl->baseEnd - cl->baseStart;
if (numBases == 1)
{
safef(str, sizeof(str), "%lu", (unsigned long int) cl->baseStart);
strWidth = mgFontStringWidth(settings->font, str);
vgText(vg, rightX - strWidth, y, MG_BLACK, settings->font, str);
continue;
}
unsigned long tickBases = (unsigned long) ceil(0.5 * basesInScale(settings->totalBases));
int tickPixels = round(tickBases/settings->totalBases * (double) chromWidth);
if (tickBases == 0)
continue;
int buffer = 10;
unsigned long firstX = floor(cl->baseStart / tickBases) * tickBases;
for (x = firstX; x <= cl->baseEnd; x += tickBases)
{
int tickX = round((x - cl->baseStart + 1) / settings->totalBases * (double) chromWidth)
+ cl->pxStart;
if (tickX <= minX || tickX > rightX)
continue;
- int tickY = labelHeight/2 - settings->fontHeight/2 - 1;
+ int tickY = labelHeight/2 - settings->fontHeight/2 - 1 + yOff;
vgBox(vg, tickX, tickY, 1, settings->fontHeight+1, MG_BLACK);
safef(str, sizeof(str), "%lu", x);
strWidth = mgFontStringWidth(settings->font, str);
if ((strWidth < tickPixels - buffer) && (tickX - strWidth > minX + buffer))
vgText(vg, tickX - strWidth - 1, y, MG_BLACK, settings->font, str);
}
}
}
vgUnclip(vg);
}
char *genomeLabelsGif(struct sqlConnection *conn, struct settings *settings)
/* Create genome label GIF file and HT that includes it. */
{
if (!settings)
return NULL;
struct hvGfx *vg;
int width = settings->width;
int height = settings->height;
if (width * height == 0)
return NULL;
struct tempName md5Tn;
char *strToHash = cartSettingsString(bbPrefix, "genomeLabelsGif");
trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
off_t size = fileSize(md5Tn.forCgi);
if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG)
{
vg = hvGfxOpenGif(width, height, md5Tn.forCgi, FALSE);
prepareImage(vg->vg, settings);
drawChromLabels(vg->vg, settings, settings->layoutList);
- drawIdeogram(vg, settings, settings->layoutList, DEFAULT_LABEL_HEIGHT);
+ drawIdeogram(vg, settings, settings->layoutList, 0);
hvGfxClose(&vg);
}
char *filename = replaceChars(md5Tn.forHtml, "..", "");
return filename;
}
double basesInScaleAndString(double numBases, char **scaleStr)
{
double bases = basesInScale(numBases);
double divisor = 1.0;
char *suffix = NULL;
if (bases == 1.0)
{
divisor = 1.0;
suffix = "base";
}
else if (bases < 10.0)
{
divisor = 1.0;
suffix = "bases";
}
else if (bases < 1000.0)
{
divisor = 1.0;
suffix = "bp";
}
else if (bases < 1000000.0)
{
divisor = 1000.0;
suffix = "Kb";
}
else if (bases < 1000000000.0)
{
divisor = 1000000.0;
suffix = "Mb";
}
else if (bases < 1000000000000.0)
{
divisor = 1000000000.0;
suffix = "Gb";
}
char str[128];
safef(str, sizeof(str), "%d %s", (int) floor(bases/divisor), suffix);
*scaleStr = cloneString(str);
return bases;
}
void drawGenomeScale(struct vGfx *vg, struct settings *settings,
int width, int height, double numBases)
{
char *scaleStr;
double scaleBases = basesInScaleAndString(numBases, &scaleStr);
if (!scaleStr)
return;
vgSetClip(vg, settings->dataOffsetX, 0, settings->dataWidth, height);
int midY = height / 2;
int buffer = 3;
int scaleWidth = round(((double) width) * scaleBases / numBases);
-int startX = round((double)width / 2.0) - scaleWidth / 2;
+int startX = round((double)width / 2.0) - scaleWidth / 2 + settings->dataOffsetX;
int dashTop = (height - settings->fontHeight)/2;
/* draw line on left/right side */
vgBox(vg, startX, dashTop, 1, settings->fontHeight, MG_BLACK);
vgBox(vg, startX+scaleWidth, dashTop, 1, settings->fontHeight, MG_BLACK);
/* draw line connecting */
vgBox(vg, startX, midY, scaleWidth, 1, MG_BLACK);
int textY = midY - settings->fontHeight/2;
vgTextRight(vg, 0, textY, startX - buffer, settings->fontHeight,
MG_BLACK, settings->font, scaleStr);
vgUnclip(vg);
}
char *genomeScaleGif(struct settings *settings)
/* Create genome GIF file and HT that includes it. */
{
if (!settings || settings->totalBases < 1)
return NULL;
int totalW = settings->width;
int totalH = settings->height;
struct tempName md5Tn;
char *strToHash = cartSettingsString(bbPrefix, "genomeScaleGif");
trashDirMD5File(&md5Tn, "hgh", ".gif", strToHash);
off_t size = fileSize(md5Tn.forCgi);
if (!fileExists(md5Tn.forCgi) || (size == 0) || DEBUG)
{
struct hvGfx *vg = hvGfxOpenGif(totalW, totalH, md5Tn.forCgi, FALSE);
prepareImage(vg->vg, settings);
drawGenomeScale(vg->vg, settings, totalW, totalH, settings->totalBases);
hvGfxClose(&vg);
}
char *filename = replaceChars(md5Tn.forHtml, "..", "");
return filename;
}