Mercurial > repos > davidmurphy > codonlogo
annotate weblogolib/__init__.py @ 15:981eb8c3a756 default tip
Uploaded
| author | davidmurphy |
|---|---|
| date | Sat, 31 Mar 2012 16:07:07 -0400 |
| parents | 778f03497adb |
| children |
| rev | line source |
|---|---|
| 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' : | |
| 14 | 983 prior = weight * equiprobable_distribution(len(alphabet)) |
| 984 return prior | |
| 985 | |
| 986 if comp.lower() == 'escherichiacoli' : | |
| 10 | 987 if(altype=="codonsT"): |
| 988 composition="{'CTT': 0.7616, 'ATG': 1.5872, 'ACA': 0.4096, 'ACG': 0.736, 'ATC': 1.1648, 'AAC': 1.5615999999999999, 'ATA': 0.2368, 'AGG': 0.1024, 'CCT': 0.5376000000000001, 'ACT': 0.512, 'AGC': 1.0624, 'AAG': 0.7744, 'AGA': 0.0896, 'CAT': 1.0112, 'AAT': 1.4016, 'ATT': 1.952, 'CTG': 3.0016, 'CTA': 0.3392, 'CTC': 0.672, 'CAC': 0.8383999999999999, 'AAA': 2.1248, 'CCG': 1.7087999999999999, 'AGT': 0.4608, 'CCA': 0.4224, 'CAA': 0.7744, 'CCC': 0.4096, 'TAT': 1.0752000000000002, 'GGT': 1.3632, 'TGT': 0.37760000000000005, 'CGA': 0.2752, 'CAG': 1.7728, 'TCT': 0.3648, 'GAT': 2.4255999999999998, 'CGG': 0.26239999999999997, 'TTT': 1.2608, 'TGC': 0.512, 'GGG': 0.5504, 'TAG': 1e-06, 'GGA': 0.5888, 'TAA': 0.1152, 'GGC': 2.1376, 'TAC': 0.9344, 'TTC': 0.96, 'TCG': 0.512, 'TTA': 0.9728, 'TTG': 0.7616, 'TCC': 0.352, 'ACC': 1.4592, 'TCA': 0.4992, 'GCA': 1.3504, 'GTA': 0.736, 'GCC': 2.0224, 'GTC': 0.7487999999999999, 'GCG': 2.464, 'GTG': 1.6896, 'GAG': 1.1776, 'GTT': 1.0752000000000002, 'GCT': 0.6848, 'TGA': 0.064, 'GAC': 1.312, 'CGT': 1.3504, 'TGG': 0.6848, 'GAA': 2.7968, 'CGC': 1.664}" | |
| 989 else: | |
| 990 composition="{'CUU': 0.7616, 'AUG': 1.5872, 'ACA': 0.4096, 'ACG': 0.736, 'AUC': 1.1648, 'AAC': 1.5615999999999999, 'AUA': 0.2368, 'AGG': 0.1024, 'CCU': 0.5376000000000001, 'ACU': 0.512, 'AGC': 1.0624, 'AAG': 0.7744, 'AGA': 0.0896, 'CAU': 1.0112, 'AAU': 1.4016, 'AUU': 1.952, 'CUG': 3.0016, 'CUA': 0.3392, 'CUC': 0.672, 'CAC': 0.8383999999999999, 'AAA': 2.1248, 'CCG': 1.7087999999999999, 'AGU': 0.4608, 'CCA': 0.4224, 'CAA': 0.7744, 'CCC': 0.4096, 'UAU': 1.0752000000000002, 'GGU': 1.3632, 'UGU': 0.37760000000000005, 'CGA': 0.2752, 'CAG': 1.7728, 'UCU': 0.3648, 'GAU': 2.4255999999999998, 'CGG': 0.26239999999999997, 'UUU': 1.2608, 'UGC': 0.512, 'GGG': 0.5504, 'UAG': 1e-06, 'GGA': 0.5888, 'UAA': 0.1152, 'GGC': 2.1376, 'UAC': 0.9344, 'UUC': 0.96, 'UCG': 0.512, 'UUA': 0.9728, 'UUG': 0.7616, 'UCC': 0.352, 'ACC': 1.4592, 'UCA': 0.4992, 'GCA': 1.3504, 'GUA': 0.736, 'GCC': 2.0224, 'GUC': 0.7487999999999999, 'GCG': 2.464, 'GUG': 1.6896, 'GAG': 1.1776, 'GUU': 1.0752000000000002, 'GCU': 0.6848, 'UGA': 0.064, 'GAC': 1.312, 'CGU': 1.3504, 'UGG': 0.6848, 'GAA': 2.7968, 'CGC': 1.664}" | |
|
9
f3462128e87c
Minor alterations to the galaxy interface with some better examples and error messages added.
davidmurphy
parents:
8
diff
changeset
|
991 elif comp.lower() == 'homosapiens' : |
| 10 | 992 if(altype=="codonsT"): |
| 993 composition="{'CTT': 0.8448, 'ATG': 1.408, 'ACA': 0.9663999999999999, 'ACG': 0.39039999999999997, 'ATC': 1.3312, 'AAC': 1.2224000000000002, 'ATA': 0.48, 'AGG': 0.768, 'CCT': 1.12, 'ACT': 0.8383999999999999, 'AGC': 1.248, 'AAG': 2.0416, 'AGA': 0.7807999999999999, 'CAT': 0.6976, 'AAT': 1.088, 'ATT': 1.024, 'CTG': 2.5344, 'CTA': 0.4608, 'CTC': 1.2544000000000002, 'CAC': 0.9663999999999999, 'AAA': 1.5615999999999999, 'CCG': 0.44160000000000005, 'AGT': 0.7744, 'CCA': 1.0816, 'CAA': 0.7872, 'CCC': 1.2672, 'TAT': 0.7807999999999999, 'GGT': 0.6912, 'TGT': 0.6784, 'CGA': 0.3968, 'CAG': 2.1888, 'TCT': 0.9728, 'GAT': 1.3952, 'CGG': 0.7296, 'TTT': 1.1264, 'TGC': 0.8064, 'GGG': 1.056, 'TAG': 0.0512, 'GGA': 1.056, 'TAA': 0.064, 'GGC': 1.4208, 'TAC': 0.9792000000000001, 'TTC': 1.2992000000000001, 'TCG': 0.2816, 'TTA': 0.4928, 'TTG': 0.8256, 'TCC': 1.1328, 'ACC': 1.2096, 'TCA': 0.7807999999999999, 'GCA': 1.0112, 'GTA': 0.45439999999999997, 'GCC': 1.7728, 'GTC': 0.928, 'GCG': 0.4736, 'GTG': 1.7984, 'GAG': 2.5344, 'GTT': 0.704, 'GCT': 1.1776, 'TGA': 0.1024, 'GAC': 1.6064, 'CGT': 0.288, 'TGG': 0.8448, 'GAA': 1.856, 'CGC': 0.6656}" | |
| 994 else: | |
| 995 composition="{'CUU': 0.8448, 'AUG': 1.408, 'ACA': 0.9663999999999999, 'ACG': 0.39039999999999997, 'AUC': 1.3312, 'AAC': 1.2224000000000002, 'AUA': 0.48, 'AGG': 0.768, 'CCU': 1.12, 'ACU': 0.8383999999999999, 'AGC': 1.248, 'AAG': 2.0416, 'AGA': 0.7807999999999999, 'CAU': 0.6976, 'AAU': 1.088, 'AUU': 1.024, 'CUG': 2.5344, 'CUA': 0.4608, 'CUC': 1.2544000000000002, 'CAC': 0.9663999999999999, 'AAA': 1.5615999999999999, 'CCG': 0.44160000000000005, 'AGU': 0.7744, 'CCA': 1.0816, 'CAA': 0.7872, 'CCC': 1.2672, 'UAU': 0.7807999999999999, 'GGU': 0.6912, 'UGU': 0.6784, 'CGA': 0.3968, 'CAG': 2.1888, 'UCU': 0.9728, 'GAU': 1.3952, 'CGG': 0.7296, 'UUU': 1.1264, 'UGC': 0.8064, 'GGG': 1.056, 'UAG': 0.0512, 'GGA': 1.056, 'UAA': 0.064, 'GGC': 1.4208, 'UAC': 0.9792000000000001, 'UUC': 1.2992000000000001, 'UCG': 0.2816, 'UUA': 0.4928, 'UUG': 0.8256, 'UCC': 1.1328, 'ACC': 1.2096, 'UCA': 0.7807999999999999, 'GCA': 1.0112, 'GUA': 0.45439999999999997, 'GCC': 1.7728, 'GUC': 0.928, 'GCG': 0.4736, 'GUG': 1.7984, 'GAG': 2.5344, 'GUU': 0.704, 'GCU': 1.1776, 'UGA': 0.1024, 'GAC': 1.6064, 'CGU': 0.288, 'UGG': 0.8448, 'GAA': 1.856, 'CGC': 0.6656}" | |
| 996 elif comp.lower() == 'saccharomycescerevisiae' : | |
| 997 if(altype=="codonsT"): | |
| 998 composition="{'CTT': 0.7872, 'ATG': 1.3376, 'ACA': 1.1392, 'ACG': 0.512, 'ATC': 1.1008, 'AAC': 1.5872, 'ATA': 1.1392, 'AGG': 0.5888, 'CCT': 0.864, 'ACT': 1.2992000000000001, 'AGC': 0.6272000000000001, 'AAG': 1.9712, 'AGA': 1.3632, 'CAT': 0.8704, 'AAT': 2.2848, 'ATT': 1.9264000000000001, 'CTG': 0.672, 'CTA': 0.8576, 'CTC': 0.3456, 'CAC': 0.4992, 'AAA': 2.6816, 'CCG': 0.3392, 'AGT': 0.9087999999999999, 'CCA': 1.1712, 'CAA': 1.7472, 'CCC': 0.4352, 'TAT': 1.2032, 'GGT': 1.5295999999999998, 'TGT': 0.5184, 'CGA': 0.192, 'CAG': 0.7744, 'TCT': 1.504, 'GAT': 2.4064, 'CGG': 0.1088, 'TTT': 1.6704, 'TGC': 0.3072, 'GGG': 0.384, 'TAG': 0.032, 'GGA': 0.6976, 'TAA': 0.0704, 'GGC': 0.6272000000000001, 'TAC': 0.9472, 'TTC': 1.1776, 'TCG': 0.5504, 'TTA': 1.6767999999999998, 'TTG': 1.7408, 'TCC': 0.9087999999999999, 'ACC': 0.8128, 'TCA': 1.1967999999999999, 'GCA': 1.0368, 'GTA': 0.7552000000000001, 'GCC': 0.8064, 'GTC': 0.7552000000000001, 'GCG': 0.3968, 'GTG': 0.6912, 'GAG': 1.2288, 'GTT': 1.4144, 'GCT': 1.3568, 'TGA': 0.0448, 'GAC': 1.2928, 'CGT': 0.4096, 'TGG': 0.6656, 'GAA': 2.9184, 'CGC': 0.1664}" | |
| 999 else: | |
| 1000 composition="{'CUU': 0.7872, 'AUG': 1.3376, 'ACA': 1.1392, 'ACG': 0.512, 'AUC': 1.1008, 'AAC': 1.5872, 'AUA': 1.1392, 'AGG': 0.5888, 'CCU': 0.864, 'ACU': 1.2992000000000001, 'AGC': 0.6272000000000001, 'AAG': 1.9712, 'AGA': 1.3632, 'CAU': 0.8704, 'AAU': 2.2848, 'AUU': 1.9264000000000001, 'CUG': 0.672, 'CUA': 0.8576, 'CUC': 0.3456, 'CAC': 0.4992, 'AAA': 2.6816, 'CCG': 0.3392, 'AGU': 0.9087999999999999, 'CCA': 1.1712, 'CAA': 1.7472, 'CCC': 0.4352, 'UAU': 1.2032, 'GGU': 1.5295999999999998, 'UGU': 0.5184, 'CGA': 0.192, 'CAG': 0.7744, 'UCU': 1.504, 'GAU': 2.4064, 'CGG': 0.1088, 'UUU': 1.6704, 'UGC': 0.3072, 'GGG': 0.384, 'UAG': 0.032, 'GGA': 0.6976, 'UAA': 0.0704, 'GGC': 0.6272000000000001, 'UAC': 0.9472, 'UUC': 1.1776, 'UCG': 0.5504, 'UUA': 1.6767999999999998, 'UUG': 1.7408, 'UCC': 0.9087999999999999, 'ACC': 0.8128, 'UCA': 1.1967999999999999, 'GCA': 1.0368, 'GUA': 0.7552000000000001, 'GCC': 0.8064, 'GUC': 0.7552000000000001, 'GCG': 0.3968, 'GUG': 0.6912, 'GAG': 1.2288, 'GUU': 1.4144, 'GCU': 1.3568, 'UGA': 0.0448, 'GAC': 1.2928, 'CGU': 0.4096, 'UGG': 0.6656, 'GAA': 2.9184, 'CGC': 0.1664}" | |
| 14 | 1001 |
| 1002 | |
| 1003 if comp.lower() == 'auto' or comp.lower() == 'automatic': | |
| 4 | 1004 if alphabet == unambiguous_protein_alphabet : |
| 1005 prior = weight * asarray(aa_composition, float64) | |
| 1006 else : | |
| 1007 prior = weight * equiprobable_distribution(len(alphabet)) | |
| 1008 elif comp in std_percentCG : | |
| 1009 prior = weight * base_distribution(std_percentCG[comp]) | |
| 13 | 1010 return prior |
| 4 | 1011 elif comp[-1] == '%' : |
| 1012 prior = weight * base_distribution( float(comp[:-1])) | |
| 13 | 1013 return prior |
| 4 | 1014 elif isfloat(comp) : |
| 1015 prior = weight * base_distribution( float(comp)*100. ) | |
| 13 | 1016 return prior |
| 14 | 1017 elif composition[0] == '{' and composition[-1] == '}' : |
| 4 | 1018 explicit = composition[1: -1] |
| 1019 explicit = explicit.replace(',',' ').replace("'", ' ').replace('"',' ').replace(':', ' ').split() | |
| 1020 | |
| 1021 if len(explicit) != len(alphabet)*2 : | |
| 1022 #print explicit | |
| 1023 raise ValueError("Explicit prior does not match length of alphabet") | |
| 1024 prior = - ones(len(alphabet), float64) | |
| 1025 try : | |
| 1026 for r in range(len(explicit)/2) : | |
| 1027 letter = explicit[r*2] | |
| 1028 index = alphabet.index(letter) | |
| 1029 value = float(explicit[r*2 +1]) | |
| 1030 prior[index] = value | |
| 1031 except ValueError : | |
| 1032 raise ValueError("Cannot parse explicit composition") | |
| 1033 | |
| 1034 if any(prior==-1.) : | |
| 1035 raise ValueError("Explicit prior does not match alphabet") | |
| 1036 prior/= sum(prior) | |
| 1037 prior *= weight | |
| 1038 | |
| 1039 | |
| 1040 else : | |
| 1041 raise ValueError("Unknown or malformed composition: %s"%composition) | |
| 1042 if len(prior) != len(alphabet) : | |
| 1043 raise ValueError( | |
| 1044 "The sequence alphabet and composition are incompatible.") | |
| 1045 | |
| 1046 return prior | |
| 1047 | |
| 1048 | |
| 1049 def base_distribution(percentCG) : | |
| 1050 A = (1. - (percentCG/100.))/2. | |
| 1051 C = (percentCG/100.)/2. | |
| 1052 G = (percentCG/100.)/2. | |
| 1053 T = (1. - (percentCG/100))/2. | |
| 1054 return asarray((A,C,G,T), float64) | |
| 1055 | |
| 1056 def equiprobable_distribution( length) : | |
| 1057 return ones( (length), float64) /length | |
| 1058 | |
| 1059 | |
| 1060 | |
| 1061 | |
| 1062 def read_seq_data(fin, input_parser=seq_io.read,alphabet=None, ignore_lower_case=False, max_file_size=0): | |
| 1063 # TODO: Document this method and enviroment variable | |
| 1064 max_file_size =int(os.environ.get("WEBLOGO_MAX_FILE_SIZE", max_file_size)) | |
| 1065 | |
| 1066 # If max_file_size is set, or if fin==stdin (which is non-seekable), we | |
| 1067 # read the data and replace fin with a StringIO object. | |
| 1068 | |
| 1069 if(max_file_size>0) : | |
| 1070 data = fin.read(max_file_size) | |
| 1071 | |
| 1072 more_data = fin.read(2) | |
| 1073 if more_data != "" : | |
| 1074 raise IOError("File exceeds maximum allowed size: %d bytes" % max_file_size) | |
| 1075 | |
| 1076 fin = StringIO(data) | |
| 1077 elif fin == sys.stdin: | |
| 1078 fin = StringIO(fin.read()) | |
| 1079 | |
| 1080 seqs = input_parser(fin) | |
| 1081 | |
| 1082 if seqs is None or len(seqs) ==0 : | |
| 1083 raise ValueError("Please provide a multiple sequence alignment") | |
| 1084 | |
| 1085 if ignore_lower_case : | |
| 1086 # Case is significant. Do not count lower case letters. | |
| 1087 for i,s in enumerate(seqs) : | |
| 1088 seqs[i] = s.mask() | |
| 1089 | |
| 1090 global altype | |
| 1091 if(altype=="codonsT" or altype=="codonsU"): | |
| 1092 if 'T' in seqs[0] or 't' in seqs[0]: | |
| 1093 altype="codonsT" | |
| 1094 if 'U' in seqs[0] or 'u' in seqs[0]: | |
| 1095 altype="codonsU" | |
| 1096 global offset | |
| 1097 global isreversed | |
| 1098 seq2=[""]*len(seqs) | |
| 1099 seq2 = [] | |
| 1100 for i in xrange(len(seqs)): | |
| 1101 seq2.append([]) | |
| 1102 if(offset%6>2): | |
| 1103 for x in range(0,len(seqs)): | |
| 1104 backs=seqs[x][::-1] | |
| 1105 for y in range(0,len(backs)): | |
| 1106 seq2[x].append(str(backs[y])) | |
| 1107 | |
| 1108 if(altype=="codonsU"): | |
| 1109 for x in range(0,len(seq2)): | |
| 1110 for y in range(0,len(seq2[x])): | |
| 1111 if(cmp(seq2[x][y],'G')==0): | |
| 1112 seq2[x][y]="C" | |
| 1113 elif(cmp(seq2[x][y],'A')==0): | |
| 1114 seq2[x][y]='U' | |
| 1115 elif(cmp(seq2[x][y],'U')==0): | |
| 1116 seq2[x][y]='A' | |
| 1117 elif(cmp(seq2[x][y],'C')==0): | |
| 1118 seq2[x][y]='G' | |
| 1119 if(altype=="codonsT"): | |
| 1120 for x in range(0,len(seq2)): | |
| 1121 for y in range(0,len(seq2[x])): | |
| 1122 if(cmp(seq2[x][y],'G')==0): | |
| 1123 seq2[x][y]='C' | |
| 1124 elif(cmp(seq2[x][y],'A')==0): | |
| 1125 seq2[x][y]='T' | |
| 1126 elif(cmp(seq2[x][y],'T')==0): | |
| 1127 seq2[x][y]='A' | |
| 1128 elif(cmp(seq2[x][y],'C')==0): | |
| 1129 seq2[x][y]='G' | |
| 1130 offset=offset%3 | |
| 1131 isreversed=True | |
| 1132 for x in range(0,len(seqs)): | |
| 1133 seqs[x]=Seq("".join(seq2[x])) | |
| 1134 | |
| 1135 | |
| 1136 # Add alphabet to seqs. | |
| 1137 if alphabet : | |
| 1138 seqs.alphabet = alphabet | |
| 1139 else : | |
| 1140 seqs.alphabet = which_alphabet(seqs) | |
| 1141 | |
| 1142 return seqs | |
| 1143 | |
| 1144 | |
| 1145 | |
| 1146 #TODO: Move to seq_io? | |
| 1147 # Would have to check that no symbol outside of full alphabet? | |
| 1148 def which_alphabet(seqs) : | |
| 1149 """ Returns the most appropriate unambiguous protien, rna or dna alphabet | |
| 1150 for a Seq or SeqList. | |
| 1151 """ | |
| 1152 alphabets = (unambiguous_protein_alphabet, | |
| 1153 unambiguous_rna_alphabet, | |
| 1154 unambiguous_dna_alphabet | |
| 1155 ) | |
| 1156 # Heuristic | |
| 1157 # Count occurances of each letter. Downweight longer alphabet. | |
| 1158 #for x in seqs: | |
| 1159 | |
| 1160 | |
| 1161 if( altype=="codonsU"): | |
| 1162 return codon_alphabetU | |
| 1163 if( altype=="codonsT"): | |
| 1164 return codon_alphabetT | |
| 1165 else: | |
| 1166 score = [1.0*asarray(seqs.tally(a)).sum()/sqrt(len(a)) for a in alphabets] | |
| 1167 #print score | |
| 1168 best = argmax(score) # Ties go to last in list. | |
| 1169 a = alphabets[best] | |
| 1170 return a | |
| 1171 | |
| 1172 | |
| 1173 | |
| 1174 class LogoData(object) : | |
| 1175 """The data needed to generate a sequence logo. | |
| 1176 | |
| 1177 - alphabet | |
| 1178 - length | |
| 1179 - counts -- An array of character counts | |
| 1180 - entropy -- The relative entropy of each column | |
| 1181 - entropy_interval -- entropy confidence interval | |
| 1182 """ | |
| 1183 | |
| 1184 def __init__(self, length=None, alphabet = None, counts =None, | |
| 1185 entropy =None, entropy_interval = None, weight=None) : | |
| 1186 """Creates a new LogoData object""" | |
| 1187 self.length = length | |
| 1188 self.alphabet = alphabet | |
| 1189 self.counts = counts | |
| 1190 self.entropy = entropy | |
| 1191 self.entropy_interval = entropy_interval | |
| 1192 self.weight = weight | |
| 1193 | |
| 1194 | |
| 1195 | |
| 1196 #@classmethod | |
| 1197 def from_counts(cls, alphabet, counts, prior= None): | |
| 1198 | |
| 1199 """Build a logodata object from counts.""" | |
| 1200 seq_length, A = counts.shape | |
| 1201 | |
| 1202 if prior is not None: prior = array(prior, float64) | |
| 1203 if prior is None : | |
| 1204 R = log(A) | |
| 1205 ent = zeros( seq_length, float64) | |
| 1206 entropy_interval = None | |
| 1207 for i in range (0, seq_length) : | |
| 1208 C = sum(counts[i]) | |
| 1209 #FIXME: fixup corebio.moremath.entropy()? | |
| 1210 if C == 0 : | |
| 1211 ent[i] = 0.0 | |
| 1212 else : | |
| 1213 ent[i] = R - entropy(counts[i]) | |
| 1214 else : | |
| 1215 ent = zeros( seq_length, float64) | |
| 1216 entropy_interval = zeros( (seq_length,2) , float64) | |
| 1217 | |
| 1218 R = log(A) | |
| 1219 | |
| 1220 for i in range (0, seq_length) : | |
| 1221 alpha = array(counts[i] , float64) | |
| 1222 alpha += prior | |
| 1223 | |
| 1224 posterior = Dirichlet(alpha) | |
| 1225 ent[i] = posterior.mean_relative_entropy(prior/sum(prior)) | |
| 1226 entropy_interval[i][0], entropy_interval[i][1] = \ | |
| 1227 posterior.interval_relative_entropy(prior/sum(prior), 0.95) | |
| 1228 weight = array( na.sum(counts,axis=1) , float) | |
| 1229 | |
| 1230 weight /= max(weight) | |
| 1231 | |
| 1232 return cls(seq_length, alphabet, counts, ent, entropy_interval, weight) | |
| 1233 from_counts = classmethod(from_counts) | |
| 1234 | |
| 1235 | |
| 1236 #@classmethod | |
| 1237 def from_seqs(cls, seqs, prior= None): | |
| 1238 | |
| 1239 | |
| 1240 alphabet=seqs.alphabet | |
| 1241 | |
| 1242 #get the offset and if it's greater than 2 swap bases to get the reverse compliment and reverse. | |
| 1243 for x in range(0,len(seqs)): | |
| 1244 seqs[x]=seqs[x].upper() | |
| 1245 counter=0 | |
| 1246 | |
| 1247 | |
| 1248 """Build a 2D array from a SeqList, a list of sequences.""" | |
| 1249 # --- VALIDATE DATA --- | |
| 1250 # check that there is at least one sequence of length 1 | |
| 1251 if len(seqs)==0 or len(seqs[0]) ==0: | |
| 1252 raise ValueError("No sequence data found.") | |
| 1253 sys.exit(0) | |
| 1254 # Check sequence lengths | |
| 1255 seq_length = len(seqs[0]) | |
| 1256 for i,s in enumerate(seqs) : | |
| 1257 #print i,s, len(s) | |
| 1258 if seq_length != len(s) : | |
| 1259 raise ArgumentError( | |
| 1260 "Sequence number %d differs in length from the previous sequences" % (i+1) ,'sequences') | |
| 1261 sys.exit(0) | |
| 1262 | |
| 1263 if(altype=="codonsT" or altype=="codonsU"): | |
| 1264 x = [[0]*64 for x in xrange(seq_length/3)] | |
| 1265 counter=offset | |
| 1266 | |
| 1267 while counter+offset<seq_length: | |
| 1268 for i in range(0,len(seqs)): | |
| 1269 if len(str(seqs[i][(counter):(counter+3)]))==3 and len(seqs[i][(counter):(counter+3)].strip("GATUC"))==0 : | |
| 1270 if(str(seqs[i][(counter):(counter+3)]) in alphabet): | |
| 1271 x[counter/3][ (alphabet.index(str(seqs[i][(counter):(counter+3)]))) ]+=1 | |
| 8 | 1272 elif show_warnings: |
| 1273 if len(seqs[i][(counter):(counter+3)].strip("GATUC"))==1 or len(seqs[i][(counter):(counter+3)].strip("GATUC"))==2 : | |
| 1274 print >>sys.stderr, 'Warning:Incomplete or non GATUC codon detected:', seqs[i][(counter):(counter+3)] | |
| 1275 print >>sys.stderr, 'Position:',counter | |
|
9
f3462128e87c
Minor alterations to the galaxy interface with some better examples and error messages added.
davidmurphy
parents:
8
diff
changeset
|
1276 print >>sys.stderr, 'Sequence:',(i+1) |
| 8 | 1277 print >>sys.stderr, 'This will be treated as ---' |
| 1278 | |
| 1279 | |
| 4 | 1280 counter=counter+3 |
| 1281 counts=asarray(x) | |
| 1282 else: | |
| 1283 counts = asarray(seqs.tally()) | |
| 1284 | |
| 1285 return cls.from_counts(alphabet, counts, prior) | |
| 1286 from_seqs = classmethod(from_seqs) | |
| 1287 | |
| 1288 | |
| 1289 | |
| 1290 def __str__(self) : | |
| 1291 out = StringIO() | |
| 1292 print >>out, '## LogoData' | |
| 1293 print >>out, '# First column is position number, couting from zero' | |
| 1294 print >>out, '# Subsequent columns are raw symbol counts' | |
| 1295 print >>out, '# Entropy is mean entropy measured in nats.' | |
| 1296 print >>out, '# Low and High are the 95% confidence limits.' | |
| 1297 print >>out, '# Weight is the fraction of non-gap symbols in the column.' | |
| 1298 print >>out, '#\t' | |
| 1299 print >>out, '#\t', | |
| 1300 for a in self.alphabet : | |
| 1301 print >>out, a, '\t', | |
| 1302 print >>out, 'Entropy\tLow\tHigh\tWeight' | |
| 1303 | |
| 1304 for i in range(self.length) : | |
| 1305 print >>out, i, '\t', | |
| 1306 for c in self.counts[i] : print >>out, c, '\t', | |
| 1307 print >>out, self.entropy[i], '\t', | |
| 1308 if self.entropy_interval is not None: | |
| 1309 print >>out, self.entropy_interval[i][0], '\t', | |
| 1310 print >>out, self.entropy_interval[i][1], '\t', | |
| 1311 else : | |
| 1312 print >>out, '\t','\t', | |
| 1313 if self.weight is not None : | |
| 1314 print >>out, self.weight[i], | |
| 1315 print >>out, '' | |
| 1316 print >>out, '# End LogoData' | |
| 1317 | |
| 1318 return out.getvalue() | |
| 1319 | |
| 1320 # ====================== Main: Parse Command line ============================= | |
| 1321 def main(): | |
| 1322 """CodonLogo command line interface """ | |
| 1323 | |
| 1324 # ------ Parse Command line ------ | |
| 1325 parser = _build_option_parser() | |
| 1326 (opts, args) = parser.parse_args(sys.argv[1:]) | |
| 1327 if args : parser.error("Unparsable arguments: %s " % args) | |
| 1328 | |
| 1329 if opts.serve: | |
| 1330 httpd_serve_forever(opts.port) # Never returns? | |
| 1331 sys.exit(0) | |
| 1332 | |
| 1333 | |
| 1334 # ------ Create Logo ------ | |
| 1335 try: | |
| 1336 data = _build_logodata(opts) | |
| 1337 | |
| 1338 | |
| 1339 format = _build_logoformat(data, opts) | |
| 1340 | |
| 1341 | |
| 1342 formatter = opts.formatter | |
| 1343 formatter(data, format, opts.fout) | |
| 1344 | |
| 1345 except ValueError, err : | |
| 1346 print >>sys.stderr, 'Error:', err | |
| 1347 sys.exit(2) | |
| 1348 except KeyboardInterrupt, err: | |
| 1349 sys.exit(0) | |
| 1350 # End main() | |
| 1351 | |
| 1352 | |
| 1353 def httpd_serve_forever(port=8080) : | |
| 1354 """ Start a webserver on a local port.""" | |
| 1355 import BaseHTTPServer | |
| 1356 import CGIHTTPServer | |
| 1357 | |
| 1358 class __HTTPRequestHandler(CGIHTTPServer.CGIHTTPRequestHandler): | |
| 1359 def is_cgi(self) : | |
| 1360 if self.path == "/create.cgi": | |
| 1361 self.cgi_info = '', 'create.cgi' | |
| 1362 return True | |
| 1363 return False | |
| 1364 | |
| 1365 # Add current directory to PYTHONPATH. This is | |
| 1366 # so that we can run the standalone server | |
| 1367 # without having to run the install script. | |
| 1368 pythonpath = os.getenv("PYTHONPATH", '') | |
| 1369 pythonpath += ":" + os.path.abspath(sys.path[0]).split()[0] | |
| 1370 os.environ["PYTHONPATH"] = pythonpath | |
| 1371 | |
| 1372 htdocs = resource_filename(__name__, 'htdocs', __file__) | |
| 1373 os.chdir(htdocs) | |
| 1374 | |
| 1375 HandlerClass = __HTTPRequestHandler | |
| 1376 ServerClass = BaseHTTPServer.HTTPServer | |
| 1377 httpd = ServerClass(('', port), HandlerClass) | |
| 1378 print "Serving HTTP on localhost:%d ..." % port | |
| 1379 | |
| 1380 try : | |
| 1381 httpd.serve_forever() | |
| 1382 except KeyboardInterrupt: | |
| 1383 sys.exit(0) | |
| 1384 # end httpd_serve_forever() | |
| 1385 | |
| 1386 def read_priors(finp, alphabet ,max_file_size=0): | |
| 1387 | |
| 1388 max_file_size =int(os.environ.get("WEBLOGO_MAX_FILE_SIZE", max_file_size)) | |
| 1389 if(max_file_size>0) : | |
| 1390 data = finp.read(max_file_size) | |
| 1391 more_data = finp.read(2) | |
| 1392 if more_data != "" : | |
| 1393 raise IOError("File exceeds maximum allowed size: %d bytes" % max_file_size) | |
| 1394 finp = StringIO(data) | |
| 1395 priordict={} | |
| 1396 while 1: | |
| 1397 line = finp.readline() | |
| 1398 if not line: | |
| 1399 break | |
| 1400 line = line.split() | |
| 8 | 1401 |
| 1402 if(altype=="codonsT"): | |
| 1403 priordict[line[0].upper().replace("U", "T")]=(float(line[1])/1000)*64 | |
| 1404 else: | |
| 1405 priordict[line[0].upper().replace("T", "U")]=(float(line[1])/1000)*64 | |
|
9
f3462128e87c
Minor alterations to the galaxy interface with some better examples and error messages added.
davidmurphy
parents:
8
diff
changeset
|
1406 if priordict[line[0].upper().replace("U", "T")] == 0: |
|
f3462128e87c
Minor alterations to the galaxy interface with some better examples and error messages added.
davidmurphy
parents:
8
diff
changeset
|
1407 priordict[line[0].upper().replace("U", "T")] = 0.000001 |
| 8 | 1408 return priordict |
| 4 | 1409 |
| 1410 def _build_logodata(options) : | |
| 1411 global offset | |
| 1412 offset=options.frame | |
| 8 | 1413 global show_warnings |
| 1414 show_warnings = options.strict | |
| 4 | 1415 options.alphabet = None |
| 1416 options.ignore_lower_case = False | |
| 1417 #options.default_color = Color.by_name("black") | |
| 1418 options.color_scheme=None | |
| 1419 #options.colors=[] | |
| 1420 options.show_ends=False | |
| 1421 seqs = read_seq_data(options.fin, | |
| 1422 options.input_parser.read, | |
| 1423 alphabet=options.alphabet, | |
| 8 | 1424 ignore_lower_case = options.ignore_lower_case) |
| 4 | 1425 if(options.priorfile!=None): |
| 1426 if(altype=="CodonsT"): | |
| 1427 options.composition= str(read_priors(options.priorfile,codon_alphabetT)) | |
| 1428 options.alphabet = codon_alphabetT | |
| 1429 else: | |
| 1430 options.composition= str(read_priors(options.priorfile,codon_alphabetU)) | |
| 1431 options.alphabet = codon_alphabetU | |
| 1432 | |
| 1433 prior = parse_prior( options.composition,seqs.alphabet, options.weight) | |
| 1434 data = LogoData.from_seqs(seqs, prior) | |
| 1435 return data | |
| 1436 | |
| 1437 | |
| 1438 def _build_logoformat( logodata, opts) : | |
| 1439 """ Extract and process relevant option values and return a | |
| 1440 LogoFormat object.""" | |
| 1441 | |
| 1442 args = {} | |
| 1443 direct_from_opts = [ | |
| 1444 "stacks_per_line", | |
| 1445 "logo_title", | |
| 1446 "yaxis_label", | |
| 1447 "show_xaxis", | |
| 8 | 1448 "strict", |
| 4 | 1449 "show_yaxis", |
| 1450 "xaxis_label", | |
| 1451 "show_ends", | |
| 1452 "fineprint", | |
| 1453 "show_errorbars", | |
| 1454 "show_boxes", | |
| 1455 "yaxis_tic_interval", | |
| 1456 "resolution", | |
| 1457 "alphabet", | |
| 1458 "debug", | |
| 1459 "show_ends", | |
| 1460 "default_color", | |
| 1461 #"show_color_key", | |
| 1462 "color_scheme", | |
| 1463 "unit_name", | |
| 1464 "logo_label", | |
| 1465 "yaxis_scale", | |
| 1466 "first_index", | |
| 1467 "logo_start", | |
| 1468 "logo_end", | |
| 1469 "scale_width", | |
| 1470 "frame", | |
| 1471 ] | |
| 1472 | |
| 1473 for k in direct_from_opts: | |
| 1474 args[k] = opts.__dict__[k] | |
| 1475 logo_size = copy.copy(opts.__dict__['logo_size']) | |
| 1476 size_from_opts = ["stack_width", "stack_height"] | |
| 1477 for k in size_from_opts : | |
| 1478 length = getattr(opts, k) | |
| 1479 if length : setattr( logo_size, k, length ) | |
| 1480 args["size"] = logo_size | |
| 1481 | |
| 1482 global col | |
| 1483 | |
| 1484 if opts.colors: | |
| 1485 color_scheme = ColorScheme() | |
| 1486 for color, symbols, desc in opts.colors: | |
| 1487 try : | |
| 1488 #c = Color.from_string(color) | |
| 1489 color_scheme.groups.append( ColorGroup(symbols, color, desc) ) | |
| 1490 #print >> sys.stderr, color | |
| 1491 | |
| 1492 col.append( ColorGroup(symbols, color, desc) ) | |
| 1493 | |
| 1494 except ValueError : | |
| 1495 raise ValueError( | |
| 1496 "error: option --color: invalid value: '%s'" % color ) | |
| 1497 if(altype!="codonsU" and altype!="codonsT") : | |
| 1498 args["color_scheme"] = color_scheme | |
| 1499 | |
| 1500 #cf = colorscheme.format_color(col[0]) | |
| 1501 #col.append( " ("+group.symbols+") " + cf ) | |
| 1502 | |
| 1503 logooptions = LogoOptions() | |
| 1504 for a, v in args.iteritems() : | |
| 1505 setattr(logooptions,a,v) | |
| 1506 | |
| 1507 | |
| 1508 | |
| 1509 theformat = LogoFormat(logodata, logooptions ) | |
| 1510 | |
| 1511 return theformat | |
| 1512 | |
| 1513 | |
| 1514 | |
| 1515 | |
| 1516 | |
| 1517 | |
| 1518 # ========================== OPTIONS ========================== | |
| 1519 def _build_option_parser() : | |
| 1520 defaults = LogoOptions() | |
| 1521 parser = DeOptionParser(usage="%prog [options] < sequence_data.fa > sequence_logo.eps", | |
| 1522 description = description, | |
| 1523 version = __version__ , | |
| 1524 add_verbose_options = False | |
| 1525 ) | |
| 1526 | |
| 1527 io_grp = OptionGroup(parser, "Input/Output Options",) | |
| 1528 data_grp = OptionGroup(parser, "Logo Data Options",) | |
| 1529 format_grp = OptionGroup(parser, "Logo Format Options", | |
| 1530 "These options control the format and display of the logo.") | |
| 1531 color_grp = OptionGroup(parser, "Color Options", | |
| 1532 "Colors can be specified using CSS2 syntax. e.g. 'red', '#FF0000', etc.") | |
| 1533 advanced_grp = OptionGroup(parser, "Advanced Format Options", | |
| 1534 "These options provide fine control over the display of the logo. ") | |
| 1535 server_grp = OptionGroup(parser, "CodonLogo Server", | |
| 1536 "Run a standalone webserver on a local port.") | |
| 1537 | |
| 1538 | |
| 1539 parser.add_option_group(io_grp) | |
| 1540 parser.add_option_group(data_grp) | |
| 1541 parser.add_option_group(format_grp) | |
| 1542 parser.add_option_group(color_grp) | |
| 1543 parser.add_option_group(advanced_grp) | |
| 1544 parser.add_option_group(server_grp) | |
| 1545 | |
| 1546 # ========================== IO OPTIONS ========================== | |
| 1547 | |
| 1548 | |
| 1549 | |
| 1550 io_grp.add_option( "-f", "--fin", | |
| 1551 dest="fin", | |
| 1552 action="store", | |
| 1553 type="file_in", | |
| 1554 default=sys.stdin, | |
| 1555 help="Sequence input file (default: stdin)", | |
| 1556 metavar="FILENAME") | |
| 1557 | |
| 1558 io_grp.add_option( "-R", "--prior", | |
| 1559 dest="priorfile", | |
| 1560 action="store", | |
| 1561 type="file_in", | |
| 1562 help="A file with 64 codons and their prior probabilities, one per line, each codon followed by a space and it's probability.", | |
| 1563 metavar="FILENAME") | |
| 1564 | |
| 1565 io_grp.add_option("", "--fin-format", | |
| 1566 dest="input_parser", | |
| 1567 action="store", type ="dict", | |
| 1568 default = seq_io, | |
| 1569 choices = seq_io.format_names(), | |
| 1570 help="Multiple sequence alignment format: (%s)" % | |
| 1571 ', '.join([ f.names[0] for f in seq_io.formats]), | |
| 1572 metavar="FORMAT") | |
| 1573 | |
| 1574 io_grp.add_option("-o", "--fout", dest="fout", | |
| 1575 type="file_out", | |
| 1576 default=sys.stdout, | |
| 1577 help="Output file (default: stdout)", | |
| 1578 metavar="FILENAME") | |
| 1579 | |
| 1580 | |
| 1581 | |
| 1582 io_grp.add_option( "-F", "--format", | |
| 1583 dest="formatter", | |
| 1584 action="store", | |
| 1585 type="dict", | |
| 1586 choices = formatters, | |
| 1587 metavar= "FORMAT", | |
| 1588 help="Format of output: eps (default), png, png_print, pdf, jpeg, txt", | |
| 1589 default = default_formatter) | |
| 1590 | |
| 1591 | |
| 1592 # ========================== Data OPTIONS ========================== | |
| 1593 | |
| 1594 data_grp.add_option("-m", "--frame", | |
| 1595 dest="frame", | |
| 1596 action="store", | |
| 1597 type="int", | |
| 1598 default=0, | |
| 1599 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.", | |
| 1600 metavar="NUMBER") | |
| 1601 | |
| 1602 #data_grp.add_option("-T", "--type", | |
| 1603 #dest="altype", | |
| 1604 #action="store", | |
| 1605 #type="boolean", | |
| 1606 #default=True, | |
| 1607 #help="Generate a codon logo rather than a sequence logo (default: True)", | |
| 1608 #metavar="YES/NO") | |
| 1609 | |
| 1610 | |
| 1611 | |
| 1612 #data_grp.add_option( "-A", "--sequence-type", | |
| 1613 #dest="alphabet", | |
| 1614 #action="store", | |
| 1615 #type="dict", | |
| 1616 #choices = std_alphabets, | |
| 1617 #help="The type of sequence data: 'protein', 'rna' or 'dna'.", | |
| 1618 #metavar="TYPE") | |
| 1619 | |
| 1620 #data_grp.add_option( "-a", "--alphabet", | |
| 1621 #dest="alphabet", | |
| 1622 #action="store", | |
| 1623 #help="The set of symbols to count, e.g. 'AGTC'. " | |
| 1624 #"All characters not in the alphabet are ignored. " | |
| 1625 #"If neither the alphabet nor sequence-type are specified then codonlogo will examine the input data and make an educated guess. " | |
| 1626 #"See also --sequence-type, --ignore-lower-case" ) | |
| 1627 | |
| 1628 # FIXME Add test? | |
| 1629 #data_grp.add_option( "", "--ignore-lower-case", | |
| 1630 #dest="ignore_lower_case", | |
| 1631 #action="store", | |
| 1632 #type = "boolean", | |
| 1633 #default=False, | |
| 1634 #metavar = "YES/NO", | |
| 1635 #help="Disregard lower case letters and only count upper case letters in sequences? (Default: No)" | |
| 1636 #) | |
| 1637 | |
| 1638 data_grp.add_option( "-U", "--units", | |
| 1639 dest="unit_name", | |
| 1640 action="store", | |
| 1641 choices = std_units.keys(), | |
| 1642 type="choice", | |
| 1643 default = defaults.unit_name, | |
| 1644 help="A unit of entropy ('bits' (default), 'nats', 'digits'), or a unit of free energy ('kT', 'kJ/mol', 'kcal/mol'), or 'probability' for probabilities", | |
| 1645 metavar = "NUMBER") | |
| 1646 | |
| 1647 | |
| 1648 data_grp.add_option( "", "--composition", | |
| 1649 dest="composition", | |
| 1650 action="store", | |
| 1651 type="string", | |
| 13 | 1652 default = "equiprobable", |
| 1653 help="The expected composition of the sequences: 'equiprobable'(default), 'none' (Do not perform any compositional adjustment), or 'escherichiacoli' 'homosapiens' 'saccharomycescerevisiae' for ecoli, human and SC codon frequencies.", | |
| 4 | 1654 metavar="COMP.") |
| 1655 | |
| 1656 data_grp.add_option( "", "--weight", | |
| 1657 dest="weight", | |
| 1658 action="store", | |
| 1659 type="float", | |
| 1660 default = None, | |
| 1661 help="The weight of prior data. Default: total pseudocounts equal to the number of monomer types.", | |
| 1662 metavar="NUMBER") | |
| 1663 | |
| 1664 data_grp.add_option( "-i", "--first-index", | |
| 1665 dest="first_index", | |
| 1666 action="store", | |
| 1667 type="int", | |
| 1668 default = 1, | |
| 1669 help="Index of first position in sequence data (default: 1)", | |
| 1670 metavar="INDEX") | |
| 1671 | |
| 1672 data_grp.add_option( "-l", "--lower", | |
| 1673 dest="logo_start", | |
| 1674 action="store", | |
| 1675 type="int", | |
| 1676 help="Lower bound of sequence to display", | |
| 1677 metavar="INDEX") | |
| 1678 | |
| 1679 data_grp.add_option( "-u", "--upper", | |
| 1680 dest="logo_end", | |
| 1681 action="store", | |
| 1682 type="int", | |
| 1683 help="Upper bound of sequence to display", | |
| 1684 metavar="INDEX") | |
| 8 | 1685 |
| 1686 data_grp.add_option( "-G", "--strict", | |
| 1687 dest="strict", | |
| 1688 action="store", | |
| 1689 type="boolean", | |
| 1690 help="Issue warnings if partial codons are encountered. Default: %default",default = defaults.strict, | |
| 1691 metavar="True/False") | |
| 4 | 1692 # ========================== FORMAT OPTIONS ========================== |
| 1693 | |
| 1694 format_grp.add_option( "-s", "--size", | |
| 1695 dest="logo_size", | |
| 1696 action="store", | |
| 1697 type ="dict", | |
| 1698 choices = std_sizes, | |
| 1699 metavar = "LOGOSIZE", | |
| 1700 default = defaults.size, | |
| 1701 help="Specify a standard logo size (small, medium (default), large)" ) | |
| 1702 | |
| 1703 | |
| 1704 | |
| 1705 format_grp.add_option( "-n", "--stacks-per-line", | |
| 1706 dest="stacks_per_line", | |
| 1707 action="store", | |
| 1708 type="int", | |
| 1709 help="Maximum number of logo stacks per logo line. (default: %default)", | |
| 1710 default = defaults.stacks_per_line, | |
| 1711 metavar="COUNT") | |
| 1712 | |
| 1713 format_grp.add_option( "-t", "--title", | |
| 1714 dest="logo_title", | |
| 1715 action="store", | |
| 1716 type="string", | |
| 1717 help="Logo title text.", | |
| 1718 default = defaults.logo_title, | |
| 1719 metavar="TEXT") | |
| 1720 | |
| 1721 format_grp.add_option( "", "--label", | |
| 1722 dest="logo_label", | |
| 1723 action="store", | |
| 1724 type="string", | |
| 1725 help="A figure label, e.g. '2a'", | |
| 1726 default = defaults.logo_label, | |
| 1727 metavar="TEXT") | |
| 1728 | |
| 1729 format_grp.add_option( "-X", "--show-xaxis", | |
| 1730 action="store", | |
| 1731 type = "boolean", | |
| 1732 default= defaults.show_xaxis, | |
| 1733 metavar = "YES/NO", | |
| 1734 help="Display sequence numbers along x-axis? (default: %default)") | |
| 1735 | |
| 1736 format_grp.add_option( "-x", "--xlabel", | |
| 1737 dest="xaxis_label", | |
| 1738 action="store", | |
| 1739 type="string", | |
| 1740 default = defaults.xaxis_label, | |
| 1741 help="X-axis label", | |
| 1742 metavar="TEXT") | |
| 1743 | |
| 1744 format_grp.add_option( "-S", "--yaxis", | |
| 1745 dest="yaxis_scale", | |
| 1746 action="store", | |
| 1747 type="float", | |
| 1748 help="Height of yaxis in units. (Default: Maximum value with uninformative prior.)", | |
| 1749 metavar = "UNIT") | |
| 1750 | |
| 1751 format_grp.add_option( "-Y", "--show-yaxis", | |
| 1752 action="store", | |
| 1753 type = "boolean", | |
| 1754 dest = "show_yaxis", | |
| 1755 default= defaults.show_yaxis, | |
| 1756 metavar = "YES/NO", | |
| 1757 help="Display entropy scale along y-axis? (default: %default)") | |
| 1758 | |
| 1759 format_grp.add_option( "-y", "--ylabel", | |
| 1760 dest="yaxis_label", | |
| 1761 action="store", | |
| 1762 type="string", | |
| 1763 help="Y-axis label (default depends on plot type and units)", | |
| 1764 metavar="TEXT") | |
| 1765 | |
| 1766 #format_grp.add_option( "-E", "--show-ends", | |
| 1767 #action="store", | |
| 1768 #type = "boolean", | |
| 1769 #default= defaults.show_ends, | |
| 1770 #metavar = "YES/NO", | |
| 1771 #help="Label the ends of the sequence? (default: %default)") | |
| 1772 | |
| 1773 format_grp.add_option( "-P", "--fineprint", | |
| 1774 dest="fineprint", | |
| 1775 action="store", | |
| 1776 type="string", | |
| 1777 default= defaults.fineprint, | |
| 1778 help="The fine print (default: Codonlogo version)", | |
| 1779 metavar="TEXT") | |
| 1780 | |
| 1781 format_grp.add_option( "", "--ticmarks", | |
| 1782 dest="yaxis_tic_interval", | |
| 1783 action="store", | |
| 1784 type="float", | |
| 1785 default= defaults.yaxis_tic_interval, | |
| 1786 help="Distance between ticmarks (default: %default)", | |
| 1787 metavar = "NUMBER") | |
| 1788 | |
| 1789 | |
| 1790 format_grp.add_option( "", "--errorbars", | |
| 1791 dest = "show_errorbars", | |
| 1792 action="store", | |
| 1793 type = "boolean", | |
| 1794 default= defaults.show_errorbars, | |
| 1795 metavar = "YES/NO", | |
| 1796 help="Display error bars? (default: %default)") | |
| 1797 | |
| 1798 | |
| 1799 | |
| 1800 # ========================== Color OPTIONS ========================== | |
| 1801 # TODO: Future Feature | |
| 1802 # color_grp.add_option( "-K", "--color-key", | |
| 1803 # dest= "show_color_key", | |
| 1804 # action="store", | |
| 1805 # type = "boolean", | |
| 1806 # default= defaults.show_color_key, | |
| 1807 # metavar = "YES/NO", | |
| 1808 # help="Display a color key (default: %default)") | |
| 1809 | |
| 1810 | |
| 1811 #color_scheme_choices = std_color_schemes.keys() | |
| 1812 #color_scheme_choices.sort() | |
| 1813 #color_grp.add_option( "-c", "--color-scheme", | |
| 1814 #dest="color_scheme", | |
| 1815 #action="store", | |
| 1816 #type ="dict", | |
| 1817 #choices = std_color_schemes, | |
| 1818 #metavar = "SCHEME", | |
| 1819 #default = None, # Auto | |
| 1820 #help="Specify a standard color scheme (%s)" % \ | |
| 1821 #", ".join(color_scheme_choices) ) | |
| 1822 | |
| 1823 color_grp.add_option( "-C", "--color", | |
| 1824 dest="colors", | |
| 1825 action="append", | |
| 1826 metavar="COLOR SYMBOLS DESCRIPTION ", | |
| 1827 nargs = 3, | |
| 1828 default=[], | |
| 1829 help="Specify symbol colors, e.g. --color black AG 'Purine' --color red TC 'Pyrimidine' ") | |
| 1830 | |
| 1831 color_grp.add_option( "", "--default-color", | |
| 1832 dest="default_color", | |
| 1833 action="store", | |
| 1834 metavar="COLOR", | |
| 1835 default= defaults.default_color, | |
| 1836 help="Symbol color if not otherwise specified.") | |
| 1837 | |
| 1838 # ========================== Advanced options ========================= | |
| 1839 | |
| 1840 advanced_grp.add_option( "-W", "--stack-width", | |
| 1841 dest="stack_width", | |
| 1842 action="store", | |
| 1843 type="float", | |
| 1844 default= None, | |
| 1845 help="Width of a logo stack (default: %s)"% defaults.size.stack_width, | |
| 1846 metavar="POINTS" ) | |
| 1847 | |
| 1848 advanced_grp.add_option( "-H", "--stack-height", | |
| 1849 dest="stack_height", | |
| 1850 action="store", | |
| 1851 type="float", | |
| 1852 default= None, | |
| 1853 help="Height of a logo stack (default: %s)"%defaults.size.stack_height, | |
| 1854 metavar="POINTS" ) | |
| 1855 | |
| 1856 advanced_grp.add_option( "", "--box", | |
| 1857 dest="show_boxes", | |
| 1858 action="store", | |
| 1859 type = "boolean", | |
| 1860 default=False, | |
| 1861 metavar = "YES/NO", | |
| 1862 help="Draw boxes around symbols? (default: no)") | |
| 1863 | |
| 1864 advanced_grp.add_option( "", "--resolution", | |
| 1865 dest="resolution", | |
| 1866 action="store", | |
| 1867 type="float", | |
| 1868 default=96, | |
| 1869 help="Bitmap resolution in dots per inch (DPI). (default: 96 DPI, except png_print, 600 DPI) Low resolution bitmaps (DPI<300) are antialiased.", | |
| 1870 metavar="DPI") | |
| 1871 | |
| 1872 advanced_grp.add_option( "", "--scale-width", | |
| 1873 dest="scale_width", | |
| 1874 action="store", | |
| 1875 type = "boolean", | |
| 1876 default= True, | |
| 1877 metavar = "YES/NO", | |
| 1878 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)") | |
| 1879 | |
| 1880 advanced_grp.add_option( "", "--debug", | |
| 1881 action="store", | |
| 1882 type = "boolean", | |
| 1883 default= defaults.debug, | |
| 1884 metavar = "YES/NO", | |
| 1885 help="Output additional diagnostic information. (default: %default)") | |
| 1886 | |
| 1887 | |
| 1888 # ========================== Server options ========================= | |
| 1889 server_grp.add_option( "", "--serve", | |
| 1890 dest="serve", | |
| 1891 action="store_true", | |
| 1892 default= False, | |
| 1893 help="Start a standalone CodonLogo server for creating sequence logos.") | |
| 1894 | |
| 1895 server_grp.add_option( "", "--port", | |
| 1896 dest="port", | |
| 1897 action="store", | |
| 1898 type="int", | |
| 1899 default= 8080, | |
| 1900 help="Listen to this local port. (Default: %default)", | |
| 1901 metavar="PORT") | |
| 1902 | |
| 1903 return parser | |
| 1904 | |
| 1905 # END _build_option_parser | |
| 1906 | |
| 1907 | |
| 1908 ############################################################## | |
| 1909 | |
| 1910 class Dirichlet(object) : | |
| 1911 """The Dirichlet probability distribution. The Dirichlet is a continuous | |
| 1912 multivariate probability distribution across non-negative unit length | |
| 1913 vectors. In other words, the Dirichlet is a probability distribution of | |
| 1914 probability distributions. It is conjugate to the multinomial | |
| 1915 distribution and is widely used in Bayesian statistics. | |
| 1916 | |
| 1917 The Dirichlet probability distribution of order K-1 is | |
| 1918 | |
| 1919 p(theta_1,...,theta_K) d theta_1 ... d theta_K = | |
| 1920 (1/Z) prod_i=1,K theta_i^{alpha_i - 1} delta(1 -sum_i=1,K theta_i) | |
| 1921 | |
| 1922 The normalization factor Z can be expressed in terms of gamma functions: | |
| 1923 | |
| 1924 Z = {prod_i=1,K Gamma(alpha_i)} / {Gamma( sum_i=1,K alpha_i)} | |
| 1925 | |
| 1926 The K constants, alpha_1,...,alpha_K, must be positive. The K parameters, | |
| 1927 theta_1,...,theta_K are nonnegative and sum to 1. | |
| 1928 | |
| 1929 Status: | |
| 1930 Alpha | |
| 1931 """ | |
| 1932 __slots__ = 'alpha', '_total', '_mean', | |
| 1933 | |
| 1934 | |
| 1935 | |
| 1936 | |
| 1937 def __init__(self, alpha) : | |
| 1938 """ | |
| 1939 Args: | |
| 1940 - alpha -- The parameters of the Dirichlet prior distribution. | |
| 1941 A vector of non-negative real numbers. | |
| 1942 """ | |
| 1943 # TODO: Check that alphas are positive | |
| 1944 #TODO : what if alpha's not one dimensional? | |
| 1945 self.alpha = asarray(alpha, float64) | |
| 1946 | |
| 1947 self._total = sum(alpha) | |
| 1948 self._mean = None | |
| 1949 | |
| 1950 | |
| 1951 def sample(self) : | |
| 1952 """Return a randomly generated probability vector. | |
| 1953 | |
| 1954 Random samples are generated by sampling K values from gamma | |
| 1955 distributions with parameters a=\alpha_i, b=1, and renormalizing. | |
| 1956 | |
| 1957 Ref: | |
| 1958 A.M. Law, W.D. Kelton, Simulation Modeling and Analysis (1991). | |
| 1959 Authors: | |
| 1960 Gavin E. Crooks <gec@compbio.berkeley.edu> (2002) | |
| 1961 """ | |
| 1962 alpha = self.alpha | |
| 1963 K = len(alpha) | |
| 1964 theta = zeros( (K,), float64) | |
| 1965 | |
| 1966 for k in range(K): | |
| 1967 theta[k] = random.gammavariate(alpha[k], 1.0) | |
| 1968 theta /= sum(theta) | |
| 1969 | |
| 1970 return theta | |
| 1971 | |
| 1972 def mean(self) : | |
| 1973 if self._mean ==None: | |
| 1974 self._mean = self.alpha / self._total | |
| 1975 return self._mean | |
| 1976 | |
| 1977 def covariance(self) : | |
| 1978 alpha = self.alpha | |
| 1979 A = sum(alpha) | |
| 1980 #A2 = A * A | |
| 1981 K = len(alpha) | |
| 1982 cv = zeros( (K,K), float64) | |
| 1983 | |
| 1984 for i in range(K) : | |
| 1985 cv[i,i] = alpha[i] * (1. - alpha[i]/A) / (A * (A+1.) ) | |
| 1986 | |
| 1987 for i in range(K) : | |
| 1988 for j in range(i+1,K) : | |
| 1989 v = - alpha[i] * alpha[j] / (A * A * (A+1.) ) | |
| 1990 cv[i,j] = v | |
| 1991 cv[j,i] = v | |
| 1992 return cv | |
| 1993 | |
| 1994 def mean_x(self, x) : | |
| 1995 x = asarray(x, float64) | |
| 1996 if shape(x) != shape(self.alpha) : | |
| 1997 raise ValueError("Argument must be same dimension as Dirichlet") | |
| 1998 return sum( x * self.mean()) | |
| 1999 | |
| 2000 def variance_x(self, x) : | |
| 2001 x = asarray(x, float64) | |
| 2002 if shape(x) != shape(self.alpha) : | |
| 2003 raise ValueError("Argument must be same dimension as Dirichlet") | |
| 2004 | |
| 2005 cv = self.covariance() | |
| 2006 var = na.dot(na.dot(na.transpose( x), cv), x) | |
| 2007 return var | |
| 2008 | |
| 2009 | |
| 2010 def mean_entropy(self) : | |
| 2011 """Calculate the average entropy of probabilities sampled | |
| 2012 from this Dirichlet distribution. | |
| 2013 | |
| 2014 Returns: | |
| 2015 The average entropy. | |
| 2016 | |
| 2017 Ref: | |
| 2018 Wolpert & Wolf, PRE 53:6841-6854 (1996) Theorem 7 | |
| 2019 (Warning: this paper contains typos.) | |
| 2020 Status: | |
| 2021 Alpha | |
| 2022 Authors: | |
| 2023 GEC 2005 | |
| 2024 | |
| 2025 """ | |
| 2026 # TODO: Optimize | |
| 2027 alpha = self.alpha | |
| 2028 A = float(sum(alpha)) | |
| 2029 ent = 0.0 | |
| 2030 for a in alpha: | |
| 2031 if a>0 : ent += - 1.0 * a * digamma( 1.0+a) # FIXME: Check | |
| 2032 ent /= A | |
| 2033 ent += digamma(A+1.0) | |
| 2034 return ent | |
| 2035 | |
| 2036 | |
| 2037 | |
| 2038 def variance_entropy(self): | |
| 2039 """Calculate the variance of the Dirichlet entropy. | |
| 2040 | |
| 2041 Ref: | |
| 2042 Wolpert & Wolf, PRE 53:6841-6854 (1996) Theorem 8 | |
| 2043 (Warning: this paper contains typos.) | |
| 2044 """ | |
| 2045 alpha = self.alpha | |
| 2046 A = float(sum(alpha)) | |
| 2047 A2 = A * (A+1) | |
| 2048 L = len(alpha) | |
| 2049 | |
| 2050 dg1 = zeros( (L) , float64) | |
| 2051 dg2 = zeros( (L) , float64) | |
| 2052 tg2 = zeros( (L) , float64) | |
| 2053 | |
| 2054 for i in range(L) : | |
| 2055 dg1[i] = digamma(alpha[i] + 1.0) | |
| 2056 dg2[i] = digamma(alpha[i] + 2.0) | |
| 2057 tg2[i] = trigamma(alpha[i] + 2.0) | |
| 2058 | |
| 2059 dg_Ap2 = digamma( A+2. ) | |
| 2060 tg_Ap2 = trigamma( A+2. ) | |
| 2061 | |
| 2062 mean = self.mean_entropy() | |
| 2063 var = 0.0 | |
| 2064 | |
| 2065 for i in range(L) : | |
| 2066 for j in range(L) : | |
| 2067 if i != j : | |
| 2068 var += ( | |
| 2069 ( dg1[i] - dg_Ap2 ) * (dg1[j] - dg_Ap2 ) - tg_Ap2 | |
| 2070 ) * (alpha[i] * alpha[j] ) / A2 | |
| 2071 else : | |
| 2072 var += ( | |
| 2073 ( dg2[i] - dg_Ap2 ) **2 + ( tg2[i] - tg_Ap2 ) | |
| 2074 ) * ( alpha[i] * (alpha[i]+1.) ) / A2 | |
| 2075 | |
| 2076 var -= mean**2 | |
| 2077 return var | |
| 2078 | |
| 2079 | |
| 2080 | |
| 2081 def mean_relative_entropy(self, pvec) : | |
| 2082 ln_p = na.log(pvec) | |
| 2083 return - self.mean_x(ln_p) - self.mean_entropy() | |
| 2084 | |
| 2085 | |
| 2086 def variance_relative_entropy(self, pvec) : | |
| 2087 ln_p = na.log(pvec) | |
| 2088 return self.variance_x(ln_p) + self.variance_entropy() | |
| 2089 | |
| 2090 | |
| 2091 def interval_relative_entropy(self, pvec, frac) : | |
| 2092 mean = self.mean_relative_entropy(pvec) | |
| 2093 variance = self.variance_relative_entropy(pvec) | |
| 2094 # If the variance is small, use the standard 95% | |
| 2095 # confidence interval: mean +/- 1.96 * sd | |
| 2096 if variance< 0.1 : | |
| 2097 sd = sqrt(variance) | |
| 2098 return max(0.0, mean - sd*1.96), mean + sd*1.96 | |
| 2099 sd = sqrt(variance) | |
| 2100 return max(0.0, mean - sd*1.96), mean + sd*1.96 | |
| 2101 | |
| 2102 g = gamma.from_mean_variance(mean, variance) | |
| 2103 low_limit = g.inverse_cdf( (1.-frac)/2.) | |
| 2104 high_limit = g.inverse_cdf( 1. - (1.-frac)/2. ) | |
| 2105 | |
| 2106 return low_limit, high_limit | |
| 2107 | |
| 2108 | |
| 2109 # Standard python voodoo for CLI | |
| 2110 if __name__ == "__main__": | |
| 2111 ## Code Profiling. Uncomment these lines | |
| 2112 #import hotshot, hotshot.stats | |
| 2113 #prof = hotshot.Profile("stones.prof") | |
| 2114 #prof.runcall(main) | |
| 2115 #prof.close() | |
| 2116 #stats = hotshot.stats.load("stones.prof") | |
| 2117 #stats.strip_dirs() | |
| 2118 #stats.sort_stats('cumulative', 'calls') | |
| 2119 #stats.print_stats(40) | |
| 2120 #sys.exit() | |
| 2121 | |
| 2122 main() | |
| 2123 | |
| 2124 | |
| 2125 | |
| 2126 | |
| 2127 | |
| 2128 |