Bar chart tool for use with taxonomy tools. Requires gnuplot and gnuplot.py. Still needs interface testing

This commit is contained in:
Anton Nekrutenko
2008-03-10 21:57:32 +00:00
parent d4c676d085
commit 85cebb59b6
4 changed files with 453 additions and 83 deletions
Binary file not shown.

After

Width:  |  Height:  |  Size: 13 KiB

+144
View File
@@ -0,0 +1,144 @@
#!/usr/bin/env python2.4
"""
histogram_gnuplot.py <datafile> <xtic column> <column_list> <title> <ylabel> <yrange_min> <yrange_max> <grath_file>
a generic histogram builder based on gnuplot backend
data_file - tab delimited file with data
xtic_column - column containing labels for x ticks [integer, 0 means no ticks]
column_list - comma separated list of columns to plot
title - title for the entire histrogram
ylabel - y axis label
yrange_max - minimal value at the y axis (integer)
yrange_max - maximal value at the y_axis (integer)
to set yrange to auto assign 0 to yrange_mmin and yrange_max
graph_file - file to write histogram image to
pdf_size - as X,Y pair in inches (e.g., 11,8 or 8,11 etc.)
This tool required gnuplot and gnuplot.py
anton nekrutenko | anton@bx.psu.edu
"""
from Numeric import *
import Gnuplot, Gnuplot.funcutils
import sys, string, tempfile, os
def stop_err(msg):
sys.stderr.write(msg)
sys.exit()
def main(tmpFileName):
skipped_lines_count = 0
skipped_lines_index = []
gf = open(tmpFileName, 'w')
try:
in_file = open( sys.argv[1], 'r' )
xtic = int( sys.argv[2] )
col_list = string.split( sys.argv[3],"," )
title = 'set title "' + sys.argv[4] + '"'
ylabel = 'set ylabel "' + sys.argv[5] + '"'
ymin = sys.argv[6]
ymax = sys.argv[7]
img_file = sys.argv[8]
pdf_size = sys.argv[9]
except:
stop_err("Check arguments\n")
try:
int( col_list[0] )
except:
stop_err('You forgot to set columns for plotting\n')
for i, line in enumerate( in_file ):
valid = True
line = line.rstrip('\r\n')
if line and not line.startswith( '#' ):
row = []
try:
fields = line.split( '\t' )
for col in col_list:
row.append( str( float( fields[int(col)-1] ) ) )
except:
valid = False
skipped_lines_count += 1
skipped_lines_index.append(i)
else:
valid = False
skipped_lines_count += 1
skipped_lines_index.append(i)
if valid and xtic > 0:
row.append( fields[xtic-1] )
elif valid and xtic == 0:
row.append( str( i ) )
if valid:
gf.write( '\t'.join( row ) )
gf.write( '\n' )
if skipped_lines_count < i:
#prepare 'using' clause of plot statement
g_plot_command = ' ';
#set the first column
if xtic > 0:
g_plot_command = "'%s' using 1:xticlabels(%s) ti 'Column %s', " % ( tmpFileName, str( len( row ) ), col_list[0] )
#g_plot_command = "'%s' using 1:xticlabels(%s), " % ( tmpFileName, str( len( row ) ) )
else:
g_plot_command = "'%s' using 1, " % ( tmpFileName )
#set subsequent columns
for i in range(1,len(col_list)):
g_plot_command += "'%s' using %s t 'Column %s', " % ( tmpFileName, str(i+1), col_list[i] )
#g_plot_command += "'%s' using %s, " % ( tmpFileName, str(i+1) )
g_plot_command = g_plot_command.rstrip(', ')
yrange = 'set yrange [' + ymin + ":" + ymax + ']'
try:
g = Gnuplot.Gnuplot()
g('reset')
g('set boxwidth 0.9 absolute')
g('set style fill solid 1.00 border -1')
g('set style histogram clustered gap 5 title offset character 0, 0, 0')
g('set xtics border in scale 1,0.5 nomirror rotate by 90 offset character 0, 0, 0')
g('set key invert reverse Left outside')
if xtic == 0: g('unset xtics')
g(title)
g(ylabel)
g_term = 'set terminal pdf size ' + pdf_size
g(g_term)
g_out = 'set output "' + img_file + '"'
if ymin != ymax:
g(yrange)
g(g_out)
g('set style data histograms')
g.plot(g_plot_command)
except:
stop_err("Gnuplot error: Data cannot be plotted")
else:
sys.stderr.write('Columns %s of your dataset do not contain valid numeric data' %sys.argv[3] )
if skipped_lines_count > 0:
sys.stderr.write('You dataset contains %d invalid lines starting with line#%d\n' % ( skipped_lines_count, skipped_lines_index[0] ) )
if __name__ == "__main__":
gp_data_file = tempfile.NamedTemporaryFile('w')
#gp_f, gp_data_file = tempfile.mkstemp(suffix="gp", text=True)
main(gp_data_file.name)
+51
View File
@@ -0,0 +1,51 @@
<tool id="barchart_gnuplot" name="Bar chart">
<description>for multiple columns</description>
<command interpreter="python2.4">
#if $xtic.userSpecified == "Yes": #bar_chart.py $input $xtic.xticColumn $colList "$title" "$ylabel" $ymin $ymax $out_file1 "$pdf_size"
#else: #bar_chart.py $input 0 $colList "$title" "$ylabel" $ymin $ymax $out_file1 "$pdf_size"
#end if
</command>
<inputs>
<param name="input" type="data" format="tabular" label="Dataset" help="Query missing? See TIP below"/>
<conditional name="xtic">
<param name="userSpecified" type="select" label="Use X Tick labels?" help="see example below">
<option value="Yes">Yes</option>
<option value="No">No</option>
</param>
<when value="Yes">
<param name="xticColumn" type="data_column" data_ref="input" numerical="False" label="Use this column for X Tick labels" />
</when>
<when value="No">
</when>
</conditional>
<param name="colList" label="Numerical columns" type="data_column" numerical="True" multiple="True" data_ref="input" help="Multi-select list - hold the appropriate key while clicking to select multiple columns" />
<param name="title" type="text" size="30" value="Bar Chart" label="Plot title"/>
<param name="ylabel" type="text" size="30" value="V1" label="Label for Y axis"/>
<param name="ymin" type="integer" size="4" value="0" label="Minimal value on Y axis" help="set to 0 for autoscaling"/>
<param name="ymax" type="integer" size="4" value="0" label="Maximal value on Y axis" help="set to 0 for autoscaling"/>
<param name="pdf_size" type="select" label="Choose chart size (inches)" help="inch = 2.54 cm">
<option value="11,8">Normal: 11 by 8</option>
<option value="5,3">Small: 5 by 3</option>
<option value="17,11">Large: 17 by 11</option>
<option value="8,11">Normal Flipped: 8 by 11</option>
<option value="3,5">Small Flipped: 3 by 5</option>
<option value="11,17">Large Flipped: 11 by 17</option>
</param>
</inputs>
<outputs>
<data format="pdf" name="out_file1" />
</outputs>
<help>
**What it does**
This tool builds a bar chart on one or more columns. Suppose you have dataset like this one::
Gene1 10 15␍ Gene2 20 14␍ Gene3 67 45␍ Gene4 55 12
Graphing columns 2 and 3 while using column 1 for X Tick Labels will produce the following plot:
.. image:: ../static/images/bar_chart.png
</help>
</tool>
+258 -83
View File
@@ -10,6 +10,7 @@
#define DEFAULT_STRING_ALLOC 16L
#define NUMBER_OF_FIELDS 24
#define NUMBER_OF_TAX_FIELDS 22
#define AVL_THRESHOLD 64
char *rankLabels [NUMBER_OF_TAX_FIELDS] =
{"root" ,
@@ -44,6 +45,8 @@ void check_pointer (void*);
char validTaxonNameChar[256];
long currentLineID = 1;
/*---------------------------------------------------------------------------------------------------- */
struct avl_table * idTagAVL = NULL,
@@ -70,6 +73,24 @@ struct vector
vaLength;
};
/*---------------------------------------------------------------------------------------------------- */
void reportError (char * theMessage)
{
fprintf (stderr, "\nERROR: %s\n", theMessage);
exit (1);
}
/*---------------------------------------------------------------------------------------------------- */
void reportErrorLine (char * theMessage, long lineID)
{
fprintf (stderr, "SKIPPED line %d: %s\n", lineID, theMessage);
}
/*---------------------------------------------------------------------------------------------------- */
struct bufferedString *allocateNewString (void)
@@ -238,9 +259,11 @@ struct treeNode
struct treeNode * parent;
long startIndex,
hitCount,
length;
length
,beenhere;
struct vector * children;
struct avl_table* cachedChildren;
} *globalTreeRoot;
@@ -251,77 +274,201 @@ struct treeNode * allocateNewTreeNode (void)
struct treeNode *newN = (struct treeNode*)malloc (sizeof (struct treeNode));
check_pointer (newN);
newN->parent = NULL;
newN->startIndex = 0;
newN->length = 0;
newN->hitCount = 0;
newN->children = allocateNewVector();
newN->startIndex = 0;
newN->length = 0;
newN->hitCount = 0;
newN->beenhere = 0;
newN->cachedChildren = NULL;
newN->children = allocateNewVector();
check_pointer (newN->children);
return newN;
}
/*---------------------------------------------------------------------------------------------------- */
struct treeNode * addAChild (struct treeNode* p, struct treeNode * c)
int compare_tree_nodes (const void *avl_a, const void *avl_b, void * xtra)
{
struct treeNode * c2 = NULL;
long i = 0;
for (; i<p->children->vLength; i++)
if (((struct treeNode**)p->children->vData)[i]->startIndex == c->startIndex)
break;
if (p->children->vLength == i)
{
if (c->parent && c->parent != p)
{
c2 = allocateNewTreeNode();
c2->startIndex = c->startIndex;
c2->length = c->length;
c2->hitCount = 1;
c = c2;
}
appendValueToVector (p->children, (long)c);
c->parent = p;
}
return c;
long t1 = ((struct treeNode*)avl_a)->startIndex,
t2 = ((struct treeNode*)avl_b)->startIndex;
if (t1 < t2)
return -1;
if (t1 > t2)
return 1;
return 0;
}
/*---------------------------------------------------------------------------------------------------- */
void destroyTreeNode (struct treeNode* n)
{
if (n->cachedChildren)
avl_destroy (n->cachedChildren, NULL);
free (n->children);
free (n);
}
/*---------------------------------------------------------------------------------------------------- */
void printTreeNode (FILE* f, struct treeNode* n)
{
long i = 0;
fprintf (f, "Node name: ");
if (n->startIndex>=0)
{
for (;i < n->length; i++)
fputc (globalNameBuffer->sData[n->startIndex + i], f);
}
else
fprintf (f, " empty node");
fprintf (f, "\nHit count %d (backup %d)\n", n->hitCount, n->beenhere);
}
/*---------------------------------------------------------------------------------------------------- */
void traverseCheck (struct treeNode* n, char first)
{
long i = 0;
if (!first && n->startIndex >= 0 && n->beenhere != n->hitCount)
{
printTreeNode (stderr, n);
//reportError ("DEATH AND DECAY, BIZNATCH!\n");
}
for (i = 0; i<n->children->vLength; i++)
traverseCheck (((struct treeNode**)n->children->vData)[i], first);
if (first)
{
if (n->children->vLength == 0)
n->beenhere = n->hitCount;
if (n->parent)
n->parent->beenhere += n->beenhere;
}
}
/*---------------------------------------------------------------------------------------------------- */
void traverseTree (FILE *summaryFile, struct treeNode* n, long maxDepth, long currentDepth)
{
long i = 0;
char c;
if (currentDepth <= maxDepth)
{
if (n->children->vLength && maxDepth > currentDepth)
{
fputc ('(',summaryFile);
for (; i<n->children->vLength; i++)
if (n->startIndex < 0)
{
traverseTree (summaryFile, ((struct treeNode**)n->children->vData)[i], maxDepth, currentDepth+1);
if (i<n->children->vLength-1)
fputc (',',summaryFile);
for (i=0; i>n->startIndex && currentDepth - i < maxDepth; i--)
fputc ('(', summaryFile);
if (i == n->startIndex)
i = 0;
}
else
fputc ('(',summaryFile);
if (i==0)
for (; i<n->children->vLength; i++)
{
traverseTree (summaryFile, ((struct treeNode**)n->children->vData)[i], maxDepth, currentDepth+1);
if (i<n->children->vLength-1)
fputc (',',summaryFile);
}
if (n->startIndex < 0)
for (i=0; i>n->startIndex+1 && currentDepth - i < maxDepth; i--)
fprintf (summaryFile,")n:%d", n->hitCount);
fputc (')',summaryFile);
}
for (i=n->startIndex; i<n->startIndex+n->length;i++)
if (n->startIndex>= 0)
{
c = globalNameBuffer->sData[i];
if (validTaxonNameChar [c])
fputc (c,summaryFile);
else
fputc ('_',summaryFile);
for (i=n->startIndex; i<n->startIndex+n->length;i++)
{
c = globalNameBuffer->sData[i];
if (validTaxonNameChar [c])
fputc (c,summaryFile);
else
fputc ('_',summaryFile);
}
fprintf (summaryFile,":%d", n->hitCount);
}
fprintf (summaryFile,":%d", n->hitCount);
else
fprintf (summaryFile, "n:%d", n->hitCount);
}
}
/*---------------------------------------------------------------------------------------------------- */
struct treeNode * addAChild (struct treeNode* p, struct treeNode * c, char killIfD)
{
struct treeNode * c2 = NULL;
long i = 0;
char addNode = 0;
if (p->cachedChildren)
{
if ((c2 = ((struct treeNode*)avl_find (p->cachedChildren, c))) == NULL)
addNode = 1;
/*else
{
if (c2->startIndex > globalNameBuffer->sLength || c2->startIndex < 0)
printf ("Reject node add %x %d %d\n", c2, c2->startIndex, c->startIndex);
}*/
}
else
{
for (; i<p->children->vLength; i++)
if (((struct treeNode**)p->children->vData)[i]->startIndex == c->startIndex)
break;
addNode = (p->children->vLength == i);
}
if (addNode)
{
if (c->parent && c->parent != p)
{
c2 = allocateNewTreeNode();
c2->startIndex = c->startIndex;
c2->length = c->length;
c2->hitCount = 1;
//if (p->cachedChildren)
// fprintf (stdout, "%x %x\n", p, c->parent);
c = c2;
}
c->parent = p;
appendValueToVector (p->children, (long)c);
if (p->children->vLength > AVL_THRESHOLD)
if (p->cachedChildren == NULL)
{
//fprintf (stdout, "Switch %d\n", currentLineID);
// switch over the the avl representation of nodes
p->cachedChildren = avl_create (compare_tree_nodes, NULL, NULL);
for (i=0; i<p->children->vLength; i++)
avl_probe(p->cachedChildren, ((struct node**)p->children->vData)[i]);
}
else
avl_probe (p->cachedChildren, c);
}
else
{
if (killIfD)
destroyTreeNode (c);
if (c2)
return c2;
else
c = ((struct treeNode**)p->children->vData)[i];
}
return c;
}
/*---------------------------------------------------------------------------------------------------- */
/*---------------------------------------------------------------------------------------------------- */
@@ -359,6 +506,7 @@ int compare_id_tags (const void *avl_a, const void *avl_b, void * xtra)
return 0;
}
/*---------------------------------------------------------------------------------------------------- */
/*---------------------------------------------------------------------------------------------------- */
@@ -473,22 +621,6 @@ void destroy_string (struct bufferedString* aStr)
}
/*---------------------------------------------------------------------------------------------------- */
void reportError (char * theMessage)
{
fprintf (stderr, "\nERROR: %s\n", theMessage);
exit (1);
}
/*---------------------------------------------------------------------------------------------------- */
void reportErrorLine (char * theMessage, long lineID)
{
fprintf (stderr, "SKIPPED line %d: %s\n", lineID, theMessage);
}
/*---------------------------------------------------------------------------------------------------- */
int main (int argc, const char * argv[])
@@ -507,23 +639,24 @@ int main (int argc, const char * argv[])
struct treeNode *currentParent = NULL,
*currentNode;
char automatonState = 0,
currentField = 0,
currentChar = 0,
showEmptyNodes = 0;
long currentLineID = 1,
expectedFields = NUMBER_OF_FIELDS,
long expectedFields = NUMBER_OF_FIELDS,
indexer,
indexer2,
indexer3,
maxTreeLevel = 0;
maxTreeLevel = 0,
setEOF = 0,
nRunCounter = 0;
FILE * treeFile = NULL,
* summaryFile = NULL,
* inFile = NULL;
globalNameBuffer = allocateNewString();
globalTreeRoot = allocateNewTreeNode();
currentBuffers = (struct bufferedString**)malloc (expectedFields*sizeof (struct bufferedString*));
@@ -571,7 +704,7 @@ int main (int argc, const char * argv[])
currentChar = fgetc(inFile);
currentField = 0;
while (!feof(inFile))
while (setEOF < 2)
{
switch (automatonState)
{
@@ -591,33 +724,35 @@ int main (int argc, const char * argv[])
break;
case 1: /* reading sequence ID */
if (isalnum(currentChar))
appendCharacterToString(currentBuffers[currentField],toupper(currentChar));
if (currentChar == '\t')
{
automatonState = 2;
continue;
}
else
if (currentChar == '\t')
{
automatonState = 2;
continue;
}
else
{
if (currentChar == '\n' || currentChar == '\r')
{
reportErrorLine ("Expected a tab following the gid",currentLineID);
automatonState = 6;
continue;
}
break;
else
appendCharacterToString(currentBuffers[currentField],toupper(currentChar));
}
break;
case 2: /* looking for a \t or a \n|\r*/
if (currentChar == '\t')
{
automatonState = 3;
currentField ++;
if (currentField == expectedFields)
{
reportErrorLine ("Too many fields",currentLineID);
automatonState = 6;
continue;
}
if (currentField < expectedFields)
automatonState = 3;
//reportErrorLine ("Too many fields",currentLineID);
//automatonState = 6;
//continue;
}
else
if (currentChar == '\n' || currentChar == '\r')
@@ -628,14 +763,19 @@ int main (int argc, const char * argv[])
aTag->taxID = atoi (currentBuffers[1]->sData);
aTag->hit_count = 1;
aTag2 = *(struct storedIDTag**)avl_probe(idTagAVL, aTag);
//fprintf (stdout, "%d %d %x %x\n", currentLineID,aTag->taxID,aTag,aTag2);
if (aTag == aTag2) // new taxID
{
// process fields
currentParent = globalTreeRoot;
nRunCounter = 0;
for (indexer = NUMBER_OF_FIELDS-NUMBER_OF_TAX_FIELDS; indexer < NUMBER_OF_FIELDS; indexer++)
{
indexer2 = strlen(currentBuffers[indexer]->sData);
if ((currentBuffers[indexer])->sData[0] != 'n' || indexer2 > 1 || showEmptyNodes && indexer2 > 0) // not 'n'
//fprintf (stdout, "%d %d\n", indexer, nRunCounter,indexer2);
if ((currentBuffers[indexer])->sData[0] == 'n' && indexer2 == 1)
nRunCounter++;
else
{
indexer3 = globalNameBuffer->sLength;
appendRangeToString (globalNameBuffer,currentBuffers[indexer],0,indexer2-1);
@@ -654,6 +794,7 @@ int main (int argc, const char * argv[])
if (sTag == sTag3) // new node
{
//fprintf (stderr, "Add node level %d %d\n",indexer-(NUMBER_OF_FIELDS-NUMBER_OF_TAX_FIELDS), sTag2->startIndex);
sTag3 = sTag;
sTag3->refNode = allocateNewTreeNode();
sTag3->refNode->length = sTag2->length;
sTag3->refNode->startIndex = sTag2->startIndex;
@@ -665,14 +806,38 @@ int main (int argc, const char * argv[])
sTag3->refNode->hitCount++;
if (showEmptyNodes && nRunCounter>0)
{
currentNode = allocateNewTreeNode();
currentNode->startIndex = -nRunCounter;
currentNode->hitCount = 1;
if (currentParent)
currentParent = addAChild (currentParent,currentNode,1);
else
reportError ("Attempting to attach an empty node a null parent.");
}
if (currentParent)
currentParent = addAChild (currentParent,sTag3->refNode);
currentParent = addAChild (currentParent,sTag3->refNode,0);
else
currentParent = sTag3->refNode;
nRunCounter = 0;
}
}
if (nRunCounter>0)
{
currentNode = allocateNewTreeNode();
currentNode->startIndex = -nRunCounter;
currentNode->hitCount = 1;
if (currentParent)
currentParent = addAChild (currentParent,currentNode,1);
else
reportError ("Attempting to attach an empty node a null parent.");
}
aTag->tNode = currentParent;
aTag = allocateIDTag();
}
else
{
@@ -683,17 +848,20 @@ int main (int argc, const char * argv[])
currentParent->hitCount++;
currentParent = currentParent->parent;
}
//printf ("Duplicate tag %d %d\n", aTag2->hit_count, aTag2->taxID);
}
//traverseCheck (globalTreeRoot);
automatonState = 5;
continue;
}
else
{
reportErrorLine ("Expected a tab following a field",currentLineID);
automatonState = 6;
continue;
if (currentField < expectedFields)
{
reportErrorLine ("Expected a tab following a field",currentLineID);
automatonState = 6;
continue;
}
}
break;
@@ -746,8 +914,13 @@ int main (int argc, const char * argv[])
break;
}
currentChar = fgetc(inFile);
if (feof (inFile))
{
setEOF ++;
currentChar = '\n';
}
}
fprintf (stderr, "Read %d unique taxIDs\n", avl_count (idTagAVL));
fclose (inFile);
@@ -781,6 +954,8 @@ int main (int argc, const char * argv[])
traverseTree (treeFile,((struct treeNode**)globalTreeRoot->children->vData)[0],(maxTreeLevel>0)?(maxTreeLevel+1):0xfffffL,0);
fclose (treeFile);
traverseCheck (((struct treeNode**)globalTreeRoot->children->vData)[0], 1);
traverseCheck (((struct treeNode**)globalTreeRoot->children->vData)[0], 0);
return 0;
}