Python RNA类代码示例

def main():
    """ for sequence string, calculate mfe structure, mfe, pf, base pair probability matrix, plot structures and bppms, calculate accessibilities"""
    print 'name mfe_low mfe_high pf_low pf_high RRS_acces_low RRS_acces_high AUG_acces_low AUG_acces_high'
    for seq_file in SeqIO.parse(sys.stdin, 'fasta'):
         sequ = str(seq_file.seq)
         fc_low = RNA.fold_compound(sequ, MODEL_LOW_TEMPERATURE)
         fc_high = RNA.fold_compound(sequ, MODEL_HIGH_TEMPERATURE)
         struct_low, mfe_low = fc_low.mfe()
         struct_high, mfe_high = fc_high.mfe()
         pfstruct_low, pf_low = fc_low.pf()
         pfstruct_high, pf_high = fc_high.pf()
         bppm_low = fc_low.bpp()
         bppm_high = fc_high.bpp()
         
         plot_2bppms(bppm_low, bppm_high, seq_file.id)
         RNA.PS_rna_plot(sequ, struct_low, '{:s}_low_ss.ps'.format(str2filename(seq_file.id)))    
         RNA.PS_rna_plot(sequ, struct_high, '{:s}_high_ss.ps'.format(str2filename(seq_file.id)))
    
         constr1, constr2 = seqconstraints(sequ,RRS,START,SPACER)
         RRS_acces_low = accessibility(sequ,MODEL_LOW_TEMPERATURE,constr1,pf_low)
         RRS_acces_high = accessibility(sequ,MODEL_HIGH_TEMPERATURE,constr1,pf_high)
         AUG_acces_low = accessibility(sequ,MODEL_LOW_TEMPERATURE,constr2,pf_low)
         AUG_acces_high = accessibility(sequ,MODEL_HIGH_TEMPERATURE,constr2,pf_high)

         print seq_file.id, mfe_low, mfe_high, pf_low, pf_high, RRS_acces_low, RRS_acces_high, AUG_acces_low, AUG_acces_high
    print versions_used()

def local_search(start_seq_, target_structs_, seq_constraint_,
                 context_front=None, context_back=None):
    global start_seq
    global seq_constraint
    global target_structs

    rna.check_struct_seq_match(target_structs_[0], start_seq_)
    rna.check_struct_seq_match(target_structs_[1], start_seq_)

    start_seq = start_seq_
    target_structs = target_structs_
    seq_constraint = seq_constraint_

    # # TODO: has to be checked
    # preset_dangles = RNA.dangles
    # if preset_dangles != 0:
    #     RNA.dangles = 1

    if (SEARCH_STRATEGY == SearchStrategy.adaptive_walk or
        SEARCH_STRATEGY == SearchStrategy.stochastic_local_search):
        seq, cost, steps = local_search_sls_pf()
    elif SEARCH_STRATEGY == SearchStrategy.full_local_search:
        seq, cost, steps = local_search_fls_pf()
    else:
        raise ValueError("Specified search strategy not valid.")

    eval_seq_container.reset()
    vienna_rna.free_pf_arrays()
    vienna_rna.free_arrays()

    # RNA.dangles = preset_dangles
    return seq, cost, steps

def getBPPM(sequence, structure = "", bppm_cutoff = 0.00001):
    """
        Requires ViennaRNAtools Python module
        Returns the base pair probability matrix using Vienna pf_fold, get_pr and free_pf_arrays functions.
        returns upper triangular matrix, whose entries exceed a threshold
    """
    bppm = {}
    
    
     #'--noPS', '-d 2', t, P
     
     
     
    if structure != "":
        RNA.cvar.fold_constrained = 1
    else:
        RNA.cvar.fold_constrained = 0
    #print "Before", structure
    RNA.pf_fold(sequence, structure)
    #print "After", structure
    seq_len = len(sequence)+1
    for i in xrange(1, seq_len):
        for j in xrange(1, seq_len):
            if i<j:
                bpp = RNA.get_pr(i,j)
                if bpp > bppm_cutoff:
                    bppm[str(i) + "_" + str(j)] = bpp
                else:
                    bppm[str(i) + "_" + str(j)] = 0
    RNA.free_pf_arrays()
    #print bppm
    #exit(1)
    return bppm

def mfe_bp_distance(S, G, masked=None):
    """This function takes an RNA sequence S, a secondary structure G,
    and returns de base pairs distance between the mfe structure (a mask
    for the folding can be provided as an optional argument) of S
    and G
    """
    Sec_struct = RNA.fold(S)[0]
    return RNA.bp_distance(Sec_struct, G)

def prep_sec2(seq_five, seq_three, seq_apta, shift, rand):
    (take, dump) = RNA.fold(seq_five + "N" + "G" * 100 + "N" * 4 + "C" * 100 + "N" + seq_three)
    (take1, dump) = RNA.fold(seq_apta)
    seq = (
        len(take.split("." + "(" * 100 + "." * 4 + ")" * 100 + ".", 1)[0]) * "."
        + "(" * rand
        + take1
        + ")" * (rand)
        + "." * shift
        + len(take.split("." + "(" * 100 + "." * 4 + ")" * 100 + ".", 1)[1]) * "."
    )
    return seq

def fold_probability(S, G=None):
    """Given a sequence S a secondary structure G (default mfe), we compute 
    the partition function of S given G as a constraint. The output
    is a triple (A,B,C) where A is the annotated partition folding,
    B is the energie of the ensemble A, and C a dictionary having as keys 
    a pair of positions and as value the probability of having the pair.
    """
    struct, energy = RNA.pf_fold(S, G) #Compute the partition function
    dict_probabilities = {}
    for left, right in ((x,y) for x in range(len(S)) for y in range(len(S))
                       if x < y):
        dict_probabilities[left,right] =RNA.get_pr(left + 1,right +1)
    return (struct, energy, dict_probabilities)

