Commits

James Taylor  committed 89b5be5

bx.motif: Drawing sequence logos from count/frequency matrices in EPS format.

  • Participants
  • Parent commits 726b490

Comments (0)

Files changed (2)

File lib/bx/motif/logo/__init__.py

+import pkg_resources
+from StringIO import StringIO
+from string import Template
+from numpy import *
+
+PAD = 2
+
+# Colors from rgb.txt, 
+DNA_DEFAULT_COLORS = {
+    'A': "0.00 1.00 0.00", # green
+    'C': "0.00 0.00 1.00", # blue
+    'G': "1.00 0.65 0.00", # orange red
+    'T': "1.00 0.00 0.00"  # red
+}
+
+# Template is adapted from Jim Kent's lib/dnaMotif.pss to support aritrary
+# alphabets.
+TEMPLATE = "template.ps"
+
+def freqs_to_heights( matrix ):
+    """
+    Calculate logo height using the method of:
+    
+    Schneider TD, Stephens RM. "Sequence logos: a new way to display consensus 
+    sequences." Nucleic Acids Res. 1990 Oct 25;18(20):6097-100.
+    """
+    # Columns are sequence positions, rows are symbol counts/frequencies
+    f = matrix.values.transpose()
+    n, m = f.shape
+    # Ensure normalized
+    f = f / sum( f, axis=0 )
+    # Shannon entropy (the where replaces 0 with 1 so that '0 log 0 == 0')
+    H = - sum( f * log2( where( f, f, 1 ) ), axis=0 )
+    # Height
+    return transpose( f * ( log2( n ) - H ) )
+    
+def logo( matrix, base_width, height, colors=DNA_DEFAULT_COLORS ):
+    """
+    Return an EPS document containing a sequence logo for matrix where each
+    bases is shown as a column of `base_width` points and the total logo
+    height is `height` points. If `colors` is provided it is a mapping from
+    characters to rgb color strings. 
+    """
+    rval = StringIO()
+    # Read header ans substitute in width / height
+    header = Template( pkg_resources.resource_string( __name__, "template.ps" ) )
+    rval.write( header.substitute( bounding_box_width = ceil( base_width * matrix.width ) + PAD,
+                                   bounding_box_height = ceil( height ) + PAD ) )
+    # Determine heights
+    heights = freqs_to_heights( matrix )
+    height_scale = height / log2( len( matrix.alphabet ) )
+    # Draw each "row" of the matrix
+    for i, row in enumerate( heights ):
+        x = ( i * base_width )
+        y = 0
+        for j, base_height in enumerate( row ):
+            char = matrix.alphabet[j]
+            page_height = height_scale * base_height
+            # print matrix.alphabet[j], base_height, height_scale, page_height
+            if page_height > 1:
+                # Draw letter
+                rval.write(  "%s setrgbcolor\n" % colors.get( char, '0 0 0' ) )
+                rval.write( "%3.2f " % x )
+                rval.write( "%3.2f " % y )
+                rval.write( "%3.2f " % ( x + base_width ) )
+                rval.write( "%3.2f " % ( y + page_height ) )
+                rval.write( "(%s) textInBox\n" % char )
+            y += page_height
+    rval.write( "showpage" )
+    return rval.getvalue()

File lib/bx/motif/logo/template.ps

+%!PS-Adobe-3.1 EPSF-3.0
+%%BoundingBox: 0 0 ${bounding_box_width} ${bounding_box_height}
+
+/logoFont {
+/Helvetica-Bold findfont
+10 scalefont
+setfont
+} def
+
+/textBounds {
+% Figure out bounding box of string in current font.  Usage:
+% call: text letterSize 
+% sets: tbX1 tbY1 tbX2 tbY2 tbW tbH
+% The bounds are relative to the current position
+gsave
+newpath
+0 0 moveto
+true charpath
+flattenpath
+pathbbox
+grestore
+/tbY2 exch def
+/tbX2 exch def
+/tbY1 exch def
+/tbX1 exch def
+/tbW tbX2 tbX1 sub def
+/tbH tbY2 tbY1 sub def
+} def
+
+/textInBox {
+% Draw text so that it fits inside of box.  Usage:
+%   x1 y1 x2 y2 text textInBox
+
+% Copy parameters from variables to stack and save context
+/tibText exch def
+/tibY2 exch def
+/tibX2 exch def
+/tibY1 exch def
+/tibX1 exch def
+gsave
+
+% move to x1/y1 adjusted for text offset
+tibText textBounds
+tibX1 tbX1 sub tibY1 tbY1 sub moveto
+
+% Set scaling
+/tibW tibX2 tibX1 sub def
+/tibH tibY2 tibY1 sub def
+tibW tbW div tibH tbH div scale
+
+% draw and return
+tibText show
+grestore
+} def
+
+/aColor { 0 0.7 0 setrgbcolor } def
+/cColor { 0 0.5 0.7 setrgbcolor } def
+/gColor { 0.8 0.5 0 setrgbcolor } def
+/tColor { 0.9 0 0 setrgbcolor } def
+
+logoFont
+