rendered paste body#!/usr/bin/env pythonimport re, Bio.SeqIO, bisect, copydef split_stemloop(seq): cache = {} def loop_pen(n): return n < 4 and 1000 or n - 2 def pair_score(n1, n2): if '-' in n1+n2: return 6.0 elif 'A' in n1+n2 and 'T' in n1+n2: return -0.9 elif 'C' in n1+n2 and 'G' in n1+n2: return -2.3 elif 'G' in n1+n2 and 'T' in n1+n2: return 1.3 else: return 3.5 def hps(i, j): if j <= i: return 1000,'' try: return cache[(i,j)] except KeyError: nogap = hps(i+1, j-1) lgap = hps(i+1, j) rgap = hps(i, j-1) cache[(i,j)] = min((loop_pen(j - i + 1),','+seq[i:j+1]+','), (pair_score(seq[i], seq[j]) + nogap[0], seq[i]+nogap[1]+seq[j]), (pair_score('-', seq[j]) + lgap[0], '-'+lgap[1]+seq[j]), (pair_score(seq[i], '-') + rgap[0], seq[i]+rgap[1]+'-')) return cache[(i,j)] score, seqs = min([hps(i,j) for i in xrange(len(seq)-1) for j in xrange(len(seq)-1,i,-1)]) return seqs.split(',') + [score]def tail_score(tail): if len(tail) != 15: raise Exception('tail length != 15') sum = 0.0 prev = 1.0 for c in tail: if c == 'T': prev *= 0.9 else: prev *= 0.6 sum += prev return -sumdef deltaG(seq): dict = {'AA':-1,'AT':-.88,'TA':-.58,'CA':-1.45,'TG':-1.45,'GT':-1.44,'AC':-1.44, 'CT':-1.28,'AG':-1.28,'GA':-1.3,'TC':-1.3,'CG':-2.17,'GC':-2.24,'GG':-1.84, 'CC':-1.84,'TT':-1} dG = 0 #3.4 for i in xrange(len(seq)-1): dG += dict[seq[i:i+2]] return dGdef hybridization(tail): i = tail.find('T', 0, 2) if i == -1: spacer = tail[:2] i = 2 else: spacer = tail[:i] proximal = tail[i:i+5] distal = tail[i+5:i+9] extra = tail[i+9:i+12] return deltaG(spacer) + deltaG(proximal) + 0.5*deltaG(distal) + 0.01*deltaG(extra) - 3.4class Gene: def __init__(self, d): self.__dict__.update(d)def parsePTT(line): m = re.match(r"(?P<left>\d+)\.\."\ "(?P<right>\d+)\s+"\ "(?P<strand>-|\+)\s+"\ "(?P<length>\d+)\s+"\ "(?P<pid>\d+)\s+"\ "(?P<gene>\S+)\s+"\ "(?P<synonym>\S+)\s+"\ "(?P<code>\S+)\s+"\ "(?P<cog>\S+)\s+"\ "(?P<product>.+)", line) if m: g = m.groupdict() for name in ['right','left','length']: g[name] = int(g[name]) g['strand'] = g['strand'] == '-' and -1 or 1 return Gene(g) else: return Noneclass DefaultDict(dict): """Dictionary with a default value for unknown keys.""" def __init__(self, default): self.default = default def __getitem__(self, key): if key in self: return self.get(key) return self.setdefault(key, copy.deepcopy(self.default))def readPTT(filename): #indexes = DefaultDict(set()) forward = DefaultDict([]) backward = DefaultDict([]) for line in open(filename): g = parsePTT(line) if g: if g.strand == 1: dir = forward else: dir = backward for i in [g.left, g.right]: dir[i].append(g) #indexes[g.left].add(g) #indexes[g.right].add(g) #active = [] #last = set() #for i in sorted(indexes): # active.append((i, last.union(indexes[i]))) # last = last.symmetric_difference(indexes[i]) #return active print len(forward) print len(backward) return [sorted(dir.iteritems()) for dir in [forward, backward]]def gene_distance(genes, t): pos = t.left + (t.right - t.left)/2 if t.direction == 1: i = bisect.bisect(genes[0], (pos,0)) - 1 if i == -1: return 100000 (stop, gene) = [(g.right,g.gene) for g in genes[0][i][1]][0] return (pos-stop, gene) else: i = bisect.bisect(genes[1], (pos,0)) if i == len(genes[1]): return 1000000 (stop, gene) = [(g.left,g.gene) for g in genes[1][i][1]][0] return (stop-pos, gene)class Terminator: def __init__(self): self.direction = '' self.left = 0 self.right = 0 self.left_context = '' self.left_stem = '' self.loop = '' self.right_stem = '' self.right_context = '' self.gene_distance = 0 self.hp_energy = 0 self.tail_energy = 0 self.hybridization = 0 def importTTHP(self, defline): fields = ['stem_len', 'loop_len', 'gaps', 'hp_energy', 'tail_energy', 'conf', 'dir', 'left_stem_base', 'right_stem_base', 'left_stem_top', 'right_stem_top', 'string'] o = dict(zip(fields, eval(defline))) (self.left_context, self.left_stem, self.loop, self.right_stem, self.right_context) = o['string'].split() self.direction = {-1:'reverse', 1:'forward'}[o['dir']] self.hp_energy = o['hp_energy'] self.tail_energy = o['tail_energy'] self.left = o['left_stem_base'] self.right = o['right_stem_base'] return self def parseRegulon(self, line, seq, genes): m = re.match(r"(?P<id>\w+)\s+"\ "(?P<left>\d+)\s+"\ "(?P<right>\d+)\s+"\ "(?P<strand>forward|reverse)\s+"\ "(?P<sequence>\w+)\s+"\ "(?P<tu_name>\S+)", line) if m: print m.group('id'), m.group('left'), m.group('right'), print m.group('strand'), m.group('sequence'), m.group('tu_name') (left, right) = [int(m.group(x)) - 1 for x in ['left', 'right']] self.left, self.right = left, right #FIXME: possible out of seq access (self.left_context, stemloop, self.right_context) = [x.tostring() for x in [seq[left-15:left], seq[left:right+1], seq[right+1:right+1+15]]] #FIXME: check consistency #(l,stemloop,r) = re.match(r"([agct]+)([AGCT]+)([agct]+)", m.group('sequence')).groups() #self.left_context = l.upper() #self.right_context = r.upper() (self.left_stem, self.loop, self.right_stem, self.hp_energy) = split_stemloop(stemloop) self.direction = m.group('strand') == 'forward' and 1 or -1 self.tail_energy = tail_score(self.right_context) self.hybridization = hybridization(self.right_context) self.gene_distance = gene_distance(genes, self) return self def read_regulonDB(): seq = list(Bio.SeqIO.parse(open('NC_000913.fna'), 'fasta'))[0].seq genes = readPTT('NC_000913.ptt') return filter(None, [Terminator().parseRegulon(line, seq, genes) for line in open('TerminatorSet.txt')])terms = read_regulonDB()for t in terms: print t.__dict__