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