basic_sec_struct.py

class sec_struct:
    def __init__(self,ss):
        self.struct = ss
        self.length = len(self.struct)
        self.energy = 0.0
        self.gen_struct_info()
        self.accuracy = None
        self.native = None

    def gen_struct_info(self):
        self.compute_pairs()
        self.num_pairs = len(self.pairs)
        self.density = (2.0*self.num_pairs)/self.length

    def is_canonical( self,x, y ):
        # GU is the wobble pair.
        canonical = ['AU', 'CG', 'GC', 'GU', 'UA', 'UG']
        #canonical = ['AU', 'CG', 'GC', 'UA']
        return x.upper()+y.upper() in canonical

    def get_pairs(self):
        return self.pairs[:]

    def get_non_canonical_pairs(self, seq ):
        if self.length!=seq.get_length():
            raise ValueError("Secondary struct and sequence length mismatch.")
        canonical = lambda ij: not self.is_canonical(seq[ij[0]-1],seq[ij[1]-1])
        return filter(canonical, self.pairs)

    def get_canonical_pairs(self, seq ):
        if self.length!=seq.get_length():
            raise ValueError("Secondary struct and sequence length mismatch.")
        canonical = lambda ij: self.is_canonical(seq[ij[0]-1],seq[ij[1]-1])
        return filter(canonical, self.pairs)

    def canonical_struct(self, seq ):
        n = len(self.struct)
        canonical_struct = list('.'*n)

        for i,j in self.get_canonical_pairs( seq ):
            canonical_struct[i-1]='('
            canonical_struct[j-1]=')'
        return secondary_struct(''.join(canonical_struct))

    def get_num_pairs(self):
        return self.num_pairs

    def get_dot_bracket(self):
        return self.struct

    def get_energy(self):
        return self.energy

    def get_length(self):
        return self.length

    def get_num_canonical_pairs(self, seq):
        return len(self.get_canonical_pairs(seq))

    def get_density(self):
        return self.density

    def get_canonical_density(self,seq):
        return (2.0*self.get_num_canonical_pairs(seq))/self.length

    def compute_pairs(self):
        sec_struct = self.struct
        stack = []
        pairs = []
        canon_pairs = []
        for i,c in enumerate(sec_struct):
            if c == '(':
                stack.append(i)
            elif c == ')':
                j = stack.pop()
                pairs.append((j+1,i+1))
        self.pairs = sorted(pairs)


    def relative_accuracy(self, native):
        if self.accuracy is None or native!=self.native:
            self._compute_accuracy(native)
        return self.accuracy

    def _compute_accuracy(self,native,seq=None):
        # We treat self as the predicted structure and it is going to be
        # compared against native.

        predicted_pairs = set(self.get_pairs() if seq is None
                           else self.get_canonical_pairs(seq))
        native_pairs = set(native.get_pairs() if seq is None
                              else native.get_canonical_pairs(seq))

        self.num_canon_pairs = len(predicted_pairs)
        self.tp = len(native_pairs.intersection(predicted_pairs))
        self.fn = len(native_pairs.difference(predicted_pairs))
        self.fp = len(predicted_pairs.difference(native_pairs))

        denom = 2.0*self.tp+self.fn+self.fp
        self.accuracy = 2*self.tp/denom

        return self.accuracy

    def __eq__(self,other):
        if isinstance(other, secondary_struct):
            return self.struct==other.struct
        return NotImplemented

    def __ne__(self,other):
        result = self.__eq__(other)
        if result is NotImplemented:
            return result
        return not result

    def __getitem__(self,i):
        return self.struct[i]

    def __str__(self):
        return self.struct