1 /*
2 Last changed Time-stamp: <2006-10-13 22:34:37 xtof>
3 $Id: Energy.cpp,v 1.12 2007/11/03 16:45:58 Kinwalker Exp $
4 */
5 #include <cstdlib>
6 #include <iostream>
7 #include <string>
8 #include <cstring>
9 #include "Energy.h"
10
11
12 extern "C" {
13 #include "fold.h"
14 #include "fold_vars.h"
15 #include "utils.h"
16 #include "pair_mat.h"
17 }
18 extern short * S;
19 extern short * S1;
20 extern short * pair_table;
21
22 static float
23 (*EnergyModel)(std::string sequence, std::string structure) = NULL;
24
25
26
27
28 double
FastEvalEnergy(std::string sequence)29 FastEvalEnergy(std::string sequence){
30 return energy_of_struct_pt(const_cast<char*>(sequence.c_str()), pair_table, S, S1)/100.;
31 //= energy_of_struct(sequence.c_str(), structure.c_str());
32 }
33 //operation: einfuege, wegnehmen
34 //jetzige gesamtenergue
35 //welche loops betroffen
36 //inline int LoopEnergy(int n1, int n2, int type, int type_2,int si1, int sj1, int sp1, int sq1);
37
38
39 // void PartialPath::RemoveFromPairTableWithLog(int first,int second,int removeCatalogIdx){
40 // // cout<<"RemoveFromPairTableWithLog idx "+Str(removeCatalogIdx)<<endl;
41 // currentPartialPathTrail.push_back(removeCatalogIdx);
42 // RemoveFromPairTable(first,second);
43 // }
44
45 // void PartialPath::AddToPairTableWithLog(int first,int second,int addCatalogIdx){
46 // // cout<<"AddToPairTableWithLog idx "+Str(addCatalogIdx)<<endl;
47 // currentPartialPathTrail.push_back(addCatalogIdx+BP_ADD_CONST);
48 // AddToPairTable(first,second);
49 // }
50
51
52
53 // void PartialPath::AddToPairTable(std::pair<int,int> bp){
54 // // cout<<"AddToPairTable"+PrintBasePair(bp) <<std::endl;
55 // pair_table[bp.first]=bp.second;
56 // pair_table[bp.second]=bp.first;
57 // }
58
59
60 // void PartialPath::RemoveFromPairTable(std::pair<int,int> bp){
61 // // cout<<"RemoveFromPairTable" +PrintBasePair(bp)<<std::endl;
62 // pair_table[bp.first]=0;
63 // pair_table[bp.second]=0;
64 // }
65
66
67
68
69 double
EvalEnergy(std::string sequence,std::string structure)70 EvalEnergy(std::string sequence, std::string structure)
71 {
72 bool debug=false;
73 if(debug){
74 std::cout<<"EvalEnergy"<<std::endl;
75 std::cout<<sequence<<std::endl;
76 std::cout<<structure<<std::endl;
77 }
78 return (EnergyModel(sequence, structure));
79 }
80
81
82 void
InitializeEnergyModel(OptionS * OptS,std::string sequence)83 InitializeEnergyModel(OptionS* OptS, std::string sequence)
84 {
85 //hard code to ViennaRNA as engies are being minimized in lots of places
86 InitializeViennaRNA(sequence,OptS->dangle,OptS->transcribed);
87 EnergyModel = FullEnergyModel;
88
89 // if (OptS->Emodel == 'M') {
90 // InitializeViennaRNA(sequence,OptS->dangle,OptS->transcribed);
91 // EnergyModel = FullEnergyModel;
92 // }
93 // else
94 // EnergyModel = NJEnergyModel;
95 }
96
97
98 bool
CanPair(char i,char j)99 CanPair(char i, char j)
100 {
101 return (pair[encode_char(i)][encode_char(j)] > 0);
102 }
103
104 /**
105 * Numerically encodes the sequence
106 */
107 std::vector<int>
EncodeSequence(std::string sequence)108 EncodeSequence(std::string sequence)
109 {
110 std::vector<int> encoded;
111
112 for (std::string::const_iterator base = sequence.begin();
113 base != sequence.end();
114 base++)
115 encoded.push_back(static_cast<int>(encode_char(*base)));
116
117 return (encoded);
118 }
119
120 /**
121 *
122 */
123
124
125
126 void
InitializeViennaRNA(std::string sequence,int dangle,int transcribed)127 InitializeViennaRNA(std::string sequence,int dangle,int transcribed)
128 {
129 int length=sequence.size();
130 S = new short [length+2];
131 S1 = new short [length+2];//NULL;//new short [length+2];
132 pair_table = new short [length+2];
133 // pair_table = make_pair_table(const_cast<char*>(std::string(length,'.').c_str()));
134 //(short *) MG_space(sizeof(short)*(length+1));
135 MakePairTable(const_cast<char*>(std::string(length,'.').c_str()));
136 make_pair_matrix();
137 S[0] = S1[0] = pair_table[0]=transcribed;
138 for (int i=0; i< length; i++) {
139 S[i+1] = encode_char(sequence[i]);
140 S1[i+1] = alias[S[i+1]];
141 }
142
143 // std::cout<<"Dangle "<<dangle<<std::endl;
144 dangles = dangle;//2;
145 // no_closingGU=1; no effect
146
147 initialize_fold(length);
148 update_fold_params();
149 // read_parameter_file("/home/mescalin/mgeis/kinwalker/noTermAU.par"); no effect
150 }
151
152 /**
153 * Calculates the energy of a nucleic acid using the energy model as
154 * implemented in the Vienna RNA package
155 */
156 float
FullEnergyModel(std::string sequence,std::string structure)157 FullEnergyModel(std::string sequence, std::string structure)
158 {
159 // initialize_fold(sequence.length());
160
161 float energy = energy_of_struct(sequence.c_str(), structure.c_str());
162 // free_arrays();
163
164 return (energy);
165 }
166
167 /**
168 * Calculates the energy of a nucleic acid using the Nussinov-Jacobson Model
169 */
170 float
NJEnergyModel(std::string sequence,std::string structure)171 NJEnergyModel(std::string sequence, std::string structure)
172 {
173 int energy = 0;
174 int length = structure.length();
175 std::vector<int> sq = EncodeSequence(sequence);
176 std::vector<int> ptbl = MakePairTable(structure);
177
178 for (int i=0; i<length; i++) {
179
180 if ( ptbl[i] == -1 ) continue; // unpaird positions
181 if ( ptbl[i] < i ) continue; // positions with j<i
182 energy += NJ_energy(sq[i], sq[ptbl[i]]);
183 }
184
185 return (static_cast<float>(energy));
186 }
187
188 /**
189 * Base pair energies rules of the Nussinov-Jacobson Model
190 */
191 int
NJ_energy(int i,int j)192 NJ_energy(int i, int j)
193 {
194 int A, U, G, C;
195 A = encode_char('A');
196 U = encode_char('U');
197 G = encode_char('G');
198 C = encode_char('C');
199
200 if ( (i == G && j == C) || (i == C && j == G) ) return (3); /* GC CG */
201 if ( (i == A && j == U) || (i == U && j == A) ) return (2); /* AU UA */
202 if ( (i == G && j == U) || (i == U && j == G) ) return (1); /* GU UG */
203 return (0); /* unknown base pair */
204 }
205
206 /**
207 * Calculates the energies for a list of structures
208 * [Structure] => [(Enegy, Structure)]
209 */
210 std::vector<std::pair<float,std::string> >
EnergyEvalStructureList(std::string sequence,std::vector<std::string> strList)211 EnergyEvalStructureList(std::string sequence, std::vector<std::string> strList)
212 {
213 std::vector<std::pair<float,std::string> > tupelList;
214
215 for (std::vector<std::string>::const_iterator structure = strList.begin();
216 structure != strList.end();
217 structure++) {
218
219 float energy = FullEnergyModel(sequence, *structure);
220 tupelList.push_back(make_pair(energy, *structure));
221 }
222
223 return (tupelList);
224 }
225
226 /**
227 * Number loops and assign each sequence position its loop-number
228 */
229 std::vector<int>
MakeLoopIndex(std::string structure)230 MakeLoopIndex(std::string structure)
231 {
232 int hx, l, nl;
233 int length;
234 std::vector<int> stack;
235 std::vector<int> loop;
236
237 length = structure.size();
238 stack = std::vector<int>(length+1, -1);
239 loop = std::vector<int>(length+2, -1);
240
241 hx = l = nl = 0;
242 for (int i=0; i<length; i++) {
243 if ( structure[i] == '(' ) {
244 nl++; l = nl;
245 stack[hx++] = i;
246 }
247 loop[i] = l;
248 if ( structure[i] ==')' ) {
249 --hx;
250 if ( hx > 0 )
251 l = loop[stack[hx-1]]; // index of enclosing loop
252 else l = 0; // external loop has index 0
253 if ( hx < 0 ) {
254 Cout("MakeLoopIndex MakePairTable(): "+structure+" unbalanced brackets");
255 exit(EXIT_FAILURE);
256 }
257 // Fatal("MakePairTable(): %s\n%s\n", "unbalanced brackets", structure.c_str());
258 }
259 }
260
261 return (loop);
262 }
263
264 #if 0
265
266 /**
267 * let table[i]=j if (i,j) is pair, -1 if i unpaired indices start at
268 * 0 in this version!
269 */
270 std::vector<int>
271 MakePairTable(std::string structure)
272 {
273 short* pt = NULL;
274 std::vector<int> tbl;
275
276 pt = make_pair_table(structure.c_str());
277 for(int i = 1; i<=pt[0]; i++) {
278 int value = (pt[i]==0) ? -1 : static_cast<int>(pt[i]-1);
279 tbl.push_back(value);
280 }
281 free(pt);
282 return (tbl);
283 }
284
285 #else
286
287 /**
288 * let table[i]=j if (i,j) is pair, -1 if i unpaired indices start at
289 * 0 in this version!
290 */
291 std::vector<int>
MakePairTable(std::string structure)292 MakePairTable(std::string structure)
293 {
294 //Cout("MakePairTable(): "+structure+"\n");
295 int length = structure.size();
296 std::vector<int> stk = std::vector<int>();
297 std::vector<int> tbl(length, -1);
298
299 for (int i=0,j=-1; i<length; i++) {
300 if ( structure[i] == '(' )
301 stk.push_back(i);
302 else if ( structure[i] == ')' ) {
303 j = stk[stk.size()-1];
304 stk.pop_back();
305 if ( j < 0 ) {
306 Cout("MakePairTable(): "+structure+" unbalanced brackets");
307 std::string s=std::string();
308 for(int r=0;r<=tbl[0];r++){
309 s+=Str(tbl[r])+" ";
310 }
311 Cout(s);
312 exit(EXIT_FAILURE);
313 }
314 // if ( j < 0 ) Fatal("MakePairTable(): %s\n%s\n",
315 // "unbalanced brackets",
316 // structure.c_str());
317
318 tbl[i] = j;
319 tbl[j] = i;
320 }
321 }
322
323 if (! stk.empty() ) {
324 Cout("MakePairTable(): "+structure+" unbalanced brackets stk not empty");
325 exit(EXIT_FAILURE);
326 }
327 //Fatal("MakePairTable(): %s\n%s\n", "too few closing brackets", structure.c_str());
328
329 return (tbl);
330 }
331
332 #endif
333
334 std::vector<std::pair<int,int> >
MakeBasePairList(std::string structure)335 MakeBasePairList (std::string structure)
336 {
337 std::vector<int> pTbl = MakePairTable(structure);
338 std::vector<std::pair<int,int> > pLst = std::vector<std::pair<int,int> > ();
339
340 int i = 0;
341 for (std::vector<int>::iterator j = pTbl.begin();
342 j != pTbl.end();
343 j++, i++) {
344
345 // skip unpaired positions and entries with pTbl[i] < i
346 if (*j == -1 || *j < i) continue;
347
348 pLst.push_back(std::make_pair(i, *j));
349 }
350
351 return (pLst);
352 }
353
354 std::vector<std::pair<int,int> >
MakeBasePairList1(std::string structure)355 MakeBasePairList1 (std::string structure)
356 {
357 std::vector<int> pTbl = MakePairTable(structure);
358 std::vector<std::pair<int,int> > pLst = std::vector<std::pair<int,int> > ();
359
360 int i = 0;
361 for (std::vector<int>::iterator j = pTbl.begin();
362 j != pTbl.end();
363 j++, i++) {
364
365 // skip unpaired positions and entries with pTbl[i] < i
366 if (*j == -1 || *j < i) continue;
367
368 pLst.push_back(std::make_pair(i+1, *j+1));
369 }
370
371 return (pLst);
372 }
373
374
375
376 // End of file
377