Mercurial > repos > davidmurphy > codonlogo

comparison weblogolib/__init__.py @ 4:4d47ab2b7bcc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidmurphy
date Fri, 13 Jan 2012 07:18:19 -0500
parents
children 5149eb3a89c2
comparison
equal deleted inserted replaced
3:09d2dac9ef73 4:4d47ab2b7bcc
1 #!/usr/bin/env python
2
3 # -------------------------------- WebLogo --------------------------------
4
5 # Copyright (c) 2003-2004 The Regents of the University of California.
6 # Copyright (c) 2005 Gavin E. Crooks
7 # Copyright (c) 2006, The Regents of the University of California, through print
8 # Lawrence Berkeley National Laboratory (subject to receipt of any required
9 # approvals from the U.S. Dept. of Energy). All rights reserved.
10
11 # This software is distributed under the new BSD Open Source License.
12 # <http://www.opensource.org/licenses/bsd-license.html>
13 #
14 # Redistribution and use in source and binary forms, with or without
15 # modification, are permitted provided that the following conditions are met:
16 #
17 # (1) Redistributions of source code must retain the above copyright notice,
18 # this list of conditions and the following disclaimer.
19 #
20 # (2) Redistributions in binary form must reproduce the above copyright
21 # notice, this list of conditions and the following disclaimer in the
22 # documentation and or other materials provided with the distribution.
23 #
24 # (3) Neither the name of the University of California, Lawrence Berkeley
25 # National Laboratory, U.S. Dept. of Energy nor the names of its contributors
26 # may be used to endorse or promote products derived from this software
27 # without specific prior written permission.
28 #
29 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
30 # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
31 # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
32 # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
33 # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
34 # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
35 # SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
36 # INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
37 # CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
38 # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
39 # POSSIBILITY OF SUCH DAMAGE.
40
41 # Replicates README.txt
42
43 """
44 WebLogo (http://code.google.com/p/weblogo/) is a tool for creating sequence
45 logos from biological sequence alignments. It can be run on the command line,
46 as a standalone webserver, as a CGI webapp, or as a python library.
47
48 The main WebLogo webserver is located at http://bespoke.lbl.gov/weblogo/
49
50
51 Codonlogo is based on Weblogo.
52
53 Please consult the manual for installation instructions and more information:
54 (Also located in the weblogolib/htdocs subdirectory.)
55
56 For help on the command line interface run
57 ./codonlogo --help
58
59 To build a simple logo run
60 ./codonlogo < cap.fa > logo0.eps
61
62 To run as a standalone webserver at localhost:8080
63 ./codonlogo --server
64
65 To create a logo in python code:
66 >>> from weblogolib import *
67 >>> fin = open('cap.fa')
68 >>> seqs = read_seq_data(fin)
69 >>> data = LogoData.from_seqs(seqs)
70 >>> options = LogoOptions()
71 >>> options.title = "A Logo Title"
72 >>> format = LogoFormat(data, options)
73 >>> fout = open('cap.eps', 'w')
74 >>> eps_formatter( data, format, fout)
75
76
77 -- Distribution and Modification --
78 This package is distributed under the new BSD Open Source License.
79 Please see the LICENSE.txt file for details on copyright and licensing.
80 The WebLogo source code can be downloaded from http://code.google.com/p/weblogo/
81 WebLogo requires Python 2.3, 2.4 or 2.5, the corebio python toolkit for computational
82 biology (http://code.google.com/p/corebio), and the python array package
83 'numpy' (http://www.scipy.org/Download)
84
85 # -------------------------------- CodonLogo --------------------------------
86
87
88 """
89
90 from math import *
91 import random
92 from itertools import izip, count
93 import sys
94 import copy
95 import os
96 from itertools import product
97 from datetime import datetime
98 from StringIO import StringIO
99
100
101
102 from corebio.data import rna_letters, dna_letters, amino_acid_letters
103 import random
104
105 # python2.3 compatability
106 from corebio._future import Template
107 from corebio._future.subprocess import *
108 from corebio._future import resource_string, resource_filename
109
110
111 from math import log, sqrt
112
113 # Avoid 'from numpy import *' since numpy has lots of names defined
114 from numpy import array, asarray, float64, ones, zeros, int32,all,any, shape
115 import numpy as na
116
117 from corebio.utils.deoptparse import DeOptionParser
118 from optparse import OptionGroup
119
120 from color import *
121 from colorscheme import *
122 from corebio.seq import Alphabet, Seq, SeqList
123 from corebio import seq_io
124 from corebio.utils import isfloat, find_command, ArgumentError
125 from corebio.moremath import *
126 from corebio.data import amino_acid_composition
127 from corebio.seq import unambiguous_rna_alphabet, unambiguous_dna_alphabet, unambiguous_protein_alphabet
128
129 codon_alphabetU=['AAA', 'AAC', 'AAG', 'AAU', 'ACA', 'ACC', 'ACG', 'ACU', 'AGA', 'AGC', 'AGG', 'AGU', 'AUA', 'AUC', 'AUG', 'AUU', 'CAA', 'CAC', 'CAG', 'CAU', 'CCA', 'CCC', 'CCG', 'CCU', 'CGA', 'CGC', 'CGG', 'CGU', 'CUA', 'CUC', 'CUG', 'CUU', 'GAA', 'GAC', 'GAG', 'GAU', 'GCA', 'GCC', 'GCG', 'GCU', 'GGA', 'GGC', 'GGG', 'GGU', 'GUA', 'GUC', 'GUG', 'GUU', 'UAA', 'UAC', 'UAG', 'UAU', 'UCA', 'UCC', 'UCG', 'UCU', 'UGA', 'UGC', 'UGG', 'UGU', 'UUA', 'UUC', 'UUG', 'UUU']
130 codon_alphabetT=['AAA', 'AAC', 'AAG', 'AAT', 'ACA', 'ACC', 'ACG', 'ACT', 'AGA', 'AGC', 'AGG', 'AGT', 'ATA', 'ATC', 'ATG', 'ATT', 'CAA', 'CAC', 'CAG', 'CAT', 'CCA', 'CCC', 'CCG', 'CCT', 'CGA', 'CGC', 'CGG', 'CGT', 'CTA', 'CTC', 'CTG', 'CTT', 'GAA', 'GAC', 'GAG', 'GAT', 'GCA', 'GCC', 'GCG', 'GCT', 'GGA', 'GGC', 'GGG', 'GGT', 'GTA', 'GTC', 'GTG', 'GTT', 'TAA', 'TAC', 'TAG', 'TAT', 'TCA', 'TCC', 'TCG', 'TCT', 'TGA', 'TGC', 'TGG', 'TGT', 'TTA', 'TTC', 'TTG', 'TTT']
131
132 altype="codonsT"
133 offset=0
134 isreversed=False
135 col=[]
136
137
138
139 __all__ = ['LogoSize',
140 'LogoOptions',
141 'description',
142 '__version__',
143 'LogoFormat',
144 'LogoData',
145 'Dirichlet',
146 'GhostscriptAPI',
147 'std_color_schemes',
148 'default_color_schemes',
149 'classic',
150 'std_units',
151 'std_sizes',
152 'std_alphabets',
153 'std_percentCG',
154 'pdf_formatter',
155 'jpeg_formatter',
156 'png_formatter',
157 'png_print_formatter',
158 'txt_formatter',
159 'eps_formatter',
160 'formatters',
161 'default_formatter',
162 'base_distribution',
163 'equiprobable_distribution',
164 'read_seq_data',
165 'which_alphabet',
166 'color',
167 'colorscheme',
168 ]
169
170 description = "Create sequence logos from biological sequence alignments."
171
172 __version__ = "1.0"
173
174 # These keywords are subsituted by subversion.
175 # The date and revision will only tell the truth after a branch or tag,
176 # since different files in trunk will have been changed at different times
177 release_date ="$Date: 2011-09-17 16:30:00 -0700 (Tue, 14 Oct 2008) $".split()[1]
178 release_build = "$Revision: 53 $".split()[1]
179 release_description = "CodonLogo %s (%s)" % (__version__, release_date)
180
181
182
183 def cgi(htdocs_directory) :
184 import weblogolib._cgi
185 weblogolib._cgi.main(htdocs_directory)
186
187 class GhostscriptAPI(object) :
188 """Interface to the command line program Ghostscript ('gs')"""
189
190 formats = ('png', 'pdf', 'jpeg')
191
192 def __init__(self, path=None) :
193 try:
194 command = find_command('gs', path=path)
195 except EnvironmentError:
196 try:
197 command = find_command('gswin32c.exe', path=path)
198 except EnvironmentError:
199 raise EnvironmentError("Could not find Ghostscript on path."
200 " There should be either a gs executable or a gswin32c.exe on your system's path")
201
202 self.command = command
203
204 def version(self) :
205 args = [self.command, '--version']
206 try :
207 p = Popen(args, stdout=PIPE)
208 (out,err) = p.communicate()
209 except OSError :
210 raise RuntimeError("Cannot communicate with ghostscript.")
211 return out.strip()
212
213 def convert(self, format, fin, fout, width, height, resolution=300) :
214 device_map = { 'png':'png16m', 'pdf':'pdfwrite', 'jpeg':'jpeg'}
215
216 try :
217 device = device_map[format]
218 except KeyError:
219 raise ValueError("Unsupported format.")
220
221 args = [self.command,
222 "-sDEVICE=%s" % device,
223 "-dPDFSETTINGS=/printer",
224 #"-q", # Quite: Do not dump messages to stdout.
225 "-sstdout=%stderr", # Redirect messages and errors to stderr
226 "-sOutputFile=-", # Stdout
227 "-dDEVICEWIDTHPOINTS=%s" % str(width),
228 "-dDEVICEHEIGHTPOINTS=%s" % str(height),
229 "-dSAFER", # For added security
230 "-dNOPAUSE",]
231
232 if device != 'pdf' :
233 args.append("-r%s" % str(resolution) )
234 if resolution < 300 : # Antialias if resolution is Less than 300 DPI
235 args.append("-dGraphicsAlphaBits=4")
236 args.append("-dTextAlphaBits=4")
237 args.append("-dAlignToPixels=0")
238
239 args.append("-") # Read from stdin. Must be last argument.
240
241 error_msg = "Unrecoverable error : Ghostscript conversion failed " \
242 "(Invalid postscript?). %s" % " ".join(args)
243
244 source = fin.read()
245
246 try :
247 p = Popen(args, stdin=PIPE, stdout = PIPE, stderr= PIPE)
248 (out,err) = p.communicate(source)
249 except OSError :
250 raise RuntimeError(error_msg)
251
252 if p.returncode != 0 :
253 error_msg += '\nReturn code: %i\n' % p.returncode
254 if err is not None : error_msg += err
255 raise RuntimeError(error_msg)
256
257 print >>fout, out
258 # end class Ghostscript
259
260
261 aa_composition = [ amino_acid_composition[_k] for _k in
262 unambiguous_protein_alphabet]
263
264
265
266 # ------ DATA ------
267
268 classic = ColorScheme([
269 ColorGroup("G", "orange" ),
270 ColorGroup("TU", "red"),
271 ColorGroup("C", "blue"),
272 ColorGroup("A", "green")
273 ] )
274
275 std_color_schemes = {"auto": None, # Depends on sequence type
276 "monochrome": monochrome,
277 "base pairing": base_pairing,
278 "classic": classic,
279 "hydrophobicity" : hydrophobicity,
280 "chemistry" : chemistry,
281 "charge" : charge,
282 }#
283
284 default_color_schemes = {
285 unambiguous_protein_alphabet: hydrophobicity,
286 unambiguous_rna_alphabet: base_pairing,
287 unambiguous_dna_alphabet: base_pairing
288 #codon_alphabet:codonsU
289 }
290
291
292 std_units = {
293 "bits" : 1./log(2),
294 "nats" : 1.,
295 "digits" : 1./log(10),
296 "kT" : 1.,
297 "kJ/mol" : 8.314472 *298.15 /1000.,
298 "kcal/mol": 1.987 *298.15 /1000.,
299 "probability" : None,
300 }
301
302 class LogoSize(object) :
303 def __init__(self, stack_width, stack_height) :
304 self.stack_width = stack_width
305 self.stack_height = stack_height
306
307 def __repr__(self):
308 return stdrepr(self)
309
310 # The base stack width is set equal to 9pt Courier.
311 # (Courier has a width equal to 3/5 of the point size.)
312 # Check that can get 80 characters in journal page @small
313 # 40 chacaters in a journal column
314 std_sizes = {
315 "small" : LogoSize( stack_width = 10, stack_height = 10*1*5),
316 "medium" : LogoSize( stack_width = 10*2, stack_height = 10*2*5),
317 "large" : LogoSize( stack_width = 10*3, stack_height = 10*3*5),
318 }
319
320
321 std_alphabets = {
322 'protein': unambiguous_protein_alphabet,
323 'rna': unambiguous_rna_alphabet,
324 'dna': unambiguous_dna_alphabet,
325 'codonsU':codon_alphabetU,
326 'codonsT':codon_alphabetT
327 }
328
329 std_percentCG = {
330 'H. sapiens' : 40.,
331 'E. coli' : 50.5,
332 'S. cerevisiae' : 38.,
333 'C. elegans' : 36.,
334 'D. melanogaster': 43.,
335 'M. musculus' : 42.,
336 'T. thermophilus' : 69.4,
337 }
338
339 # Thermus thermophilus: Henne A, Bruggemann H, Raasch C, Wiezer A, Hartsch T,
340 # Liesegang H, Johann A, Lienard T, Gohl O, Martinez-Arias R, Jacobi C,
341 # Starkuviene V, Schlenczeck S, Dencker S, Huber R, Klenk HP, Kramer W,
342 # Merkl R, Gottschalk G, Fritz HJ: The genome sequence of the extreme
343 # thermophile Thermus thermophilus.
344 # Nat Biotechnol 2004, 22:547-53
345
346 def stdrepr(obj) :
347 attr = vars(obj).items()
348
349
350 attr.sort()
351 args = []
352 for a in attr :
353 if a[0][0]=='_' : continue
354 args.append( '%s=%s' % ( a[0], repr(a[1])) )
355 args = ',\n'.join(args).replace('\n', '\n ')
356 return '%s(\n %s\n)' % (obj.__class__.__name__, args)
357
358
359 class LogoOptions(object) :
360 """ A container for all logo formating options. Not all of these
361 are directly accessible through the CLI or web interfaces.
362
363 To display LogoOption defaults:
364 >>> from weblogolib import *
365 >>> LogoOptions()
366
367
368 Attributes:
369 o alphabet
370 o creator_text -- Embedded as comment in figures.
371 o logo_title
372 o logo_label
373 o stacks_per_line
374 o unit_name
375 o show_yaxis
376 o yaxis_label -- Default depends on other settings.
377 o yaxis_tic_interval
378 o yaxis_minor_tic_ratio
379 o yaxis_scale
380 o show_xaxis
381 o xaxis_label
382 o xaxis_tic_interval
383 o rotate_numbers
384 o number_interval
385 o show_ends
386 o show_fineprint
387 o fineprint
388 o show_boxes
389 o shrink_fraction
390 o show_errorbars
391 o errorbar_fraction
392 o errorbar_width_fraction
393 o errorbar_gray
394 o resolution -- Dots per inch
395 o default_color
396 o color_scheme
397 o debug
398 o logo_margin
399 o stroke_width
400 o tic_length
401 o size
402 o stack_margin
403 o pad_right
404 o small_fontsize
405 o fontsize
406 o title_fontsize
407 o number_fontsize
408 o text_font
409 o logo_font
410 o title_font
411 o first_index
412 o logo_start
413 o logo_end
414 o scale_width
415 """
416
417 def __init__(self, **kwargs) :
418 """ Create a new LogoOptions instance.
419
420 >>> L = LogoOptions(logo_title = "Some Title String")
421 >>> L.show_yaxis = False
422 >>> repr(L)
423 """
424
425 self.creator_text = release_description,
426 self.alphabet = None
427
428 self.logo_title = ""
429 self.logo_label = ""
430 self.stacks_per_line = 20
431
432 self.unit_name = "bits"
433
434 self.show_yaxis = True
435 # yaxis_lable default depends on other settings. See LogoFormat
436 self.yaxis_label = None
437 self.yaxis_tic_interval = 1.
438 self.yaxis_minor_tic_ratio = 5
439 self.yaxis_scale = None
440
441 self.show_xaxis = True
442 self.xaxis_label = ""
443 self.xaxis_tic_interval =1
444 self.rotate_numbers = False
445 self.number_interval = 5
446 self.show_ends = False
447
448 self.show_fineprint = True
449 self.fineprint = "CodonLogo "+__version__
450
451 self.show_boxes = False
452 self.shrink_fraction = 0.5
453
454 self.show_errorbars = True
455 self.altype = True
456
457 self.errorbar_fraction = 0.90
458 self.errorbar_width_fraction = 0.25
459 self.errorbar_gray = 0.50
460
461 self.resolution = 96. # Dots per inch
462
463 self.default_color = Color.by_name("black")
464 self.color_scheme = None
465 #self.show_color_key = False # NOT yet implemented
466
467 self.debug = False
468
469 self.logo_margin = 2
470 self.stroke_width = 0.5
471 self.tic_length = 5
472
473 self.size = std_sizes["medium"]
474
475 self.stack_margin = 0.5
476 self.pad_right = False
477
478 self.small_fontsize = 6
479 self.fontsize = 10
480 self.title_fontsize = 12
481 self.number_fontsize = 8
482
483 self.text_font = "ArialMT"
484 self.logo_font = "Arial-BoldMT"
485 self.title_font = "ArialMT"
486
487 self.first_index = 1
488 self.logo_start = None
489 self.logo_end=None
490
491 # Scale width of characters proportional to gaps
492 self.scale_width = True
493
494 from corebio.utils import update
495 update(self, **kwargs)
496
497 def __repr__(self) :
498 attr = vars(self).items()
499 attr.sort()
500 args = []
501 for a in attr :
502 if a[0][0]=='_' : continue
503 args.append( '%s=%s' % ( a[0], repr(a[1])) )
504 args = ',\n'.join(args).replace('\n', '\n ')
505 return '%s(\n %s\n)' % (self.__class__.__name__, args)
506 # End class LogoOptions
507
508
509
510
511 class LogoFormat(LogoOptions) :
512 """ Specifies the format of the logo. Requires a LogoData and LogoOptions
513 objects.
514
515 >>> data = LogoData.from_seqs(seqs )
516 >>> options = LogoOptions()
517 >>> options.title = "A Logo Title"
518 >>> format = LogoFormat(data, options)
519 """
520 # TODO: Raise ArgumentErrors instead of ValueError and document
521 def __init__(self, data, options= None) :
522
523 LogoOptions.__init__(self)
524 #global offset
525 if options is not None :
526 self.__dict__.update(options.__dict__)
527
528 #offset=options.frame
529
530 self.alphabet = data.alphabet
531 self.seqlen = data.length
532 self.altype = True
533 self.show_title = False
534 self.show_xaxis_label = False
535 self.yaxis_minor_tic_interval = None
536 self.lines_per_logo = None
537 self.char_width = None
538 self.line_margin_left = None
539 self.line_margin_right = None
540 self.line_margin_bottom = None
541 self.line_margin_top = None
542 self.title_height = None
543 self.xaxis_label_height = None
544 self.line_height = None
545 self.line_width = None
546 self.logo_height = None
547 self.logo_width = None
548 self.creation_date = None
549 self.end_type = None
550
551 if self.stacks_per_line< 1 :
552 raise ArgumentError("Stacks per line should be greater than zero.",
553 "stacks_per_line" )
554
555 if self.size.stack_height<=0.0 :
556 raise ArgumentError(
557 "Stack height must be greater than zero.", "stack_height")
558 if (self.small_fontsize <= 0 or self.fontsize <=0 or
559 self.title_fontsize<=0 ):
560 raise ValueError("Font sizes must be positive.")
561
562 if self.errorbar_fraction<0.0 or self.errorbar_fraction>1.0 :
563 raise ValueError(
564 "The visible fraction of the error bar must be between zero and one.")
565
566 if self.yaxis_tic_interval<=0.0 :
567 raise ArgumentError( "The yaxis tic interval cannot be negative.",
568 'yaxis_tic_interval')
569
570 if self.size.stack_width <= 0.0 :
571 raise ValueError(
572 "The width of a stack should be a positive number.")
573
574 if self.yaxis_minor_tic_interval and \
575 self.yaxis_minor_tic_interval<=0.0 :
576 raise ValueError("Distances cannot be negative.")
577
578 if self.xaxis_tic_interval<=0 :
579 raise ValueError("Tic interval must be greater than zero.")
580
581 if self.number_interval<=0 :
582 raise ValueError("Invalid interval between numbers.")
583
584 if self.shrink_fraction<0.0 or self.shrink_fraction>1.0 :
585 raise ValueError("Invalid shrink fraction.")
586
587 if self.stack_margin<=0.0 :
588 raise ValueError("Invalid stack margin." )
589
590 if self.logo_margin<=0.0 :
591 raise ValueError("Invalid logo margin." )
592
593 if self.stroke_width<=0.0 :
594 raise ValueError("Invalid stroke width.")
595
596 if self.tic_length<=0.0 :
597 raise ValueError("Invalid tic length.")
598
599 # FIXME: More validation
600
601 # Inclusive upper and lower bounds
602 # FIXME: Validate here. Move from eps_formatter
603 if self.logo_start is None: self.logo_start = self.first_index
604
605 if self.logo_end is None :
606 self.logo_end = self.seqlen + self.first_index -1
607
608 self.total_stacks = self.logo_end - self.logo_start +1
609
610 if self.logo_start - self.first_index <0 :
611 raise ArgumentError(
612 "Logo range extends before start of available sequence.",
613 'logo_range')
614
615 if self.logo_end - self.first_index >= self.seqlen :
616 raise ArgumentError(
617 "Logo range extends beyond end of available sequence.",
618 'logo_range')
619
620 if self.logo_title : self.show_title = True
621 if not self.fineprint : self.show_fineprint = False
622 if self.xaxis_label : self.show_xaxis_label = True
623
624 if self.yaxis_label is None :
625 self.yaxis_label = self.unit_name
626
627 if self.yaxis_label :
628 self.show_yaxis_label = True
629 else :
630 self.show_yaxis_label = False
631 self.show_ends = False
632
633 if not self.yaxis_scale :
634 conversion_factor = std_units[self.unit_name]
635 if conversion_factor :
636 self.yaxis_scale=log(len(self.alphabet))*conversion_factor
637 #self.yaxis_scale=max(data.entropy)*conversion_factor
638 #marker# this is where I handle the max height. needs revision.
639 else :
640 self.yaxis_scale=1.0 # probability units
641
642 if self.yaxis_scale<=0.0 :
643 raise ValueError(('yaxis_scale', "Invalid yaxis scale"))
644 if self.yaxis_tic_interval >= self.yaxis_scale:
645 self.yaxis_tic_interval /= 2.
646
647 self.yaxis_minor_tic_interval \
648 = float(self.yaxis_tic_interval)/self.yaxis_minor_tic_ratio
649
650 if self.color_scheme is None :
651 #if self.alphabet in default_color_schemes :
652 #self.color_scheme = default_color_schemes[self.alphabet]
653 #else :
654 self.color_scheme = codonsT
655 #else:
656 #for color, symbols, desc in options.colors:
657 #try :
658 #self.color_scheme.append( ColorGroup(symbols, color, desc) )
659 #print >>sys.stderr, color_scheme.groups[2]
660 #except ValueError :
661 #raise ValueError(
662 #"error: option --color: invalid value: '%s'" % color )
663
664
665 self.lines_per_logo = 1+ ( (self.total_stacks-1) / self.stacks_per_line)
666
667 if self.lines_per_logo==1 and not self.pad_right:
668 self.stacks_per_line = min(self.stacks_per_line, self.total_stacks)
669
670 self.char_width = self.size.stack_width - 2* self.stack_margin
671
672
673 if self.show_yaxis :
674 self.line_margin_left = self.fontsize * 3.0
675 else :
676 self.line_margin_left = 0
677
678 if self.show_ends :
679 self.line_margin_right = self.fontsize *1.5
680 else :
681 self.line_margin_right = self.fontsize
682
683 if self.show_xaxis :
684 if self.rotate_numbers :
685 self.line_margin_bottom = self.number_fontsize *2.5
686 else:
687 self.line_margin_bottom = self.number_fontsize *1.5
688 else :
689 self.line_margin_bottom = 4
690
691 self.line_margin_top = 4
692
693 if self.show_title :
694 self.title_height = self.title_fontsize
695 else :
696 self.title_height = 0
697
698 self.xaxis_label_height =0.
699 if self.show_xaxis_label :
700 self.xaxis_label_height += self.fontsize
701 if self.show_fineprint :
702 self.xaxis_label_height += self.small_fontsize
703
704 self.line_height = (self.size.stack_height + self.line_margin_top +
705 self.line_margin_bottom )
706 self.line_width = (self.size.stack_width*self.stacks_per_line +
707 self.line_margin_left + self.line_margin_right )
708
709 self.logo_height = int(2*self.logo_margin + self.title_height \
710 + self.xaxis_label_height + self.line_height*self.lines_per_logo)
711 self.logo_width = int(2*self.logo_margin + self.line_width )
712
713
714 self.creation_date = datetime.now().isoformat(' ')
715
716 end_type = '-'
717 end_types = {
718 unambiguous_protein_alphabet: 'p',
719 unambiguous_rna_alphabet: '-',
720 unambiguous_dna_alphabet: 'd'
721 }
722 if self.show_ends and self.alphabet in end_types:
723 end_type = end_types[self.alphabet]
724 self.end_type = end_type
725 # End __init__
726 # End class LogoFormat
727
728
729
730 # ------ Logo Formaters ------
731 # Each formatter is a function f(LogoData, LogoFormat, output file).
732 # that draws a represntation of the logo into the given file.
733 # The main graphical formatter is eps_formatter. A mapping 'formatters'
734 # containing all available formatters is located after the formatter
735 # definitions.
736
737 def pdf_formatter(data, format, fout) :
738 """ Generate a logo in PDF format."""
739
740 feps = StringIO()
741 eps_formatter(data, format, feps)
742 feps.seek(0)
743
744 gs = GhostscriptAPI()
745 gs.convert('pdf', feps, fout, format.logo_width, format.logo_height)
746
747
748 def _bitmap_formatter(data, format, fout, device) :
749 feps = StringIO()
750 eps_formatter(data, format, feps)
751 feps.seek(0)
752
753 gs = GhostscriptAPI()
754 gs.convert(device, feps, fout,
755 format.logo_width, format.logo_height, format.resolution)
756
757
758 def jpeg_formatter(data, format, fout) :
759 """ Generate a logo in JPEG format."""
760 _bitmap_formatter(data, format, fout, device="jpeg")
761
762
763 def png_formatter(data, format, fout) :
764 """ Generate a logo in PNG format."""
765
766 _bitmap_formatter(data, format, fout, device="png")
767
768
769 def png_print_formatter(data, format, fout) :
770 """ Generate a logo in PNG format with print quality (600 DPI) resolution."""
771 format.resolution = 600
772 _bitmap_formatter(data, format, fout, device="png")
773
774
775 def txt_formatter( logodata, format, fout) :
776 """ Create a text representation of the logo data.
777 """
778 print >>fout, str(logodata)
779
780
781
782
783 def eps_formatter( logodata, format, fout) :
784 """ Generate a logo in Encapsulated Postscript (EPS)"""
785
786 subsitutions = {}
787 from_format =[
788 "creation_date", "logo_width", "logo_height",
789 "lines_per_logo", "line_width", "line_height",
790 "line_margin_right","line_margin_left", "line_margin_bottom",
791 "line_margin_top", "title_height", "xaxis_label_height",
792 "creator_text", "logo_title", "logo_margin",
793 "stroke_width", "tic_length",
794 "stacks_per_line", "stack_margin",
795 "yaxis_label", "yaxis_tic_interval", "yaxis_minor_tic_interval",
796 "xaxis_label", "xaxis_tic_interval", "number_interval",
797 "fineprint", "shrink_fraction", "errorbar_fraction",
798 "errorbar_width_fraction",
799 "errorbar_gray", "small_fontsize", "fontsize",
800 "title_fontsize", "number_fontsize", "text_font",
801 "logo_font", "title_font",
802 "logo_label", "yaxis_scale", "end_type",
803 "debug", "show_title", "show_xaxis",
804 "show_xaxis_label", "show_yaxis", "show_yaxis_label",
805 "show_boxes", "show_errorbars", "show_fineprint",
806 "rotate_numbers", "show_ends", "altype",
807
808 ]
809
810 for s in from_format :
811 subsitutions[s] = getattr(format,s)
812
813
814 from_format_size = ["stack_height", "stack_width"]
815 for s in from_format_size :
816 subsitutions[s] = getattr(format.size,s)
817
818 subsitutions["shrink"] = str(format.show_boxes).lower()
819
820
821 # --------- COLORS --------------
822 def format_color(color):
823 return " ".join( ("[",str(color.red) , str(color.green),
824 str(color.blue), "]"))
825
826 subsitutions["default_color"] = format_color(format.default_color)
827 global col
828 colors = []
829
830 if altype=="codonsT" or altype=="codonsU":
831 for group in col:
832 cf = format_color(group.color)
833 colors.append( " ("+group.symbols+") " + cf )
834 for group in format.color_scheme.groups :
835 cf = format_color(group.color)
836
837 colors.append( " ("+group.symbols+") " + cf )
838 #print >>sys.stderr,opts.colors
839 #print >>sys.stderr,logodata.options
840 #print >>sys.stderr, group.symbols
841 #print >>sys.stderr, cf
842
843
844
845 else:
846 for group in format.color_scheme.groups :
847 cf = format_color(group.color)
848 for s in group.symbols :
849 colors.append( " ("+s+") " + cf )
850
851 subsitutions["color_dict"] = "\n".join(colors)
852 data = []
853
854 # Unit conversion. 'None' for probability units
855 conv_factor = std_units[format.unit_name]
856
857 data.append("StartLine")
858
859
860 seq_from = format.logo_start- format.first_index
861 seq_to = format.logo_end - format.first_index +1
862
863 # seq_index : zero based index into sequence data
864 # logo_index : User visible coordinate, first_index based
865 # stack_index : zero based index of visible stacks
866 for seq_index in range(seq_from, seq_to) :
867 logo_index = seq_index + format.first_index
868 stack_index = seq_index - seq_from
869
870 if stack_index!=0 and (stack_index % format.stacks_per_line) ==0 :
871 data.append("")
872 data.append("EndLine")
873 data.append("StartLine")
874 data.append("")
875
876
877 if isreversed :
878 data.append("(%d) StartStack" % logo_index)
879 else :
880 data.append("(%d) StartStack" % logo_index)
881
882
883
884 if conv_factor:
885 stack_height = logodata.entropy[seq_index] * std_units[format.unit_name]
886 else :
887 stack_height = 1.0 # Probability
888
889 # if logodata.entropy_interval is not None and conv_factor:
890 # Draw Error bars
891 # low, high = logodata.entropy_interval[seq_index]
892 # center = logodata.entropy[seq_index]
893
894
895 # down = (center - low) * conv_factor
896 # up = (high - center) * conv_factor
897 # data.append(" %f %f %f DrawErrorbarFirst" % (down, up, stack_height) )
898
899 s = zip(logodata.counts[seq_index], logodata.alphabet)
900 def mycmp( c1, c2 ) :
901 # Sort by frequency. If equal frequency then reverse alphabetic
902 if c1[0] == c2[0] : return cmp(c2[1], c1[1])
903 return cmp(c1[0], c2[0])
904
905 s.sort(mycmp)
906
907 C = float(sum(logodata.counts[seq_index]))
908 if C > 0.0 :
909 fraction_width = 1.0
910 if format.scale_width :
911 fraction_width = logodata.weight[seq_index]
912 for c in s:
913 data.append(" %f %f (%s) ShowSymbol" % (fraction_width, c[0]*stack_height/C, c[1]) )
914
915 # Draw error bar on top of logo. Replaced by DrawErrorbarFirst above.
916 if logodata.entropy_interval is not None and conv_factor:
917 low, high = logodata.entropy_interval[seq_index]
918 center = logodata.entropy[seq_index]
919
920 down = (center - low) * conv_factor
921 up = (high - center) * conv_factor
922 data.append(" %f %f DrawErrorbar" % (down, up) )
923
924 data.append("EndStack")
925 data.append("")
926
927 data.append("EndLine")
928 subsitutions["logo_data"] = "\n".join(data)
929
930
931 # Create and output logo
932 template = resource_string( __name__, 'template.eps', __file__)
933 logo = Template(template).substitute(subsitutions)
934 print >>fout, logo
935
936
937 # map between output format names and logo
938 formatters = {
939 'eps': eps_formatter,
940 'pdf': pdf_formatter,
941 'png': png_formatter,
942 'png_print' : png_print_formatter,
943 'jpeg' : jpeg_formatter,
944 'txt' : txt_formatter,
945 }
946
947 default_formatter = eps_formatter
948
949
950
951
952
953 def parse_prior(composition, alphabet, weight=None) :
954 """ Parse a description of the expected monomer distribution of a sequence.
955
956 Valid compositions:
957
958 - None or 'none' : No composition sepecified
959 - 'auto' or 'automatic': Use the typical average distribution
960 for proteins and an equiprobable distribution for
961 everything else.
962 - 'equiprobable' : All monomers have the same probability.
963 - a percentage, e.g. '45%' or a fraction '0.45':
964 The fraction of CG bases for nucleotide alphabets
965 - a species name, e.g. 'E. coli', 'H. sapiens' :
966 Use the average CG percentage for the specie's
967 genome.
968 - An explicit distribution, e.g. {'A':10, 'C':40, 'G':40, 'T':10}
969 """
970
971
972 if composition is None: return None
973 comp = composition.strip()
974
975 if comp.lower() == 'none': return None
976
977 if weight is None and alphabet is not None: weight = float(len(alphabet))
978 if comp.lower() == 'equiprobable' :
979 prior = weight * equiprobable_distribution(len(alphabet))
980
981
982 elif comp.lower() == 'auto' or comp.lower() == 'automatic':
983 if alphabet == unambiguous_protein_alphabet :
984 prior = weight * asarray(aa_composition, float64)
985 else :
986 prior = weight * equiprobable_distribution(len(alphabet))
987 elif comp in std_percentCG :
988 prior = weight * base_distribution(std_percentCG[comp])
989
990 elif comp[-1] == '%' :
991 prior = weight * base_distribution( float(comp[:-1]))
992
993 elif isfloat(comp) :
994 prior = weight * base_distribution( float(comp)*100. )
995
996 elif composition[0] == '{' and composition[-1] == '}' :
997 explicit = composition[1: -1]
998 explicit = explicit.replace(',',' ').replace("'", ' ').replace('"',' ').replace(':', ' ').split()
999
1000 if len(explicit) != len(alphabet)*2 :
1001 #print explicit
1002 raise ValueError("Explicit prior does not match length of alphabet")
1003 prior = - ones(len(alphabet), float64)
1004 try :
1005 for r in range(len(explicit)/2) :
1006 letter = explicit[r*2]
1007 index = alphabet.index(letter)
1008 value = float(explicit[r*2 +1])
1009 prior[index] = value
1010 except ValueError :
1011 raise ValueError("Cannot parse explicit composition")
1012
1013 if any(prior==-1.) :
1014 raise ValueError("Explicit prior does not match alphabet")
1015 prior/= sum(prior)
1016 prior *= weight
1017
1018
1019 else :
1020 raise ValueError("Unknown or malformed composition: %s"%composition)
1021 if len(prior) != len(alphabet) :
1022 raise ValueError(
1023 "The sequence alphabet and composition are incompatible.")
1024
1025 return prior
1026
1027
1028 def base_distribution(percentCG) :
1029 A = (1. - (percentCG/100.))/2.
1030 C = (percentCG/100.)/2.
1031 G = (percentCG/100.)/2.
1032 T = (1. - (percentCG/100))/2.
1033 return asarray((A,C,G,T), float64)
1034
1035 def equiprobable_distribution( length) :
1036 return ones( (length), float64) /length
1037
1038
1039
1040
1041 def read_seq_data(fin, input_parser=seq_io.read,alphabet=None, ignore_lower_case=False, max_file_size=0):
1042 # TODO: Document this method and enviroment variable
1043 max_file_size =int(os.environ.get("WEBLOGO_MAX_FILE_SIZE", max_file_size))
1044
1045 # If max_file_size is set, or if fin==stdin (which is non-seekable), we
1046 # read the data and replace fin with a StringIO object.
1047
1048 if(max_file_size>0) :
1049 data = fin.read(max_file_size)
1050
1051 more_data = fin.read(2)
1052 if more_data != "" :
1053 raise IOError("File exceeds maximum allowed size: %d bytes" % max_file_size)
1054
1055 fin = StringIO(data)
1056 elif fin == sys.stdin:
1057 fin = StringIO(fin.read())
1058
1059 seqs = input_parser(fin)
1060
1061 if seqs is None or len(seqs) ==0 :
1062 raise ValueError("Please provide a multiple sequence alignment")
1063
1064 if ignore_lower_case :
1065 # Case is significant. Do not count lower case letters.
1066 for i,s in enumerate(seqs) :
1067 seqs[i] = s.mask()
1068
1069 global altype
1070 if(altype=="codonsT" or altype=="codonsU"):
1071 if 'T' in seqs[0] or 't' in seqs[0]:
1072 altype="codonsT"
1073 if 'U' in seqs[0] or 'u' in seqs[0]:
1074 altype="codonsU"
1075 global offset
1076 global isreversed
1077 seq2=[""]*len(seqs)
1078 seq2 = []
1079 for i in xrange(len(seqs)):
1080 seq2.append([])
1081 if(offset%6>2):
1082 for x in range(0,len(seqs)):
1083 backs=seqs[x][::-1]
1084 for y in range(0,len(backs)):
1085 seq2[x].append(str(backs[y]))
1086
1087 if(altype=="codonsU"):
1088 for x in range(0,len(seq2)):
1089 for y in range(0,len(seq2[x])):
1090 if(cmp(seq2[x][y],'G')==0):
1091 seq2[x][y]="C"
1092 elif(cmp(seq2[x][y],'A')==0):
1093 seq2[x][y]='U'
1094 elif(cmp(seq2[x][y],'U')==0):
1095 seq2[x][y]='A'
1096 elif(cmp(seq2[x][y],'C')==0):
1097 seq2[x][y]='G'
1098 if(altype=="codonsT"):
1099 for x in range(0,len(seq2)):
1100 for y in range(0,len(seq2[x])):
1101 if(cmp(seq2[x][y],'G')==0):
1102 seq2[x][y]='C'
1103 elif(cmp(seq2[x][y],'A')==0):
1104 seq2[x][y]='T'
1105 elif(cmp(seq2[x][y],'T')==0):
1106 seq2[x][y]='A'
1107 elif(cmp(seq2[x][y],'C')==0):
1108 seq2[x][y]='G'
1109 offset=offset%3
1110 isreversed=True
1111 for x in range(0,len(seqs)):
1112 seqs[x]=Seq("".join(seq2[x]))
1113
1114
1115 # Add alphabet to seqs.
1116 if alphabet :
1117 seqs.alphabet = alphabet
1118 else :
1119 seqs.alphabet = which_alphabet(seqs)
1120
1121 return seqs
1122
1123
1124
1125 #TODO: Move to seq_io?
1126 # Would have to check that no symbol outside of full alphabet?
1127 def which_alphabet(seqs) :
1128 """ Returns the most appropriate unambiguous protien, rna or dna alphabet
1129 for a Seq or SeqList.
1130 """
1131 alphabets = (unambiguous_protein_alphabet,
1132 unambiguous_rna_alphabet,
1133 unambiguous_dna_alphabet
1134 )
1135 # Heuristic
1136 # Count occurances of each letter. Downweight longer alphabet.
1137 #for x in seqs:
1138
1139
1140 if( altype=="codonsU"):
1141 return codon_alphabetU
1142 if( altype=="codonsT"):
1143 return codon_alphabetT
1144 else:
1145 score = [1.0*asarray(seqs.tally(a)).sum()/sqrt(len(a)) for a in alphabets]
1146 #print score
1147 best = argmax(score) # Ties go to last in list.
1148 a = alphabets[best]
1149 return a
1150
1151
1152
1153 class LogoData(object) :
1154 """The data needed to generate a sequence logo.
1155
1156 - alphabet
1157 - length
1158 - counts -- An array of character counts
1159 - entropy -- The relative entropy of each column
1160 - entropy_interval -- entropy confidence interval
1161 """
1162
1163 def __init__(self, length=None, alphabet = None, counts =None,
1164 entropy =None, entropy_interval = None, weight=None) :
1165 """Creates a new LogoData object"""
1166 self.length = length
1167 self.alphabet = alphabet
1168 self.counts = counts
1169 self.entropy = entropy
1170 self.entropy_interval = entropy_interval
1171 self.weight = weight
1172
1173
1174
1175 #@classmethod
1176 def from_counts(cls, alphabet, counts, prior= None):
1177
1178 """Build a logodata object from counts."""
1179 seq_length, A = counts.shape
1180
1181 if prior is not None: prior = array(prior, float64)
1182 if prior is None :
1183 R = log(A)
1184 ent = zeros( seq_length, float64)
1185 entropy_interval = None
1186 for i in range (0, seq_length) :
1187 C = sum(counts[i])
1188 #FIXME: fixup corebio.moremath.entropy()?
1189 if C == 0 :
1190 ent[i] = 0.0
1191 else :
1192 ent[i] = R - entropy(counts[i])
1193 else :
1194 ent = zeros( seq_length, float64)
1195 entropy_interval = zeros( (seq_length,2) , float64)
1196
1197 R = log(A)
1198
1199 for i in range (0, seq_length) :
1200 alpha = array(counts[i] , float64)
1201 alpha += prior
1202
1203 posterior = Dirichlet(alpha)
1204 ent[i] = posterior.mean_relative_entropy(prior/sum(prior))
1205 entropy_interval[i][0], entropy_interval[i][1] = \
1206 posterior.interval_relative_entropy(prior/sum(prior), 0.95)
1207 weight = array( na.sum(counts,axis=1) , float)
1208
1209 weight /= max(weight)
1210
1211 return cls(seq_length, alphabet, counts, ent, entropy_interval, weight)
1212 from_counts = classmethod(from_counts)
1213
1214
1215 #@classmethod
1216 def from_seqs(cls, seqs, prior= None):
1217
1218
1219 alphabet=seqs.alphabet
1220
1221 #get the offset and if it's greater than 2 swap bases to get the reverse compliment and reverse.
1222 for x in range(0,len(seqs)):
1223 seqs[x]=seqs[x].upper()
1224 counter=0
1225
1226
1227 """Build a 2D array from a SeqList, a list of sequences."""
1228 # --- VALIDATE DATA ---
1229 # check that there is at least one sequence of length 1
1230 if len(seqs)==0 or len(seqs[0]) ==0:
1231 raise ValueError("No sequence data found.")
1232 sys.exit(0)
1233 # Check sequence lengths
1234 seq_length = len(seqs[0])
1235 for i,s in enumerate(seqs) :
1236 #print i,s, len(s)
1237 if seq_length != len(s) :
1238 raise ArgumentError(
1239 "Sequence number %d differs in length from the previous sequences" % (i+1) ,'sequences')
1240 sys.exit(0)
1241
1242 if(altype=="codonsT" or altype=="codonsU"):
1243 x = [[0]*64 for x in xrange(seq_length/3)]
1244 counter=offset
1245
1246 while counter+offset<seq_length:
1247 for i in range(0,len(seqs)):
1248 if len(str(seqs[i][(counter):(counter+3)]))==3 and len(seqs[i][(counter):(counter+3)].strip("GATUC"))==0 :
1249 if(str(seqs[i][(counter):(counter+3)]) in alphabet):
1250 x[counter/3][ (alphabet.index(str(seqs[i][(counter):(counter+3)]))) ]+=1
1251 counter=counter+3
1252 counts=asarray(x)
1253 else:
1254 counts = asarray(seqs.tally())
1255
1256 return cls.from_counts(alphabet, counts, prior)
1257 from_seqs = classmethod(from_seqs)
1258
1259
1260
1261 def __str__(self) :
1262 out = StringIO()
1263 print >>out, '## LogoData'
1264 print >>out, '# First column is position number, couting from zero'
1265 print >>out, '# Subsequent columns are raw symbol counts'
1266 print >>out, '# Entropy is mean entropy measured in nats.'
1267 print >>out, '# Low and High are the 95% confidence limits.'
1268 print >>out, '# Weight is the fraction of non-gap symbols in the column.'
1269 print >>out, '#\t'
1270 print >>out, '#\t',
1271 for a in self.alphabet :
1272 print >>out, a, '\t',
1273 print >>out, 'Entropy\tLow\tHigh\tWeight'
1274
1275 for i in range(self.length) :
1276 print >>out, i, '\t',
1277 for c in self.counts[i] : print >>out, c, '\t',
1278 print >>out, self.entropy[i], '\t',
1279 if self.entropy_interval is not None:
1280 print >>out, self.entropy_interval[i][0], '\t',
1281 print >>out, self.entropy_interval[i][1], '\t',
1282 else :
1283 print >>out, '\t','\t',
1284 if self.weight is not None :
1285 print >>out, self.weight[i],
1286 print >>out, ''
1287 print >>out, '# End LogoData'
1288
1289 return out.getvalue()
1290
1291 # ====================== Main: Parse Command line =============================
1292 def main():
1293 """CodonLogo command line interface """
1294
1295 # ------ Parse Command line ------
1296 parser = _build_option_parser()
1297 (opts, args) = parser.parse_args(sys.argv[1:])
1298 if args : parser.error("Unparsable arguments: %s " % args)
1299
1300 if opts.serve:
1301 httpd_serve_forever(opts.port) # Never returns?
1302 sys.exit(0)
1303
1304
1305 # ------ Create Logo ------
1306 try:
1307 data = _build_logodata(opts)
1308
1309
1310 format = _build_logoformat(data, opts)
1311
1312
1313 formatter = opts.formatter
1314 formatter(data, format, opts.fout)
1315
1316 except ValueError, err :
1317 print >>sys.stderr, 'Error:', err
1318 sys.exit(2)
1319 except KeyboardInterrupt, err:
1320 sys.exit(0)
1321 # End main()
1322
1323
1324 def httpd_serve_forever(port=8080) :
1325 """ Start a webserver on a local port."""
1326 import BaseHTTPServer
1327 import CGIHTTPServer
1328
1329 class __HTTPRequestHandler(CGIHTTPServer.CGIHTTPRequestHandler):
1330 def is_cgi(self) :
1331 if self.path == "/create.cgi":
1332 self.cgi_info = '', 'create.cgi'
1333 return True
1334 return False
1335
1336 # Add current directory to PYTHONPATH. This is
1337 # so that we can run the standalone server
1338 # without having to run the install script.
1339 pythonpath = os.getenv("PYTHONPATH", '')
1340 pythonpath += ":" + os.path.abspath(sys.path[0]).split()[0]
1341 os.environ["PYTHONPATH"] = pythonpath
1342
1343 htdocs = resource_filename(__name__, 'htdocs', __file__)
1344 os.chdir(htdocs)
1345
1346 HandlerClass = __HTTPRequestHandler
1347 ServerClass = BaseHTTPServer.HTTPServer
1348 httpd = ServerClass(('', port), HandlerClass)
1349 print "Serving HTTP on localhost:%d ..." % port
1350
1351 try :
1352 httpd.serve_forever()
1353 except KeyboardInterrupt:
1354 sys.exit(0)
1355 # end httpd_serve_forever()
1356
1357 def read_priors(finp, alphabet ,max_file_size=0):
1358
1359 max_file_size =int(os.environ.get("WEBLOGO_MAX_FILE_SIZE", max_file_size))
1360 if(max_file_size>0) :
1361 data = finp.read(max_file_size)
1362 more_data = finp.read(2)
1363 if more_data != "" :
1364 raise IOError("File exceeds maximum allowed size: %d bytes" % max_file_size)
1365 finp = StringIO(data)
1366 priordict={}
1367 while 1:
1368 line = finp.readline()
1369 if not line:
1370 break
1371 line = line.split()
1372 priordict[line[0]]=line[1]
1373 return priordict
1374
1375 def _build_logodata(options) :
1376 global offset
1377 offset=options.frame
1378 options.alphabet = None
1379 options.ignore_lower_case = False
1380 #options.default_color = Color.by_name("black")
1381 options.color_scheme=None
1382 #options.colors=[]
1383 options.show_ends=False
1384 seqs = read_seq_data(options.fin,
1385 options.input_parser.read,
1386 alphabet=options.alphabet,
1387 ignore_lower_case = options.ignore_lower_case)
1388 if(options.priorfile!=None):
1389 if(altype=="CodonsT"):
1390 options.composition= str(read_priors(options.priorfile,codon_alphabetT))
1391 options.alphabet = codon_alphabetT
1392 else:
1393 options.composition= str(read_priors(options.priorfile,codon_alphabetU))
1394 options.alphabet = codon_alphabetU
1395
1396 prior = parse_prior( options.composition,seqs.alphabet, options.weight)
1397 data = LogoData.from_seqs(seqs, prior)
1398 return data
1399
1400
1401 def _build_logoformat( logodata, opts) :
1402 """ Extract and process relevant option values and return a
1403 LogoFormat object."""
1404
1405 args = {}
1406 direct_from_opts = [
1407 "stacks_per_line",
1408 "logo_title",
1409 "yaxis_label",
1410 "show_xaxis",
1411 "show_yaxis",
1412 "xaxis_label",
1413 "show_ends",
1414 "fineprint",
1415 "show_errorbars",
1416 "show_boxes",
1417 "yaxis_tic_interval",
1418 "resolution",
1419 "alphabet",
1420 "debug",
1421 "show_ends",
1422 "default_color",
1423 #"show_color_key",
1424 "color_scheme",
1425 "unit_name",
1426 "logo_label",
1427 "yaxis_scale",
1428 "first_index",
1429 "logo_start",
1430 "logo_end",
1431 "scale_width",
1432 "frame",
1433 ]
1434
1435 for k in direct_from_opts:
1436 args[k] = opts.__dict__[k]
1437 logo_size = copy.copy(opts.__dict__['logo_size'])
1438 size_from_opts = ["stack_width", "stack_height"]
1439 for k in size_from_opts :
1440 length = getattr(opts, k)
1441 if length : setattr( logo_size, k, length )
1442 args["size"] = logo_size
1443
1444 global col
1445
1446 if opts.colors:
1447 color_scheme = ColorScheme()
1448 for color, symbols, desc in opts.colors:
1449 try :
1450 #c = Color.from_string(color)
1451 color_scheme.groups.append( ColorGroup(symbols, color, desc) )
1452 #print >> sys.stderr, color
1453
1454 col.append( ColorGroup(symbols, color, desc) )
1455
1456 except ValueError :
1457 raise ValueError(
1458 "error: option --color: invalid value: '%s'" % color )
1459 if(altype!="codonsU" and altype!="codonsT") :
1460 args["color_scheme"] = color_scheme
1461
1462 #cf = colorscheme.format_color(col[0])
1463 #col.append( " ("+group.symbols+") " + cf )
1464
1465 logooptions = LogoOptions()
1466 for a, v in args.iteritems() :
1467 setattr(logooptions,a,v)
1468
1469
1470
1471 theformat = LogoFormat(logodata, logooptions )
1472
1473 return theformat
1474
1475
1476
1477
1478
1479
1480 # ========================== OPTIONS ==========================
1481 def _build_option_parser() :
1482 defaults = LogoOptions()
1483 parser = DeOptionParser(usage="%prog [options] < sequence_data.fa > sequence_logo.eps",
1484 description = description,
1485 version = __version__ ,
1486 add_verbose_options = False
1487 )
1488
1489 io_grp = OptionGroup(parser, "Input/Output Options",)
1490 data_grp = OptionGroup(parser, "Logo Data Options",)
1491 format_grp = OptionGroup(parser, "Logo Format Options",
1492 "These options control the format and display of the logo.")
1493 color_grp = OptionGroup(parser, "Color Options",
1494 "Colors can be specified using CSS2 syntax. e.g. 'red', '#FF0000', etc.")
1495 advanced_grp = OptionGroup(parser, "Advanced Format Options",
1496 "These options provide fine control over the display of the logo. ")
1497 server_grp = OptionGroup(parser, "CodonLogo Server",
1498 "Run a standalone webserver on a local port.")
1499
1500
1501 parser.add_option_group(io_grp)
1502 parser.add_option_group(data_grp)
1503 parser.add_option_group(format_grp)
1504 parser.add_option_group(color_grp)
1505 parser.add_option_group(advanced_grp)
1506 parser.add_option_group(server_grp)
1507
1508 # ========================== IO OPTIONS ==========================
1509
1510
1511
1512 io_grp.add_option( "-f", "--fin",
1513 dest="fin",
1514 action="store",
1515 type="file_in",
1516 default=sys.stdin,
1517 help="Sequence input file (default: stdin)",
1518 metavar="FILENAME")
1519
1520 io_grp.add_option( "-R", "--prior",
1521 dest="priorfile",
1522 action="store",
1523 type="file_in",
1524 help="A file with 64 codons and their prior probabilities, one per line, each codon followed by a space and it's probability.",
1525 metavar="FILENAME")
1526
1527 io_grp.add_option("", "--fin-format",
1528 dest="input_parser",
1529 action="store", type ="dict",
1530 default = seq_io,
1531 choices = seq_io.format_names(),
1532 help="Multiple sequence alignment format: (%s)" %
1533 ', '.join([ f.names[0] for f in seq_io.formats]),
1534 metavar="FORMAT")
1535
1536 io_grp.add_option("-o", "--fout", dest="fout",
1537 type="file_out",
1538 default=sys.stdout,
1539 help="Output file (default: stdout)",
1540 metavar="FILENAME")
1541
1542
1543
1544 io_grp.add_option( "-F", "--format",
1545 dest="formatter",
1546 action="store",
1547 type="dict",
1548 choices = formatters,
1549 metavar= "FORMAT",
1550 help="Format of output: eps (default), png, png_print, pdf, jpeg, txt",
1551 default = default_formatter)
1552
1553
1554 # ========================== Data OPTIONS ==========================
1555
1556 data_grp.add_option("-m", "--frame",
1557 dest="frame",
1558 action="store",
1559 type="int",
1560 default=0,
1561 help="Offset of the reading frame you wish to look in (default: 0) 3-6 indicate the reverse complement in which case codonlogo will read from the last symbol in the sequences backwards and replace each base with it's complement.",
1562 metavar="NUMBER")
1563
1564 #data_grp.add_option("-T", "--type",
1565 #dest="altype",
1566 #action="store",
1567 #type="boolean",
1568 #default=True,
1569 #help="Generate a codon logo rather than a sequence logo (default: True)",
1570 #metavar="YES/NO")
1571
1572
1573
1574 #data_grp.add_option( "-A", "--sequence-type",
1575 #dest="alphabet",
1576 #action="store",
1577 #type="dict",
1578 #choices = std_alphabets,
1579 #help="The type of sequence data: 'protein', 'rna' or 'dna'.",
1580 #metavar="TYPE")
1581
1582 #data_grp.add_option( "-a", "--alphabet",
1583 #dest="alphabet",
1584 #action="store",
1585 #help="The set of symbols to count, e.g. 'AGTC'. "
1586 #"All characters not in the alphabet are ignored. "
1587 #"If neither the alphabet nor sequence-type are specified then codonlogo will examine the input data and make an educated guess. "
1588 #"See also --sequence-type, --ignore-lower-case" )
1589
1590 # FIXME Add test?
1591 #data_grp.add_option( "", "--ignore-lower-case",
1592 #dest="ignore_lower_case",
1593 #action="store",
1594 #type = "boolean",
1595 #default=False,
1596 #metavar = "YES/NO",
1597 #help="Disregard lower case letters and only count upper case letters in sequences? (Default: No)"
1598 #)
1599
1600 data_grp.add_option( "-U", "--units",
1601 dest="unit_name",
1602 action="store",
1603 choices = std_units.keys(),
1604 type="choice",
1605 default = defaults.unit_name,
1606 help="A unit of entropy ('bits' (default), 'nats', 'digits'), or a unit of free energy ('kT', 'kJ/mol', 'kcal/mol'), or 'probability' for probabilities",
1607 metavar = "NUMBER")
1608
1609
1610 data_grp.add_option( "", "--composition",
1611 dest="composition",
1612 action="store",
1613 type="string",
1614 default = "auto",
1615 help="The expected composition of the sequences: 'auto' (default), 'equiprobable', 'none' (Do not perform any compositional adjustment), ",
1616 metavar="COMP.")
1617
1618 data_grp.add_option( "", "--weight",
1619 dest="weight",
1620 action="store",
1621 type="float",
1622 default = None,
1623 help="The weight of prior data. Default: total pseudocounts equal to the number of monomer types.",
1624 metavar="NUMBER")
1625
1626 data_grp.add_option( "-i", "--first-index",
1627 dest="first_index",
1628 action="store",
1629 type="int",
1630 default = 1,
1631 help="Index of first position in sequence data (default: 1)",
1632 metavar="INDEX")
1633
1634 data_grp.add_option( "-l", "--lower",
1635 dest="logo_start",
1636 action="store",
1637 type="int",
1638 help="Lower bound of sequence to display",
1639 metavar="INDEX")
1640
1641 data_grp.add_option( "-u", "--upper",
1642 dest="logo_end",
1643 action="store",
1644 type="int",
1645 help="Upper bound of sequence to display",
1646 metavar="INDEX")
1647
1648 # ========================== FORMAT OPTIONS ==========================
1649
1650 format_grp.add_option( "-s", "--size",
1651 dest="logo_size",
1652 action="store",
1653 type ="dict",
1654 choices = std_sizes,
1655 metavar = "LOGOSIZE",
1656 default = defaults.size,
1657 help="Specify a standard logo size (small, medium (default), large)" )
1658
1659
1660
1661 format_grp.add_option( "-n", "--stacks-per-line",
1662 dest="stacks_per_line",
1663 action="store",
1664 type="int",
1665 help="Maximum number of logo stacks per logo line. (default: %default)",
1666 default = defaults.stacks_per_line,
1667 metavar="COUNT")
1668
1669 format_grp.add_option( "-t", "--title",
1670 dest="logo_title",
1671 action="store",
1672 type="string",
1673 help="Logo title text.",
1674 default = defaults.logo_title,
1675 metavar="TEXT")
1676
1677 format_grp.add_option( "", "--label",
1678 dest="logo_label",
1679 action="store",
1680 type="string",
1681 help="A figure label, e.g. '2a'",
1682 default = defaults.logo_label,
1683 metavar="TEXT")
1684
1685 format_grp.add_option( "-X", "--show-xaxis",
1686 action="store",
1687 type = "boolean",
1688 default= defaults.show_xaxis,
1689 metavar = "YES/NO",
1690 help="Display sequence numbers along x-axis? (default: %default)")
1691
1692 format_grp.add_option( "-x", "--xlabel",
1693 dest="xaxis_label",
1694 action="store",
1695 type="string",
1696 default = defaults.xaxis_label,
1697 help="X-axis label",
1698 metavar="TEXT")
1699
1700 format_grp.add_option( "-S", "--yaxis",
1701 dest="yaxis_scale",
1702 action="store",
1703 type="float",
1704 help="Height of yaxis in units. (Default: Maximum value with uninformative prior.)",
1705 metavar = "UNIT")
1706
1707 format_grp.add_option( "-Y", "--show-yaxis",
1708 action="store",
1709 type = "boolean",
1710 dest = "show_yaxis",
1711 default= defaults.show_yaxis,
1712 metavar = "YES/NO",
1713 help="Display entropy scale along y-axis? (default: %default)")
1714
1715 format_grp.add_option( "-y", "--ylabel",
1716 dest="yaxis_label",
1717 action="store",
1718 type="string",
1719 help="Y-axis label (default depends on plot type and units)",
1720 metavar="TEXT")
1721
1722 #format_grp.add_option( "-E", "--show-ends",
1723 #action="store",
1724 #type = "boolean",
1725 #default= defaults.show_ends,
1726 #metavar = "YES/NO",
1727 #help="Label the ends of the sequence? (default: %default)")
1728
1729 format_grp.add_option( "-P", "--fineprint",
1730 dest="fineprint",
1731 action="store",
1732 type="string",
1733 default= defaults.fineprint,
1734 help="The fine print (default: Codonlogo version)",
1735 metavar="TEXT")
1736
1737 format_grp.add_option( "", "--ticmarks",
1738 dest="yaxis_tic_interval",
1739 action="store",
1740 type="float",
1741 default= defaults.yaxis_tic_interval,
1742 help="Distance between ticmarks (default: %default)",
1743 metavar = "NUMBER")
1744
1745
1746 format_grp.add_option( "", "--errorbars",
1747 dest = "show_errorbars",
1748 action="store",
1749 type = "boolean",
1750 default= defaults.show_errorbars,
1751 metavar = "YES/NO",
1752 help="Display error bars? (default: %default)")
1753
1754
1755
1756 # ========================== Color OPTIONS ==========================
1757 # TODO: Future Feature
1758 # color_grp.add_option( "-K", "--color-key",
1759 # dest= "show_color_key",
1760 # action="store",
1761 # type = "boolean",
1762 # default= defaults.show_color_key,
1763 # metavar = "YES/NO",
1764 # help="Display a color key (default: %default)")
1765
1766
1767 #color_scheme_choices = std_color_schemes.keys()
1768 #color_scheme_choices.sort()
1769 #color_grp.add_option( "-c", "--color-scheme",
1770 #dest="color_scheme",
1771 #action="store",
1772 #type ="dict",
1773 #choices = std_color_schemes,
1774 #metavar = "SCHEME",
1775 #default = None, # Auto
1776 #help="Specify a standard color scheme (%s)" % \
1777 #", ".join(color_scheme_choices) )
1778
1779 color_grp.add_option( "-C", "--color",
1780 dest="colors",
1781 action="append",
1782 metavar="COLOR SYMBOLS DESCRIPTION ",
1783 nargs = 3,
1784 default=[],
1785 help="Specify symbol colors, e.g. --color black AG 'Purine' --color red TC 'Pyrimidine' ")
1786
1787 color_grp.add_option( "", "--default-color",
1788 dest="default_color",
1789 action="store",
1790 metavar="COLOR",
1791 default= defaults.default_color,
1792 help="Symbol color if not otherwise specified.")
1793
1794 # ========================== Advanced options =========================
1795
1796 advanced_grp.add_option( "-W", "--stack-width",
1797 dest="stack_width",
1798 action="store",
1799 type="float",
1800 default= None,
1801 help="Width of a logo stack (default: %s)"% defaults.size.stack_width,
1802 metavar="POINTS" )
1803
1804 advanced_grp.add_option( "-H", "--stack-height",
1805 dest="stack_height",
1806 action="store",
1807 type="float",
1808 default= None,
1809 help="Height of a logo stack (default: %s)"%defaults.size.stack_height,
1810 metavar="POINTS" )
1811
1812 advanced_grp.add_option( "", "--box",
1813 dest="show_boxes",
1814 action="store",
1815 type = "boolean",
1816 default=False,
1817 metavar = "YES/NO",
1818 help="Draw boxes around symbols? (default: no)")
1819
1820 advanced_grp.add_option( "", "--resolution",
1821 dest="resolution",
1822 action="store",
1823 type="float",
1824 default=96,
1825 help="Bitmap resolution in dots per inch (DPI). (default: 96 DPI, except png_print, 600 DPI) Low resolution bitmaps (DPI<300) are antialiased.",
1826 metavar="DPI")
1827
1828 advanced_grp.add_option( "", "--scale-width",
1829 dest="scale_width",
1830 action="store",
1831 type = "boolean",
1832 default= True,
1833 metavar = "YES/NO",
1834 help="Scale the visible stack width by the fraction of symbols in the column? (i.e. columns with many gaps of unknowns are narrow.) (default: yes)")
1835
1836 advanced_grp.add_option( "", "--debug",
1837 action="store",
1838 type = "boolean",
1839 default= defaults.debug,
1840 metavar = "YES/NO",
1841 help="Output additional diagnostic information. (default: %default)")
1842
1843
1844 # ========================== Server options =========================
1845 server_grp.add_option( "", "--serve",
1846 dest="serve",
1847 action="store_true",
1848 default= False,
1849 help="Start a standalone CodonLogo server for creating sequence logos.")
1850
1851 server_grp.add_option( "", "--port",
1852 dest="port",
1853 action="store",
1854 type="int",
1855 default= 8080,
1856 help="Listen to this local port. (Default: %default)",
1857 metavar="PORT")
1858
1859 return parser
1860
1861 # END _build_option_parser
1862
1863
1864 ##############################################################
1865
1866 class Dirichlet(object) :
1867 """The Dirichlet probability distribution. The Dirichlet is a continuous
1868 multivariate probability distribution across non-negative unit length
1869 vectors. In other words, the Dirichlet is a probability distribution of
1870 probability distributions. It is conjugate to the multinomial
1871 distribution and is widely used in Bayesian statistics.
1872
1873 The Dirichlet probability distribution of order K-1 is
1874
1875 p(theta_1,...,theta_K) d theta_1 ... d theta_K =
1876 (1/Z) prod_i=1,K theta_i^{alpha_i - 1} delta(1 -sum_i=1,K theta_i)
1877
1878 The normalization factor Z can be expressed in terms of gamma functions:
1879
1880 Z = {prod_i=1,K Gamma(alpha_i)} / {Gamma( sum_i=1,K alpha_i)}
1881
1882 The K constants, alpha_1,...,alpha_K, must be positive. The K parameters,
1883 theta_1,...,theta_K are nonnegative and sum to 1.
1884
1885 Status:
1886 Alpha
1887 """
1888 __slots__ = 'alpha', '_total', '_mean',
1889
1890
1891
1892
1893 def __init__(self, alpha) :
1894 """
1895 Args:
1896 - alpha -- The parameters of the Dirichlet prior distribution.
1897 A vector of non-negative real numbers.
1898 """
1899 # TODO: Check that alphas are positive
1900 #TODO : what if alpha's not one dimensional?
1901 self.alpha = asarray(alpha, float64)
1902
1903 self._total = sum(alpha)
1904 self._mean = None
1905
1906
1907 def sample(self) :
1908 """Return a randomly generated probability vector.
1909
1910 Random samples are generated by sampling K values from gamma
1911 distributions with parameters a=\alpha_i, b=1, and renormalizing.
1912
1913 Ref:
1914 A.M. Law, W.D. Kelton, Simulation Modeling and Analysis (1991).
1915 Authors:
1916 Gavin E. Crooks <gec@compbio.berkeley.edu> (2002)
1917 """
1918 alpha = self.alpha
1919 K = len(alpha)
1920 theta = zeros( (K,), float64)
1921
1922 for k in range(K):
1923 theta[k] = random.gammavariate(alpha[k], 1.0)
1924 theta /= sum(theta)
1925
1926 return theta
1927
1928 def mean(self) :
1929 if self._mean ==None:
1930 self._mean = self.alpha / self._total
1931 return self._mean
1932
1933 def covariance(self) :
1934 alpha = self.alpha
1935 A = sum(alpha)
1936 #A2 = A * A
1937 K = len(alpha)
1938 cv = zeros( (K,K), float64)
1939
1940 for i in range(K) :
1941 cv[i,i] = alpha[i] * (1. - alpha[i]/A) / (A * (A+1.) )
1942
1943 for i in range(K) :
1944 for j in range(i+1,K) :
1945 v = - alpha[i] * alpha[j] / (A * A * (A+1.) )
1946 cv[i,j] = v
1947 cv[j,i] = v
1948 return cv
1949
1950 def mean_x(self, x) :
1951 x = asarray(x, float64)
1952 if shape(x) != shape(self.alpha) :
1953 raise ValueError("Argument must be same dimension as Dirichlet")
1954 return sum( x * self.mean())
1955
1956 def variance_x(self, x) :
1957 x = asarray(x, float64)
1958 if shape(x) != shape(self.alpha) :
1959 raise ValueError("Argument must be same dimension as Dirichlet")
1960
1961 cv = self.covariance()
1962 var = na.dot(na.dot(na.transpose( x), cv), x)
1963 return var
1964
1965
1966 def mean_entropy(self) :
1967 """Calculate the average entropy of probabilities sampled
1968 from this Dirichlet distribution.
1969
1970 Returns:
1971 The average entropy.
1972
1973 Ref:
1974 Wolpert & Wolf, PRE 53:6841-6854 (1996) Theorem 7
1975 (Warning: this paper contains typos.)
1976 Status:
1977 Alpha
1978 Authors:
1979 GEC 2005
1980
1981 """
1982 # TODO: Optimize
1983 alpha = self.alpha
1984 A = float(sum(alpha))
1985 ent = 0.0
1986 for a in alpha:
1987 if a>0 : ent += - 1.0 * a * digamma( 1.0+a) # FIXME: Check
1988 ent /= A
1989 ent += digamma(A+1.0)
1990 return ent
1991
1992
1993
1994 def variance_entropy(self):
1995 """Calculate the variance of the Dirichlet entropy.
1996
1997 Ref:
1998 Wolpert & Wolf, PRE 53:6841-6854 (1996) Theorem 8
1999 (Warning: this paper contains typos.)
2000 """
2001 alpha = self.alpha
2002 A = float(sum(alpha))
2003 A2 = A * (A+1)
2004 L = len(alpha)
2005
2006 dg1 = zeros( (L) , float64)
2007 dg2 = zeros( (L) , float64)
2008 tg2 = zeros( (L) , float64)
2009
2010 for i in range(L) :
2011 dg1[i] = digamma(alpha[i] + 1.0)
2012 dg2[i] = digamma(alpha[i] + 2.0)
2013 tg2[i] = trigamma(alpha[i] + 2.0)
2014
2015 dg_Ap2 = digamma( A+2. )
2016 tg_Ap2 = trigamma( A+2. )
2017
2018 mean = self.mean_entropy()
2019 var = 0.0
2020
2021 for i in range(L) :
2022 for j in range(L) :
2023 if i != j :
2024 var += (
2025 ( dg1[i] - dg_Ap2 ) * (dg1[j] - dg_Ap2 ) - tg_Ap2
2026 ) * (alpha[i] * alpha[j] ) / A2
2027 else :
2028 var += (
2029 ( dg2[i] - dg_Ap2 ) **2 + ( tg2[i] - tg_Ap2 )
2030 ) * ( alpha[i] * (alpha[i]+1.) ) / A2
2031
2032 var -= mean**2
2033 return var
2034
2035
2036
2037 def mean_relative_entropy(self, pvec) :
2038 ln_p = na.log(pvec)
2039 return - self.mean_x(ln_p) - self.mean_entropy()
2040
2041
2042 def variance_relative_entropy(self, pvec) :
2043 ln_p = na.log(pvec)
2044 return self.variance_x(ln_p) + self.variance_entropy()
2045
2046
2047 def interval_relative_entropy(self, pvec, frac) :
2048 mean = self.mean_relative_entropy(pvec)
2049 variance = self.variance_relative_entropy(pvec)
2050 # If the variance is small, use the standard 95%
2051 # confidence interval: mean +/- 1.96 * sd
2052 if variance< 0.1 :
2053 sd = sqrt(variance)
2054 return max(0.0, mean - sd*1.96), mean + sd*1.96
2055 sd = sqrt(variance)
2056 return max(0.0, mean - sd*1.96), mean + sd*1.96
2057
2058 g = gamma.from_mean_variance(mean, variance)
2059 low_limit = g.inverse_cdf( (1.-frac)/2.)
2060 high_limit = g.inverse_cdf( 1. - (1.-frac)/2. )
2061
2062 return low_limit, high_limit
2063
2064
2065 # Standard python voodoo for CLI
2066 if __name__ == "__main__":
2067 ## Code Profiling. Uncomment these lines
2068 #import hotshot, hotshot.stats
2069 #prof = hotshot.Profile("stones.prof")
2070 #prof.runcall(main)
2071 #prof.close()
2072 #stats = hotshot.stats.load("stones.prof")
2073 #stats.strip_dirs()
2074 #stats.sort_stats('cumulative', 'calls')
2075 #stats.print_stats(40)
2076 #sys.exit()
2077
2078 main()
2079
2080
2081
2082
2083
2084