MFE Prediction (flat interface)
import RNA
seq = "GAGUAGUGGAACCAGGCUAUGUUUGUGACUCGCAGACUAACA"
(ss, mfe) = RNA.fold(seq)
print("{}\n{} [ {:6.2f} ]".format(seq, ss, mfe))
MFE Prediction (object oriented interface)
import RNA;
sequence = "CGCAGGGAUACCCGCG"
fc = RNA.fold_compound(sequence)
(ss, mfe) = fc.mfe()
print("{} [ {:6.2f} ]".format(ss, mfe))
Suboptimal Structure Prediction
import RNA
sequence = "GGGGAAAACCCC"
RNA.cvar.uniq_ML = 1
subopt_data = { 'counter' : 1, 'sequence' : sequence }
def print_subopt_result(structure, energy, data):
if not structure == None:
print(">subopt {:d}".format(data['counter']))
print("{}\n{} [{:6.2f}]".format(data['sequence'], structure, energy))
data['counter'] = data['counter'] + 1
a = RNA.fold_compound(sequence)
a.subopt_cb(500, print_subopt_result, subopt_data);
Boltzmann Sampling (a.k.a. Probabilistic Backtracing)
import RNA
sequence = "UGGGAAUAGUCUCUUCCGAGUCUCGCGGGCGACGGGCGAUCUUCGAAAGUGGAAUCCGUACUUAUACCGCCUGUGCGGACUACUAUCCUGACCACAUAGU"
def store_structure(s, data):
"""
A simple callback function that stores
a structure sample into a list
"""
if s:
data.append(s)
"""
First we prepare a fold_compound object
"""
md = RNA.md()
md.uniq_ML = 1
fc = RNA.fold_compound(sequence, md)
(ss, mfe) = fc.mfe()
fc.exp_params_rescale(mfe)
fc.pf()
"""
Now we are ready to perform Boltzmann sampling
"""
print("{}".format(fc.pbacktrack5(10)))
print("{}".format(fc.pbacktrack5(50)))
for s in fc.pbacktrack5(20, 10):
print("{}".format(s))
for s in fc.pbacktrack5(100, 50):
print("{}".format(s))
print("{}".format(fc.pbacktrack()))
for s in fc.pbacktrack(100):
print("{}".format(s))
for s in fc.pbacktrack(100, RNA.PBACKTRACK_NON_REDUNDANT):
print("{}".format(s))
num_samples = 500
iterations = 15
d = None
s_list = []
for i in range(0, iterations):
d, ss = fc.pbacktrack(num_samples, d, RNA.PBACKTRACK_NON_REDUNDANT)
s_list = s_list + list(ss)
for s in s_list:
print("{}".format(s))
ss = []
i = fc.pbacktrack5(100, 50, store_structure, ss)
for s in ss:
print("{}".format(s))
ss = list()
i = fc.pbacktrack(100, store_structure, ss)
for s in ss:
print("{}".format(s))
ss = list()
i = fc.pbacktrack(100, store_structure, ss, RNA.PBACKTRACK_NON_REDUNDANT)
for s in ss:
print("{}".format(s))
ss = []
d = None
for i in range(0, iterations):
d, i = fc.pbacktrack(num_samples, store_structure, ss, d, RNA.PBACKTRACK_NON_REDUNDANT)
for s in ss:
print("{}".format(s))
print("{}".format(fc.pbacktrack_sub(10, 50)))
for s in fc.pbacktrack_sub(100, 10, 50):
print("{}".format(s))
for s in fc.pbacktrack_sub(100, 10, 50, RNA.PBACKTRACK_NON_REDUNDANT):
print("{}".format(s))
RNAfold -p –MEA equivalent
import RNA
seq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGAACACGGUGGAUAUGUUCGACGACAGGAUCGGCGCACUA"
fc = RNA.fold_compound(seq)
(mfe_struct, mfe) = fc.mfe()
fc.exp_params_rescale(mfe)
(pp, pf) = fc.pf()
(centroid_struct, dist) = fc.centroid()
centroid_en = fc.eval_structure(centroid_struct)
(MEA_struct, MEA) = fc.MEA()
MEA_en = fc.eval_structure(MEA_struct)
print("{}\n{} ({:6.2f})".format(seq, mfe_struct, mfe))
print("{} [{:6.2f}]".format(pp, pf))
print("{} {{{:6.2f} d={:.2f}}}".format(centroid_struct, centroid_en, dist))
print("{} {{{:6.2f} MEA={:.2f}}}".format(MEA_struct, MEA_en, MEA))
print(" frequency of mfe structure in ensemble {:g}; ensemble diversity {:-6.2f}".format(fc.pr_structure(mfe_struct), fc.mean_bp_distance()))
Fun with Soft Constraints
import RNA
seq1 = "CUCGUCGCCUUAAUCCAGUGCGGGCGCUAGACAUCUAGUUAUCGCCGCAA"
RNA.cvar.dangles = 0
mm_data = { 'dummy': RNA.fold_compound(seq1), 'params': RNA.param() }
revert_NN = {
RNA.DECOMP_PAIR_HP: lambda i, j, k, l, f, p: - f.eval_hp_loop(i, j) - 100,
RNA.DECOMP_PAIR_IL: lambda i, j, k, l, f, p: - f.eval_int_loop(i, j, k, l) - 100,
RNA.DECOMP_PAIR_ML: lambda i, j, k, l, f, p: - p.MLclosing - p.MLintern[0] - (j - i - k + l - 2) * p.MLbase - 100,
RNA.DECOMP_ML_ML_STEM: lambda i, j, k, l, f, p: - p.MLintern[0] - (l - k - 1) * p.MLbase,
RNA.DECOMP_ML_STEM: lambda i, j, k, l, f, p: - p.MLintern[0] - (j - i - k + l) * p.MLbase,
RNA.DECOMP_ML_ML: lambda i, j, k, l, f, p: - (j - i - k + l) * p.MLbase,
RNA.DECOMP_ML_ML_ML: lambda i, j, k, l, f, p: 0,
RNA.DECOMP_ML_UP: lambda i, j, k, l, f, p: - (j - i + 1) * p.MLbase,
RNA.DECOMP_EXT_STEM: lambda i, j, k, l, f, p: - f.eval_ext_stem(k, l),
RNA.DECOMP_EXT_EXT: lambda i, j, k, l, f, p: 0,
RNA.DECOMP_EXT_STEM_EXT: lambda i, j, k, l, f, p: - f.eval_ext_stem(i, k),
RNA.DECOMP_EXT_EXT_STEM: lambda i, j, k, l, f, p: - f.eval_ext_stem(l, j),
}
def MaximumMatching(i, j, k, l, d, data):
return revert_NN[d](i, j, k, l, data['dummy'], data['params'])
fc = RNA.fold_compound(seq1)
fc.sc_add_f(MaximumMatching)
fc.sc_add_data(mm_data, None)
(s, mm) = fc.mfe()
print("{}\n{} (MM: {:d})".format(seq1, s, int(-mm)))