Mercurial > repos > xuebing > sharplabtool

comparison mytools/sequence.py @ 7:f0dc65e7f6c0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Fri, 09 Mar 2012 19:59:07 -0500
parents
children
comparison
equal deleted inserted replaced
6:094f377db4e5 7:f0dc65e7f6c0
1 #!@WHICHPYTHON@
2
3 import copy, string, sys
4
5 #------------------ Alphabet -------------------
6
7 class Alphabet(object):
8 """Biological alphabet class.
9 This defines the set of symbols from which various objects can be built, e.g. sequences and motifs.
10 The symbol set is immutable and accessed as a tuple.
11 symstr: symbols in alphabet as either a tuple or string
12 complement: dictionary defining letters and their complement
13 """
14 def __init__(self, symstr, complement = None):
15 """Construct an alphabet from a string or tuple of characters.
16 Lower case characters will be converted to upper case.
17 An optional mapping for complements may be provided.
18 Example:
19 >>> alpha = sequence.Alphabet('ACGTttga', {'A':'C', 'G':'T'})
20 >>> alpha.getSymbols()
21 will construct the DNA alphabet and output:
22 ('A', 'C', 'G', 'T')
23 """
24 symlst = []
25 for s in [str(sym).upper()[0] for sym in symstr]:
26 if not s in symlst:
27 symlst.append(s)
28 self.symbols = tuple(symlst)
29 if complement != None:
30 # expand the mapping and check for contradictions
31 cmap = {}
32 for s in self.symbols:
33 c = complement.get(s, None)
34 if c != None:
35 if s in cmap and cmap[s] != c:
36 raise RuntimeError("Alphabet complement map "
37 "contains contradictory mapping")
38 cmap[s] = c
39 cmap[c] = s
40 # replace mapping with indicies
41 cimap = {}
42 for idx in range (len(self.symbols)):
43 s = self.symbols[idx]
44 if s in cmap:
45 cimap[cmap[s]] = idx
46 # create tuple
47 cidxlst = []
48 for idx in range (len(self.symbols)):
49 cidxlst.append(cimap.get(self.symbols[idx], None))
50 self.complements = tuple(cidxlst)
51 else:
52 self.complements = None
53
54 def getSymbols(self):
55 """Retrieve a tuple with all symbols, immutable membership and order"""
56 return self.symbols
57
58 def getComplements(self):
59 """Retrieve a tuple with all complement indicies, immutable"""
60 return self.complements
61
62 def isValidSymbol(self, sym):
63 """Check if the symbol is a member of alphabet"""
64 return any([s==sym for s in self.symbols])
65
66 def getIndex(self, sym):
67 """Retrieve the index of the symbol (immutable)"""
68 for idx in range (len(self.symbols)):
69 if self.symbols[idx] == sym:
70 return idx
71 raise RuntimeError("Symbol " + sym + " does not exist in alphabet")
72
73 def isComplementable(self):
74 return self.complements != None
75
76 def getComplement(self, sym):
77 """Retrieve the complement of the symbol (immutable)"""
78 return self.symbols[self.complements[self.getIndex(sym)]];
79
80 def isValidString(self, symstr):
81 """Check if the string contains only symbols that belong to the alphabet"""
82 found = True
83 for sym in symstr:
84 if self.isValidSymbol(sym) == False:
85 return False
86 return True
87
88 def getLen(self):
89 """Retrieve the number of symbols in (the length of) the alphabet"""
90 return len(self.symbols)
91
92 # pre-defined alphabets that can be specified by their name
93 predefAlphabets = [
94 ("DNA" , Alphabet('ACGT', {'A':'T', 'G':'C'})),
95 ("RNA" , Alphabet('ACGU')),
96 ("Extended DNA" , Alphabet('ACGTYRN')),
97 ("Protein" , Alphabet('ACDEFGHIKLMNPQRSTVWY')),
98 ("Extended Protein" , Alphabet('ACDEFGHIKLMNPQRSTVWYX')),
99 ("TM Labels" , Alphabet('MIO'))
100 ]
101
102 def getAlphabet(name):
103 """Retrieve a pre-defined alphabet by name.
104 Currently, "Protein", "DNA", "RNA", "Extended DNA", "Extended Protein" and "TM Labels" are available.
105 Example:
106 >>> alpha = sequence.getAlphabet('Protein')
107 >>> alpha.getSymbols()
108 will retrieve the 20 amino acid alphabet and output the tuple:
109 ('A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y')
110 """
111 for (xname, xalpha) in predefAlphabets:
112 if xname == name:
113 return xalpha
114 return None
115
116 #------------------ Sequence -------------------
117
118 class Sequence(object):
119 """Biological sequence class. Sequence data is immutable.
120
121 data: the sequence data as a tuple or string
122 alpha: the alphabet from which symbols are taken
123 name: the sequence name, if any
124 info: can contain additional sequence information apart from the name
125 """
126 def __init__(self, sequence, alpha = None, name = "", seqinfo = ""):
127 """Create a sequence with sequence data.
128 Specifying the alphabet is optional, so is the name and info.
129 Example:
130 >>> myseq = sequence.Sequence('MVSAKKVPAIAMSFGVSF')
131 will create a sequence with name "", and assign one of the predefined alphabets on basis of what symbols were used.
132 >>> myseq.getAlphabet().getSymbols()
133 will most likely output the standard protein alphabet:
134 ('A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y')
135 """
136 if type(sequence) is str:
137 self.data = tuple(sequence.upper())
138 elif type(sequence) is tuple:
139 self.data = sequence
140 elif type(sequence) is list:
141 self.data = tuple([s.upper() for s in sequence])
142 else:
143 raise RuntimeError("Sequence data is not specified correctly: must be string or tuple")
144 # Resolve choice of alphabet
145 validAlphabet = False
146 if alpha == None: # Alphabet is not set, attempt to set it automatically...
147 for (xname, xalpha) in predefAlphabets: # Iterate through each predefined alphabet, in order
148 if xalpha.isValidString( self.data ): # This alphabet works, go with it
149 self.alpha = alpha = xalpha
150 validAlphabet = True
151 break
152 self.name = name
153 self.info = seqinfo
154 if validAlphabet == False: # we were either unsuccessful above or the alphabet was specified so test it
155 if type(alpha) is str: # check if name is a predefined alphabet
156 for (xname, xalpha) in predefAlphabets: # Iterate through each predefined alphabet, check for name
157 if (xname == alpha):
158 alpha = xalpha
159 break
160 if type(alpha) is Alphabet: # the alphabet is specified
161 if alpha.isValidString(self.data) == False:
162 raise RuntimeError("Invalid alphabet specified: "+"".join(alpha.getSymbols())+" is not compatible with sequence '"+"".join(self.data)+"'")
163 else:
164 self.alpha = alpha
165 else:
166 raise RuntimeError("Could not identify alphabet from sequence")
167
168 #basic getters and setters for the class
169 def getName(self):
170 """Get the name of the sequence"""
171 return self.name
172 def getInfo(self):
173 """Get additional info of the sequence (e.g. from the defline in a FASTA file)"""
174 return self.info
175 def getAlphabet(self):
176 """Retrieve the alphabet that is assigned to this sequence"""
177 return self.alpha
178 def setName(self, name):
179 """Change the name of the sequence"""
180 self.name = name
181 def setAlphabet(self, alpha):
182 """Set the alphabet, throws an exception if it is not compatible with the sequence data"""
183 if type(alpha) is Alphabet:
184 if alpha.isValid( sequence ) == False:
185 raise RuntimeError( "Invalid alphabet specified" )
186 #sequence functions
187 def getSequence(self):
188 """Retrieve the sequence data (a tuple of symbols)"""
189 return self.data
190 def getString(self):
191 """Retrieve the sequence data as a string (copy of actual data)"""
192 return "".join(self.data)
193 def getLen(self):
194 """Get the length of the sequence (number of symbols)"""
195 return len(self.data)
196 def getSite(self, position, length = 1):
197 """Retrieve a site in the sequence of desired length.
198 Note that positions go from 0 to length-1, and that if the requested site
199 extends beyond those the method throws an exception.
200 """
201 if position >= 0 and position <= self.getLen() - length:
202 if length == 1:
203 return self.data[position]
204 else:
205 return self.data[position:position+length]
206 else:
207 raise RuntimeError( "Attempt to access invalid position in sequence "+self.getName() )
208
209 def nice(self):
210 """ A short description of the sequence """
211 print self.getName(), ":", self.getLen()
212
213 def readStrings(filename):
214 """ Read one or more lines of text from a file--for example an alignment.
215 Return as a list of strings.
216 filename: name of file
217 """
218 txtlist = []
219 f = open(filename)
220 for line in f.readlines():
221 txtlist.extend(line.split())
222 return txtlist
223
224 def readFASTA(filename, alpha = None):
225 """ Read one or more sequences from a file in FASTA format.
226 filename: name of file to load sequences from
227 alpha: alphabet that is used (if left unspecified, an attempt is made to identify the alphabet for each individual sequence)
228 """
229 seqlist = []
230 seqname = None
231 seqinfo = None
232 seqdata = []
233 fh = open(filename)
234 thisline = fh.readline()
235 while (thisline):
236 if (thisline[0] == '>'): # new sequence
237 if (seqname): # take care of the data that is already in the buffer before processing the new sequence
238 try:
239 seqnew = Sequence(seqdata, alpha, seqname, seqinfo)
240 seqlist.append(seqnew)
241 except RuntimeError, e:
242 print >> sys.stderr, "Warning: "+seqname+" is invalid (ignored): ", e
243 seqinfo = thisline[1:-1] # everything on the defline is "info"
244 seqname = seqinfo.split()[0] # up to first space
245 seqdata = []
246 else: # pull out the sequence data
247 cleanline = thisline.split()
248 for line in cleanline:
249 seqdata.extend(tuple(line.strip('*'))) # sometimes a line ends with an asterisk in FASTA files
250 thisline = fh.readline()
251
252 if (seqname):
253 try:
254 seqnew = Sequence(seqdata, alpha, seqname, seqinfo)
255 seqlist.append(seqnew)
256 except RuntimeError, e:
257 print >> sys.stderr, "Warning: " + seqname + " is invalid (ignored): ", e
258 else:
259 raise RuntimeError("No sequences on FASTA format found in this file")
260 fh.close()
261 return seqlist
262
263 def _writeOneFASTA(sequence, filehandle):
264 """Write one sequence in FASTA format to an already open file"""
265 filehandle.write(">" + sequence.getName()+"\n")
266 data = sequence.getSequence()
267 lines = ( sequence.getLen() - 1) / 60 + 1
268 for i in range(lines):
269 #note: python lets us get the last line (var length) free
270 #lineofseq = data[i*60 : (i+1)*60] + "\n"
271 lineofseq = "".join(data[i*60 : (i+1)*60]) + "\n"
272 filehandle.write(lineofseq)
273
274 def writeFASTA(sequence, filename):
275 """Write a list (or a single) of sequences to a file in the FASTA format"""
276 fh = open(filename, "w")
277 if isinstance(sequence, Sequence):
278 _writeOneFASTA(sequence, fh)
279 else:
280 for seq in sequence:
281 if isinstance(seq, Sequence):
282 _writeOneFASTA(seq, fh)
283 else:
284 print >> sys.stderr, "Warning: could not write " + seq.getName() + " (ignored)."
285 fh.flush()
286 fh.close()
287
288 #------------------ Distrib -------------------
289
290 class Distrib(object):
291 """Class for storing a multinomial probability distribution over the symbols in an alphabet"""
292 def __init__(self, alpha, pseudo_count = 0.0):
293 self.alpha = alpha
294 self.tot = pseudo_count * self.alpha.getLen()
295 self.cnt = [pseudo_count for _ in range( self.alpha.getLen() )]
296
297 def __deepcopy__(self, memo):
298 dup = Distrib(self.alpha)
299 dup.tot = copy.deepcopy(self.tot, memo)
300 dup.cnt = copy.deepcopy(self.cnt, memo)
301 return dup
302
303 def count(self, syms = None ):
304 """Count an observation of a symbol"""
305 if syms == None:
306 syms = self.alpha.getSymbols()
307 for sym in syms:
308 idx = self.alpha.getIndex( sym )
309 self.cnt[idx] += 1.0
310 self.tot += 1
311
312 def complement(self):
313 """Complement the counts, throw an error if this is impossible"""
314 if not self.alpha.isComplementable():
315 raise RuntimeError("Attempt to complement a Distrib "
316 "based on a non-complementable alphabet.")
317 coms = self.alpha.getComplements()
318 new_count = []
319 for idx in range(len(coms)):
320 cidx = coms[idx]
321 if cidx == None:
322 cidx = idx
323 new_count.append(self.cnt[cidx])
324 self.cnt = new_count
325 return self
326
327 def reset(self):
328 """Reset the distribution, that is, restart counting."""
329 self.tot = 0
330 self.cnt = [0.0 for _ in range( self.alpha.getLen() )]
331
332 def getFreq(self, sym = None):
333 """Determine the probability distribution from the current counts.
334 The order in which probabilities are given follow the order of the symbols in the alphabet."""
335 if self.tot > 0:
336 if sym == None:
337 freq = tuple([ y / self.tot for y in self.cnt ])
338 return freq
339 else:
340 idx = self.alpha.getIndex( sym )
341 return self.cnt[idx] / self.tot
342 return None
343
344 def pretty(self):
345 """Retrieve the probabilites for all symbols and return as a pretty table (a list of text strings)"""
346 table = ["".join(["%4s " % s for s in self.alpha.getSymbols()])]
347 table.append("".join(["%3.2f " % y for y in Distrib.getFreq(self)]))
348 return table
349
350 def getSymbols(self):
351 """Get the symbols in the alphabet in the same order as probabilities are given."""
352 return self.alpha.getSymbols()
353
354 def getAlphabet(self):
355 """Get the alphabet over which the distribution is defined."""
356 return self.alpha
357
358 #------------------ Motif (and subclasses) -------------------
359
360 class Motif(object):
361 """ Sequence motif class--defining a pattern that can be searched in sequences.
362 This class is not intended for direct use. Instead use and develop sub-classes (see below).
363 """
364 def __init__(self, alpha):
365 self.len = 0
366 self.alpha = alpha
367
368 def getLen(self):
369 """Get the length of the motif"""
370 return self.len
371
372 def getAlphabet(self):
373 """Get the alphabet that is used in the motif"""
374 return self.alpha
375
376 def isAlphabet(self, seqstr):
377 """Check if the sequence can be processed by this motif"""
378 mystr = seqstr
379 if type(seqstr) is Sequence:
380 mystr = seqstr.getString()
381 return self.getAlphabet().isValidString(mystr)
382
383 import re
384
385 class RegExp(Motif):
386 """A motif class that defines the pattern in terms of a regular expression"""
387 def __init__(self, alpha, re_string):
388 Motif.__init__(self, alpha)
389 self.pattern = re.compile(re_string)
390
391 def match(self, seq):
392 """Find matches to the motif in a specified sequence.
393 The method is a generator, hence subsequent hits can be retrieved using next().
394 The returned result is a tuple (position, match-sequence, score), where score is
395 always 1.0 since a regular expression is either true or false (not returned).
396 """
397 myseq = seq
398 if not type(seq) is Sequence:
399 myseq = Sequence(seq, self.alpha)
400 mystr = myseq.getString()
401 if not Motif.isAlphabet(self, mystr):
402 raise RuntimeError("Motif alphabet is not valid for sequence " + myseq.getName())
403 for m in re.finditer(self.pattern, mystr):
404 yield (m.start(), m.group(), 1.0)
405
406 import math, time
407
408 # Variables used by the PWM for creating an EPS file
409 _colour_def = (
410 "/black [0 0 0] def\n"
411 "/red [0.8 0 0] def\n"
412 "/green [0 0.5 0] def\n"
413 "/blue [0 0 0.8] def\n"
414 "/yellow [1 1 0] def\n"
415 "/purple [0.8 0 0.8] def\n"
416 "/magenta [1.0 0 1.0] def\n"
417 "/cyan [0 1.0 1.0] def\n"
418 "/pink [1.0 0.8 0.8] def\n"
419 "/turquoise [0.2 0.9 0.8] def\n"
420 "/orange [1 0.7 0] def\n"
421 "/lightred [0.8 0.56 0.56] def\n"
422 "/lightgreen [0.35 0.5 0.35] def\n"
423 "/lightblue [0.56 0.56 0.8] def\n"
424 "/lightyellow [1 1 0.71] def\n"
425 "/lightpurple [0.8 0.56 0.8] def\n"
426 "/lightmagenta [1.0 0.7 1.0] def\n"
427 "/lightcyan [0.7 1.0 1.0] def\n"
428 "/lightpink [1.0 0.9 0.9] def\n"
429 "/lightturquoise [0.81 0.9 0.89] def\n"
430 "/lightorange [1 0.91 0.7] def\n")
431 _colour_dict = (
432 "/fullColourDict <<\n"
433 " (G) orange\n"
434 " (T) green\n"
435 " (C) blue\n"
436 " (A) red\n"
437 " (U) green\n"
438 ">> def\n"
439 "/mutedColourDict <<\n"
440 " (G) lightorange\n"
441 " (T) lightgreen\n"
442 " (C) lightblue\n"
443 " (A) lightred\n"
444 " (U) lightgreen\n"
445 ">> def\n"
446 "/colorDict fullColourDict def\n")
447
448 _eps_defaults = {
449 'LOGOTYPE': 'NA',
450 'FONTSIZE': '12',
451 'TITLEFONTSIZE': '12',
452 'SMALLFONTSIZE': '6',
453 'TOPMARGIN': '0.9',
454 'BOTTOMMARGIN': '0.9',
455 'YAXIS': 'true',
456 'YAXISLABEL': 'bits',
457 'XAXISLABEL': '',
458 'TITLE': '',
459 'ERRORBARFRACTION': '1.0',
460 'SHOWINGBOX': 'false',
461 'BARBITS': '2.0',
462 'TICBITS': '1',
463 'COLORDEF': _colour_def,
464 'COLORDICT': _colour_dict,
465 'SHOWENDS': 'false',
466 'NUMBERING': 'true',
467 'OUTLINE': 'false',
468 }
469 class PWM(Motif):
470 """This motif subclass defines a pattern in terms of a position weight matrix.
471 An alphabet must be provided. A pseudo-count to be added to each count is
472 optional. A uniform background distribution is used by default.
473 """
474 def __init__(self, alpha):
475 Motif.__init__(self, alpha) # set alphabet of this multinomial distribution
476 self.background = Distrib(alpha) # the default background ...
477 self.background.count(alpha.getSymbols()) # ... is uniform
478 self.nsites = 0
479
480 def setFromAlignment(self, aligned, pseudo_count = 0.0):
481 """Set the probabilities in the PWM from an alignment.
482 The alignment is a list of equal-length strings (see readStrings), OR
483 a list of Sequence.
484 """
485 self.cols = -1
486 self.nsites = len(aligned)
487 seqs = []
488 # Below we create a list of Sequence from the alignment,
489 # while doing some error checking, and figure out the number of columns
490 for s in aligned:
491 # probably a text string, so we make a nameless sequence from it
492 if not type(s) is Sequence:
493 s=Sequence(s, Motif.getAlphabet(self))
494 else:
495 # it was a sequence, so we check that the alphabet in
496 # this motif will be able to process it
497 if not Motif.isAlphabet(self, s):
498 raise RuntimeError("Motif alphabet is not valid for sequence " + s.getName())
499 if self.cols == -1:
500 self.cols = s.getLen()
501 elif self.cols != s.getLen():
502 raise RuntimeError("Sequences in alignment are not of equal length")
503 seqs.append(s)
504 # The line below initializes the list of Distrib (one for each column of the alignment)
505 self.counts = [Distrib(Motif.getAlphabet(self), pseudo_count) for _ in range(self.cols)]
506 # Next, we do the counting, column by column
507 for c in range( self.cols ): # iterate through columns
508 for s in seqs: # iterate through rows
509 # determine the index of the symbol we find at this position (row, column c)
510 self.counts[c].count(s.getSite(c))
511 # Update the length
512 self.len = self.cols
513
514 def reverseComplement(self):
515 """Reverse complement the PWM"""
516 i = 0
517 j = len(self.counts)-1
518 while (i < j):
519 temp = self.counts[i];
520 self.counts[i] = self.counts[j]
521 self.counts[j] = temp
522 self.counts[i].complement()
523 self.counts[j].complement()
524 i += 1;
525 j -= 1;
526 if i == j:
527 self.counts[i].complement()
528 return self
529
530 def getNSites(self):
531 """Get the number of sites that made the PWM"""
532 return self.nsites
533
534 def setBackground(self, distrib):
535 """Set the background distribution"""
536 if not distrib.getAlphabet() == Motif.getAlphabet(self):
537 raise RuntimeError("Incompatible alphabets")
538 self.background = distrib
539
540 def getFreq(self, col = None, sym = None):
541 """Get the probabilities for all positions in the PWM (a list of Distribs)"""
542 if (col == None):
543 return [y.getFreq() for y in self.counts]
544 else:
545 return self.counts[col].getFreq(sym)
546
547 def pretty(self):
548 """Retrieve the probabilites for all positions in the PWM as a pretty table (a list of text strings)"""
549 #table = ["".join(["%8s " % s for s in self.alpha.getSymbols()])]
550 table = []
551 for row in PWM.getFreq(self):
552 table.append("".join(["%8.6f " % y for y in row]))
553 return table
554
555 def logoddsPretty(self, bkg):
556 """Retrieve the (base-2) log-odds for all positions in the PWM as a pretty table (a list of text strings)"""
557 table = []
558 for row in PWM.getFreq(self):
559 #table.append("".join(["%8.6f " % (math.log((row[i]+1e-6)/bkg[i])/math.log(2)) for i in range(len(row))]))
560 table.append("".join(["%8.6f " % (math.log((row[i])/bkg[i])/math.log(2)) for i in range(len(row))]))
561 #table.append("".join(["%8.6f " % row[i] for i in range(len(row))]))
562 return table
563
564
565 def consensus_sequence(self):
566 """
567 Get the consensus sequence corresponding to a PWM.
568 Consensus sequence is the letter in each column
569 with the highest probability.
570 """
571 consensus = ""
572 alphabet = Motif.getAlphabet(self).getSymbols()
573 for pos in range(self.cols):
574 best_letter = alphabet[0]
575 best_p = self.counts[pos].getFreq(best_letter)
576 for letter in alphabet[1:]:
577 p = self.counts[pos].getFreq(letter)
578 if p > best_p:
579 best_p = p
580 best_letter = letter
581 consensus += best_letter
582 return consensus
583
584
585 def consensus(self):
586 """
587 Get the consensus corresponding to a PWM.
588 Consensus at each column of motif is a list of
589 characters with non-zero probabilities.
590 """
591 consensus = []
592 for pos in range(self.cols):
593 matches = []
594 for letter in Motif.getAlphabet(self).getSymbols():
595 p = self.counts[pos].getFreq(letter)
596 if p > 0:
597 matches += letter
598 consensus.append(matches)
599 return consensus
600
601
602 def getScore(self, seq, start):
603 """Score this particular list of symbols using the PFM (background needs to be set separately)"""
604 sum = 0.0
605 seqdata = seq.getSequence()[start : start+self.cols]
606 for pos in range(len(seqdata)):
607 q = self.counts[pos].getFreq(seqdata[pos])
608 if q == 0:
609 q = 0.0001 # to avoid log(0) == -Infinity
610 logodds = math.log(q / self.background.getFreq(seqdata[pos]))
611 sum += logodds
612 return sum
613
614 def match(self, seq, _LOG0 = -10):
615 """Find matches to the motif in a specified sequence.
616 The method is a generator, hence subsequent hits can be retrieved using next().
617 The returned result is a tuple (position, match-sequence, score).
618 The optional parameter _LOG0 specifies a lower bound on reported logodds scores.
619 """
620 myseq = seq
621 if not type(seq) is Sequence:
622 myseq = Sequence(seq, self.alpha)
623 if not Motif.isAlphabet(self, myseq):
624 raise RuntimeError("Motif alphabet is not valid for sequence " + myseq.getName())
625 for pos in range(myseq.getLen() - self.cols):
626 score = PWM.getScore(self, myseq, pos)
627 if score > _LOG0:
628 yield (pos, "".join(myseq.getSite(pos, self.cols)), score)
629
630 def writeEPS(self, program, template_file, eps_fh,
631 timestamp = time.localtime()):
632 """Write out a DNA motif to EPS format."""
633 small_dfmt = "%d.%m.%Y %H:%M"
634 full_dfmt = "%d.%m.%Y %H:%M:%S %Z"
635 small_date = time.strftime(small_dfmt, timestamp)
636 full_date = time.strftime(full_dfmt, timestamp)
637 points_per_cm = 72.0 / 2.54
638 height = 4.5
639 width = self.getLen() * 0.8 + 2
640 width = min(30, width)
641 points_height = int(height * points_per_cm)
642 points_width = int(width * points_per_cm)
643 defaults = _eps_defaults.copy()
644 defaults['CREATOR'] = program
645 defaults['CREATIONDATE'] = full_date
646 defaults['LOGOHEIGHT'] = str(height)
647 defaults['LOGOWIDTH'] = str(width)
648 defaults['FINEPRINT'] = program + ' ' + small_date
649 defaults['CHARSPERLINE'] = str(self.getLen())
650 defaults['BOUNDINGHEIGHT'] = str(points_height)
651 defaults['BOUNDINGWIDTH'] = str(points_width)
652 defaults['LOGOLINEHEIGHT'] = str(height)
653 with open(template_file, 'r') as template_fh:
654 m_var = re.compile("\{\$([A-Z]+)\}")
655 for line in template_fh:
656 last = 0
657 match = m_var.search(line)
658 while (match):
659 if (last < match.start()):
660 prev = line[last:match.start()]
661 eps_fh.write(prev)
662 key = match.group(1)
663 if (key == "DATA"):
664 eps_fh.write("\nStartLine\n")
665 for pos in range(self.getLen()):
666 eps_fh.write("({0:d}) startstack\n".format(pos+1))
667 stack = []
668 # calculate the stack information content
669 alpha_ic = 2
670 h = 0
671 for sym in self.getAlphabet().getSymbols():
672 freq = self.getFreq(pos, sym)
673 if (freq == 0):
674 continue
675 h -= (freq * math.log(freq, 2))
676 stack_ic = alpha_ic - h
677 # calculate the heights of each symbol
678 for sym in self.getAlphabet().getSymbols():
679 freq = self.getFreq(pos, sym)
680 if (freq == 0):
681 continue
682 stack.append((freq * stack_ic, sym))
683 stack.sort();
684 # output the symbols
685 for symh, sym in stack:
686 eps_fh.write(" {0:f} ({1:s}) numchar\n".format(
687 symh, sym))
688 eps_fh.write("endstack\n\n")
689 eps_fh.write("EndLine\n")
690 elif (key in defaults):
691 eps_fh.write(defaults[key])
692 else:
693 raise RuntimeError('Unknown variable "' + key +
694 '" in EPS template')
695 last = match.end();
696 match = m_var.search(line, last)
697 if (last < len(line)):
698 eps_fh.write(line[last:])
699
700
701 #------------------ Main method -------------------
702 # Executed if you run this file from the operating system prompt, e.g.
703 # > python sequence.py
704
705 if __name__=='__main__':
706 alpha = getAlphabet('Extended DNA')
707 #seqs = readFASTA('pos.fasta')
708 seqs = []
709 aln = readStrings('tmp0')
710 #regexp = RegExp(alpha, '[AG]G.[DE]TT[AS].')
711 pwm = PWM(alpha)
712 pwm.setFromAlignment(aln)
713 for row in pwm.pretty():
714 print row
715 for s in seqs:
716 print s.getName(), s.getLen(), s.getAlphabet().getSymbols()
717 for m in regexp.match( s ):
718 print "pos: %d pat: %s %4.2f" % (m[0], m[1], m[2])
719 for m in pwm.match( s ):
720 print "pos: %d pat: %s %4.2f" % (m[0], m[1], m[2])