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
