Mercurial > repos > davidmurphy > codonlogo
comparison weblogolib/__init__.py @ 7:8d676bbd1f2d
Uploaded
author | davidmurphy |
---|---|
date | Mon, 16 Jan 2012 07:03:36 -0500 |
parents | 4d47ab2b7bcc |
children | 5149eb3a89c2 |
comparison
equal
deleted
inserted
replaced
6:4a4aca3d57c9 | 7:8d676bbd1f2d |
---|---|
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 |