def temperature_reactivity( sequence, structure, temperature1, temperature2 ):
    """Evaluate temperature-dependent difference in energy, entropy, enthalpy."""
    temperature_model1 = RNA.md()
    temperature_model2 = RNA.md()
    temperature_model1.temperature = temperature1
    temperature_model2.temperature = temperature2
    fc1 = RNA.fold_compound(sequence, temperature_model1)
    fc2 = RNA.fold_compound(sequence, temperature_model2)
    energy_of_struct1 = fc1.eval_structure(structure)
    energy_of_struct2 = fc2.eval_structure(structure)
    # normalize delta_energy, add 0.001 to prevent division by 0
    delta_energy = abs((energy_of_struct2 - energy_of_struct1) / (energy_of_struct2 + 0.001))

    return (delta_energy, energy_of_struct1, energy_of_struct2)

    def Rewards(self,k):
        #copy_unpairedposition=list(unpairedposition)
        #copy_bppused=list(bppused)
        if k > len(str_uindex)-1:
            posbasep=self.position[len(str_uindex):self.n]
            posbase=self.position[0:len(str_uindex)]
            e=list(itertools.chain(*posbasep))
            for i in range(len(a)):
                posbase.insert(b[i],e[c[i]])
            mutated_s= ''.join(map(str, posbase))
            mutated_str1=RNA.fold(mutated_s)
            mutated_str=mutated_str1[0]

            d=0.0
            g=0.0
            n=len(s)
            for i in range(len(s)):
                if mutated_str[i]!=s[i]:
                    d=d+1
            g=(n-d)/n
            if g==1.0:
                solution.append(mutated_s)
                return g
            else:
                return g




        if k <= len(str_uindex)-1:
            posbasep=self.position[len(str_uindex):self.n]
            posbase=self.position[0:len(str_uindex)]
            e=list(itertools.chain(*posbasep))
            for i in range(len(a)):
                posbase.insert(b[i],e[c[i]])
            mutated_s= ''.join(map(str, posbase))
            mutated_str1=RNA.fold(mutated_s)
            mutated_str=mutated_str1[0]
            d=0.0
            g=0.0
            n=len(s)
            for i in range(len(s)):
                if mutated_str[i]!=s[i]:
                    d=d+1
            g=(n-d)/n
            if g==1.0:
                solution.append(mutated_s)
                return g
            else:
                return g

    def __init__(self, fullseq, vrna_md):
        super(TrafoLandscape, self).__init__()

        self._full_sequence = fullseq
        self._model_details = vrna_md
        self._fold_compound = RNA.fold_compound(fullseq, vrna_md)

        # Adjust simulation parameters
        self._RT = 0.61632077549999997
        if vrna_md.temperature != 37.0:
            kelvin = 273.15 + vrna_md.temperature
            self._RT = (self._RT/310.15) * kelvin

        # Private instance variables:
        self._transcript_length = 0
        self._total_time = 0
        self._nodeid = 0

        # Default parameters:
        self._p_min = 0.01  # probability threshold
        self._fpath = 20    # findpath_search_width

        self._k0 = 2e5    # set directly
        self._dG_max = 0  # set using t_slow
        self._dG_min = 0  # set using t_fast

    def test_eval_structure_pt(self):
        print "test_eval_structure_pt\n"
        fc=RNA.fold_compound(seq1)
        energy= fc.eval_structure_pt(struct1_pt) /100; #/100 for dcal

        self.assertEqual("%6.2f" % energy, "%6.2f" % -5.60)
        print  struct1, "[%6.2f" % energy,"]\n"

 def test_centroid(self):
     print "test_centroid\n"
     fc=RNA.fold_compound(align)
     fc.pf()
     (sc,dist) = fc.centroid()
     print  sc,"\tDistance of :  %6.2f" %dist ,"\n"
     self.assertTrue(sc and dist)

def bt(i,j,k,l,d,data=None):
    """
    The backtracking callback must return a list of base pairs
    Here, the base pairs may be given in one of the three ways
    shown below:
    """

    if d == RNA.DECOMP_PAIR_HP:
        """
        1. We create a list of dictionaries with 'i' and 'j'
        keys that specify the coordinates of the base pair (i,j)
        """
        bp = { 'i' : i+1, 'j' : j-1 }

        """
        2. We create a list of tuples (i,j)
        """
        bp = (i+1, j-1)

        """
        3. We create a list of RNA::basepair objects
        """
        bp = RNA.basepair()
        bp.i = i+1
        bp.j = j-1

        return [ bp ]

    return None

def main():
    for monster in monsters:
        # calculate 1) foldcompound 2) partition function 3) base pair probability matrix in that order (!)
        foldmonster = RNA.fold_compound(str(monster))
        pfstruct, pf = foldmonster.pf()
        bppm = foldmonster.bpp()
        plot_bppm(bppm, monster.id)
    print versions_used()

def score_match(query):
    
    motif = 'CCTCCT'
    length = len(motif)
    
    score = RNA.cofold(query + '&' + motif)
        
    return score[1]

 def test_pf(self):
     print "test_pf"
     fc= RNA.fold_compound(seq1)
     (ss,gfe) = fc.pf()
     print ss, "[ %6.2f" %gfe ,"]\n"
     self.assertTrue(ss)
     bp_dis = fc.mean_bp_distance()
     print seq1 ,"\t meanBPDistance : ", bp_dis,"\n"
     self.assertTrue(bp_dis)