3 /* Last changed Time-stamp: <2007-08-26 11:59:45 ivo> */
5 compute the duplex structure of two RNA strands,
6 allowing only inter-strand base pairs.
7 see cofold() for computing hybrid structures without
21 #include "energy_par.h"
22 #include "fold_vars.h"
28 /* #include "fold.h" */
30 #include "loop_energies.h"
32 static char rcsid[] UNUSED = "$Id: duplex.c,v 1.8 2007/08/26 10:08:44 ivo Exp $";
34 #define MINPSCORE -2*UNIT
36 #define PRIVATE static
38 #define STACK_BULGE1 1 /* stacking energies for bulges of size 1 */
39 #define NEW_NINIO 1 /* new asymetry penalty */
43 PRIVATE void encode_seqs(const char *s1, const char *s2);
44 PRIVATE short *encode_seq(const char *seq);
46 PRIVATE void find_max_snoop(const char *s1, const char *s2, const int max,
47 const int alignment_length, const int* position,
48 const int delta, const int distance, const int penalty,
49 const int threshloop, const int threshLE, const int threshRE,
50 const int threshDE, const int threshTE, const int threshSE, const int threshD,
51 const int half_stem, const int max_half_stem, const int min_s2,
52 const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const char* name, const int fullStemEnergy);
54 PRIVATE void find_max_snoop_XS(const char *s1, const char *s2, const int **access_s1, const int max,
55 const int alignment_length, const int* position, const int *position_j,
56 const int delta, const int distance, const int penalty,
57 const int threshloop, const int threshLE, const int threshRE,
58 const int threshDE, const int threshTE, const int threshSE, const int threshD,
59 const int half_stem, const int max_half_stem, const int min_s2,
60 const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const char *name, const int fullStemEnergy);
66 PRIVATE char * alisnoop_backtrack(int i, int j, const char ** s2,
67 int* Duplex_El, int* Duplex_Er, int* Loop_E, int *Loop_D, int *u,
68 int *pscd, int *psct, int *pscg,
69 const int penalty, const int threshloop,
70 const int threshLE, const int threshRE, const int threshDE, const int threshD,
71 const int half_stem, const int max_half_stem,
72 const int min_s2, const int max_s2, const int min_s1,
73 const int max_s1, const int min_d1, const int min_d2,
74 const short **S1, const short **S2);
77 PRIVATE char * snoop_backtrack(int i, int j, const char* s2, int* Duplex_El, int* Duplex_Er, int* Loop_E, int *Loop_D, int *u,
78 const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
80 const int half_stem, const int max_half_stem,
81 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2);
83 PRIVATE char * snoop_backtrack_XS(int i, int j, const char* s2, int* Duplex_El, int* Duplex_Er, int* Loop_E, int *Loop_D, int *u,
84 const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
86 const int half_stem, const int max_half_stem,
87 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2);
92 PRIVATE int compare(const void *sub1, const void *sub2);
93 PRIVATE int covscore(const int *types, int n_seq);
94 PRIVATE short * aliencode_seq(const char *sequence);
96 PUBLIC int snoop_subopt_sorted=0; /* from subopt.c, default 0 */
105 #define MIN2(A, B) ((A) < (B) ? (A) : (B))
106 #define MAX2(A, B) ((A) > (B) ? (A) : (B))
108 PRIVATE paramT *P = NULL;
110 PRIVATE int **c = NULL; /* energy array, given that i-j pair */
111 PRIVATE int **r = NULL;
112 PRIVATE int **lc = NULL; /* energy array, given that i-j pair */
113 PRIVATE int **lr = NULL;
114 PRIVATE int **c_fill = NULL;
115 PRIVATE int **r_fill = NULL;
116 PRIVATE int **lpair = NULL;
119 PRIVATE short *S1 = NULL, *SS1 = NULL, *S2 = NULL, *SS2 = NULL;
120 PRIVATE short *S1_fill = NULL, *SS1_fill = NULL, *S2_fill = NULL, *SS2_fill = NULL;
121 PRIVATE int n1,n2; /* sequence lengths */
123 extern int cut_point;
125 PRIVATE int delay_free=0;
126 /*--------------------------------------------------------------------------*/
128 snoopT alisnoopfold(const char **s1, const char **s2,
129 const int penalty, const int threshloop,
130 const int threshLE, const int threshRE, const int threshDE, const int threshD,
131 const int half_stem, const int max_half_stem,
132 const int min_s2, const int max_s2, const int min_s1,
133 const int max_s1, const int min_d1, const int min_d2) {
136 int i, j, E, l1,Emin=INF, i_min=0, j_min=0;
142 folden **foldlist; folden **foldlist_XS;
143 int Duplex_El, Duplex_Er,pscd,psct,pscg;
147 short **Sali1,**Sali2;
148 int *type,*type2,*type3;
149 Duplex_El=0;Duplex_Er=0;Loop_E=0; Loop_D=0;pscd=0;psct=0;pscg=0;
150 snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS);
151 n1 = (int) strlen(s1[0]);
152 n2 = (int) strlen(s2[0]);
154 for (s=0; s1[s]!=NULL; s++);
156 for (s=0; s2[s]!=NULL; s++);
157 if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");
159 if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
160 snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
164 c = (int **) space(sizeof(int *) * (n1+1));
165 r = (int **) space(sizeof(int *) * (n1+1));
166 for (i=0; i<=n1; i++) {
167 c[i] = (int *) space(sizeof(int) * (n2+1));
168 r[i] = (int *) space(sizeof(int) * (n2+1));
169 for(j=n2; j>-1; j--){
174 Sali1 = (short **) space((n_seq+1)*sizeof(short *));
175 Sali2 = (short **) space((n_seq+1)*sizeof(short *));
176 for (s=0; s<n_seq; s++) {
177 if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
178 if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
179 Sali1[s] = aliencode_seq(s1[s]);
180 Sali2[s] = aliencode_seq(s2[s]);
182 type = (int *) space(n_seq*sizeof(int));
183 type2 = (int *) space(n_seq*sizeof(int));
184 type3 = (int *) space(n_seq*sizeof(int));
185 /* encode_seqs(s1, s2); */
186 for (i=6; i<=n1-5; i++) {
188 for (s=0; s<n_seq; s++) {
191 U = (U==(n_seq)*4?1:0);
192 for (j=n2-min_d2; j>min_d1; j--) {
193 int type4, k,l,psc,psc2,psc3;
194 for (s=0; s<n_seq; s++) {
195 type[s] = pair[Sali1[s][i]][Sali2[s][j]];
197 psc = covscore(type, n_seq);
198 for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
199 c[i][j] = (psc>=MINPSCORE) ? (n_seq*P->DuplexInit) : INF;
200 if (psc<MINPSCORE) continue;
201 if(/* pair[Sali1[i+1]][Sali2[j-1]] && */
202 U && j < max_s1 && j > min_s1 &&
203 j > n2 - max_s2 - max_half_stem &&
204 j < n2 -min_s2 -half_stem ) { /*constraint on s2 and i*/
209 for (s=0; s<n_seq; s++) {
210 type2[s]= pair[Sali1[s][i-3]][Sali2[s][k+1]];
211 type3[s]= pair[Sali1[s][i-4]][Sali2[s][k+1]];
213 psc2 = covscore(type2, n_seq);
214 psc3 = covscore(type3, n_seq);
215 if(psc2 > MINPSCORE){
216 r[i][j]=MIN2(r[i][j],c[i-3][k+1]+temp->energy);
218 if(psc3 > MINPSCORE){
219 r[i][j]=MIN2(r[i][j],c[i-4][k+1]+temp->energy);
224 /* dangle 5'SIDE relative to the mRNA */
225 for (s=0; s<n_seq; s++) {
226 c[i][j] += E_ExtLoop(type[s], Sali1[s][i-1],Sali2[s][j+1],P);
228 for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
229 for (l=j+1; l<=n2 ; l++) {
230 if (i-k+l-j>2*MAXLOOP_L-2) break;
231 if (abs(i-k-l+j) >= ASS ) continue;
232 for (E=s=0; s<n_seq; s++) {
233 type4 = pair[Sali1[s][k]][Sali2[s][l]];
234 if (type4==0) type4=7;
235 E += E_IntLoop(i-k-1, l-j-1, type4, rtype[type[s]],
236 Sali1[s][k+1], Sali2[s][l-1], Sali1[s][i-1], Sali2[s][j+1],P);
238 c[i][j] = MIN2(c[i][j], c[k][l] + E);
239 r[i][j] = MIN2(r[i][j], r[k][l] + E);
245 for (s=0; s<n_seq; s++) {
246 E+= E_ExtLoop(rtype[type[s]], Sali2[s][j-1], Sali1[s][i+1], P);
248 *** if (i<n1) E += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
249 *** if (j>1) E += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
250 *** if (type[s]>2) E += P->TerminalAU;
254 Emin=E; i_min=i; j_min=j;
259 printf("no target found under the constraints chosen\n");
260 for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
263 for(s=0; s<n_seq;s++){
267 free(Sali1); free(Sali2);
268 free(S2); free(SS1); free(SS2);free(type);free(type2);free(type3);
273 struc = alisnoop_backtrack(i_min, j_min,(const char**) s2,
274 &Duplex_El, &Duplex_Er, &Loop_E,
275 &Loop_D, &u, &pscd, &psct, &pscg,
276 penalty, threshloop, threshLE,
277 threshRE,threshDE, threshD,
278 half_stem, max_half_stem, min_s2,
279 max_s2, min_s1, max_s1, min_d1,
280 min_d2,(const short**) Sali1,(const short**) Sali2);
281 /* if (i_min<n1-5) i_min++; */
282 /* if (j_min>6 ) j_min--; */
283 l1 = strchr(struc, '&')-struc;
287 mfe.Duplex_Er = (float) Duplex_Er/100;
288 mfe.Duplex_El = (float) Duplex_El/100;
289 mfe.Loop_D = (float) Loop_D/100;
290 mfe.Loop_E = (float) Loop_E/100;
291 mfe.energy = (float) Emin/100 ;
292 /* mfe.fullStemEnergy = (float) fullStemEnergy/100; */
295 mfe.structure = struc;
296 for(s=0; s<n_seq;s++){
297 free(Sali1[s]);free(Sali2[s]);
299 free(Sali1);free(Sali2);free(type);free(type2);free(type3);
302 for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
305 free(S2); free(SS1); free(SS2);
310 PUBLIC snoopT *alisnoop_subopt(const char **s1, const char **s2, int delta, int w,
311 const int penalty, const int threshloop,
312 const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
313 const int distance, const int half_stem, const int max_half_stem,
314 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2) {
317 short **Sali1, **Sali2;
318 /* printf("%d %d\n", min_s2, max_s2); */
319 int i,j,s,n_seq, n1, n2, E, n_subopt=0, n_max;
325 int Duplex_El, Duplex_Er, Loop_E,pscd,psct,pscg;
327 Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;pscd=0;psct=0;pscg=0;
331 subopt = (snoopT *) space(n_max*sizeof(snoopT));
333 mfe = alisnoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE,threshD,
334 half_stem, max_half_stem,
335 min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
341 thresh = MIN2((int) ((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E)*100+0.1 + 410) + delta, threshTE );
342 /* subopt[n_subopt++]=mfe; */
346 for (s=0; s1[s]!=NULL; s++);
348 Sali1 = (short **) space((n_seq+1)*sizeof(short *));
349 Sali2 = (short **) space((n_seq+1)*sizeof(short *));
350 for (s=0; s<n_seq; s++) {
351 if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
352 if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
353 Sali1[s] = aliencode_seq(s1[s]);
354 Sali2[s] = aliencode_seq(s2[s]);
356 Sali1[n_seq]=NULL; Sali2[n_seq]=NULL;
357 type = (int *) space(n_seq*sizeof(int));
358 for (i=n1; i>1; i--){
359 for (j=1; j<=n2; j++) {
360 int ii,jj, Ed,psc,skip;
361 for (s=0; s<n_seq; s++) {
362 type[s] = pair[Sali2[s][j]][Sali1[s][i]];
364 psc = covscore(type, n_seq);
365 for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
366 if (psc<MINPSCORE) continue;
368 for (s=0; s<n_seq; s++) {
369 /* if (i<n1-5) Ed += P->dangle3[type[s]][Sali1[s][i+1]]; */
370 /* if (j>6) Ed += P->dangle5[type[s]][Sali2[s][j-1]]; */
371 if (type[s]>2) Ed += P->TerminalAU;
373 if (Ed>thresh) continue;
374 /* too keep output small, remove hits that are dominated by a
375 better one close (w) by. For simplicity we do test without
376 adding dangles, which is slightly inaccurate.
379 for (skip=0, ii=MAX2(i-w,1); (ii<=MIN2(i+w,n1)) && type; ii++) {
380 for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n2); jj++)
381 if (r[ii][jj]<E) {skip=1; break;}
386 struc = alisnoop_backtrack(i,j,s2, &Duplex_El,
387 &Duplex_Er, &Loop_E, &Loop_D, &u, &pscd, &psct,&pscg,
388 penalty, threshloop,threshLE,threshRE,threshDE, threshD,
389 half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,(const short int**) Sali1,(const int short **) Sali2);
391 if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
392 (Duplex_Er + Duplex_El) > threshDE ||
393 (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
394 (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
395 /* printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n" */
396 /* " Duplex_Er + Duplex_El %d threshDE %d \n" */
397 /* " Duplex_Er + Duplex_El + Loop_E %d threshTE %d \n" */
398 /* " Duplex_Er + Duplex_El + Loop_E + Loop_D %d threshSE %d \n", */
399 /* Duplex_Er , threshRE , Duplex_El ,threshLE, */
400 /* Duplex_Er + Duplex_El, threshDE, */
401 /* Duplex_Er + Duplex_El+ Loop_E , threshTE, */
402 /* Duplex_Er + Duplex_El+ Loop_E + Loop_D, threshSE); */
412 if (n_subopt+1>=n_max) {
414 subopt = (snoopT *) xrealloc(subopt, n_max*sizeof(snoopT));
417 subopt[n_subopt].i = i-5;
418 subopt[n_subopt].j = j-5;
419 subopt[n_subopt].u = u-5;
420 subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
421 subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
422 subopt[n_subopt].Loop_E = Loop_E * 0.01;
423 subopt[n_subopt].Loop_D = Loop_D * 0.01;
424 subopt[n_subopt].energy = (Duplex_Er +Duplex_El + Loop_E + Loop_D + 410) * 0.01 ;
425 subopt[n_subopt].pscd = pscd * 0.01;
426 subopt[n_subopt].psct = -psct * 0.01;
427 subopt[n_subopt++].structure = struc;
430 Duplex_Er=0; Duplex_El=0; Loop_E=0; Loop_D=0;u=0;pscd=0;psct=0;
434 for (i=0; i<=n1; i++) {free(c[i]);free(r[i]);}
436 for (s=0; s<n_seq; s++) {
437 free(Sali1[s]); free(Sali2[s]);
439 free(Sali1); free(Sali2); free(type);
441 if (snoop_subopt_sorted)
442 qsort(subopt, n_subopt, sizeof(snoopT), compare);
443 subopt[n_subopt].i =0;
444 subopt[n_subopt].j =0;
445 subopt[n_subopt].structure = NULL;
455 PRIVATE char *alisnoop_backtrack(int i, int j, const char ** snoseq, int *Duplex_El,
456 int *Duplex_Er, int *Loop_E, int *Loop_D, int *u,
457 int *pscd, int *psct, int *pscg,
458 const int penalty, const int threshloop, const int threshLE,
459 const int threshRE, const int threshDE, const int threshD, const int half_stem,
460 const int max_half_stem,
461 const int min_s2, const int max_s2, const int min_s1,
463 const int min_d1, const int min_d2,const short **Sali1, const short **Sali2) {
464 /* backtrack structure going backwards from i, and forwards from j
465 return structure in bracket notation with & as separator */
466 int k, l, *type,*type2,*type3,type4, E, traced, i0, j0,s,n_seq,psc;
467 int traced_r=0; /* flag for following backtrack in c or r */
468 char *st1, *st2, *struc;
470 n1 = (int) Sali1[0][0];
471 n2 = (int) Sali2[0][0];
473 for (s=0; Sali1[s]!=NULL; s++);
475 for (s=0; Sali2[s]!=NULL; s++);
476 if (n_seq != s) nrerror("unequal number of sequences in alibacktrack()\n");
478 st1 = (char *) space(sizeof(char)*(n1+1));
479 st2 = (char *) space(sizeof(char)*(n2+1));
480 type = (int *) space(n_seq*sizeof(int));
481 type2 = (int *) space(n_seq*sizeof(int));
482 type3 = (int *) space(n_seq*sizeof(int));
486 folden **foldlist, **foldlist_XS;
487 snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
488 i0=i; j0=j; /* MIN2(i+1,n1); j0=MAX2(j-1,1);!modified */
489 for (s=0; s<n_seq; s++) {
490 type[s] = pair[Sali1[s][i]][Sali2[s][j]];
491 if(type[s]==0) type[s] = 7;
492 *Duplex_Er += E_ExtLoop(rtype[type[s]], (j>1) ? Sali2[s][j-1] : -1, (i<n1) ? Sali1[s][i+1] : -1, P);
494 *** if (i<n1) *Duplex_Er += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
495 *** if (j>1) *Duplex_Er += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
496 *** if (type[s]>2) *Duplex_Er += P->TerminalAU;
499 while (i>0 && j<=n2-min_d2 ) {
501 E = r[i][j]; traced=0;
504 for (s=0; s<n_seq; s++) {
505 type[s] = pair[Sali1[s][i]][Sali2[s][j]];
507 psc = covscore(type,n_seq);
508 for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
511 for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
512 for (l=j+1; l<=n2 ; l++) {
514 if (i-k+l-j>2*MAXLOOP_L-2) break;
515 if (abs(i-k-l+j) >= ASS) continue;
516 for (s=LE=0; s<n_seq; s++) {
517 type4 = pair[Sali1[s][k]][Sali2[s][l]];
518 if (type4==0) type4=7;
519 LE += E_IntLoop(i-k-1, l-j-1, type4, rtype[type[s]], Sali1[s][k+1], Sali2[s][l-1], Sali1[s][i-1], Sali2[s][j+1],P);
521 if (E == r[k][l]+LE) {
532 for (s=0; s<n_seq; s++) {
535 U = (U==(n_seq)*4?1:0);
536 if(/* pair[Sali1[i+1]][Sali2[j-1]] && */ /* only U's are allowed */
537 U && j < max_s1 && j > min_s1 &&
538 j > n2 - max_s2 - max_half_stem &&
539 j < n2 -min_s2 -half_stem ) {
541 max_k = MIN2(n2-min_s2,j+max_half_stem+1);
542 min_k = MAX2(j+half_stem+1, n2-max_s2);
548 for (s=0; s<n_seq; s++) {
549 type2[s]= pair[Sali1[s][i-3]][Sali2[s][k+1]];
550 type3[s]= pair[Sali1[s][i-4]][Sali2[s][k+1]];
552 psc2 = covscore(type2, n_seq);
553 psc3 = covscore(type3, n_seq);
554 if(psc2>MINPSCORE /*&& pair[Sali1[i-4]][Sali2[k+2]]*/ ){ /* introduce structure from RNAfold */
555 if(E==c[i-3][k+1]+temp->energy){
556 *Loop_E=temp->energy;
560 /* int fix_ij=indx[k-1+1]+j+1; */
561 for(a=0; a< MISMATCH ;a++){
562 for(b=0; b< MISMATCH ; b++){
563 int ij=indx[k-1-a+1]+j+1+b;
564 if(cLoop[ij]==temp->energy) {
566 struc_loop=alisnobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1,psct);
577 if (psc3>MINPSCORE /*&& pair[Sali1[i-5]][Sali2[k+2]]*/){ /* introduce structure from RNAfold */
578 if(E==c[i-4][k+1]+temp->energy){
579 *Loop_E=temp->energy;
583 /* int fix_ij=indx[k-1+1]+j+1; */
584 for(a=0; a< MISMATCH ;a++){
585 for(b=0; b< MISMATCH ; b++){
586 int ij=indx[k-1-a+1]+j+1+b;
587 if(cLoop[ij]==temp->energy) {
589 struc_loop=alisnobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1,psct);
601 } /* while temp-> next */
606 E = c[i][j]; traced=0;
609 for (s=0; s<n_seq; s++) {
610 type[s] = pair[Sali1[s][i]][Sali2[s][j]];
612 psc = covscore(type,n_seq);
613 for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
616 if (!type) nrerror("backtrack failed in fold duplex c");
617 for (k=i-1; (i-k)<MAXLOOP_L; k--) {
618 for (l=j+1; l<=n2; l++) {
620 if (i-k+l-j>2*MAXLOOP_L-2) break;
621 if (abs(i-k-l+j) >= ASS) continue;
622 for (s=LE=0; s<n_seq; s++) {
623 type4 = pair[Sali1[s][k]][Sali2[s][l]];
624 if (type4==0) type4=7;
625 LE += E_IntLoop(i-k-1, l-j-1, type4, rtype[type[s]], Sali1[s][k+1], Sali2[s][l-1], Sali1[s][i-1], Sali2[s][j+1],P);
627 if (E == c[k][l]+LE) {
638 for (s=0; s<n_seq; s++) {
640 correction = E_ExtLoop(type[s], (i>1) ? Sali1[s][i-1] : -1, (j<n2) ? Sali2[s][j+1] : -1, P);
641 *Duplex_El += correction;
644 *** if (i>1) {E -= P->dangle5[type[s]][Sali1[s][i-1]]; *Duplex_El +=P->dangle5[type[s]][Sali1[s][i-1]];}
645 *** if (j<n2) {E -= P->dangle3[type[s]][Sali2[s][j+1]]; *Duplex_El +=P->dangle3[type[s]][Sali2[s][j+1]];}
646 *** if (type[s]>2) {E -= P->TerminalAU; *Duplex_El +=P->TerminalAU;}
649 if (E != n_seq * P->DuplexInit) {
650 nrerror("backtrack failed in fold duplex end");
656 /* struc = (char *) space(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
657 struc = (char *) space(i0-i+1+n2-1+1+2); /* declare final duplex structure */
659 struc2 = (char *) space(n2+1);
660 /* char * struct_const; */
661 for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
662 /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
663 /* char * struct_const; */
664 /* struct_const = (char *) space(sizeof(char)*(n2+1)); */
665 for (k=1; k<=n2; k++) {
666 if (!st2[k-1]) st2[k-1] = struc_loop[k-1];/* '.'; */
667 struc2[k-1] = st2[k-1];/* '.'; */
668 /* if (k>=j0 && k<=j){ */
669 /* struct_const[k-1]='x'; */
672 /* if(k<j0) {struct_const[k-1]='<';} */
673 /* if(k>j) {struct_const[k-1]='>';} */
677 /* char duplexseq_1[j0+1]; */
678 /* char duplexseq_2[n2-j+3]; */
680 char **duplexseq_1, **duplexseq_2;
681 duplexseq_1 = (char**) space((n_seq+1) * sizeof(char*));
682 duplexseq_2 = (char**) space((n_seq+1) * sizeof(char*));
683 for(s=0; s<n_seq; s++){
684 duplexseq_1[s] = (char*) space((j0)*sizeof(char)); /* modfied j0+1 */
685 duplexseq_2[s] = (char*) space((n2-j+2)*sizeof(char)); /* modified j+3 */
686 strncpy(duplexseq_1[s], snoseq[s], j0-1); /* modified j0 */
687 strcpy(duplexseq_2[s], snoseq[s] + j); /* modified j-1 */
688 duplexseq_1[s][j0-1]='\0'; /* modified j0 */
689 duplexseq_2[s][n2-j+1]='\0';/* modified j+2 */
691 duplexseq_1[n_seq]=NULL;
692 duplexseq_2[n_seq]=NULL;
694 temp=aliduplexfold((const char**)duplexseq_1, (const char**)duplexseq_2);
695 *Loop_D = MIN2(0,-410 + (int) 100 * temp.energy*n_seq);
697 int l1,ibegin, iend, jbegin, jend;
698 l1=strchr(temp.structure, '&')-temp.structure;
702 jend =temp.j+strlen(temp.structure)-l1-2-1;
703 for(k=ibegin+1; k<=iend+1; k++){
704 struc2[k-1]=temp.structure[k-ibegin-1];
706 for(k=jbegin+j; k<=jend+j; k++){
707 struc2[k-1]=temp.structure[l1+k-j-jbegin+1];
710 for(s=0; s<n_seq; s++){
711 free(duplexseq_1[s]);
712 free(duplexseq_2[s]);
714 free(duplexseq_1);free(duplexseq_2);
715 free(temp.structure);
717 strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&");
718 /* strcat(struc, st2); */
719 strncat(struc, struc2+5, strlen(struc2)-10);
722 free(st1); free(st2);
723 free(type);free(type2);free(type3);
734 void Lsnoop_subopt(const char *s1, const char *s2, int delta, int w,
735 const int penalty, const int threshloop,
736 const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
738 const int half_stem, const int max_half_stem,
739 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int alignment_length, const char* name, const int fullStemEnergy)
746 /* int nsubopt=10; */
747 n1 = (int) strlen(s1);
748 n2 = (int) strlen(s2);
750 position = (int*) space((n1+3)*sizeof(int));
753 /* int Eminj, Emin_l; */
754 int i, j; /* l1, Emin=INF, i_min=0, j_min=0; */
760 folden **foldlist, **foldlist_XS;
761 int Duplex_El, Duplex_Er;
765 Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
766 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
767 if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
768 snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
772 lc = (int **) space(sizeof(int *) * (5));
773 lr = (int **) space(sizeof(int *) * (5));
774 for (i=0; i<5; i++) {
775 lc[i] = (int *) space(sizeof(int) * (n2+1));
776 lr[i] = (int *) space(sizeof(int) * (n2+1));
777 for(j=n2; j>-1; j--){
783 for (i=1; i<=n1; i++) {
789 for (j=n2-min_d2; j>min_d1; j--) {
791 type = pair[S1[i]][S2[j]];
792 lc[idx][j] = (type) ? P->DuplexInit + 2*penalty : INF;
795 if( /*pair[S1[i+1]][S2[j-1]] && check that we have a solid base stack after the mLoop */
796 j < max_s1 && j > min_s1 &&
797 j > n2 - max_s2 - max_half_stem &&
798 j < n2 -min_s2 -half_stem && S1[i-2]==4) { /*constraint on s2 and i*/
800 max_k = MIN2(n2-min_s2,j+max_half_stem+1);
801 min_k = MAX2(j+half_stem+1, n2-max_s2);
802 for(k=min_k; k <= max_k ; k++){
803 if(mLoop[indx[k-1]+j+1] < 0){
805 if(pair[S1[i-3]][S2[k]] /*genau zwei ungepaarte nucleotiden --NU--*/
806 && mLoop[indx[k-1]+j+1] < threshloop){
807 lr[idx][j]=MIN2(lr[idx][j], lc[idx_3][k]+mLoop[indx[k-1]+j+1]);
809 else if(pair[S1[i-4]][S2[k]] && mLoop[indx[k-1]+j+1] < threshloop){/*--NUN--*/
810 lr[idx][j]=MIN2(lr[idx][j], lc[idx_4][k]+mLoop[indx[k-1]+j+1]);
814 /* dangle 5'SIDE relative to the mRNA */
815 lc[idx][j] += E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, P);
817 *** if (i>1) lc[idx][j] += P->dangle5[type][SS1[i-1]];
818 *** if (j<n2) lc[idx][j] += P->dangle3[type][SS2[j+1]];
819 *** if (type>2) lc[idx][j] += P->TerminalAU;
823 type2=pair[S1[i-1]][S2[j+1]];
825 lc[idx][j]=MIN2(lc[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+2*penalty, lc[idx][j]);
826 lr[idx][j]=MIN2(lr[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+2*penalty, lr[idx][j]);
830 type2=pair[S1[i-2]][S2[j+2]];
832 lc[idx][j]=MIN2(lc[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P)+4*penalty, lc[idx][j]);
833 lr[idx][j]=MIN2(lr[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P)+4*penalty, lr[idx][j]);
837 type2 = pair[S1[i-3]][S2[j+3]];
839 lc[idx][j]=MIN2(lc[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+6*penalty,lc[idx][j]);
840 lr[idx][j]=MIN2(lr[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+6*penalty,lr[idx][j]);
844 *** (type>2?P->TerminalAU:0)+(i<(n1)?P->dangle3[rtype[type]][SS1[i+1]]+penalty:0)+(j>1?P->dangle5[rtype[type]][SS2[j-1]]+penalty:0)
846 min_colonne=MIN2(lr[idx][j]+E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P), min_colonne);
848 position[i]=min_colonne;
849 if(max>=min_colonne){
856 free(S1); free(S2); free(SS1); free(SS2);
858 find_max_snoop(s1, s2, max, alignment_length, position, delta,
859 distance, penalty, threshloop, threshLE, threshRE, threshDE,
860 threshTE, threshSE, threshD, half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name, fullStemEnergy);
862 for (i=1; i<5; i++) {free(lc[i]);free(lr[i]);}
863 free(lc[0]);free(lr[0]);
871 void Lsnoop_subopt_list(const char *s1, const char *s2, int delta, int w,
872 const int penalty, const int threshloop,
873 const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
875 const int half_stem, const int max_half_stem,
876 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int alignment_length,const char *name,const int fullStemEnergy)
883 /* int nsubopt=10; */
884 n1 = (int) strlen(s1);
885 n2 = (int) strlen(s2);
887 position = (int*) space((n1+3)*sizeof(int));
890 /* int Eminj, Emin_l; */
891 int i, j;/* l1, Emin=INF, i_min=0, j_min=0; */
897 folden **foldlist, **foldlist_XS;
898 int Duplex_El, Duplex_Er;
902 Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
903 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
904 if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
905 snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
909 lpair = (int **) space(sizeof(int *) * (6));
910 lc = (int **) space(sizeof(int *) * (6));
911 lr = (int **) space(sizeof(int *) * (6));
912 for (i=0; i<6; i++) {
913 lc[i] = (int *) space(sizeof(int) * (n2+1));
914 lr[i] = (int *) space(sizeof(int) * (n2+1));
915 lpair[i] = (int *) space(sizeof(int) * (n2+1));
916 for(j=n2; j>-1; j--){
923 int lim_maxj=n2-min_d2 ;
926 for (i=5; i<=lim_maxi; i++) {
933 for (j=lim_maxj; j>lim_minj; j--) {
934 int type, type2;/* E, k,l; */
935 type = pair[S1[i]][S2[j]];
936 lpair[idx][j] = type;
937 lc[idx][j] = (type) ? P->DuplexInit + 2*penalty : INF;
940 if( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
941 j < max_s1 && j > min_s1 &&
942 j > n2 - max_s2 - max_half_stem &&
943 j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
945 max_k = MIN2(n2-min_s2,j+max_half_stem+1);
946 min_k = MAX2(j+half_stem+1, n2-max_s2);
951 /* if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour */
952 if(lpair[idx_3][k+1] /*&& lpair[idx_4][k+2]*/){
953 lr[idx][j]=MIN2(lr[idx][j], lc[idx_3][k+1]+temp->energy);/*--NU--*/
955 /*else*/ if(lpair[idx_4][k+1]){/*--NUN--*/
956 lr[idx][j]=MIN2(lr[idx][j], lc[idx_4][k+1]+temp->energy);
962 /* dangle 5'SIDE relative to the mRNA */
963 lc[idx][j] += E_ExtLoop(type, SS1[i-1] , SS2[j+1] , P);
965 *** lc[idx][j] += P->dangle5[type][SS1[i-1]];
966 *** lc[idx][j] += P->dangle3[type][SS2[j+1]];
967 *** if (type>2) lc[idx][j] += P->TerminalAU;
969 /* if(j<n2 && i>1){ */
970 /* type2=pair[S1[i-1]][S2[j+1]]; */
971 type2=lpair[idx_1][j+1];
973 lc[idx][j]=MIN2(lc[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+2*penalty, lc[idx][j]);
974 lr[idx][j]=MIN2(lr[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+2*penalty, lr[idx][j]);
977 /* if(j<n2-1 && i>2){ */
978 /* type2=pair[S1[i-2]][S2[j+2]]; */
979 type2=lpair[idx_2][j+2];
981 lc[idx][j]=MIN2(lc[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P), lc[idx][j]);
982 lr[idx][j]=MIN2(lr[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P), lr[idx][j]);
985 /* if(j<n2-2 && i>3){ */
986 /* type2 = pair[S1[i-3]][S2[j+3]]; */
987 type2 =lpair[idx_3][j+3];
989 lc[idx][j]=MIN2(lc[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+6*penalty,lc[idx][j]);
990 lr[idx][j]=MIN2(lr[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+6*penalty,lr[idx][j]);
993 /* min_colonne=MIN2(lr[idx][j]+(type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]], min_colonne); */
995 bla=lr[idx][j]+E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1], P)+2*penalty;
996 min_colonne=MIN2(bla, min_colonne);
998 position[i]=min_colonne;
999 if(max>=min_colonne){
1006 free(S1); free(S2); free(SS1); free(SS2);
1008 find_max_snoop(s1, s2, max, alignment_length, position,
1009 delta, distance, penalty, threshloop,
1010 threshLE, threshRE, threshDE, threshTE, threshSE, threshD,
1011 half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name, fullStemEnergy);
1013 for (i=1; i<6; i++) {free(lc[i]);free(lr[i]);free(lpair[i]);}
1014 free(lc[0]);free(lr[0]);free(lpair[0]);
1015 free(lc);free(lr);free(lpair);
1020 PRIVATE void find_max_snoop(const char *s1, const char *s2,const int max, const int alignment_length, const int* position, const int delta,
1021 const int distance, const int penalty, const int threshloop, const int threshLE, const int threshRE,
1022 const int threshDE, const int threshTE, const int threshSE, const int threshD,
1023 const int half_stem, const int max_half_stem, const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const char* name, const int fullStemEnergy)
1027 int threshold = MIN2(threshTE , max + delta );
1028 /* printf("threshTE %d max %d\n", threshTE, max); */
1029 /* #pragma omp parallel for */
1030 /* for(pos=n1+1;pos>distance;pos--){ */
1033 if(position[pos]<(threshold)){
1035 search_range=distance+1;
1036 while(--search_range){
1037 if(position[pos-search_range]<=position[pos-temp_min]){
1038 temp_min=search_range;
1042 int begin=MAX2(6, pos-alignment_length+1);
1043 char *s3 = (char*) space(sizeof(char)*(pos-begin+3+12));
1044 strcpy(s3, "NNNNN");
1045 strncat(s3, (s1+begin-1), pos-begin+2);
1046 strcat(s3,"NNNNN\0");
1047 /* printf("%s s3\n", s3); */
1049 test = snoopfold(s3, s2, penalty, threshloop, threshLE, threshRE, threshDE, threshD, half_stem, max_half_stem, min_s2, max_s2, min_s1,
1050 max_s1, min_d1, min_d2,fullStemEnergy);
1051 if(test.energy==INF){
1055 if(test.Duplex_El > threshLE * 0.01 || test.Duplex_Er > threshRE * 0.01 ||
1056 test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
1057 (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE*0.01) {
1058 free(test.structure);free(s3);
1062 l1 = strchr(test.structure, '&')-test.structure;
1065 if(test.i > strlen(s3)-10){
1073 char *target_struct = (char*) space(sizeof(char) * (strlen(test.structure)+1));
1074 strncpy(target_struct, test.structure+shift, l1);
1075 strncat(target_struct, test.structure + (strchr(test.structure, '&')-
1076 test.structure), strlen(test.structure) - (strchr(test.structure, '&')-
1078 strcat(target_struct,"\0");
1080 target = (char *) space(l1+1);
1081 strncpy(target, (s3+test.i+5-l1), l1);
1084 s4 = (char*) space(sizeof(char) *(strlen(s2)-9));
1085 strncpy(s4, s2+5, strlen(s2)-10);
1086 s4[strlen(s2)-10]='\0';
1087 printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + 4.1 ) (%5.2f) \n%s&%s\n",
1088 target_struct,begin + test.i-5-l1,begin + test.i -6 , begin + test.u -6,
1089 test.j+1, test.j + (strrchr(test.structure,'>') - strchr(test.structure,'>'))+1 ,
1090 test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10, test.Duplex_El,
1091 test.Duplex_Er, test.Loop_E, test.Loop_D,test.fullStemEnergy, target,s4);
1096 temp_seq = (char*) space(sizeof(char)*(l1+n2-9));
1097 temp_struc = (char*) space(sizeof(char)*(l1+n2-9));
1098 strcpy(temp_seq, target);
1099 strcat(temp_seq, s4);
1100 strncpy(temp_struc, target_struct, l1);
1101 strcat(temp_struc, target_struct+l1+1);
1102 temp_seq[n2+l1-10]='\0';
1103 temp_struc[n2+l1-10]='\0';
1105 char str[16];char upos[16];
1106 strcpy(psoutput,"sno_");
1107 sprintf(str,"%d",count);
1108 strcat(psoutput,str);
1109 sprintf(upos,"%d",begin + test.u - 6);
1110 strcat(psoutput,"_u_");
1111 strcat(psoutput,upos);
1112 strcat(psoutput,"_");
1113 strcat(psoutput,name);
1114 strcat(psoutput,".ps\0");
1115 PS_rna_plot_snoop_a(temp_seq, temp_struc, psoutput, NULL, NULL);
1120 /* free(psoutput); */
1123 free(test.structure);
1124 free(target_struct);
1139 snoopT snoopfold(const char *s1, const char *s2,
1140 const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
1142 const int half_stem, const int max_half_stem,
1143 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int fullStemEnergy) {
1144 /* int Eminj, Emin_l; */
1145 int i, j, l1, Emin=INF, i_min=0, j_min=0;
1151 folden** foldlist, **foldlist_XS;
1152 int Duplex_El, Duplex_Er;
1156 Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
1157 snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
1158 n1 = (int) strlen(s1);
1159 n2 = (int) strlen(s2);
1161 if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
1162 snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
1166 c = (int **) space(sizeof(int *) * (n1+1));
1167 r = (int **) space(sizeof(int *) * (n1+1));
1168 for (i=0; i<=n1; i++) {
1169 c[i] = (int *) space(sizeof(int) * (n2+1));
1170 r[i] = (int *) space(sizeof(int) * (n2+1));
1171 for(j=n2; j>-1; j--){
1176 encode_seqs(s1, s2);
1177 for (i=6; i<=n1-5; i++) {
1178 for (j=n2-min_d2; j>min_d1; j--) {
1179 int type, type2, E, k,l;
1180 type = pair[S1[i]][S2[j]];
1181 c[i][j] = (type ) ? P->DuplexInit : INF;
1183 if(/* pair[S1[i+1]][S2[j-1]] && */
1184 j < max_s1 && j > min_s1 &&
1185 j > n2 - max_s2 - max_half_stem &&
1186 j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
1188 max_k = MIN2(n2-min_s2,j+max_half_stem);
1189 min_k = MAX2(j+half_stem, n2-max_s2);
1194 /* if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour */
1195 if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){
1196 r[i][j]=MIN2(r[i][j], c[i-3][k+1]+temp->energy);
1198 /*else*/ if(pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/ ){
1199 r[i][j]=MIN2(r[i][j], c[i-4][k+1]+temp->energy);
1205 /* dangle 5'SIDE relative to the mRNA */
1207 *** c[i][j] += P->dangle5[type][SS1[i-1]];
1208 *** c[i][j] += P->dangle3[type][SS2[j+1]];
1209 *** if (type>2) c[i][j] += P->TerminalAU;
1211 c[i][j]+=E_ExtLoop(type, SS1[i-1] , SS2[j+1], P);
1212 for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
1213 for (l=j+1; l<=n2 ; l++) {
1214 if (i-k+l-j>2*MAXLOOP_L-2) break;
1215 if (abs(i-k-l+j) >= ASS ) continue;
1216 type2 = pair[S1[k]][S2[l]];
1217 if (!type2) continue;
1218 E = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1219 SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1220 c[i][j] = MIN2(c[i][j], c[k][l]+E+(i-k+l-j)*penalty);
1221 r[i][j] = MIN2(r[i][j], r[k][l]+E+(i-k+l-j)*penalty);
1226 *** if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];
1227 *** if (j>1) E += P->dangle5[rtype[type]][SS2[j-1]];
1228 *** f (type>2) E += P->TerminalAU;
1230 E+=E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
1232 Emin=E; i_min=i; j_min=j;
1237 printf("no target found under the constraints chosen\n");
1238 for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
1241 free(S1); free(S2); free(SS1); free(SS2);
1245 struc = snoop_backtrack(i_min, j_min,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D,
1246 &u, penalty, threshloop, threshLE, threshRE,threshDE, threshD,
1247 half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
1248 /* if (i_min<n1-5) i_min++; */
1249 /* if (j_min>1 ) j_min--; */
1250 l1 = strchr(struc, '&')-struc;
1254 mfe.Duplex_Er = (float) Duplex_Er/100;
1255 mfe.Duplex_El = (float) Duplex_El/100;
1256 mfe.Loop_D = (float) Loop_D/100;
1257 mfe.Loop_E = (float) Loop_E/100;
1258 mfe.energy = (float) Emin/100 ;
1259 mfe.fullStemEnergy = (float) fullStemEnergy/100;
1260 mfe.structure = struc;
1262 for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
1265 free(S1); free(S2); free(SS1); free(SS2);
1270 PRIVATE int snoopfold_XS_fill(const char *s1, const char *s2, const int **access_s1,
1271 const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
1273 const int half_stem, const int max_half_stem,
1274 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2) {
1275 /* int Eminj, Emin_l; */
1276 int i, j, Emin=INF, i_min=0, j_min=0;
1282 folden** foldlist, **foldlist_XS;
1283 int Duplex_El, Duplex_Er;
1287 Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
1288 snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
1289 n1 = (int) strlen(s1);
1290 n2 = (int) strlen(s2);
1292 if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
1293 snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
1297 c_fill = (int **) space(sizeof(int *) * (n1+1));
1298 r_fill = (int **) space(sizeof(int *) * (n1+1));
1299 for (i=0; i<=n1; i++) {
1300 c_fill[i] = (int *) space(sizeof(int) * (n2+1));
1301 r_fill[i] = (int *) space(sizeof(int) * (n2+1));
1302 for(j=n2; j>-1; j--){
1307 encode_seqs(s1, s2);
1311 for (i=6; i<=n1-5; i++) {
1312 di[1]=access_s1[5][i] - access_s1[4][i-1];
1313 di[2]=access_s1[5][i-1] - access_s1[4][i-2] + di[1];
1314 di[3]=access_s1[5][i-2] - access_s1[4][i-3] + di[2];
1315 di[4]=access_s1[5][i-3] - access_s1[4][i-4] + di[3];
1316 di[1]=MIN2(di[1],165);
1317 di[2]=MIN2(di[2],330);
1318 di[3]=MIN2(di[3],495);
1319 di[4]=MIN2(di[4],660);
1320 for (j=n2-min_d2; j>min_d1; j--) {
1321 int type, type2, E, k,l;
1322 type = pair[S1[i]][S2[j]];
1323 c_fill[i][j] = (type ) ? P->DuplexInit : INF;
1325 if(/* pair[S1[i+1]][S2[j-1]] && */
1326 j < max_s1 && j > min_s1 &&
1327 j > n2 - max_s2 - max_half_stem &&
1328 j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
1330 max_k = MIN2(n2-min_s2,j+max_half_stem);
1331 min_k = MAX2(j+half_stem, n2-max_s2);
1336 /* if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour */
1337 if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){
1338 r_fill[i][j]=MIN2(r_fill[i][j], c_fill[i-3][k+1]+temp->energy+ di[3]);
1340 /*else*/ if(pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/ ){
1341 r_fill[i][j]=MIN2(r_fill[i][j], c_fill[i-4][k+1]+temp->energy + di[4]);
1347 /* dangle 5'SIDE relative to the mRNA */
1349 *** c_fill[i][j] += P->dangle5[type][SS1[i-1]];
1350 *** c_fill[i][j] += P->dangle3[type][SS2[j+1]];
1351 *** if (type>2) c_fill[i][j] += P->TerminalAU;
1353 c_fill[i][j]+= E_ExtLoop(type, SS1[i-1], SS2[j+1],P);
1354 for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
1355 for (l=j+1; l<=n2 ; l++) {
1356 if (i-k+l-j>2*MAXLOOP_L-2) break;
1357 if (abs(i-k-l+j) >= ASS ) continue;
1358 type2 = pair[S1[k]][S2[l]];
1359 if (!type2) continue;
1360 E = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1361 SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1362 c_fill[i][j] = MIN2(c_fill[i][j], c_fill[k][l]+E+di[i-k]);
1363 r_fill[i][j] = MIN2(r_fill[i][j], r_fill[k][l]+E+di[i-k]);
1368 *** if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];
1369 *** if (j>1) E += P->dangle5[rtype[type]][SS2[j-1]];
1370 *** if (type>2) E += P->TerminalAU;
1372 E+= E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
1374 Emin=E; i_min=i; j_min=j;
1383 PUBLIC snoopT *snoop_subopt(const char *s1, const char *s2, int delta, int w,
1384 const int penalty, const int threshloop,
1385 const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
1386 const int distance, const int half_stem, const int max_half_stem,
1387 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int fullStemEnergy) {
1391 /* printf("%d %d\n", min_s2, max_s2); */
1392 int i,j, n1, n2, E, n_subopt=0, n_max;
1398 int Duplex_El, Duplex_Er, Loop_E;
1400 Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;
1404 subopt = (snoopT *) space(n_max*sizeof(snoopT));
1406 mfe = snoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE,threshD,
1407 half_stem, max_half_stem,
1408 min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, fullStemEnergy);
1417 thresh = MIN2((int) ((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E)*100+0.1 + 410) + delta, threshTE );
1418 /* subopt[n_subopt++]=mfe; */
1419 free(mfe.structure);
1421 n1 = strlen(s1); n2=strlen(s2);
1422 for (i=n1; i>0; i--) {
1423 for (j=1; j<=n2; j++) {
1425 type = pair[S2[j]][S1[i]];
1426 if (!type) continue;
1429 *** if (i<n1) Ed += P->dangle3[type][SS1[i+1]];
1430 *** if (j>1) Ed += P->dangle5[type][SS2[j-1]];
1431 *** if (type>2) Ed += P->TerminalAU;
1433 Ed+= E_ExtLoop(type, (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
1434 if (Ed>thresh) continue;
1435 /* too keep output small, remove hits that are dominated by a
1436 better one close (w) by. For simplicity we do test without
1437 adding dangles, which is slightly inaccurate.
1440 /* for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n1)) && type; ii++) { */
1441 /* for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n2); jj++) */
1442 /* if (r[ii][jj]<E) {type=0; break;} */
1444 if (!type) continue;
1446 struc = snoop_backtrack(i,j,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D, &u, penalty, threshloop,threshLE,threshRE,threshDE,threshD,
1447 half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
1448 if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
1449 (Duplex_Er + Duplex_El) > threshDE ||
1450 (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
1451 (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
1452 /* printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n" */
1453 /* " Duplex_Er + Duplex_El %d threshDE %d \n" */
1454 /* " Duplex_Er + Duplex_El + Loop_E %d threshTE %d \n" */
1455 /* " Duplex_Er + Duplex_El + Loop_E + Loop_D %d threshSE %d \n", */
1456 /* Duplex_Er , threshRE , Duplex_El ,threshLE, */
1457 /* Duplex_Er + Duplex_El, threshDE, */
1458 /* Duplex_Er + Duplex_El+ Loop_E , threshTE, */
1459 /* Duplex_Er + Duplex_El+ Loop_E + Loop_D, threshSE); */
1469 if (n_subopt+1>=n_max) {
1471 subopt = (snoopT *) xrealloc(subopt, n_max*sizeof(snoopT));
1473 subopt[n_subopt].i = i-5;
1474 subopt[n_subopt].j = j-5;
1475 subopt[n_subopt].u = u-5;
1476 subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
1477 subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
1478 subopt[n_subopt].Loop_E = Loop_E * 0.01;
1479 subopt[n_subopt].Loop_D = Loop_D * 0.01;
1480 subopt[n_subopt].energy = (Duplex_Er +Duplex_El + Loop_E + Loop_D + 410) * 0.01 ;
1481 subopt[n_subopt].fullStemEnergy = (float) fullStemEnergy * 0.01;
1482 subopt[n_subopt++].structure = struc;
1484 Duplex_Er=0; Duplex_El=0; Loop_E=0; Loop_D=0;u=0;
1488 for (i=0; i<=n1; i++) {free(c[i]);free(r[i]);}
1490 free(S1); free(S2); free(SS1); free(SS2);
1493 if (snoop_subopt_sorted) qsort(subopt, n_subopt, sizeof(snoopT), compare);
1494 subopt[n_subopt].i =0;
1495 subopt[n_subopt].j =0;
1496 subopt[n_subopt].structure = NULL;
1500 PUBLIC void snoop_subopt_XS(const char *s1, const char *s2, const int **access_s1, int delta, int w,
1501 const int penalty, const int threshloop,
1502 const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
1503 const int distance, const int half_stem, const int max_half_stem,
1504 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int alignment_length, const char *name, const int fullStemEnergy) {
1508 /* printf("%d %d\n", min_s2, max_s2); */
1515 int Duplex_El, Duplex_Er, Loop_E;
1517 Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;
1522 int Emin = snoopfold_XS_fill(s1, s2, access_s1,penalty, threshloop, threshLE, threshRE, threshDE,threshD,
1523 half_stem, max_half_stem,
1524 min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
1528 thresh = MIN2(-100, threshTE +alignment_length*30);
1529 /* n1=strlen(s1); */
1530 /* n2=strlen(s2); */
1534 S1_fill = (short*)space(sizeof(short)*(n3+2));
1535 S2_fill = (short*)space(sizeof(short)*(n4+2));
1536 SS1_fill = (short*)space(sizeof(short)*(n3+1));
1537 SS2_fill = (short*)space(sizeof(short)*(n4+1));
1538 memcpy(S1_fill, S1, sizeof(short)*n3+2);
1539 memcpy(S2_fill, S2, sizeof(short)*n4+2);
1540 memcpy(SS1_fill, SS1, sizeof(short)*n3+1);
1541 memcpy(SS2_fill, SS2, sizeof(short)*n4+1);
1542 free(S1);free(S2);free(SS1);free(SS2);
1544 for (i=n3-5; i>0; i--) {
1545 for (j=1; j<=n4; j++) {
1547 type = pair[S2_fill[j]][S1_fill[i]];
1548 if (!type) continue;
1549 E = Ed = r_fill[i][j];
1551 ***if (i<n3) Ed += P->dangle3[type][SS1_fill[i+1]];
1552 ***if (j>1) Ed += P->dangle5[type][SS2_fill[j-1]];
1553 ***if (type>2) Ed += P->TerminalAU;
1555 Ed+=E_ExtLoop(type, (j > 1) ? SS2[j-1] : -1, (i<n3) ? SS1[i+1] : -1, P);
1556 if (Ed>thresh) continue;
1558 /* to keep output small, remove hits that are dominated by a
1559 better one close (w) by. For simplicity we do test without
1560 adding dangles, which is slightly inaccurate.
1563 /* for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n3-5)) && type; ii++) { */
1564 /* for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n4-5); jj++) */
1565 /* if (r_fill[ii][jj]<E) {type=0; break;} */
1568 if (!type) continue;
1569 int begin=MAX2(5, i-alignment_length);
1570 int end =MIN2(n3-5, i-1);
1571 char *s3 = (char*) space(sizeof(char)*(end-begin+2)+5);
1572 strncpy(s3, (s1+begin), end - begin +1);
1573 strcat(s3,"NNNNN\0");
1574 int n5 = strlen(s3);
1575 snoopT test = snoopfold_XS(s3, s2, access_s1, i, j ,penalty,
1576 threshloop, threshLE, threshRE,
1577 threshDE, threshD, half_stem,
1578 max_half_stem, min_s2, max_s2, min_s1,
1579 max_s1, min_d1, min_d2,fullStemEnergy);
1580 if(test.energy==INF){
1584 if( test.Duplex_El > threshLE * 0.01 ||test.Duplex_Er > threshRE * 0.01 ||
1585 test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
1586 (test.Duplex_Er + test.Duplex_El + test.Loop_E) > threshTE*0.01 || (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE*0.01)
1588 free(test.structure);free(s3);
1592 s4 = (char*) space(sizeof(char) *(n4-9));
1593 strncpy(s4, s2+5, n4-10);
1596 char *s5 = space(sizeof(char) * n5-test.i+2-5);
1597 strncpy(s5,s3+test.i-1,n5-test.i+1-5);
1598 s5[n5-test.i+1-5]='\0';
1599 float dE = ((float) (access_s1[n5-test.i+1-5][i]))*0.01;
1600 printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + %5.2f + 4.10) (%5.2f)\n%s&%s\n" ,
1601 test.structure, i - (n5 - test.i) ,i - 5, i - (n5 - test.u ),
1602 j-5, j-5 + (strrchr(test.structure,'>') - strchr(test.structure,'>')),
1603 test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10+dE, test.Duplex_El,
1604 test.Duplex_Er, test.Loop_E, test.Loop_D,dE , test.fullStemEnergy, s5,s4);
1606 int begin_t, end_t, begin_q, end_q, and, pipe,k;
1611 end_t=n5-test.i+ 1-5;
1613 pipe=test.u -test.i + 1;
1614 cut_point = end_t +1 ;
1615 char *catseq, *catstruct;/* *fname; */
1616 catseq = (char*) space(n5 + end_q -begin_q +2);
1617 catstruct = (char*) space(n5 + end_q -begin_q +2);
1619 strncpy(catstruct, test.structure, end_t);
1621 strncat(catstruct, test.structure+end_t+1, end_q-begin_q+1);
1622 catstruct[end_t - begin_t + end_q-begin_q+2]='\0';
1623 catseq[end_t - begin_t + end_q-begin_q+2]='\0';
1624 int *relative_access;
1625 relative_access = space(sizeof(int)*strlen(s5));
1626 relative_access[0] = access_s1[1][i - (n5 - test.i) + 5];
1627 for(k=1;k<strlen(s5);k++){
1628 relative_access[k] = access_s1[k+1][i - (n5 - test.i) + k + 5] - access_s1[k][i - (n5 - test.i) + k + 4];
1630 char str[16];char upos[16];
1631 strcpy(psoutput,"sno_XS_");
1632 sprintf(str,"%d",count);
1633 strcat(psoutput,str);
1634 sprintf(upos,"%d",i - (n5 - test.u ));
1635 strcat(psoutput,"_u_");
1636 strcat(psoutput,upos);
1637 strcat(psoutput,"_");
1638 strcat(psoutput,name);
1639 strcat(psoutput,".ps\0");
1640 PS_rna_plot_snoop_a(catseq, catstruct, psoutput,relative_access,NULL);
1641 free(catseq);free(catstruct);free(relative_access);
1644 free(s3);free(s4);free(s5);free(test.structure);
1647 for (i=0; i<=n3; i++) {free(c_fill[i]);free(r_fill[i]);}
1648 free(c_fill);free(r_fill);
1649 free(S1_fill); free(S2_fill); free(SS1_fill); free(SS2_fill);
1656 PRIVATE char *snoop_backtrack(int i, int j, const char* snoseq,
1657 int *Duplex_El, int *Duplex_Er,
1658 int *Loop_E, int *Loop_D, int *u,
1659 const int penalty, const int threshloop,
1660 const int threshLE, const int threshRE, const int threshDE, const int threshD,
1661 const int half_stem, const int max_half_stem,
1662 const int min_s2, const int max_s2, const int min_s1,
1663 const int max_s1, const int min_d1, const int min_d2) {
1664 /* backtrack structure going backwards from i, and forwards from j
1665 return structure in bracket notation with & as separator */
1666 int k, l, type, type2, E, traced, i0, j0;
1667 int traced_r=0; /* flag for following backtrack in c or r */
1668 char *st1, *st2, *struc;
1671 st1 = (char *) space(sizeof(char)*(n1+1));
1672 st2 = (char *) space(sizeof(char)*(n2+1));
1676 folden **foldlist, **foldlist_XS;
1677 type=pair[S1[i]][S2[j]];
1678 snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
1681 *** if (i<n1) *Duplex_Er += P->dangle3[rtype[type]][SS1[i+1]];
1682 *** if (j>1) *Duplex_Er += P->dangle5[rtype[type]][SS2[j-1]];
1683 *** if (type>2) *Duplex_Er += P->TerminalAU;
1685 *Duplex_Er += E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
1686 while (i>0 && j<=n2-min_d2 ) {
1688 E = r[i][j]; traced=0;
1691 type = pair[S1[i]][S2[j]];
1692 if (!type) nrerror("backtrack failed in fold duplex r");
1693 for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
1694 for (l=j+1; l<=n2 ; l++) {
1696 if (i-k+l-j>2*MAXLOOP_L-2) break;
1697 if (abs(i-k-l+j) >= ASS) continue;
1699 type2 = pair[S1[k]][S2[l]];
1700 if (!type2) continue;
1701 LE = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1702 SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1703 if (E == r[k][l]+LE+(i-k+l-j)*penalty) {
1713 if(/* pair[S1[i+1]][S2[j-1]] && */
1714 j < max_s1 && j > min_s1 &&
1715 j > n2 - max_s2 - max_half_stem &&
1716 j < n2 -min_s2 -half_stem &&
1719 max_k = MIN2(n2-min_s2,j+max_half_stem+1);
1720 min_k = MAX2(j+half_stem+1, n2-max_s2);
1725 if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){ /* introduce structure from RNAfold */
1726 if(E==c[i-3][k+1]+temp->energy){
1727 *Loop_E=temp->energy;
1731 /* int fix_ij=indx[k-1+1]+j+1; */
1732 for(a=0; a< MISMATCH ;a++){
1733 for(b=0; b< MISMATCH ; b++){
1734 int ij=indx[k-1-a+1]+j+1+b;
1735 if(cLoop[ij]==temp->energy) {
1736 struc_loop=snobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1);
1747 /*else*/ if (pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/){ /* introduce structure from RNAfold */
1748 if(E==c[i-4][k+1]+temp->energy){
1749 *Loop_E=temp->energy;
1753 /* int fix_ij=indx[k-1+1]+j+1; */
1754 for(a=0; a< MISMATCH ;a++){
1755 for(b=0; b< MISMATCH ; b++){
1756 int ij=indx[k-1-a+1]+j+1+b;
1757 if(cLoop[ij]==temp->energy) {
1758 struc_loop=snobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1);
1770 } /* while temp-> next */
1775 E = c[i][j]; traced=0;
1778 type = pair[S1[i]][S2[j]];
1779 if (!type) nrerror("backtrack failed in fold duplex c");
1780 for (k=i-1; (i-k)<MAXLOOP_L; k--) {
1781 for (l=j+1; l<=n2; l++) {
1783 if (i-k+l-j>2*MAXLOOP_L-2) break;
1784 if (abs(i-k-l+j) >= ASS) continue;
1785 type2 = pair[S1[k]][S2[l]];
1786 if (!type2) continue;
1787 LE = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
1788 SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
1789 if (E == c[k][l]+LE+(i-k+l-j)*penalty) {
1801 correction = E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, P);
1803 *Duplex_El+=correction;
1805 *** if (i>1) {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];}
1806 *** if (j<n2) {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];}
1807 *** if (type>2) {E -= P->TerminalAU; *Duplex_El +=P->TerminalAU;}
1809 if (E != P->DuplexInit) {
1810 nrerror("backtrack failed in fold duplex end");
1815 /* if (j<n2) j++; */
1816 /* struc = (char *) space(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
1817 struc = (char *) space(i0-i+1+n2-1+1+2); /* declare final duplex structure */
1819 struc2 = (char *) space(n2+1);
1820 /* char * struct_const; */
1821 for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
1822 /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
1823 /* char * struct_const; */
1824 /* struct_const = (char *) space(sizeof(char)*(n2+1)); */
1825 for (k=1; k<=n2; k++) {
1826 if (!st2[k-1]) st2[k-1] = struc_loop[k-1];/* '.'; */
1827 struc2[k-1] = st2[k-1];/* '.'; */
1828 /* if (k>=j0 && k<=j){ */
1829 /* struct_const[k-1]='x'; */
1832 /* if(k<j0) {struct_const[k-1]='<';} */
1833 /* if(k>j) {struct_const[k-1]='>';} */
1836 char duplexseq_1[j0];
1837 char duplexseq_2[n2-j+2];
1839 strncpy(duplexseq_1, snoseq, j0-1);
1840 strcpy(duplexseq_2, snoseq + j);
1841 duplexseq_1[j0-1]='\0';
1842 duplexseq_2[n2-j+1]='\0';
1844 temp=duplexfold(duplexseq_1, duplexseq_2);
1845 *Loop_D = MIN2(0,-410 + (int) 100 * temp.energy);
1847 int l1,ibegin, iend, jbegin, jend;
1848 l1=strchr(temp.structure, '&')-temp.structure;
1852 jend =temp.j+strlen(temp.structure)-l1-2-1;
1853 for(k=ibegin+1; k<=iend+1; k++){
1854 struc2[k-1]=temp.structure[k-ibegin-1];
1856 for(k=jbegin+j; k<=jend+j; k++){
1857 struc2[k-1]=temp.structure[l1+k-j-jbegin+1];
1860 free(temp.structure);
1862 strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&");
1863 /* strcat(struc, st2); */
1864 strncat(struc, struc2+5, strlen(struc2)-10);
1867 free(st1); free(st2);
1868 /* free_arrays(); */
1872 void Lsnoop_subopt_list_XS(const char *s1, const char *s2, const int **access_s1, int delta, int w,
1873 const int penalty, const int threshloop,
1874 const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
1876 const int half_stem, const int max_half_stem,
1877 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int alignment_length, const char *name, const int fullStemEnergy)
1880 int min_colonne=INF;
1884 /* int nsubopt=10; */
1885 n1 = (int) strlen(s1);
1886 n2 = (int) strlen(s2);
1891 position = (int*) space((n1+3)*sizeof(int));
1892 position_j = (int*) space((n1+3)*sizeof(int));
1894 /* int Eminj, Emin_l; */
1895 int i, j;/* l1, Emin=INF, i_min=0, j_min=0; */
1901 folden **foldlist, **foldlist_XS;
1902 int Duplex_El, Duplex_Er;
1906 Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
1907 snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
1908 if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
1909 snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
1913 lpair = (int **) space(sizeof(int *) * (6));
1914 lc = (int **) space(sizeof(int *) * (6));
1915 lr = (int **) space(sizeof(int *) * (6));
1916 for (i=0; i<6; i++) {
1917 lc[i] = (int *) space(sizeof(int) * (n2+1));
1918 lr[i] = (int *) space(sizeof(int) * (n2+1));
1919 lpair[i] = (int *) space(sizeof(int) * (n2+1));
1920 for(j=n2; j>-1; j--){
1926 encode_seqs(s1, s2);
1927 int lim_maxj=n2-min_d2 ;
1928 int lim_minj=min_d1;
1930 for (i=5; i<=lim_maxi; i++) {
1936 int di1, di2, di3, di4;
1937 di1 = access_s1[5][i] - access_s1[4][i-1];
1938 di2 =access_s1[5][i-1] - access_s1[4][i-2] + di1;
1939 di3 = access_s1[5][i-2] - access_s1[4][i-3] + di2;
1940 di4 = access_s1[5][i-3] - access_s1[4][i-4] + di3;
1945 for (j=lim_maxj; j>lim_minj; j--) {
1946 int type, type2;/* E, k,l; */
1947 type = pair[S1[i]][S2[j]];
1948 lpair[idx][j] = type;
1949 lc[idx][j] = (type) ? P->DuplexInit + access_s1[1][i] : INF;
1952 if( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
1953 j < max_s1 && j > min_s1 &&
1954 j > n2 - max_s2 - max_half_stem &&
1955 j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
1957 max_k = MIN2(n2-min_s2,j+max_half_stem+1);
1958 min_k = MAX2(j+half_stem+1, n2-max_s2);
1963 /* if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour */
1964 if(lpair[idx_3][k+1] && lc[idx_3][k+1] /*+di3*/ < 411 /*&& lpair[idx_4][k+2]*/){ /* remove second condition */
1965 lr[idx][j]=MIN2(lr[idx][j], di3 + lc[idx_3][k+1]+temp->energy);/*--NU--*/
1967 /*else*/ if(lpair[idx_4][k+1] && /*di4 +*/ lc[idx_4][k+1] < 411 ){/*--NUN--*/ /* remove second condition */
1968 lr[idx][j]=MIN2(lr[idx][j], di4 + lc[idx_4][k+1]+temp->energy);
1974 /* dangle 5'SIDE relative to the mRNA */
1976 *** lc[idx][j] += P->dangle5[type][SS1[i-1]];
1977 *** lc[idx][j] += P->dangle3[type][SS2[j+1]];
1978 *** if (type>2) lc[idx][j] += P->TerminalAU;
1980 lc[idx][j]+=E_ExtLoop(type, SS1[i-1] , SS2[j+1] , P);
1981 /* if(j<n2 && i>1){ */
1982 /* type2=pair[S1[i-1]][S2[j+1]]; */
1983 type2=lpair[idx_1][j+1];
1985 lc[idx][j]=MIN2(lc[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+di1, lc[idx][j]);
1986 lr[idx][j]=MIN2(lr[idx_1][j+1]+E_IntLoop(0,0,type2, rtype[type],SS1[i], SS2[j], SS1[i-1], SS2[j+1], P)+di1, lr[idx][j]);
1988 type2=lpair[idx_2][j+2];
1990 lc[idx][j]=MIN2(lc[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P)+di2, lc[idx][j]);
1991 lr[idx][j]=MIN2(lr[idx_2][j+2]+E_IntLoop(1,1,type2, rtype[type],SS1[i-1], SS2[j+1], SS1[i-1], SS2[j+1], P)+di2, lr[idx][j]);
1994 type2 =lpair[idx_3][j+3];
1996 lc[idx][j]=MIN2(lc[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+di3,lc[idx][j]);
1997 lr[idx][j]=MIN2(lr[idx_3][j+3]+E_IntLoop(2,2,type2, rtype[type],SS1[i-2], SS2[j+2], SS1[i-1], SS2[j+1], P)+di3,lr[idx][j]);
2003 bla=lr[idx][j] + E_ExtLoop(rtype[type], SS2[j-1], SS1[i+1] , P);
2005 *** (type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]];
2007 min_colonne=MIN2(bla, min_colonne);
2008 if(temp2>min_colonne){
2012 position[i]=min_colonne;
2013 if(max>=min_colonne){
2016 max_pos_j=min_j_colonne;
2018 position_j[i]=min_j_colonne;
2021 free(S1); free(S2); free(SS1); free(SS2);
2023 if(max<threshTE + 30*alignment_length){
2024 find_max_snoop_XS(s1, s2, access_s1,max,alignment_length, position, position_j,
2025 delta, distance, penalty, threshloop,
2026 threshLE, threshRE, threshDE, threshTE, threshSE, threshD,
2027 half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name,fullStemEnergy);
2029 for (i=1; i<6; i++) {free(lc[i]);free(lr[i]);free(lpair[i]);}
2030 free(lc[0]);free(lr[0]);free(lpair[0]);
2031 free(lc);free(lr);free(lpair);
2032 free(position);free(position_j);
2036 PRIVATE void find_max_snoop_XS(const char *s1, const char *s2, const int **access_s1,
2037 const int max, const int alignment_length,
2038 const int* position, const int* position_j, const int delta,
2039 const int distance, const int penalty, const int threshloop, const int threshLE, const int threshRE,
2040 const int threshDE, const int threshTE, const int threshSE, const int threshD,
2041 const int half_stem, const int max_half_stem, const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const char *name, const int fullStemEnergy){
2047 int threshold = MIN2(threshTE + alignment_length*30, -100);
2048 /* printf("threshTE %d max %d\n", threshTE, max); */
2049 /* #pragma omp parallel for */
2050 /* for(pos=n1+1;pos>distance;pos--){ */
2053 if(position[pos]<(threshold)){
2055 search_range=distance+1;
2056 while(--search_range){
2057 if(position[pos-search_range]<=position[pos-temp_min]){
2058 temp_min=search_range;
2062 max_pos_j=position_j[pos];
2063 int begin=MAX2(5, pos-alignment_length);
2064 int end =MIN2(n3-5, pos-1);
2065 char *s3 = (char*) space(sizeof(char)*(end-begin+2)+5);
2066 strncpy(s3, (s1+begin), end - begin +1);
2067 strcat(s3,"NNNNN\0");
2069 int n5 = strlen(s3);
2071 test = snoopfold_XS(s3, s2, access_s1, pos, max_pos_j ,penalty,
2072 threshloop, threshLE, threshRE,
2073 threshDE, threshD, half_stem,
2074 max_half_stem, min_s2, max_s2, min_s1,
2075 max_s1, min_d1, min_d2, fullStemEnergy);
2076 if(test.energy==INF){
2080 if( test.Duplex_El > threshLE * 0.01 ||test.Duplex_Er > threshRE * 0.01 ||
2081 test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
2082 (test.Duplex_Er + test.Duplex_El + test.Loop_E) > threshTE*0.01 || (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE*0.01) {
2083 free(test.structure);free(s3);
2088 s4 = (char*) space(sizeof(char) *(n4-9));
2089 strncpy(s4, s2+5, n4-10);
2092 char *s5 = space(sizeof(char) * n5-test.i+2-5);
2093 strncpy(s5,s3+test.i-1,n5-test.i+1-5);
2094 s5[n5-test.i+1-5]='\0';
2095 float dE = ((float) (access_s1[n5-test.i+1-5][pos]))*0.01;
2096 printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + %5.2f + 4.10) (%5.2f)\n%s&%s\n" ,
2097 test.structure, pos - (n5 - test.i) ,pos - 5, pos - (n5 - test.u ),
2098 max_pos_j-5, max_pos_j-5 + (strrchr(test.structure,'>') - strchr(test.structure,'>')),
2099 test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10+dE, test.Duplex_El,
2100 test.Duplex_Er, test.Loop_E, test.Loop_D,dE ,test.fullStemEnergy, s5,s4);
2102 int begin_t, end_t, begin_q, end_q, and, pipe, i;
2107 end_t=n5-test.i+ 1-5;
2109 pipe=test.u -test.i + 1;
2110 cut_point = end_t +1 ;
2111 char *catseq, *catstruct;/* *fname; */
2112 catseq = (char*) space(n5 + end_q -begin_q +2);
2113 catstruct = (char*) space(n5 + end_q -begin_q +2);
2115 strncpy(catstruct, test.structure, end_t);
2117 strncat(catstruct, test.structure+end_t+1, end_q-begin_q+1);
2118 catstruct[end_t - begin_t + end_q-begin_q+2]='\0';
2119 catseq[end_t - begin_t + end_q-begin_q+2]='\0';
2120 int *relative_access;
2121 relative_access = space(sizeof(int)*strlen(s5));
2123 relative_access[0] = access_s1[1][pos - (n5 - test.i) + 5];
2124 for(i=1;i<strlen(s5);i++){
2125 relative_access[i] = access_s1[i+1][pos - (n5 - test.i) + i + 5] - access_s1[i][pos - (n5 - test.i) + i + 4];
2129 strcpy(psoutput,"sno_XS_");
2130 sprintf(str,"%d",count);
2131 strcat(psoutput,str);
2132 sprintf(upos,"%d",pos - (n5 - test.u ));
2133 strcat(psoutput,"_u_");
2134 strcat(psoutput,upos);
2135 strcat(psoutput,"_");
2136 strcat(psoutput,name);
2137 strcat(psoutput,".ps\0");
2138 PS_rna_plot_snoop_a(catseq, catstruct, psoutput,relative_access,NULL);
2139 free(catseq);free(catstruct);free(relative_access);
2142 free(s3);free(s4);free(s5);free(test.structure);
2147 snoopT snoopfold_XS(const char *s1, const char *s2, const int **access_s1, const int pos_i, const int pos_j,
2148 const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
2150 const int half_stem, const int max_half_stem,
2151 const int min_s2, const int max_s2, const int min_s1, const int max_s1, const int min_d1, const int min_d2, const int fullStemEnergy) {
2152 /* int Eminj, Emin_l; */
2153 int a,b,i, j, Emin=INF, a_min=0, b_min=0;
2159 folden** foldlist, **foldlist_XS;
2160 int Duplex_El, Duplex_Er;
2165 Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
2166 snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
2167 n1 = (int) strlen(s1);
2168 n2 = (int) strlen(s2);
2170 if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
2171 snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
2175 c = (int **) space(sizeof(int *) * (n1+1));
2176 r = (int **) space(sizeof(int *) * (n1+1));
2177 for (i=0; i<=n1; i++) {
2178 c[i] = (int *) space(sizeof(int) * (n2+1));
2179 r[i] = (int *) space(sizeof(int) * (n2+1));
2180 for(j=n2; j>-1; j--){
2185 encode_seqs(s1, s2);
2188 /* printf("tar: %s\nsno: %s\n ", s1, s2); */
2189 /* printf("pos_i %d pos_j %d\n", pos_i, pos_j); */
2190 /* printf("type %d n1 %d n2 %d S1[n1] %d S2[n2] %d", pair[S1[i]][S2[j]], n1, n2, S1[n1], S2[n2]); */
2191 int type, type2, E, p,q;
2192 r[i][j] = P->DuplexInit;
2193 /* r[i][j] += P->dangle3[rtype[type]][SS1[i+1]] + P->dangle5[rtype[type]][SS2[j-1]]; */
2195 if(pair[S1[i]][S2[j]]>2) r[i][j]+=P->TerminalAU;
2196 for(a=i-1; a>0;a--){ /* i-1 */
2198 for(b=j+1; b<=n2-min_d2;b++){ /* j+1 */
2200 type = pair[S1[a]][S2[b]];
2204 temp=foldlist_XS[b-1];
2207 if(pair[S1[a+3]][S2[k-1]] && k< max_s1 && k > min_s1 && k > n2 - max_s2 - max_half_stem && k < n2 -min_s2 -half_stem /*&& r[a+3][k-1] + access_s1[i-(a+3)+1][pos_i] < 411*/) { /* remove last condition last condition is to check if the interaction is stable enough */
2208 c[a][b]=MIN2(c[a][b], r[a+3][k-1]+temp->energy);
2215 temp=foldlist_XS[b-1];
2218 if(pair[S1[a+4]][S2[k-1]] && k< max_s1 && k > min_s1 && k > n2 - max_s2 - max_half_stem && k < n2 -min_s2 -half_stem /*&& r[a+4][k-1] + access_s1[i-(a+4)+1][pos_i] < 411 */ ) { /* remove last condition */
2219 c[a][b]=MIN2(c[a][b], r[a+4][k-1]+temp->energy);
2224 for(p=a+1; p<n1 && (p-a) < MAXLOOP_L; p++) { /* p < n1 */
2225 for (q=b-1; q > 1 ; q--) { /* q > 1 */
2226 if (p-a+b-q>2*MAXLOOP_L-2) break;
2227 if (abs((p-a)-(b-q)) >= ASS ) continue;
2228 type2 = pair[S1[p]][S2[q]];
2229 if (!type2) continue;
2230 E = E_IntLoop(p-a-1, b-q-1, type2, rtype[type],SS1[a+1], SS2[b-1], SS1[p-1], SS2[q+1],P);
2231 c[a][b] = MIN2(c[a][b], c[p][q]+E);
2232 r[a][b] = MIN2(r[a][b], r[p][q]+E);
2236 if (type>2) E += P->TerminalAU;
2237 /* E +=P->dangle5[rtype[type]][SS1[i+1]]; */
2238 /* E +=P->dangle5[rtype[type]][SS2[j-1]]; */
2239 E+=access_s1[i-a+1][pos_i];
2241 Emin=E; a_min=a; b_min=b;
2246 printf("no target found under the constraints chosen\n");
2247 for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
2250 free(S1); free(S2); free(SS1); free(SS2);
2254 type2=pair[S1[a_min]][S2[b_min]];
2255 if(type2>2) Emin +=P->TerminalAU;
2256 mfe.energy = ((float) (Emin))/100;
2257 struc = snoop_backtrack_XS(a_min, b_min,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D,
2258 &u, penalty, threshloop, threshLE, threshRE,threshDE, threshD,
2259 half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
2264 mfe.Duplex_Er = (float) Duplex_Er/100;
2265 mfe.Duplex_El = (float) Duplex_El/100;
2266 mfe.Loop_D = (float) Loop_D/100;
2267 mfe.Loop_E = (float) Loop_E/100;
2268 mfe.energy = (float) Emin/100 ;
2269 mfe.fullStemEnergy = (float) fullStemEnergy/100;
2270 mfe.structure = struc;
2274 PRIVATE char *snoop_backtrack_XS(int i, int j, const char* snoseq,
2275 int *Duplex_El, int *Duplex_Er,
2276 int *Loop_E, int *Loop_D, int *u,
2277 const int penalty, const int threshloop,
2278 const int threshLE, const int threshRE, const int threshDE, const int threshD,
2279 const int half_stem, const int max_half_stem,
2280 const int min_s2, const int max_s2, const int min_s1,
2281 const int max_s1, const int min_d1, const int min_d2) {
2282 /* backtrack structure going backwards from i, and forwards from j
2283 return structure in bracket notation with & as separator */
2284 int k, l, type, type2, E, traced, i0, j0;
2285 int traced_c=0; /* flag for following backtrack in c or r */
2286 char *st1, *st2, *struc;
2289 st1 = (char *) space(sizeof(char)*(n1+1));
2290 st2 = (char *) space(sizeof(char)*(n2+1));
2294 folden **foldlist, **foldlist_XS;
2295 type=pair[S1[i]][S2[j]];
2296 snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
2298 /* i0=MAX2(i,1); j0=MIN2(j+1,n2); */
2299 while (i<=n1 && j>=1 ) {
2301 E = c[i][j]; traced=0;
2304 type = pair[S1[i]][S2[j]];
2305 if (!type) nrerror("backtrack failed in fold duplex c");
2306 for (k=i+1; k>0 && (k-i)<MAXLOOP_L; k++) {
2307 for (l=j-1; l>=1 ; l--) {
2309 if (k-i+j-l>2*MAXLOOP_L-2) break;
2310 if (abs(k-i-j+l) >= ASS) continue;
2311 type2 = pair[S1[k]][S2[l]];
2312 if (!type2) continue;
2313 LE = E_IntLoop(k-i-1, j-l-1, type2, rtype[type],
2314 SS1[i+1], SS2[j-1], SS1[k-1], SS2[l+1],P);
2315 if (E == c[k][l]+LE) {
2327 temp=foldlist_XS[j-1];
2330 if(pair[S1[i+3]][S2[k-1]] && k< max_s1 && k > min_s1 && k > n2 - max_s2 - max_half_stem && k < n2 -min_s2 -half_stem ) {
2331 if(E==r[i+3][k-1]+temp->energy){
2332 *Loop_E=temp->energy;
2337 for(a=0; a< MISMATCH ;a++){
2338 for(b=0; b< MISMATCH ; b++){
2339 int ij=indx[j-1-a]+k+b;
2340 if(cLoop[ij]==temp->energy) {
2341 struc_loop=snobacktrack_fold_from_pair(snoseq, k+b, j-1-a);
2355 if (S1[i+2]==4){ /* introduce structure from RNAfold */
2357 temp=foldlist_XS[j-1];
2360 if(pair[S1[i+4]][S2[k-1]] && k< max_s1 && k > min_s1 && k > n2 - max_s2 - max_half_stem && k < n2 -min_s2 -half_stem ) {
2361 if(E==r[i+4][k-1]+temp->energy){
2362 *Loop_E=temp->energy;
2364 st1[i+1]=st1[i+3]='.';
2367 for(a=0; a< MISMATCH ;a++){
2368 for(b=0; b< MISMATCH ; b++){
2369 int ij=indx[j-1-a]+k+b;
2370 if(cLoop[ij]==temp->energy) {
2371 struc_loop=snobacktrack_fold_from_pair(snoseq, k+b, j-a-1);
2388 E = r[i][j]; traced=0;
2391 type = pair[S1[i]][S2[j]];
2392 if (!type) nrerror("backtrack failed in fold duplex r");
2393 for (k=i+1; k>0 && (k-i)<MAXLOOP_L; k++) {
2394 for (l=j-1; l>=1 ; l--) {
2396 if (k-i+j-l>2*MAXLOOP_L-2) break;
2397 if (abs(k-i-j+l) >= ASS) continue;
2398 type2 = pair[S1[k]][S2[l]];
2399 if (!type2) continue;
2400 LE = E_IntLoop(k-i-1, j-l-1, type2, rtype[type],
2401 SS1[i+1], SS2[j-1], SS1[k-1], SS2[l+1],P);
2402 if (E == r[k][l]+LE) {
2413 /* if (i>1) {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];} */
2414 /* if (j<n2) {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];} */
2415 if (type>2) {E -= P->TerminalAU; *Duplex_Er +=P->TerminalAU;}
2416 if (E != P->DuplexInit) {
2417 nrerror("backtrack failed in fold duplex end");
2423 /* struc = (char *) space(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
2424 struc = (char *) space(i-i0+1+n2); /* declare final duplex structure */
2426 struc2 = (char *) space(n2+1);
2427 /* char * struct_const; */
2429 for (k=MIN2(i0,1); k<=i; k++) if (!st1[k-1]) st1[k-1] = '.';
2430 /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
2431 /* char * struct_const; */
2432 /* struct_const = (char *) space(sizeof(char)*(n2+1)); */
2433 for (k=1; k<=n2; k++) {
2434 if (!st2[k-1]) st2[k-1] = struc_loop[k-1];/* '.'; */
2435 struc2[k-1] = st2[k-1];/* '.'; */
2436 /* if (k>=j0 && k<=j){ */
2437 /* struct_const[k-1]='x'; */
2440 /* if(k<j0) {struct_const[k-1]='<';} */
2441 /* if(k>j) {struct_const[k-1]='>';} */
2444 char duplexseq_1[j];
2445 char duplexseq_2[n2-j0+2];
2447 strncpy(duplexseq_1, snoseq, j-1);
2448 strcpy(duplexseq_2, snoseq + j0);
2449 duplexseq_1[j-1]='\0';
2450 duplexseq_2[n2-j0+1]='\0';
2452 temp=duplexfold(duplexseq_1, duplexseq_2);
2453 *Loop_D = MIN2(0,-410 + (int) 100 * temp.energy);
2455 int l1,ibegin, iend, jbegin, jend;
2456 l1=strchr(temp.structure, '&')-temp.structure;
2460 jend =temp.j+strlen(temp.structure)-l1-2-1;
2461 for(k=ibegin+1; k<=iend+1; k++){
2462 struc2[k-1]=temp.structure[k-ibegin-1];
2464 for(k=jbegin+j0; k<=jend+j0; k++){
2465 struc2[k-1]=temp.structure[l1+k-j0-jbegin+1];
2468 free(temp.structure);
2470 strcpy(struc, st1+MAX2(i0,1)); strcat(struc, "&");
2471 /* strcat(struc, st2); */
2472 strncat(struc, struc2+5, strlen(struc2)-10);
2475 free(st1); free(st2);
2477 for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
2480 free(S1);free(S2);free(SS1);free(SS2);
2481 /* free_arrays(); */
2485 PRIVATE int covscore(const int *types, int n_seq) {
2486 /* calculate co-variance bonus for a pair depending on */
2487 /* compensatory/consistent mutations and incompatible seqs */
2488 /* should be 0 for conserved pairs, >0 for good pairs */
2489 #define NONE -10000 /* score for forbidden pairs */
2490 int k,l,s,score, pscore;
2491 int dm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
2492 {0,0,2,2,1,2,2} /* CG */,
2493 {0,2,0,1,2,2,2} /* GC */,
2494 {0,2,1,0,2,1,2} /* GU */,
2495 {0,1,2,2,0,2,1} /* UG */,
2496 {0,2,2,1,2,0,2} /* AU */,
2497 {0,2,2,2,1,2,0} /* UA */};
2499 int pfreq[8]={0,0,0,0,0,0,0,0};
2500 for (s=0; s<n_seq; s++)
2503 if (pfreq[0]*2>n_seq) return NONE;
2504 for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
2505 for (l=k+1; l<=6; l++)
2506 /* scores for replacements between pairtypes */
2507 /* consistent or compensatory mutations score 1 or 2 */
2508 score += pfreq[k]*pfreq[l]*dm[k][l];
2510 /* counter examples score -1, gap-gap scores -0.25 */
2512 ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
2516 /*---------------------------------------------------------------------------*/
2518 PRIVATE short * aliencode_seq(const char *sequence) {
2521 l = strlen(sequence);
2522 Stemp = (short *) space(sizeof(short)*(l+2));
2523 Stemp[0] = (short) l;
2525 /* make numerical encoding of sequence */
2526 for (i=1; i<=l; i++)
2527 Stemp[i]= (short) encode_char(toupper(sequence[i-1]));
2529 /* for circular folding add first base at position n+1 */
2530 /* Stemp[l+1] = Stemp[1]; */
2535 PRIVATE short * encode_seq(const char *sequence) {
2538 l = strlen(sequence);
2539 extern double nc_fact;
2540 S = (short *) space(sizeof(short)*(l+2));
2543 /* make numerical encoding of sequence */
2544 for (i=1; i<=l; i++)
2545 S[i]= (short) encode_char(toupper(sequence[i-1]));
2546 /* for circular folding add first base at position n+1 */
2552 PRIVATE void encode_seqs(const char *s1, const char *s2) {
2556 S1 = encode_seq(s1);
2557 SS1= (short *) space(sizeof(short)*(l+1));
2558 /* SS1 exists only for the special X K and I bases and energy_set!=0 */
2560 for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
2561 SS1[i] = alias[S1[i]]; /* for mismatches of nostandard bases */
2565 S2 = encode_seq(s2);
2566 SS2= (short *) space(sizeof(short)*(l+1));
2567 /* SS2 exists only for the special X K and I bases and energy_set!=0 */
2569 for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
2570 SS2[i] = alias[S2[i]]; /* for mismatches of nostandard bases */
2574 PRIVATE int compare(const void *sub1, const void *sub2) {
2576 if (((snoopT *) sub1)->energy > ((snoopT *) sub2)->energy)
2578 if (((snoopT *) sub1)->energy < ((snoopT *) sub2)->energy)
2580 d = ((snoopT *) sub1)->i - ((snoopT *) sub2)->i;
2582 return ((snoopT *) sub1)->j - ((snoopT *) sub2)->j;