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