--- /dev/null
+
+
+/* Last changed Time-stamp: <2007-08-26 11:59:45 ivo> */
+/*
+ compute the duplex structure of two RNA strands,
+ allowing only inter-strand base pairs.
+ see cofold() for computing hybrid structures without
+ restriction.
+
+ Ivo Hofacker
+ Vienna RNA package
+*/
+
+#include <config.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <ctype.h>
+#include <string.h>
+#include "utils.h"
+#include "energy_par.h"
+#include "fold_vars.h"
+#include "snofold.h"
+#include "pair_mat.h"
+#include "params.h"
+#include "snoop.h"
+#include "PS_dot.h"
+/* #include "fold.h" */
+#include "duplex.h"
+#include "loop_energies.h"
+/*@unused@*/
+static char rcsid[] UNUSED = "$Id: duplex.c,v 1.8 2007/08/26 10:08:44 ivo Exp $";
+#define UNIT 100
+#define MINPSCORE -2*UNIT
+#define PUBLIC
+#define PRIVATE static
+
+#define STACK_BULGE1 1 /* stacking energies for bulges of size 1 */
+#define NEW_NINIO 1 /* new asymetry penalty */
+
+
+
+PRIVATE void encode_seqs(const char *s1, const char *s2);
+PRIVATE short *encode_seq(const char *seq);
+
+PRIVATE void find_max_snoop(const char *s1, const char *s2, const int max,
+ const int alignment_length, const int* position,
+ const int delta, const int distance, const int penalty,
+ const int threshloop, const int threshLE, const int threshRE,
+ const int threshDE, const int threshTE, const int threshSE, const int threshD,
+ 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);
+
+PRIVATE void find_max_snoop_XS(const char *s1, const char *s2, const int **access_s1, const int max,
+ const int alignment_length, const int* position, const int *position_j,
+ const int delta, const int distance, const int penalty,
+ const int threshloop, const int threshLE, const int threshRE,
+ const int threshDE, const int threshTE, const int threshSE, const int threshD,
+ 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);
+
+
+
+
+
+PRIVATE char * alisnoop_backtrack(int i, int j, const char ** s2,
+ int* Duplex_El, int* Duplex_Er, int* Loop_E, int *Loop_D, int *u,
+ int *pscd, int *psct, int *pscg,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshD,
+ 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 short **S1, const short **S2);
+
+
+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,
+ const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
+ const int threshD,
+ 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);
+
+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,
+ const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
+ const int threshD,
+ 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);
+
+
+
+
+PRIVATE int compare(const void *sub1, const void *sub2);
+PRIVATE int covscore(const int *types, int n_seq);
+PRIVATE short * aliencode_seq(const char *sequence);
+
+PUBLIC int snoop_subopt_sorted=0; /* from subopt.c, default 0 */
+
+
+/*@unused@*/
+
+
+
+
+#define MAXLOOP_L 3
+#define MIN2(A, B) ((A) < (B) ? (A) : (B))
+#define MAX2(A, B) ((A) > (B) ? (A) : (B))
+#define ASS 1
+PRIVATE paramT *P = NULL;
+
+PRIVATE int **c = NULL; /* energy array, given that i-j pair */
+PRIVATE int **r = NULL;
+PRIVATE int **lc = NULL; /* energy array, given that i-j pair */
+PRIVATE int **lr = NULL;
+PRIVATE int **c_fill = NULL;
+PRIVATE int **r_fill = NULL;
+PRIVATE int **lpair = NULL;
+
+
+PRIVATE short *S1 = NULL, *SS1 = NULL, *S2 = NULL, *SS2 = NULL;
+PRIVATE short *S1_fill = NULL, *SS1_fill = NULL, *S2_fill = NULL, *SS2_fill = NULL;
+PRIVATE int n1,n2; /* sequence lengths */
+
+extern int cut_point;
+
+PRIVATE int delay_free=0;
+/*--------------------------------------------------------------------------*/
+
+snoopT alisnoopfold(const char **s1, const char **s2,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshD,
+ 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) {
+
+ int s,n_seq;
+ int i, j, E, l1,Emin=INF, i_min=0, j_min=0;
+ char *struc;
+ snoopT mfe;
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden **foldlist; folden **foldlist_XS;
+ int Duplex_El, Duplex_Er,pscd,psct,pscg;
+ int Loop_D;
+ int u;
+ int Loop_E;
+ short **Sali1,**Sali2;
+ int *type,*type2,*type3;
+ Duplex_El=0;Duplex_Er=0;Loop_E=0; Loop_D=0;pscd=0;psct=0;pscg=0;
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS);
+ n1 = (int) strlen(s1[0]);
+ n2 = (int) strlen(s2[0]);
+
+ for (s=0; s1[s]!=NULL; s++);
+ n_seq = s;
+ for (s=0; s2[s]!=NULL; s++);
+ if (n_seq != s) nrerror("unequal number of sequences in aliduplexfold()\n");
+
+ if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
+ snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
+ make_pair_matrix();
+ }
+
+ c = (int **) space(sizeof(int *) * (n1+1));
+ r = (int **) space(sizeof(int *) * (n1+1));
+ for (i=0; i<=n1; i++) {
+ c[i] = (int *) space(sizeof(int) * (n2+1));
+ r[i] = (int *) space(sizeof(int) * (n2+1));
+ for(j=n2; j>-1; j--){
+ c[i][j]=INF;
+ r[i][j]=INF;
+ }
+ }
+ Sali1 = (short **) space((n_seq+1)*sizeof(short *));
+ Sali2 = (short **) space((n_seq+1)*sizeof(short *));
+ for (s=0; s<n_seq; s++) {
+ if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
+ if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
+ Sali1[s] = aliencode_seq(s1[s]);
+ Sali2[s] = aliencode_seq(s2[s]);
+ }
+ type = (int *) space(n_seq*sizeof(int));
+ type2 = (int *) space(n_seq*sizeof(int));
+ type3 = (int *) space(n_seq*sizeof(int));
+ /* encode_seqs(s1, s2); */
+ for (i=6; i<=n1-5; i++) {
+ int U; U=0;
+ for (s=0; s<n_seq; s++) {
+ U+=Sali1[s][i-2];
+ }
+ U = (U==(n_seq)*4?1:0);
+ for (j=n2-min_d2; j>min_d1; j--) {
+ int type4, k,l,psc,psc2,psc3;
+ for (s=0; s<n_seq; s++) {
+ type[s] = pair[Sali1[s][i]][Sali2[s][j]];
+ }
+ psc = covscore(type, n_seq);
+ for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
+ c[i][j] = (psc>=MINPSCORE) ? (n_seq*P->DuplexInit) : INF;
+ if (psc<MINPSCORE) continue;
+ if(/* pair[Sali1[i+1]][Sali2[j-1]] && */
+ U && j < max_s1 && j > min_s1 &&
+ j > n2 - max_s2 - max_half_stem &&
+ j < n2 -min_s2 -half_stem ) { /*constraint on s2 and i*/
+ folden *temp;
+ temp=foldlist[j+1];
+ while(temp->next){
+ int k = temp->k;
+ for (s=0; s<n_seq; s++) {
+ type2[s]= pair[Sali1[s][i-3]][Sali2[s][k+1]];
+ type3[s]= pair[Sali1[s][i-4]][Sali2[s][k+1]];
+ }
+ psc2 = covscore(type2, n_seq);
+ psc3 = covscore(type3, n_seq);
+ if(psc2 > MINPSCORE){
+ r[i][j]=MIN2(r[i][j],c[i-3][k+1]+temp->energy);
+ }
+ if(psc3 > MINPSCORE){
+ r[i][j]=MIN2(r[i][j],c[i-4][k+1]+temp->energy);
+ }
+ temp=temp->next;
+ }
+ }
+ /* dangle 5'SIDE relative to the mRNA */
+ for (s=0; s<n_seq; s++) {
+ c[i][j] += E_ExtLoop(type[s], Sali1[s][i-1],Sali2[s][j+1],P);
+ }
+ for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
+ for (l=j+1; l<=n2 ; l++) {
+ if (i-k+l-j>2*MAXLOOP_L-2) break;
+ if (abs(i-k-l+j) >= ASS ) continue;
+ for (E=s=0; s<n_seq; s++) {
+ type4 = pair[Sali1[s][k]][Sali2[s][l]];
+ if (type4==0) type4=7;
+ E += 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);
+ }
+ c[i][j] = MIN2(c[i][j], c[k][l] + E);
+ r[i][j] = MIN2(r[i][j], r[k][l] + E);
+ }
+ }
+ c[i][j]-=psc;
+ r[i][j]-=psc;
+ E = r[i][j];
+ for (s=0; s<n_seq; s++) {
+ E+= E_ExtLoop(rtype[type[s]], Sali2[s][j-1], Sali1[s][i+1], P);
+ /**
+ *** if (i<n1) E += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
+ *** if (j>1) E += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
+ *** if (type[s]>2) E += P->TerminalAU;
+ **/
+ }
+ if (E<Emin) {
+ Emin=E; i_min=i; j_min=j;
+ }
+ }
+ }
+ if(Emin > 0){
+ printf("no target found under the constraints chosen\n");
+ for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
+ free(c);
+ free(r);
+ for(s=0; s<n_seq;s++){
+ free(Sali1[s]);
+ free(Sali2[s]);
+ }
+ free(Sali1); free(Sali2);
+ free(S2); free(SS1); free(SS2);free(type);free(type2);free(type3);
+ mfe.energy=INF;
+ mfe.structure=NULL;
+ return mfe;
+ }
+ struc = alisnoop_backtrack(i_min, j_min,(const char**) s2,
+ &Duplex_El, &Duplex_Er, &Loop_E,
+ &Loop_D, &u, &pscd, &psct, &pscg,
+ penalty, threshloop, threshLE,
+ threshRE,threshDE, threshD,
+ half_stem, max_half_stem, min_s2,
+ max_s2, min_s1, max_s1, min_d1,
+ min_d2,(const short**) Sali1,(const short**) Sali2);
+ /* if (i_min<n1-5) i_min++; */
+ /* if (j_min>6 ) j_min--; */
+ l1 = strchr(struc, '&')-struc;
+ mfe.i = i_min-5;
+ mfe.j = j_min-5;
+ mfe.u = u -5;
+ mfe.Duplex_Er = (float) Duplex_Er/100;
+ mfe.Duplex_El = (float) Duplex_El/100;
+ mfe.Loop_D = (float) Loop_D/100;
+ mfe.Loop_E = (float) Loop_E/100;
+ mfe.energy = (float) Emin/100 ;
+ /* mfe.fullStemEnergy = (float) fullStemEnergy/100; */
+ mfe.pscd = pscd;
+ mfe.psct = psct;
+ mfe.structure = struc;
+ for(s=0; s<n_seq;s++){
+ free(Sali1[s]);free(Sali2[s]);
+ }
+ free(Sali1);free(Sali2);free(type);free(type2);free(type3);
+
+ if (!delay_free) {
+ for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
+ free(c);
+ free(r);
+ free(S2); free(SS1); free(SS2);
+ }
+ return mfe;
+}
+
+PUBLIC snoopT *alisnoop_subopt(const char **s1, const char **s2, int delta, int w,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
+ const int distance, 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) {
+
+
+ short **Sali1, **Sali2;
+ /* printf("%d %d\n", min_s2, max_s2); */
+ int i,j,s,n_seq, n1, n2, E, n_subopt=0, n_max;
+ char *struc;
+ snoopT mfe;
+ snoopT *subopt;
+ int thresh;
+ int *type;
+ int Duplex_El, Duplex_Er, Loop_E,pscd,psct,pscg;
+ int Loop_D;
+ Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;pscd=0;psct=0;pscg=0;
+ int u;
+ u=0;
+ n_max=16;
+ subopt = (snoopT *) space(n_max*sizeof(snoopT));
+ delay_free=1;
+ mfe = alisnoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE,threshD,
+ half_stem, max_half_stem,
+ min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
+ if(mfe.energy > 0){
+ free(subopt);
+ delay_free=0;
+ return NULL;
+ }
+ thresh = MIN2((int) ((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E)*100+0.1 + 410) + delta, threshTE );
+ /* subopt[n_subopt++]=mfe; */
+ free(mfe.structure);
+ n1 = strlen(s1[0]);
+ n2 = strlen(s2[0]);
+ for (s=0; s1[s]!=NULL; s++);
+ n_seq = s;
+ Sali1 = (short **) space((n_seq+1)*sizeof(short *));
+ Sali2 = (short **) space((n_seq+1)*sizeof(short *));
+ for (s=0; s<n_seq; s++) {
+ if (strlen(s1[s]) != n1) nrerror("uneqal seqence lengths");
+ if (strlen(s2[s]) != n2) nrerror("uneqal seqence lengths");
+ Sali1[s] = aliencode_seq(s1[s]);
+ Sali2[s] = aliencode_seq(s2[s]);
+ }
+ Sali1[n_seq]=NULL; Sali2[n_seq]=NULL;
+ type = (int *) space(n_seq*sizeof(int));
+ for (i=n1; i>1; i--){
+ for (j=1; j<=n2; j++) {
+ int ii,jj, Ed,psc,skip;
+ for (s=0; s<n_seq; s++) {
+ type[s] = pair[Sali2[s][j]][Sali1[s][i]];
+ }
+ psc = covscore(type, n_seq);
+ for (s=0; s<n_seq; s++) if (type[s]==0) type[s]=7;
+ if (psc<MINPSCORE) continue;
+ E = Ed = r[i][j];
+ for (s=0; s<n_seq; s++) {
+ /* if (i<n1-5) Ed += P->dangle3[type[s]][Sali1[s][i+1]]; */
+ /* if (j>6) Ed += P->dangle5[type[s]][Sali2[s][j-1]]; */
+ if (type[s]>2) Ed += P->TerminalAU;
+ }
+ if (Ed>thresh) continue;
+ /* too keep output small, remove hits that are dominated by a
+ better one close (w) by. For simplicity we do test without
+ adding dangles, which is slightly inaccurate.
+ */
+ w=1;
+ for (skip=0, ii=MAX2(i-w,1); (ii<=MIN2(i+w,n1)) && type; ii++) {
+ for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n2); jj++)
+ if (r[ii][jj]<E) {skip=1; break;}
+ }
+ if (skip){continue;}
+ psct=0;
+ pscg=0;
+ struc = alisnoop_backtrack(i,j,s2, &Duplex_El,
+ &Duplex_Er, &Loop_E, &Loop_D, &u, &pscd, &psct,&pscg,
+ penalty, threshloop,threshLE,threshRE,threshDE, threshD,
+ 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);
+
+ if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
+ (Duplex_Er + Duplex_El) > threshDE ||
+ (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
+ (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
+ /* printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n" */
+ /* " Duplex_Er + Duplex_El %d threshDE %d \n" */
+ /* " Duplex_Er + Duplex_El + Loop_E %d threshTE %d \n" */
+ /* " Duplex_Er + Duplex_El + Loop_E + Loop_D %d threshSE %d \n", */
+ /* Duplex_Er , threshRE , Duplex_El ,threshLE, */
+ /* Duplex_Er + Duplex_El, threshDE, */
+ /* Duplex_Er + Duplex_El+ Loop_E , threshTE, */
+ /* Duplex_Er + Duplex_El+ Loop_E + Loop_D, threshSE); */
+ Duplex_Er=0;
+ Duplex_El=0;
+ Loop_E = 0;
+ Loop_D = 0;
+ u=0,
+ free(struc);
+ continue;
+ }
+
+ if (n_subopt+1>=n_max) {
+ n_max *= 2;
+ subopt = (snoopT *) xrealloc(subopt, n_max*sizeof(snoopT));
+ }
+
+ subopt[n_subopt].i = i-5;
+ subopt[n_subopt].j = j-5;
+ subopt[n_subopt].u = u-5;
+ subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
+ subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
+ subopt[n_subopt].Loop_E = Loop_E * 0.01;
+ subopt[n_subopt].Loop_D = Loop_D * 0.01;
+ subopt[n_subopt].energy = (Duplex_Er +Duplex_El + Loop_E + Loop_D + 410) * 0.01 ;
+ subopt[n_subopt].pscd = pscd * 0.01;
+ subopt[n_subopt].psct = -psct * 0.01;
+ subopt[n_subopt++].structure = struc;
+
+ /* i=u; */
+ Duplex_Er=0; Duplex_El=0; Loop_E=0; Loop_D=0;u=0;pscd=0;psct=0;
+ }
+ }
+
+ for (i=0; i<=n1; i++) {free(c[i]);free(r[i]);}
+ free(c);free(r);
+ for (s=0; s<n_seq; s++) {
+ free(Sali1[s]); free(Sali2[s]);
+ }
+ free(Sali1); free(Sali2); free(type);
+
+ if (snoop_subopt_sorted)
+ qsort(subopt, n_subopt, sizeof(snoopT), compare);
+ subopt[n_subopt].i =0;
+ subopt[n_subopt].j =0;
+ subopt[n_subopt].structure = NULL;
+ return subopt;
+}
+
+
+
+
+
+
+
+PRIVATE char *alisnoop_backtrack(int i, int j, const char ** snoseq, int *Duplex_El,
+ int *Duplex_Er, int *Loop_E, int *Loop_D, int *u,
+ int *pscd, int *psct, int *pscg,
+ const int penalty, const int threshloop, const int threshLE,
+ const int threshRE, const int threshDE, const int threshD, 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 short **Sali1, const short **Sali2) {
+ /* backtrack structure going backwards from i, and forwards from j
+ return structure in bracket notation with & as separator */
+ int k, l, *type,*type2,*type3,type4, E, traced, i0, j0,s,n_seq,psc;
+ int traced_r=0; /* flag for following backtrack in c or r */
+ char *st1, *st2, *struc;
+ char *struc_loop;
+ n1 = (int) Sali1[0][0];
+ n2 = (int) Sali2[0][0];
+
+ for (s=0; Sali1[s]!=NULL; s++);
+ n_seq = s;
+ for (s=0; Sali2[s]!=NULL; s++);
+ if (n_seq != s) nrerror("unequal number of sequences in alibacktrack()\n");
+
+ st1 = (char *) space(sizeof(char)*(n1+1));
+ st2 = (char *) space(sizeof(char)*(n2+1));
+ type = (int *) space(n_seq*sizeof(int));
+ type2 = (int *) space(n_seq*sizeof(int));
+ type3 = (int *) space(n_seq*sizeof(int));
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden **foldlist, **foldlist_XS;
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
+ i0=i; j0=j; /* MIN2(i+1,n1); j0=MAX2(j-1,1);!modified */
+ for (s=0; s<n_seq; s++) {
+ type[s] = pair[Sali1[s][i]][Sali2[s][j]];
+ if(type[s]==0) type[s] = 7;
+ *Duplex_Er += E_ExtLoop(rtype[type[s]], (j>1) ? Sali2[s][j-1] : -1, (i<n1) ? Sali1[s][i+1] : -1, P);
+ /**
+ *** if (i<n1) *Duplex_Er += P->dangle3[rtype[type[s]]][Sali1[s][i+1]];
+ *** if (j>1) *Duplex_Er += P->dangle5[rtype[type[s]]][Sali2[s][j-1]];
+ *** if (type[s]>2) *Duplex_Er += P->TerminalAU;
+ **/
+ }
+ while (i>0 && j<=n2-min_d2 ) {
+ if(!traced_r) {
+ E = r[i][j]; traced=0;
+ st1[i-1] = '<';
+ st2[j-1] = '>';
+ for (s=0; s<n_seq; s++) {
+ type[s] = pair[Sali1[s][i]][Sali2[s][j]];
+ }
+ psc = covscore(type,n_seq);
+ for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
+ E += psc;
+ *pscd +=psc;
+ for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
+ for (l=j+1; l<=n2 ; l++) {
+ int LE;
+ if (i-k+l-j>2*MAXLOOP_L-2) break;
+ if (abs(i-k-l+j) >= ASS) continue;
+ for (s=LE=0; s<n_seq; s++) {
+ type4 = pair[Sali1[s][k]][Sali2[s][l]];
+ if (type4==0) type4=7;
+ 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);
+ }
+ if (E == r[k][l]+LE) {
+ traced=1;
+ i=k; j=l;
+ *Duplex_Er+=LE;
+ break;
+ }
+ }
+ if (traced) break;
+ }
+ if(!traced){
+ int U=0;
+ for (s=0; s<n_seq; s++) {
+ U+=Sali1[s][i-2];
+ }
+ U = (U==(n_seq)*4?1:0);
+ if(/* pair[Sali1[i+1]][Sali2[j-1]] && */ /* only U's are allowed */
+ U && j < max_s1 && j > min_s1 &&
+ j > n2 - max_s2 - max_half_stem &&
+ j < n2 -min_s2 -half_stem ) {
+ int min_k, max_k;
+ max_k = MIN2(n2-min_s2,j+max_half_stem+1);
+ min_k = MAX2(j+half_stem+1, n2-max_s2);
+ folden * temp;
+ temp=foldlist[j+1];
+ while(temp->next) {
+ int psc2, psc3;
+ int k = temp->k;
+ for (s=0; s<n_seq; s++) {
+ type2[s]= pair[Sali1[s][i-3]][Sali2[s][k+1]];
+ type3[s]= pair[Sali1[s][i-4]][Sali2[s][k+1]];
+ }
+ psc2 = covscore(type2, n_seq);
+ psc3 = covscore(type3, n_seq);
+ if(psc2>MINPSCORE /*&& pair[Sali1[i-4]][Sali2[k+2]]*/ ){ /* introduce structure from RNAfold */
+ if(E==c[i-3][k+1]+temp->energy){
+ *Loop_E=temp->energy;
+ st1[i-3]= '|';
+ *u=i-2;
+ int a,b;
+ /* int fix_ij=indx[k-1+1]+j+1; */
+ for(a=0; a< MISMATCH ;a++){
+ for(b=0; b< MISMATCH ; b++){
+ int ij=indx[k-1-a+1]+j+1+b;
+ if(cLoop[ij]==temp->energy) {
+ /* int bla; */
+ struc_loop=alisnobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1,psct);
+ a=INF; b=INF;
+ }
+ }
+ }
+ traced=1;
+ traced_r=1;
+ i=i-3;j=k+1;
+ break;
+ }
+ }
+ if (psc3>MINPSCORE /*&& pair[Sali1[i-5]][Sali2[k+2]]*/){ /* introduce structure from RNAfold */
+ if(E==c[i-4][k+1]+temp->energy){
+ *Loop_E=temp->energy;
+ st1[i-3]= '|';
+ *u=i-2;
+ int a,b;
+ /* int fix_ij=indx[k-1+1]+j+1; */
+ for(a=0; a< MISMATCH ;a++){
+ for(b=0; b< MISMATCH ; b++){
+ int ij=indx[k-1-a+1]+j+1+b;
+ if(cLoop[ij]==temp->energy) {
+ /* int bla; */
+ struc_loop=alisnobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1,psct);
+ a=INF; b=INF;
+ }
+ }
+ }
+ traced=1;
+ traced_r=1;
+ i=i-4;j=k+1;
+ break;
+ }
+ } /* else if */
+ temp=temp->next;
+ } /* while temp-> next */
+ } /* test on j */
+ }/* traced? */
+ }/* traced_r? */
+ else{
+ E = c[i][j]; traced=0;
+ st1[i-1] = '<';
+ st2[j-1] = '>';
+ for (s=0; s<n_seq; s++) {
+ type[s] = pair[Sali1[s][i]][Sali2[s][j]];
+ }
+ psc = covscore(type,n_seq);
+ for (s=0; s<n_seq; s++) if (type[s]==0) type[s] = 7;
+ E += psc;
+ *pscd+=psc;
+ if (!type) nrerror("backtrack failed in fold duplex c");
+ for (k=i-1; (i-k)<MAXLOOP_L; k--) {
+ for (l=j+1; l<=n2; l++) {
+ int LE;
+ if (i-k+l-j>2*MAXLOOP_L-2) break;
+ if (abs(i-k-l+j) >= ASS) continue;
+ for (s=LE=0; s<n_seq; s++) {
+ type4 = pair[Sali1[s][k]][Sali2[s][l]];
+ if (type4==0) type4=7;
+ 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);
+ }
+ if (E == c[k][l]+LE) {
+ traced=1;
+ i=k; j=l;
+ *Duplex_El+=LE;
+ break;
+ }
+ }
+ if (traced) break;
+ }
+ }
+ if (!traced) {
+ for (s=0; s<n_seq; s++) {
+ int correction;
+ correction = E_ExtLoop(type[s], (i>1) ? Sali1[s][i-1] : -1, (j<n2) ? Sali2[s][j+1] : -1, P);
+ *Duplex_El += correction;
+ E -= correction;
+ /**
+ *** if (i>1) {E -= P->dangle5[type[s]][Sali1[s][i-1]]; *Duplex_El +=P->dangle5[type[s]][Sali1[s][i-1]];}
+ *** if (j<n2) {E -= P->dangle3[type[s]][Sali2[s][j+1]]; *Duplex_El +=P->dangle3[type[s]][Sali2[s][j+1]];}
+ *** if (type[s]>2) {E -= P->TerminalAU; *Duplex_El +=P->TerminalAU;}
+ **/
+ }
+ if (E != n_seq * P->DuplexInit) {
+ nrerror("backtrack failed in fold duplex end");
+ } else break;
+ }
+ }
+/* if (i>1) i--; */
+/* if (j<n2) j++; */
+ /* struc = (char *) space(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
+ struc = (char *) space(i0-i+1+n2-1+1+2); /* declare final duplex structure */
+ char * struc2;
+ struc2 = (char *) space(n2+1);
+ /* char * struct_const; */
+ for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
+ /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
+ /* char * struct_const; */
+ /* struct_const = (char *) space(sizeof(char)*(n2+1)); */
+ for (k=1; k<=n2; k++) {
+ if (!st2[k-1]) st2[k-1] = struc_loop[k-1];/* '.'; */
+ struc2[k-1] = st2[k-1];/* '.'; */
+ /* if (k>=j0 && k<=j){ */
+ /* struct_const[k-1]='x'; */
+ /* } */
+ /* else{ */
+ /* if(k<j0) {struct_const[k-1]='<';} */
+ /* if(k>j) {struct_const[k-1]='>';} */
+ /* } */
+ }
+
+ /* char duplexseq_1[j0+1]; */
+ /* char duplexseq_2[n2-j+3]; */
+ if(j<n2){
+ char **duplexseq_1, **duplexseq_2;
+ duplexseq_1 = (char**) space((n_seq+1) * sizeof(char*));
+ duplexseq_2 = (char**) space((n_seq+1) * sizeof(char*));
+ for(s=0; s<n_seq; s++){
+ duplexseq_1[s] = (char*) space((j0)*sizeof(char)); /* modfied j0+1 */
+ duplexseq_2[s] = (char*) space((n2-j+2)*sizeof(char)); /* modified j+3 */
+ strncpy(duplexseq_1[s], snoseq[s], j0-1); /* modified j0 */
+ strcpy(duplexseq_2[s], snoseq[s] + j); /* modified j-1 */
+ duplexseq_1[s][j0-1]='\0'; /* modified j0 */
+ duplexseq_2[s][n2-j+1]='\0';/* modified j+2 */
+ }
+ duplexseq_1[n_seq]=NULL;
+ duplexseq_2[n_seq]=NULL;
+ duplexT temp;
+ temp=aliduplexfold((const char**)duplexseq_1, (const char**)duplexseq_2);
+ *Loop_D = MIN2(0,-410 + (int) 100 * temp.energy*n_seq);
+ if(*Loop_D){
+ int l1,ibegin, iend, jbegin, jend;
+ l1=strchr(temp.structure, '&')-temp.structure;
+ ibegin=temp.i-l1;
+ iend =temp.i-1;
+ jbegin=temp.j;
+ jend =temp.j+strlen(temp.structure)-l1-2-1;
+ for(k=ibegin+1; k<=iend+1; k++){
+ struc2[k-1]=temp.structure[k-ibegin-1];
+ }
+ for(k=jbegin+j; k<=jend+j; k++){
+ struc2[k-1]=temp.structure[l1+k-j-jbegin+1];
+ }
+ }
+ for(s=0; s<n_seq; s++){
+ free(duplexseq_1[s]);
+ free(duplexseq_2[s]);
+ }
+ free(duplexseq_1);free(duplexseq_2);
+ free(temp.structure);
+ }
+ strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&");
+ /* strcat(struc, st2); */
+ strncat(struc, struc2+5, strlen(struc2)-10);
+ free(struc2);
+ free(struc_loop);
+ free(st1); free(st2);
+ free(type);free(type2);free(type3);
+ /* free_arrays(); */
+ return struc;
+}
+
+
+
+
+
+
+
+void Lsnoop_subopt(const char *s1, const char *s2, int delta, int w,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
+ const int distance,
+ 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 int alignment_length, const char* name, const int fullStemEnergy)
+{
+
+ int min_colonne=INF;
+ int max_pos;
+ int max;max=INF;
+ /* int temp; */
+ /* int nsubopt=10; */
+ n1 = (int) strlen(s1);
+ n2 = (int) strlen(s2);
+ int *position;
+ position = (int*) space((n1+3)*sizeof(int));
+
+
+ /* int Eminj, Emin_l; */
+ int i, j; /* l1, Emin=INF, i_min=0, j_min=0; */
+ /* char *struc; */
+ /* snoopT mfe; */
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden **foldlist, **foldlist_XS;
+ int Duplex_El, Duplex_Er;
+ int Loop_D;
+ /* int u; */
+ int Loop_E;
+ Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
+ if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
+ snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
+ make_pair_matrix();
+ }
+
+ lc = (int **) space(sizeof(int *) * (5));
+ lr = (int **) space(sizeof(int *) * (5));
+ for (i=0; i<5; i++) {
+ lc[i] = (int *) space(sizeof(int) * (n2+1));
+ lr[i] = (int *) space(sizeof(int) * (n2+1));
+ for(j=n2; j>-1; j--){
+ lc[i][j]=INF;
+ lr[i][j]=INF;
+ }
+ }
+ encode_seqs(s1, s2);
+ for (i=1; i<=n1; i++) {
+ int idx=i%5;
+ int idx_1=(i-1)%5;
+ int idx_2=(i-2)%5;
+ int idx_3=(i-3)%5;
+ int idx_4=(i-4)%5;
+ for (j=n2-min_d2; j>min_d1; j--) {
+ int type, type2, k;
+ type = pair[S1[i]][S2[j]];
+ lc[idx][j] = (type) ? P->DuplexInit + 2*penalty : INF;
+ lr[idx][j] = INF;
+ if(!type) continue;
+ if( /*pair[S1[i+1]][S2[j-1]] && check that we have a solid base stack after the mLoop */
+ j < max_s1 && j > min_s1 &&
+ j > n2 - max_s2 - max_half_stem &&
+ j < n2 -min_s2 -half_stem && S1[i-2]==4) { /*constraint on s2 and i*/
+ int min_k, max_k;
+ max_k = MIN2(n2-min_s2,j+max_half_stem+1);
+ min_k = MAX2(j+half_stem+1, n2-max_s2);
+ for(k=min_k; k <= max_k ; k++){
+ if(mLoop[indx[k-1]+j+1] < 0){
+ }
+ if(pair[S1[i-3]][S2[k]] /*genau zwei ungepaarte nucleotiden --NU--*/
+ && mLoop[indx[k-1]+j+1] < threshloop){
+ lr[idx][j]=MIN2(lr[idx][j], lc[idx_3][k]+mLoop[indx[k-1]+j+1]);
+ }
+ else if(pair[S1[i-4]][S2[k]] && mLoop[indx[k-1]+j+1] < threshloop){/*--NUN--*/
+ lr[idx][j]=MIN2(lr[idx][j], lc[idx_4][k]+mLoop[indx[k-1]+j+1]);
+ }
+ }
+ }
+ /* dangle 5'SIDE relative to the mRNA */
+ lc[idx][j] += E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, P);
+ /**
+ *** if (i>1) lc[idx][j] += P->dangle5[type][SS1[i-1]];
+ *** if (j<n2) lc[idx][j] += P->dangle3[type][SS2[j+1]];
+ *** if (type>2) lc[idx][j] += P->TerminalAU;
+ **/
+
+ if(j<n2 && i>1){
+ type2=pair[S1[i-1]][S2[j+1]];
+ if(type2>0){
+ 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]);
+ 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]);
+ }
+ }
+ if(j<n2-1 && i>2){
+ type2=pair[S1[i-2]][S2[j+2]];
+ if(type2>0 ){
+ 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]);
+ 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]);
+ }
+ }
+ if(j<n2-2 && i>3){
+ type2 = pair[S1[i-3]][S2[j+3]];
+ if(type2>0){
+ 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]);
+ 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]);
+ }
+ }
+ /**
+ *** (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)
+ **/
+ 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);
+ }
+ position[i]=min_colonne;
+ if(max>=min_colonne){
+ max=min_colonne;
+ max_pos=i;
+ }
+ min_colonne=INF;
+ }
+
+ free(S1); free(S2); free(SS1); free(SS2);
+ if(max<threshTE){
+ find_max_snoop(s1, s2, max, alignment_length, position, delta,
+ distance, penalty, threshloop, threshLE, threshRE, threshDE,
+ threshTE, threshSE, threshD, half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name, fullStemEnergy);
+ }
+ for (i=1; i<5; i++) {free(lc[i]);free(lr[i]);}
+ free(lc[0]);free(lr[0]);
+ free(lc);free(lr);
+ free(position);
+
+}
+
+
+
+void Lsnoop_subopt_list(const char *s1, const char *s2, int delta, int w,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
+ const int distance,
+ 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 int alignment_length,const char *name,const int fullStemEnergy)
+{
+
+ int min_colonne=INF;
+ int max_pos;
+ int max;max=INF;
+ /* int temp; */
+ /* int nsubopt=10; */
+ n1 = (int) strlen(s1);
+ n2 = (int) strlen(s2);
+ int *position;
+ position = (int*) space((n1+3)*sizeof(int));
+
+
+ /* int Eminj, Emin_l; */
+ int i, j;/* l1, Emin=INF, i_min=0, j_min=0; */
+ /* char *struc; */
+ /* snoopT mfe; */
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden **foldlist, **foldlist_XS;
+ int Duplex_El, Duplex_Er;
+ int Loop_D;
+ /* int u; */
+ int Loop_E;
+ Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
+ if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
+ snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
+ make_pair_matrix();
+ }
+
+ lpair = (int **) space(sizeof(int *) * (6));
+ lc = (int **) space(sizeof(int *) * (6));
+ lr = (int **) space(sizeof(int *) * (6));
+ for (i=0; i<6; i++) {
+ lc[i] = (int *) space(sizeof(int) * (n2+1));
+ lr[i] = (int *) space(sizeof(int) * (n2+1));
+ lpair[i] = (int *) space(sizeof(int) * (n2+1));
+ for(j=n2; j>-1; j--){
+ lc[i][j]=INF;
+ lr[i][j]=INF;
+ lpair[i][j]=0;
+ }
+ }
+ encode_seqs(s1, s2);
+ int lim_maxj=n2-min_d2 ;
+ int lim_minj=min_d1;
+ int lim_maxi=n1;
+ for (i=5; i<=lim_maxi; i++) {
+ int idx=i%5;
+ int idx_1=(i-1)%5;
+ int idx_2=(i-2)%5;
+ int idx_3=(i-3)%5;
+ int idx_4=(i-4)%5;
+
+ for (j=lim_maxj; j>lim_minj; j--) {
+ int type, type2;/* E, k,l; */
+ type = pair[S1[i]][S2[j]];
+ lpair[idx][j] = type;
+ lc[idx][j] = (type) ? P->DuplexInit + 2*penalty : INF;
+ lr[idx][j] = INF;
+ if(!type) continue;
+ if( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
+ j < max_s1 && j > min_s1 &&
+ j > n2 - max_s2 - max_half_stem &&
+ j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
+ int min_k, max_k;
+ max_k = MIN2(n2-min_s2,j+max_half_stem+1);
+ min_k = MAX2(j+half_stem+1, n2-max_s2);
+ folden * temp;
+ temp=foldlist[j+1];
+ while(temp->next){
+ int k = temp->k;
+ /* if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour */
+ if(lpair[idx_3][k+1] /*&& lpair[idx_4][k+2]*/){
+ lr[idx][j]=MIN2(lr[idx][j], lc[idx_3][k+1]+temp->energy);/*--NU--*/
+ }
+ /*else*/ if(lpair[idx_4][k+1]){/*--NUN--*/
+ lr[idx][j]=MIN2(lr[idx][j], lc[idx_4][k+1]+temp->energy);
+ }
+ /* } */
+ temp=temp->next;
+ }
+ }
+ /* dangle 5'SIDE relative to the mRNA */
+ lc[idx][j] += E_ExtLoop(type, SS1[i-1] , SS2[j+1] , P);
+ /**
+ *** lc[idx][j] += P->dangle5[type][SS1[i-1]];
+ *** lc[idx][j] += P->dangle3[type][SS2[j+1]];
+ *** if (type>2) lc[idx][j] += P->TerminalAU;
+ **/
+ /* if(j<n2 && i>1){ */
+ /* type2=pair[S1[i-1]][S2[j+1]]; */
+ type2=lpair[idx_1][j+1];
+ if(type2>0 ){
+ 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]);
+ 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]);
+ }
+ /* } */
+ /* if(j<n2-1 && i>2){ */
+ /* type2=pair[S1[i-2]][S2[j+2]]; */
+ type2=lpair[idx_2][j+2];
+ if(type2>0 ){
+ 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]);
+ 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]);
+ /* } */
+ }
+ /* if(j<n2-2 && i>3){ */
+ /* type2 = pair[S1[i-3]][S2[j+3]]; */
+ type2 =lpair[idx_3][j+3];
+ if(type2>0 ){
+ 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]);
+ 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]);
+ /* } */
+ }
+ /* 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); */
+ int bla;
+ bla=lr[idx][j]+E_ExtLoop(rtype[type], SS2[j-1] , SS1[i+1], P)+2*penalty;
+ min_colonne=MIN2(bla, min_colonne);
+ }
+ position[i]=min_colonne;
+ if(max>=min_colonne){
+ max=min_colonne;
+ max_pos=i;
+ }
+ min_colonne=INF;
+ }
+
+ free(S1); free(S2); free(SS1); free(SS2);
+ if(max<threshTE){
+ find_max_snoop(s1, s2, max, alignment_length, position,
+ delta, distance, penalty, threshloop,
+ threshLE, threshRE, threshDE, threshTE, threshSE, threshD,
+ half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name, fullStemEnergy);
+ }
+ for (i=1; i<6; i++) {free(lc[i]);free(lr[i]);free(lpair[i]);}
+ free(lc[0]);free(lr[0]);free(lpair[0]);
+ free(lc);free(lr);free(lpair);
+ free(position);
+}
+
+
+PRIVATE void find_max_snoop(const char *s1, const char *s2,const int max, const int alignment_length, const int* position, const int delta,
+ const int distance, const int penalty, const int threshloop, const int threshLE, const int threshRE,
+ const int threshDE, const int threshTE, const int threshSE, const int threshD,
+ 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)
+{
+ int count=0;
+ int pos=n1+1;
+ int threshold = MIN2(threshTE , max + delta );
+ /* printf("threshTE %d max %d\n", threshTE, max); */
+ /* #pragma omp parallel for */
+ /* for(pos=n1+1;pos>distance;pos--){ */
+ while(pos-- > 5){
+ int temp_min=0;
+ if(position[pos]<(threshold)){
+ int search_range;
+ search_range=distance+1;
+ while(--search_range){
+ if(position[pos-search_range]<=position[pos-temp_min]){
+ temp_min=search_range;
+ }
+ }
+ pos-=temp_min;
+ int begin=MAX2(6, pos-alignment_length+1);
+ char *s3 = (char*) space(sizeof(char)*(pos-begin+3+12));
+ strcpy(s3, "NNNNN");
+ strncat(s3, (s1+begin-1), pos-begin+2);
+ strcat(s3,"NNNNN\0");
+ /* printf("%s s3\n", s3); */
+ snoopT test;
+ test = snoopfold(s3, s2, penalty, threshloop, threshLE, threshRE, threshDE, threshD, half_stem, max_half_stem, min_s2, max_s2, min_s1,
+ max_s1, min_d1, min_d2,fullStemEnergy);
+ if(test.energy==INF){
+ free(s3);
+ continue;
+ }
+ if(test.Duplex_El > threshLE * 0.01 || test.Duplex_Er > threshRE * 0.01 ||
+ test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
+ (test.Duplex_Er + test.Duplex_El + test.Loop_E + test.Loop_D + 410) > threshSE*0.01) {
+ free(test.structure);free(s3);
+ continue;
+ }
+ int l1;
+ l1 = strchr(test.structure, '&')-test.structure;
+
+ int shift=0;
+ if(test.i > strlen(s3)-10){
+ test.i--;
+ l1--;
+ }
+ if(test.i-l1<0){
+ l1--;
+ shift++;
+ }
+ char *target_struct = (char*) space(sizeof(char) * (strlen(test.structure)+1));
+ strncpy(target_struct, test.structure+shift, l1);
+ strncat(target_struct, test.structure + (strchr(test.structure, '&')-
+ test.structure), strlen(test.structure) - (strchr(test.structure, '&')-
+ test.structure));
+ strcat(target_struct,"\0");
+ char *target;
+ target = (char *) space(l1+1);
+ strncpy(target, (s3+test.i+5-l1), l1);
+ target[l1]='\0';
+ char *s4;
+ s4 = (char*) space(sizeof(char) *(strlen(s2)-9));
+ strncpy(s4, s2+5, strlen(s2)-10);
+ s4[strlen(s2)-10]='\0';
+ 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",
+ target_struct,begin + test.i-5-l1,begin + test.i -6 , begin + test.u -6,
+ test.j+1, test.j + (strrchr(test.structure,'>') - strchr(test.structure,'>'))+1 ,
+ test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10, test.Duplex_El,
+ test.Duplex_Er, test.Loop_E, test.Loop_D,test.fullStemEnergy, target,s4);
+ if(name){
+ char *temp_seq;
+ char *temp_struc;
+ char psoutput[100];
+ temp_seq = (char*) space(sizeof(char)*(l1+n2-9));
+ temp_struc = (char*) space(sizeof(char)*(l1+n2-9));
+ strcpy(temp_seq, target);
+ strcat(temp_seq, s4);
+ strncpy(temp_struc, target_struct, l1);
+ strcat(temp_struc, target_struct+l1+1);
+ temp_seq[n2+l1-10]='\0';
+ temp_struc[n2+l1-10]='\0';
+ cut_point = l1+1;
+ char str[16];char upos[16];
+ strcpy(psoutput,"sno_");
+ sprintf(str,"%d",count);
+ strcat(psoutput,str);
+ sprintf(upos,"%d",begin + test.u - 6);
+ strcat(psoutput,"_u_");
+ strcat(psoutput,upos);
+ strcat(psoutput,"_");
+ strcat(psoutput,name);
+ strcat(psoutput,".ps\0");
+ PS_rna_plot_snoop_a(temp_seq, temp_struc, psoutput, NULL, NULL);
+ cut_point = -1;
+ free(temp_seq);
+ free(temp_struc);
+ count++;
+ /* free(psoutput); */
+ }
+ free(s4);
+ free(test.structure);
+ free(target_struct);
+ free(target);
+ free(s3);
+ }
+ }
+
+}
+
+
+
+
+
+
+
+
+snoopT snoopfold(const char *s1, const char *s2,
+ const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
+ const int threshD,
+ 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 int fullStemEnergy) {
+ /* int Eminj, Emin_l; */
+ int i, j, l1, Emin=INF, i_min=0, j_min=0;
+ char *struc;
+ snoopT mfe;
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden** foldlist, **foldlist_XS;
+ int Duplex_El, Duplex_Er;
+ int Loop_D;
+ int u;
+ int Loop_E;
+ Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
+ n1 = (int) strlen(s1);
+ n2 = (int) strlen(s2);
+
+ if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
+ snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
+ make_pair_matrix();
+ }
+
+ c = (int **) space(sizeof(int *) * (n1+1));
+ r = (int **) space(sizeof(int *) * (n1+1));
+ for (i=0; i<=n1; i++) {
+ c[i] = (int *) space(sizeof(int) * (n2+1));
+ r[i] = (int *) space(sizeof(int) * (n2+1));
+ for(j=n2; j>-1; j--){
+ c[i][j]=INF;
+ r[i][j]=INF;
+ }
+ }
+ encode_seqs(s1, s2);
+ for (i=6; i<=n1-5; i++) {
+ for (j=n2-min_d2; j>min_d1; j--) {
+ int type, type2, E, k,l;
+ type = pair[S1[i]][S2[j]];
+ c[i][j] = (type ) ? P->DuplexInit : INF;
+ if(!type) continue;
+ if(/* pair[S1[i+1]][S2[j-1]] && */
+ j < max_s1 && j > min_s1 &&
+ j > n2 - max_s2 - max_half_stem &&
+ j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
+ int min_k, max_k;
+ max_k = MIN2(n2-min_s2,j+max_half_stem);
+ min_k = MAX2(j+half_stem, n2-max_s2);
+ folden * temp;
+ temp=foldlist[j+1];
+ while(temp->next){
+ int k = temp->k;
+ /* if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour */
+ if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){
+ r[i][j]=MIN2(r[i][j], c[i-3][k+1]+temp->energy);
+ }
+ /*else*/ if(pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/ ){
+ r[i][j]=MIN2(r[i][j], c[i-4][k+1]+temp->energy);
+ }
+ /* } */
+ temp=temp->next;
+ }
+ }
+ /* dangle 5'SIDE relative to the mRNA */
+ /**
+ *** c[i][j] += P->dangle5[type][SS1[i-1]];
+ *** c[i][j] += P->dangle3[type][SS2[j+1]];
+ *** if (type>2) c[i][j] += P->TerminalAU;
+ **/
+ c[i][j]+=E_ExtLoop(type, SS1[i-1] , SS2[j+1], P);
+ for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
+ for (l=j+1; l<=n2 ; l++) {
+ if (i-k+l-j>2*MAXLOOP_L-2) break;
+ if (abs(i-k-l+j) >= ASS ) continue;
+ type2 = pair[S1[k]][S2[l]];
+ if (!type2) continue;
+ E = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
+ SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
+ c[i][j] = MIN2(c[i][j], c[k][l]+E+(i-k+l-j)*penalty);
+ r[i][j] = MIN2(r[i][j], r[k][l]+E+(i-k+l-j)*penalty);
+ }
+ }
+ E = r[i][j];
+ /**
+ *** if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];
+ *** if (j>1) E += P->dangle5[rtype[type]][SS2[j-1]];
+ *** f (type>2) E += P->TerminalAU;
+ **/
+ E+=E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
+ if (E<Emin) {
+ Emin=E; i_min=i; j_min=j;
+ }
+ }
+ }
+ if(Emin > 0){
+ printf("no target found under the constraints chosen\n");
+ for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
+ free(c);
+ free(r);
+ free(S1); free(S2); free(SS1); free(SS2);
+ mfe.energy=INF;
+ return mfe;
+ }
+ struc = snoop_backtrack(i_min, j_min,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D,
+ &u, penalty, threshloop, threshLE, threshRE,threshDE, threshD,
+ half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
+/* if (i_min<n1-5) i_min++; */
+/* if (j_min>1 ) j_min--; */
+ l1 = strchr(struc, '&')-struc;
+ mfe.i = i_min-5;
+ mfe.j = j_min-5;
+ mfe.u = u -5;
+ mfe.Duplex_Er = (float) Duplex_Er/100;
+ mfe.Duplex_El = (float) Duplex_El/100;
+ mfe.Loop_D = (float) Loop_D/100;
+ mfe.Loop_E = (float) Loop_E/100;
+ mfe.energy = (float) Emin/100 ;
+ mfe.fullStemEnergy = (float) fullStemEnergy/100;
+ mfe.structure = struc;
+ if (!delay_free) {
+ for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
+ free(c);
+ free(r);
+ free(S1); free(S2); free(SS1); free(SS2);
+ }
+ return mfe;
+}
+
+PRIVATE int snoopfold_XS_fill(const char *s1, const char *s2, const int **access_s1,
+ const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
+ const int threshD,
+ 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) {
+ /* int Eminj, Emin_l; */
+ int i, j, Emin=INF, i_min=0, j_min=0;
+ /* char *struc; */
+ /* snoopT mfe; */
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden** foldlist, **foldlist_XS;
+ int Duplex_El, Duplex_Er;
+ int Loop_D;
+ /* int u; */
+ int Loop_E;
+ Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
+ n1 = (int) strlen(s1);
+ n2 = (int) strlen(s2);
+
+ if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
+ snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
+ make_pair_matrix();
+ }
+
+ c_fill = (int **) space(sizeof(int *) * (n1+1));
+ r_fill = (int **) space(sizeof(int *) * (n1+1));
+ for (i=0; i<=n1; i++) {
+ c_fill[i] = (int *) space(sizeof(int) * (n2+1));
+ r_fill[i] = (int *) space(sizeof(int) * (n2+1));
+ for(j=n2; j>-1; j--){
+ c_fill[i][j]=INF;
+ r_fill[i][j]=INF;
+ }
+ }
+ encode_seqs(s1, s2);
+
+ int di[5];
+ di[0]=0;
+ for (i=6; i<=n1-5; i++) {
+ di[1]=access_s1[5][i] - access_s1[4][i-1];
+ di[2]=access_s1[5][i-1] - access_s1[4][i-2] + di[1];
+ di[3]=access_s1[5][i-2] - access_s1[4][i-3] + di[2];
+ di[4]=access_s1[5][i-3] - access_s1[4][i-4] + di[3];
+ di[1]=MIN2(di[1],165);
+ di[2]=MIN2(di[2],330);
+ di[3]=MIN2(di[3],495);
+ di[4]=MIN2(di[4],660);
+ for (j=n2-min_d2; j>min_d1; j--) {
+ int type, type2, E, k,l;
+ type = pair[S1[i]][S2[j]];
+ c_fill[i][j] = (type ) ? P->DuplexInit : INF;
+ if(!type) continue;
+ if(/* pair[S1[i+1]][S2[j-1]] && */
+ j < max_s1 && j > min_s1 &&
+ j > n2 - max_s2 - max_half_stem &&
+ j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
+ int min_k, max_k;
+ max_k = MIN2(n2-min_s2,j+max_half_stem);
+ min_k = MAX2(j+half_stem, n2-max_s2);
+ folden * temp;
+ temp=foldlist[j+1];
+ while(temp->next){
+ int k = temp->k;
+ /* if(k >= min_k-1 && k < max_k){ uncomment to recovernormal behaviour */
+ if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){
+ r_fill[i][j]=MIN2(r_fill[i][j], c_fill[i-3][k+1]+temp->energy+ di[3]);
+ }
+ /*else*/ if(pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/ ){
+ r_fill[i][j]=MIN2(r_fill[i][j], c_fill[i-4][k+1]+temp->energy + di[4]);
+ }
+ /* } */
+ temp=temp->next;
+ }
+ }
+ /* dangle 5'SIDE relative to the mRNA */
+ /**
+ *** c_fill[i][j] += P->dangle5[type][SS1[i-1]];
+ *** c_fill[i][j] += P->dangle3[type][SS2[j+1]];
+ *** if (type>2) c_fill[i][j] += P->TerminalAU;
+ **/
+ c_fill[i][j]+= E_ExtLoop(type, SS1[i-1], SS2[j+1],P);
+ for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
+ for (l=j+1; l<=n2 ; l++) {
+ if (i-k+l-j>2*MAXLOOP_L-2) break;
+ if (abs(i-k-l+j) >= ASS ) continue;
+ type2 = pair[S1[k]][S2[l]];
+ if (!type2) continue;
+ E = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
+ SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
+ c_fill[i][j] = MIN2(c_fill[i][j], c_fill[k][l]+E+di[i-k]);
+ r_fill[i][j] = MIN2(r_fill[i][j], r_fill[k][l]+E+di[i-k]);
+ }
+ }
+ E = r_fill[i][j];
+ /**
+ *** if (i<n1) E += P->dangle3[rtype[type]][SS1[i+1]];
+ *** if (j>1) E += P->dangle5[rtype[type]][SS2[j-1]];
+ *** if (type>2) E += P->TerminalAU;
+ **/
+ E+= E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
+ if (E<Emin) {
+ Emin=E; i_min=i; j_min=j;
+ }
+ }
+ }
+ return Emin;
+}
+
+
+
+PUBLIC snoopT *snoop_subopt(const char *s1, const char *s2, int delta, int w,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
+ const int distance, 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 int fullStemEnergy) {
+
+
+
+ /* printf("%d %d\n", min_s2, max_s2); */
+ int i,j, n1, n2, E, n_subopt=0, n_max;
+ char *struc;
+ snoopT mfe;
+ snoopT *subopt;
+ int thresh;
+
+ int Duplex_El, Duplex_Er, Loop_E;
+ int Loop_D;
+ Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;
+ int u;
+ u=0;
+ n_max=16;
+ subopt = (snoopT *) space(n_max*sizeof(snoopT));
+ delay_free=1;
+ mfe = snoopfold(s1, s2, penalty, threshloop, threshLE, threshRE, threshDE,threshD,
+ half_stem, max_half_stem,
+ min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, fullStemEnergy);
+
+
+
+ if(mfe.energy > 0){
+ free(subopt);
+ delay_free=0;
+ return NULL;
+ }
+ thresh = MIN2((int) ((mfe.Duplex_Er + mfe.Duplex_El + mfe.Loop_E)*100+0.1 + 410) + delta, threshTE );
+ /* subopt[n_subopt++]=mfe; */
+ free(mfe.structure);
+
+ n1 = strlen(s1); n2=strlen(s2);
+ for (i=n1; i>0; i--) {
+ for (j=1; j<=n2; j++) {
+ int type, Ed;
+ type = pair[S2[j]][S1[i]];
+ if (!type) continue;
+ E = Ed = r[i][j];
+ /**
+ *** if (i<n1) Ed += P->dangle3[type][SS1[i+1]];
+ *** if (j>1) Ed += P->dangle5[type][SS2[j-1]];
+ *** if (type>2) Ed += P->TerminalAU;
+ **/
+ Ed+= E_ExtLoop(type, (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
+ if (Ed>thresh) continue;
+ /* too keep output small, remove hits that are dominated by a
+ better one close (w) by. For simplicity we do test without
+ adding dangles, which is slightly inaccurate.
+ */
+ /* w=1; */
+/* for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n1)) && type; ii++) { */
+/* for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n2); jj++) */
+/* if (r[ii][jj]<E) {type=0; break;} */
+/* } */
+ if (!type) continue;
+
+ struc = snoop_backtrack(i,j,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D, &u, penalty, threshloop,threshLE,threshRE,threshDE,threshD,
+ half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
+ if (Duplex_Er > threshRE || Duplex_El > threshLE || Loop_D > threshD ||
+ (Duplex_Er + Duplex_El) > threshDE ||
+ (Duplex_Er + Duplex_El + Loop_E) > threshTE ||
+ (Duplex_Er + Duplex_El + Loop_E + Loop_D + 410) > threshSE) {
+ /* printf(" Duplex_Er %d threshRE %d Duplex_El %d threshLE %d \n" */
+ /* " Duplex_Er + Duplex_El %d threshDE %d \n" */
+ /* " Duplex_Er + Duplex_El + Loop_E %d threshTE %d \n" */
+ /* " Duplex_Er + Duplex_El + Loop_E + Loop_D %d threshSE %d \n", */
+ /* Duplex_Er , threshRE , Duplex_El ,threshLE, */
+ /* Duplex_Er + Duplex_El, threshDE, */
+ /* Duplex_Er + Duplex_El+ Loop_E , threshTE, */
+ /* Duplex_Er + Duplex_El+ Loop_E + Loop_D, threshSE); */
+ Duplex_Er=0;
+ Duplex_El=0;
+ Loop_E = 0;
+ Loop_D = 0;
+ u=0,
+ free(struc);
+ continue;
+ }
+
+ if (n_subopt+1>=n_max) {
+ n_max *= 2;
+ subopt = (snoopT *) xrealloc(subopt, n_max*sizeof(snoopT));
+ }
+ subopt[n_subopt].i = i-5;
+ subopt[n_subopt].j = j-5;
+ subopt[n_subopt].u = u-5;
+ subopt[n_subopt].Duplex_Er = Duplex_Er * 0.01;
+ subopt[n_subopt].Duplex_El = Duplex_El * 0.01;
+ subopt[n_subopt].Loop_E = Loop_E * 0.01;
+ subopt[n_subopt].Loop_D = Loop_D * 0.01;
+ subopt[n_subopt].energy = (Duplex_Er +Duplex_El + Loop_E + Loop_D + 410) * 0.01 ;
+ subopt[n_subopt].fullStemEnergy = (float) fullStemEnergy * 0.01;
+ subopt[n_subopt++].structure = struc;
+
+ Duplex_Er=0; Duplex_El=0; Loop_E=0; Loop_D=0;u=0;
+ }
+ }
+
+ for (i=0; i<=n1; i++) {free(c[i]);free(r[i]);}
+ free(c);free(r);
+ free(S1); free(S2); free(SS1); free(SS2);
+ delay_free=0;
+
+ if (snoop_subopt_sorted) qsort(subopt, n_subopt, sizeof(snoopT), compare);
+ subopt[n_subopt].i =0;
+ subopt[n_subopt].j =0;
+ subopt[n_subopt].structure = NULL;
+ return subopt;
+}
+
+PUBLIC void snoop_subopt_XS(const char *s1, const char *s2, const int **access_s1, int delta, int w,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshTE, const int threshSE, const int threshD,
+ const int distance, 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 int alignment_length, const char *name, const int fullStemEnergy) {
+
+
+
+ /* printf("%d %d\n", min_s2, max_s2); */
+ int i,j, E, n_max;
+ /* char *struc; */
+ /* snoopT mfe; */
+
+ int thresh;
+
+ int Duplex_El, Duplex_Er, Loop_E;
+ int Loop_D;
+ Duplex_El=0; Duplex_Er=0; Loop_E=0;Loop_D=0;
+ int u;
+ u=0;
+ n_max=16;
+ delay_free=1;
+ int Emin = snoopfold_XS_fill(s1, s2, access_s1,penalty, threshloop, threshLE, threshRE, threshDE,threshD,
+ half_stem, max_half_stem,
+ min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
+ if(Emin > 0){
+ delay_free=0;
+ }
+ thresh = MIN2(-100, threshTE +alignment_length*30);
+ /* n1=strlen(s1); */
+ /* n2=strlen(s2); */
+
+ int n3=strlen(s1);
+ int n4=strlen(s2);
+ S1_fill = (short*)space(sizeof(short)*(n3+2));
+ S2_fill = (short*)space(sizeof(short)*(n4+2));
+ SS1_fill = (short*)space(sizeof(short)*(n3+1));
+ SS2_fill = (short*)space(sizeof(short)*(n4+1));
+ memcpy(S1_fill, S1, sizeof(short)*n3+2);
+ memcpy(S2_fill, S2, sizeof(short)*n4+2);
+ memcpy(SS1_fill, SS1, sizeof(short)*n3+1);
+ memcpy(SS2_fill, SS2, sizeof(short)*n4+1);
+ free(S1);free(S2);free(SS1);free(SS2);
+ int count=0;
+ for (i=n3-5; i>0; i--) {
+ for (j=1; j<=n4; j++) {
+ int type, Ed;
+ type = pair[S2_fill[j]][S1_fill[i]];
+ if (!type) continue;
+ E = Ed = r_fill[i][j];
+ /**
+ ***if (i<n3) Ed += P->dangle3[type][SS1_fill[i+1]];
+ ***if (j>1) Ed += P->dangle5[type][SS2_fill[j-1]];
+ ***if (type>2) Ed += P->TerminalAU;
+ **/
+ Ed+=E_ExtLoop(type, (j > 1) ? SS2[j-1] : -1, (i<n3) ? SS1[i+1] : -1, P);
+ if (Ed>thresh) continue;
+
+ /* to keep output small, remove hits that are dominated by a
+ better one close (w) by. For simplicity we do test without
+ adding dangles, which is slightly inaccurate.
+ */
+/* w=10; */
+/* for (ii=MAX2(i-w,1); (ii<=MIN2(i+w,n3-5)) && type; ii++) { */
+/* for (jj=MAX2(j-w,1); jj<=MIN2(j+w,n4-5); jj++) */
+/* if (r_fill[ii][jj]<E) {type=0; break;} */
+/* } */
+/* i=ii;j=jj; */
+ if (!type) continue;
+ int begin=MAX2(5, i-alignment_length);
+ int end =MIN2(n3-5, i-1);
+ char *s3 = (char*) space(sizeof(char)*(end-begin+2)+5);
+ strncpy(s3, (s1+begin), end - begin +1);
+ strcat(s3,"NNNNN\0");
+ int n5 = strlen(s3);
+ snoopT test = snoopfold_XS(s3, s2, access_s1, i, j ,penalty,
+ threshloop, threshLE, threshRE,
+ threshDE, threshD, half_stem,
+ max_half_stem, min_s2, max_s2, min_s1,
+ max_s1, min_d1, min_d2,fullStemEnergy);
+ if(test.energy==INF){
+ free(s3);
+ continue;
+ }
+ if( test.Duplex_El > threshLE * 0.01 ||test.Duplex_Er > threshRE * 0.01 ||
+ test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
+ (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)
+ {
+ free(test.structure);free(s3);
+ continue;
+ }
+ char *s4;
+ s4 = (char*) space(sizeof(char) *(n4-9));
+ strncpy(s4, s2+5, n4-10);
+ s4[n4-10]='\0';
+
+ char *s5 = space(sizeof(char) * n5-test.i+2-5);
+ strncpy(s5,s3+test.i-1,n5-test.i+1-5);
+ s5[n5-test.i+1-5]='\0';
+ float dE = ((float) (access_s1[n5-test.i+1-5][i]))*0.01;
+ 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" ,
+ test.structure, i - (n5 - test.i) ,i - 5, i - (n5 - test.u ),
+ j-5, j-5 + (strrchr(test.structure,'>') - strchr(test.structure,'>')),
+ test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10+dE, test.Duplex_El,
+ test.Duplex_Er, test.Loop_E, test.Loop_D,dE , test.fullStemEnergy, s5,s4);
+ if(name){
+ int begin_t, end_t, begin_q, end_q, and, pipe,k;
+ char psoutput[100];
+ begin_q=0;
+ end_q=n4-10;
+ begin_t=0;
+ end_t=n5-test.i+ 1-5;
+ and=end_t+1;
+ pipe=test.u -test.i + 1;
+ cut_point = end_t +1 ;
+ char *catseq, *catstruct;/* *fname; */
+ catseq = (char*) space(n5 + end_q -begin_q +2);
+ catstruct = (char*) space(n5 + end_q -begin_q +2);
+ strcpy(catseq, s5);
+ strncpy(catstruct, test.structure, end_t);
+ strcat(catseq, s4);
+ strncat(catstruct, test.structure+end_t+1, end_q-begin_q+1);
+ catstruct[end_t - begin_t + end_q-begin_q+2]='\0';
+ catseq[end_t - begin_t + end_q-begin_q+2]='\0';
+ int *relative_access;
+ relative_access = space(sizeof(int)*strlen(s5));
+ relative_access[0] = access_s1[1][i - (n5 - test.i) + 5];
+ for(k=1;k<strlen(s5);k++){
+ relative_access[k] = access_s1[k+1][i - (n5 - test.i) + k + 5] - access_s1[k][i - (n5 - test.i) + k + 4];
+ }
+ char str[16];char upos[16];
+ strcpy(psoutput,"sno_XS_");
+ sprintf(str,"%d",count);
+ strcat(psoutput,str);
+ sprintf(upos,"%d",i - (n5 - test.u ));
+ strcat(psoutput,"_u_");
+ strcat(psoutput,upos);
+ strcat(psoutput,"_");
+ strcat(psoutput,name);
+ strcat(psoutput,".ps\0");
+ PS_rna_plot_snoop_a(catseq, catstruct, psoutput,relative_access,NULL);
+ free(catseq);free(catstruct);free(relative_access);
+ count++;
+ }
+ free(s3);free(s4);free(s5);free(test.structure);
+ }
+ }
+ for (i=0; i<=n3; i++) {free(c_fill[i]);free(r_fill[i]);}
+ free(c_fill);free(r_fill);
+ free(S1_fill); free(S2_fill); free(SS1_fill); free(SS2_fill);
+ delay_free=0;
+}
+
+
+
+
+PRIVATE char *snoop_backtrack(int i, int j, const char* snoseq,
+ int *Duplex_El, int *Duplex_Er,
+ int *Loop_E, int *Loop_D, int *u,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshD,
+ 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) {
+ /* backtrack structure going backwards from i, and forwards from j
+ return structure in bracket notation with & as separator */
+ int k, l, type, type2, E, traced, i0, j0;
+ int traced_r=0; /* flag for following backtrack in c or r */
+ char *st1, *st2, *struc;
+ char *struc_loop;
+
+ st1 = (char *) space(sizeof(char)*(n1+1));
+ st2 = (char *) space(sizeof(char)*(n2+1));
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden **foldlist, **foldlist_XS;
+ type=pair[S1[i]][S2[j]];
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
+ i0=i; j0=j;
+ /**
+ *** if (i<n1) *Duplex_Er += P->dangle3[rtype[type]][SS1[i+1]];
+ *** if (j>1) *Duplex_Er += P->dangle5[rtype[type]][SS2[j-1]];
+ *** if (type>2) *Duplex_Er += P->TerminalAU;
+ **/
+ *Duplex_Er += E_ExtLoop(rtype[type], (j > 1) ? SS2[j-1] : -1, (i<n1) ? SS1[i+1] : -1, P);
+ while (i>0 && j<=n2-min_d2 ) {
+ if(!traced_r) {
+ E = r[i][j]; traced=0;
+ st1[i-1] = '<';
+ st2[j-1] = '>';
+ type = pair[S1[i]][S2[j]];
+ if (!type) nrerror("backtrack failed in fold duplex r");
+ for (k=i-1; k>0 && (i-k)<MAXLOOP_L; k--) {
+ for (l=j+1; l<=n2 ; l++) {
+ int LE;
+ if (i-k+l-j>2*MAXLOOP_L-2) break;
+ if (abs(i-k-l+j) >= ASS) continue;
+
+ type2 = pair[S1[k]][S2[l]];
+ if (!type2) continue;
+ LE = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
+ SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
+ if (E == r[k][l]+LE+(i-k+l-j)*penalty) {
+ traced=1;
+ i=k; j=l;
+ *Duplex_Er+=LE;
+ break;
+ }
+ }
+ if (traced) break;
+ }
+ if(!traced){
+ if(/* pair[S1[i+1]][S2[j-1]] && */
+ j < max_s1 && j > min_s1 &&
+ j > n2 - max_s2 - max_half_stem &&
+ j < n2 -min_s2 -half_stem &&
+ S1[i-2]==4) {
+ int min_k, max_k;
+ max_k = MIN2(n2-min_s2,j+max_half_stem+1);
+ min_k = MAX2(j+half_stem+1, n2-max_s2);
+ folden * temp;
+ temp=foldlist[j+1];
+ while(temp->next) {
+ int k = temp->k;
+ if(pair[S1[i-3]][S2[k+1]] /*&& pair[S1[i-4]][S2[k+2]]*/ ){ /* introduce structure from RNAfold */
+ if(E==c[i-3][k+1]+temp->energy){
+ *Loop_E=temp->energy;
+ st1[i-3]= '|';
+ *u=i-2;
+ int a,b;
+ /* int fix_ij=indx[k-1+1]+j+1; */
+ for(a=0; a< MISMATCH ;a++){
+ for(b=0; b< MISMATCH ; b++){
+ int ij=indx[k-1-a+1]+j+1+b;
+ if(cLoop[ij]==temp->energy) {
+ struc_loop=snobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1);
+ a=INF; b=INF;
+ }
+ }
+ }
+ traced=1;
+ traced_r=1;
+ i=i-3;j=k+1;
+ break;
+ }
+ }
+ /*else*/ if (pair[S1[i-4]][S2[k+1]] /*&& pair[S1[i-5]][S2[k+2]]*/){ /* introduce structure from RNAfold */
+ if(E==c[i-4][k+1]+temp->energy){
+ *Loop_E=temp->energy;
+ st1[i-3]= '|';
+ *u=i-2;
+ int a,b;
+ /* int fix_ij=indx[k-1+1]+j+1; */
+ for(a=0; a< MISMATCH ;a++){
+ for(b=0; b< MISMATCH ; b++){
+ int ij=indx[k-1-a+1]+j+1+b;
+ if(cLoop[ij]==temp->energy) {
+ struc_loop=snobacktrack_fold_from_pair(snoseq, j+1+b, k-a-1+1);
+ a=INF; b=INF;
+ }
+ }
+ }
+ traced=1;
+ traced_r=1;
+ i=i-4;j=k+1;
+ break;
+ }
+ } /* else if */
+ temp=temp->next;
+ } /* while temp-> next */
+ } /* test on j */
+ }/* traced? */
+ }/* traced_r? */
+ else{
+ E = c[i][j]; traced=0;
+ st1[i-1] = '<';
+ st2[j-1] = '>';
+ type = pair[S1[i]][S2[j]];
+ if (!type) nrerror("backtrack failed in fold duplex c");
+ for (k=i-1; (i-k)<MAXLOOP_L; k--) {
+ for (l=j+1; l<=n2; l++) {
+ int LE;
+ if (i-k+l-j>2*MAXLOOP_L-2) break;
+ if (abs(i-k-l+j) >= ASS) continue;
+ type2 = pair[S1[k]][S2[l]];
+ if (!type2) continue;
+ LE = E_IntLoop(i-k-1, l-j-1, type2, rtype[type],
+ SS1[k+1], SS2[l-1], SS1[i-1], SS2[j+1],P);
+ if (E == c[k][l]+LE+(i-k+l-j)*penalty) {
+ traced=1;
+ i=k; j=l;
+ *Duplex_El+=LE;
+ break;
+ }
+ }
+ if (traced) break;
+ }
+ }
+ if (!traced) {
+ int correction;
+ correction = E_ExtLoop(type, (i>1) ? SS1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, P);
+ E-=correction;
+ *Duplex_El+=correction;
+ /**
+ *** if (i>1) {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];}
+ *** if (j<n2) {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];}
+ *** if (type>2) {E -= P->TerminalAU; *Duplex_El +=P->TerminalAU;}
+ **/
+ if (E != P->DuplexInit) {
+ nrerror("backtrack failed in fold duplex end");
+ } else break;
+ }
+ }
+/* if (i>1) i--; */
+/* if (j<n2) j++; */
+ /* struc = (char *) space(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
+ struc = (char *) space(i0-i+1+n2-1+1+2); /* declare final duplex structure */
+ char * struc2;
+ struc2 = (char *) space(n2+1);
+ /* char * struct_const; */
+ for (k=MAX2(i,1); k<=i0; k++) if (!st1[k-1]) st1[k-1] = '.';
+ /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
+ /* char * struct_const; */
+ /* struct_const = (char *) space(sizeof(char)*(n2+1)); */
+ for (k=1; k<=n2; k++) {
+ if (!st2[k-1]) st2[k-1] = struc_loop[k-1];/* '.'; */
+ struc2[k-1] = st2[k-1];/* '.'; */
+ /* if (k>=j0 && k<=j){ */
+ /* struct_const[k-1]='x'; */
+ /* } */
+ /* else{ */
+ /* if(k<j0) {struct_const[k-1]='<';} */
+ /* if(k>j) {struct_const[k-1]='>';} */
+ /* } */
+ }
+ char duplexseq_1[j0];
+ char duplexseq_2[n2-j+2];
+ if(j<n2){
+ strncpy(duplexseq_1, snoseq, j0-1);
+ strcpy(duplexseq_2, snoseq + j);
+ duplexseq_1[j0-1]='\0';
+ duplexseq_2[n2-j+1]='\0';
+ duplexT temp;
+ temp=duplexfold(duplexseq_1, duplexseq_2);
+ *Loop_D = MIN2(0,-410 + (int) 100 * temp.energy);
+ if(*Loop_D){
+ int l1,ibegin, iend, jbegin, jend;
+ l1=strchr(temp.structure, '&')-temp.structure;
+ ibegin=temp.i-l1;
+ iend =temp.i-1;
+ jbegin=temp.j;
+ jend =temp.j+strlen(temp.structure)-l1-2-1;
+ for(k=ibegin+1; k<=iend+1; k++){
+ struc2[k-1]=temp.structure[k-ibegin-1];
+ }
+ for(k=jbegin+j; k<=jend+j; k++){
+ struc2[k-1]=temp.structure[l1+k-j-jbegin+1];
+ }
+ }
+ free(temp.structure);
+ }
+ strcpy(struc, st1+MAX2(i-1,0)); strcat(struc, "&");
+ /* strcat(struc, st2); */
+ strncat(struc, struc2+5, strlen(struc2)-10);
+ free(struc2);
+ free(struc_loop);
+ free(st1); free(st2);
+ /* free_arrays(); */
+ return struc;
+}
+
+void Lsnoop_subopt_list_XS(const char *s1, const char *s2, const int **access_s1, int delta, int w,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshTE,const int threshSE,const int threshD,
+ const int distance,
+ 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 int alignment_length, const char *name, const int fullStemEnergy)
+{
+
+ int min_colonne=INF;
+ int max_pos;
+ int max;max=INF;
+ /* int temp; */
+ /* int nsubopt=10; */
+ n1 = (int) strlen(s1);
+ n2 = (int) strlen(s2);
+ int *position;
+ int *position_j;
+ int min_j_colonne;
+ int max_pos_j=INF;
+ position = (int*) space((n1+3)*sizeof(int));
+ position_j = (int*) space((n1+3)*sizeof(int));
+
+ /* int Eminj, Emin_l; */
+ int i, j;/* l1, Emin=INF, i_min=0, j_min=0; */
+ /* char *struc; */
+ /* snoopT mfe; */
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden **foldlist, **foldlist_XS;
+ int Duplex_El, Duplex_Er;
+ int Loop_D;
+ /* int u; */
+ int Loop_E;
+ Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop, &foldlist, &foldlist_XS);
+ if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
+ snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
+ make_pair_matrix();
+ }
+
+ lpair = (int **) space(sizeof(int *) * (6));
+ lc = (int **) space(sizeof(int *) * (6));
+ lr = (int **) space(sizeof(int *) * (6));
+ for (i=0; i<6; i++) {
+ lc[i] = (int *) space(sizeof(int) * (n2+1));
+ lr[i] = (int *) space(sizeof(int) * (n2+1));
+ lpair[i] = (int *) space(sizeof(int) * (n2+1));
+ for(j=n2; j>-1; j--){
+ lc[i][j]=INF;
+ lr[i][j]=INF;
+ lpair[i][j]=0;
+ }
+ }
+ encode_seqs(s1, s2);
+ int lim_maxj=n2-min_d2 ;
+ int lim_minj=min_d1;
+ int lim_maxi=n1-5;
+ for (i=5; i<=lim_maxi; i++) {
+ int idx=i%5;
+ int idx_1=(i-1)%5;
+ int idx_2=(i-2)%5;
+ int idx_3=(i-3)%5;
+ int idx_4=(i-4)%5;
+ int di1, di2, di3, di4;
+ di1 = access_s1[5][i] - access_s1[4][i-1];
+ di2 =access_s1[5][i-1] - access_s1[4][i-2] + di1;
+ di3 = access_s1[5][i-2] - access_s1[4][i-3] + di2;
+ di4 = access_s1[5][i-3] - access_s1[4][i-4] + di3;
+ di1=MIN2(di1,165);
+ di2=MIN2(di2,330);
+ di3=MIN2(di3,495);
+ di4=MIN2(di4,660);
+ for (j=lim_maxj; j>lim_minj; j--) {
+ int type, type2;/* E, k,l; */
+ type = pair[S1[i]][S2[j]];
+ lpair[idx][j] = type;
+ lc[idx][j] = (type) ? P->DuplexInit + access_s1[1][i] : INF;
+ lr[idx][j] = INF;
+ if(!type) continue;
+ if( /*pair[S1[i+1]][S2[j-1]] && Be sure it binds*/
+ j < max_s1 && j > min_s1 &&
+ j > n2 - max_s2 - max_half_stem &&
+ j < n2 -min_s2 -half_stem && S1[i-2]==4 ) { /*constraint on s2 and i*/
+ int min_k, max_k;
+ max_k = MIN2(n2-min_s2,j+max_half_stem+1);
+ min_k = MAX2(j+half_stem+1, n2-max_s2);
+ folden * temp;
+ temp=foldlist[j+1];
+ while(temp->next){
+ int k = temp->k;
+ /* if(k >= min_k-1 && k < max_k){ comment to recover normal behaviour */
+ if(lpair[idx_3][k+1] && lc[idx_3][k+1] /*+di3*/ < 411 /*&& lpair[idx_4][k+2]*/){ /* remove second condition */
+ lr[idx][j]=MIN2(lr[idx][j], di3 + lc[idx_3][k+1]+temp->energy);/*--NU--*/
+ }
+ /*else*/ if(lpair[idx_4][k+1] && /*di4 +*/ lc[idx_4][k+1] < 411 ){/*--NUN--*/ /* remove second condition */
+ lr[idx][j]=MIN2(lr[idx][j], di4 + lc[idx_4][k+1]+temp->energy);
+ }
+ /* } */
+ temp=temp->next;
+ }
+ }
+ /* dangle 5'SIDE relative to the mRNA */
+ /**
+ *** lc[idx][j] += P->dangle5[type][SS1[i-1]];
+ *** lc[idx][j] += P->dangle3[type][SS2[j+1]];
+ *** if (type>2) lc[idx][j] += P->TerminalAU;
+ **/
+ lc[idx][j]+=E_ExtLoop(type, SS1[i-1] , SS2[j+1] , P);
+ /* if(j<n2 && i>1){ */
+ /* type2=pair[S1[i-1]][S2[j+1]]; */
+ type2=lpair[idx_1][j+1];
+ if(type2>0 ){
+ 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]);
+ 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]);
+ }
+ type2=lpair[idx_2][j+2];
+ if(type2>0 ){
+ 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]);
+ 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]);
+
+ }
+ type2 =lpair[idx_3][j+3];
+ if(type2>0 ){
+ 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]);
+ 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]);
+
+ }
+ int bla;
+ int temp2;
+ temp2=min_colonne;
+ bla=lr[idx][j] + E_ExtLoop(rtype[type], SS2[j-1], SS1[i+1] , P);
+ /**
+ *** (type>2?P->TerminalAU:0)+P->dangle3[rtype[type]][SS1[i+1]]+P->dangle5[rtype[type]][SS2[j-1]];
+ **/
+ min_colonne=MIN2(bla, min_colonne);
+ if(temp2>min_colonne){
+ min_j_colonne=j;
+ }
+ }
+ position[i]=min_colonne;
+ if(max>=min_colonne){
+ max=min_colonne;
+ max_pos=i;
+ max_pos_j=min_j_colonne;
+ }
+ position_j[i]=min_j_colonne;
+ min_colonne=INF;
+ }
+ free(S1); free(S2); free(SS1); free(SS2);
+
+ if(max<threshTE + 30*alignment_length){
+ find_max_snoop_XS(s1, s2, access_s1,max,alignment_length, position, position_j,
+ delta, distance, penalty, threshloop,
+ threshLE, threshRE, threshDE, threshTE, threshSE, threshD,
+ half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,name,fullStemEnergy);
+ }
+ for (i=1; i<6; i++) {free(lc[i]);free(lr[i]);free(lpair[i]);}
+ free(lc[0]);free(lr[0]);free(lpair[0]);
+ free(lc);free(lr);free(lpair);
+ free(position);free(position_j);
+}
+
+
+PRIVATE void find_max_snoop_XS(const char *s1, const char *s2, const int **access_s1,
+ const int max, const int alignment_length,
+ const int* position, const int* position_j, const int delta,
+ const int distance, const int penalty, const int threshloop, const int threshLE, const int threshRE,
+ const int threshDE, const int threshTE, const int threshSE, const int threshD,
+ 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){
+ int count=0;
+ int n3=strlen(s1);
+ int n4=strlen(s2);
+ int pos=n1-4;
+ int max_pos_j;
+ int threshold = MIN2(threshTE + alignment_length*30, -100);
+ /* printf("threshTE %d max %d\n", threshTE, max); */
+ /* #pragma omp parallel for */
+ /* for(pos=n1+1;pos>distance;pos--){ */
+ while(pos-->5){
+ int temp_min=0;
+ if(position[pos]<(threshold)){
+ int search_range;
+ search_range=distance+1;
+ while(--search_range){
+ if(position[pos-search_range]<=position[pos-temp_min]){
+ temp_min=search_range;
+ }
+ }
+ pos-=temp_min;
+ max_pos_j=position_j[pos];
+ int begin=MAX2(5, pos-alignment_length);
+ int end =MIN2(n3-5, pos-1);
+ char *s3 = (char*) space(sizeof(char)*(end-begin+2)+5);
+ strncpy(s3, (s1+begin), end - begin +1);
+ strcat(s3,"NNNNN\0");
+
+ int n5 = strlen(s3);
+ snoopT test;
+ test = snoopfold_XS(s3, s2, access_s1, pos, max_pos_j ,penalty,
+ threshloop, threshLE, threshRE,
+ threshDE, threshD, half_stem,
+ max_half_stem, min_s2, max_s2, min_s1,
+ max_s1, min_d1, min_d2, fullStemEnergy);
+ if(test.energy==INF){
+ free(s3);
+ continue;
+ }
+ if( test.Duplex_El > threshLE * 0.01 ||test.Duplex_Er > threshRE * 0.01 ||
+ test.Loop_D > threshD * 0.01 || (test.Duplex_Er + test.Duplex_El) > threshDE * 0.01 ||
+ (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) {
+ free(test.structure);free(s3);
+ continue;
+ }
+
+ char *s4;
+ s4 = (char*) space(sizeof(char) *(n4-9));
+ strncpy(s4, s2+5, n4-10);
+ s4[n4-10]='\0';
+
+ char *s5 = space(sizeof(char) * n5-test.i+2-5);
+ strncpy(s5,s3+test.i-1,n5-test.i+1-5);
+ s5[n5-test.i+1-5]='\0';
+ float dE = ((float) (access_s1[n5-test.i+1-5][pos]))*0.01;
+ 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" ,
+ test.structure, pos - (n5 - test.i) ,pos - 5, pos - (n5 - test.u ),
+ max_pos_j-5, max_pos_j-5 + (strrchr(test.structure,'>') - strchr(test.structure,'>')),
+ test.Loop_D + test.Duplex_El + test.Duplex_Er + test.Loop_E + 4.10+dE, test.Duplex_El,
+ test.Duplex_Er, test.Loop_E, test.Loop_D,dE ,test.fullStemEnergy, s5,s4);
+ if(name){
+ int begin_t, end_t, begin_q, end_q, and, pipe, i;
+ char psoutput[100];
+ begin_q=0;
+ end_q=n4-10;
+ begin_t=0;
+ end_t=n5-test.i+ 1-5;
+ and=end_t+1;
+ pipe=test.u -test.i + 1;
+ cut_point = end_t +1 ;
+ char *catseq, *catstruct;/* *fname; */
+ catseq = (char*) space(n5 + end_q -begin_q +2);
+ catstruct = (char*) space(n5 + end_q -begin_q +2);
+ strcpy(catseq, s5);
+ strncpy(catstruct, test.structure, end_t);
+ strcat(catseq, s4);
+ strncat(catstruct, test.structure+end_t+1, end_q-begin_q+1);
+ catstruct[end_t - begin_t + end_q-begin_q+2]='\0';
+ catseq[end_t - begin_t + end_q-begin_q+2]='\0';
+ int *relative_access;
+ relative_access = space(sizeof(int)*strlen(s5));
+
+ relative_access[0] = access_s1[1][pos - (n5 - test.i) + 5];
+ for(i=1;i<strlen(s5);i++){
+ relative_access[i] = access_s1[i+1][pos - (n5 - test.i) + i + 5] - access_s1[i][pos - (n5 - test.i) + i + 4];
+ }
+ char str[16];
+ char upos[16];
+ strcpy(psoutput,"sno_XS_");
+ sprintf(str,"%d",count);
+ strcat(psoutput,str);
+ sprintf(upos,"%d",pos - (n5 - test.u ));
+ strcat(psoutput,"_u_");
+ strcat(psoutput,upos);
+ strcat(psoutput,"_");
+ strcat(psoutput,name);
+ strcat(psoutput,".ps\0");
+ PS_rna_plot_snoop_a(catseq, catstruct, psoutput,relative_access,NULL);
+ free(catseq);free(catstruct);free(relative_access);
+ count++;
+ }
+ free(s3);free(s4);free(s5);free(test.structure);
+ }
+ }
+}
+
+snoopT snoopfold_XS(const char *s1, const char *s2, const int **access_s1, const int pos_i, const int pos_j,
+ const int penalty, const int threshloop, const int threshLE, const int threshRE, const int threshDE,
+ const int threshD,
+ 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 int fullStemEnergy) {
+ /* int Eminj, Emin_l; */
+ int a,b,i, j, Emin=INF, a_min=0, b_min=0;
+ char *struc;
+ snoopT mfe;
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden** foldlist, **foldlist_XS;
+ int Duplex_El, Duplex_Er;
+ int Loop_D;
+ int u;
+ int Loop_E;
+
+ Duplex_El=0;Duplex_Er=0;Loop_E=0, Loop_D=0;
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
+ n1 = (int) strlen(s1);
+ n2 = (int) strlen(s2);
+
+ if ((!P) || (fabs(P->temperature - temperature)>1e-6)) {
+ snoupdate_fold_params(); if(P) free(P); P = scale_parameters();
+ make_pair_matrix();
+ }
+
+ c = (int **) space(sizeof(int *) * (n1+1));
+ r = (int **) space(sizeof(int *) * (n1+1));
+ for (i=0; i<=n1; i++) {
+ c[i] = (int *) space(sizeof(int) * (n2+1));
+ r[i] = (int *) space(sizeof(int) * (n2+1));
+ for(j=n2; j>-1; j--){
+ c[i][j]=INF;
+ r[i][j]=INF;
+ }
+ }
+ encode_seqs(s1, s2);
+ i=n1-5;
+ j=pos_j;
+ /* printf("tar: %s\nsno: %s\n ", s1, s2); */
+ /* printf("pos_i %d pos_j %d\n", pos_i, pos_j); */
+ /* printf("type %d n1 %d n2 %d S1[n1] %d S2[n2] %d", pair[S1[i]][S2[j]], n1, n2, S1[n1], S2[n2]); */
+ int type, type2, E, p,q;
+ r[i][j] = P->DuplexInit;
+ /* r[i][j] += P->dangle3[rtype[type]][SS1[i+1]] + P->dangle5[rtype[type]][SS2[j-1]]; */
+
+ if(pair[S1[i]][S2[j]]>2) r[i][j]+=P->TerminalAU;
+ for(a=i-1; a>0;a--){ /* i-1 */
+ r[a+1][0]=INF;
+ for(b=j+1; b<=n2-min_d2;b++){ /* j+1 */
+ r[a][b]=INF;
+ type = pair[S1[a]][S2[b]];
+ if(!type) continue;
+ if(S1[a+1]==4){
+ folden * temp;
+ temp=foldlist_XS[b-1];
+ while(temp->next) {
+ int k = temp->k;
+ 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 */
+ c[a][b]=MIN2(c[a][b], r[a+3][k-1]+temp->energy);
+ }
+ temp=temp->next;
+ }
+ }
+ if(S1[a+2]==4){
+ folden * temp;
+ temp=foldlist_XS[b-1];
+ while(temp->next){
+ int k = temp->k;
+ 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 */
+ c[a][b]=MIN2(c[a][b], r[a+4][k-1]+temp->energy);
+ }
+ temp=temp->next;
+ }
+ }
+ for(p=a+1; p<n1 && (p-a) < MAXLOOP_L; p++) { /* p < n1 */
+ for (q=b-1; q > 1 ; q--) { /* q > 1 */
+ if (p-a+b-q>2*MAXLOOP_L-2) break;
+ if (abs((p-a)-(b-q)) >= ASS ) continue;
+ type2 = pair[S1[p]][S2[q]];
+ if (!type2) continue;
+ 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);
+ c[a][b] = MIN2(c[a][b], c[p][q]+E);
+ r[a][b] = MIN2(r[a][b], r[p][q]+E);
+ }
+ }
+ E = c[a][b];
+ if (type>2) E += P->TerminalAU;
+ /* E +=P->dangle5[rtype[type]][SS1[i+1]]; */
+ /* E +=P->dangle5[rtype[type]][SS2[j-1]]; */
+ E+=access_s1[i-a+1][pos_i];
+ if (E<Emin) {
+ Emin=E; a_min=a; b_min=b;
+ }
+ }
+ }
+ if(Emin > 0){
+ printf("no target found under the constraints chosen\n");
+ for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
+ free(c);
+ free(r);
+ free(S1); free(S2); free(SS1); free(SS2);
+ mfe.energy=INF;
+ return mfe;
+ }
+ type2=pair[S1[a_min]][S2[b_min]];
+ if(type2>2) Emin +=P->TerminalAU;
+ mfe.energy = ((float) (Emin))/100;
+ struc = snoop_backtrack_XS(a_min, b_min,s2, &Duplex_El, &Duplex_Er, &Loop_E, &Loop_D,
+ &u, penalty, threshloop, threshLE, threshRE,threshDE, threshD,
+ half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2);
+
+ mfe.i = a_min;
+ mfe.j = b_min;
+ mfe.u = u;
+ mfe.Duplex_Er = (float) Duplex_Er/100;
+ mfe.Duplex_El = (float) Duplex_El/100;
+ mfe.Loop_D = (float) Loop_D/100;
+ mfe.Loop_E = (float) Loop_E/100;
+ mfe.energy = (float) Emin/100 ;
+ mfe.fullStemEnergy = (float) fullStemEnergy/100;
+ mfe.structure = struc;
+ return mfe;
+}
+
+PRIVATE char *snoop_backtrack_XS(int i, int j, const char* snoseq,
+ int *Duplex_El, int *Duplex_Er,
+ int *Loop_E, int *Loop_D, int *u,
+ const int penalty, const int threshloop,
+ const int threshLE, const int threshRE, const int threshDE, const int threshD,
+ 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) {
+ /* backtrack structure going backwards from i, and forwards from j
+ return structure in bracket notation with & as separator */
+ int k, l, type, type2, E, traced, i0, j0;
+ int traced_c=0; /* flag for following backtrack in c or r */
+ char *st1, *st2, *struc;
+ char *struc_loop;
+
+ st1 = (char *) space(sizeof(char)*(n1+1));
+ st2 = (char *) space(sizeof(char)*(n2+1));
+ int *indx;
+ int *mLoop;
+ int *cLoop;
+ folden **foldlist, **foldlist_XS;
+ type=pair[S1[i]][S2[j]];
+ snoexport_fold_arrays(&indx, &mLoop, &cLoop,&foldlist, &foldlist_XS );
+ i0=i;j0=j;
+ /* i0=MAX2(i,1); j0=MIN2(j+1,n2); */
+ while (i<=n1 && j>=1 ) {
+ if(!traced_c) {
+ E = c[i][j]; traced=0;
+ st1[i] = '<';
+ st2[j-1] = '>';
+ type = pair[S1[i]][S2[j]];
+ if (!type) nrerror("backtrack failed in fold duplex c");
+ for (k=i+1; k>0 && (k-i)<MAXLOOP_L; k++) {
+ for (l=j-1; l>=1 ; l--) {
+ int LE;
+ if (k-i+j-l>2*MAXLOOP_L-2) break;
+ if (abs(k-i-j+l) >= ASS) continue;
+ type2 = pair[S1[k]][S2[l]];
+ if (!type2) continue;
+ LE = E_IntLoop(k-i-1, j-l-1, type2, rtype[type],
+ SS1[i+1], SS2[j-1], SS1[k-1], SS2[l+1],P);
+ if (E == c[k][l]+LE) {
+ traced=1;
+ i=k; j=l;
+ *Duplex_El+=LE;
+ break;
+ }
+ }
+ if (traced) break;
+ }
+ if(!traced){
+ if(S1[i+1]==4) {
+ folden * temp;
+ temp=foldlist_XS[j-1];
+ while(temp->next) {
+ int k = temp->k;
+ 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 ) {
+ if(E==r[i+3][k-1]+temp->energy){
+ *Loop_E=temp->energy;
+ st1[i+1]= '|';
+ st1[i+2]='.';
+ *u=i+1;
+ int a,b;
+ for(a=0; a< MISMATCH ;a++){
+ for(b=0; b< MISMATCH ; b++){
+ int ij=indx[j-1-a]+k+b;
+ if(cLoop[ij]==temp->energy) {
+ struc_loop=snobacktrack_fold_from_pair(snoseq, k+b, j-1-a);
+ a=INF; b=INF;
+ }
+ }
+ }
+ traced=1;
+ traced_c=1;
+ i=i+3;j=k-1;
+ break;
+ }
+ }
+ temp=temp->next;
+ }
+ }
+ if (S1[i+2]==4){ /* introduce structure from RNAfold */
+ folden * temp;
+ temp=foldlist_XS[j-1];
+ while(temp->next) {
+ int k = temp->k;
+ 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 ) {
+ if(E==r[i+4][k-1]+temp->energy){
+ *Loop_E=temp->energy;
+ st1[i+2]= '|';
+ st1[i+1]=st1[i+3]='.';
+ *u=i+2;
+ int a,b;
+ for(a=0; a< MISMATCH ;a++){
+ for(b=0; b< MISMATCH ; b++){
+ int ij=indx[j-1-a]+k+b;
+ if(cLoop[ij]==temp->energy) {
+ struc_loop=snobacktrack_fold_from_pair(snoseq, k+b, j-a-1);
+ a=INF; b=INF;
+ }
+ }
+ }
+ traced=1;
+ traced_c=1;
+ i=i+4;j=k-1;
+ break;
+ }
+ }
+ temp=temp->next;
+ }
+ }
+ }/* traced? */
+ }/* traced_r? */
+ else{
+ E = r[i][j]; traced=0;
+ st1[i] = '<';
+ st2[j-1] = '>';
+ type = pair[S1[i]][S2[j]];
+ if (!type) nrerror("backtrack failed in fold duplex r");
+ for (k=i+1; k>0 && (k-i)<MAXLOOP_L; k++) {
+ for (l=j-1; l>=1 ; l--) {
+ int LE;
+ if (k-i+j-l>2*MAXLOOP_L-2) break;
+ if (abs(k-i-j+l) >= ASS) continue;
+ type2 = pair[S1[k]][S2[l]];
+ if (!type2) continue;
+ LE = E_IntLoop(k-i-1, j-l-1, type2, rtype[type],
+ SS1[i+1], SS2[j-1], SS1[k-1], SS2[l+1],P);
+ if (E == r[k][l]+LE) {
+ traced=1;
+ i=k; j=l;
+ *Duplex_Er+=LE;
+ break;
+ }
+ }
+ if (traced) break;
+ }
+ }
+ if (!traced) {
+/* if (i>1) {E -= P->dangle5[type][SS1[i-1]]; *Duplex_El +=P->dangle5[type][SS1[i-1]];} */
+/* if (j<n2) {E -= P->dangle3[type][SS2[j+1]]; *Duplex_El +=P->dangle3[type][SS2[j+1]];} */
+ if (type>2) {E -= P->TerminalAU; *Duplex_Er +=P->TerminalAU;}
+ if (E != P->DuplexInit) {
+ nrerror("backtrack failed in fold duplex end");
+ } else break;
+ }
+ }
+
+
+ /* struc = (char *) space(i0-i+1+j-j0+1+2); */ /* declare final duplex structure */
+ struc = (char *) space(i-i0+1+n2); /* declare final duplex structure */
+ char * struc2;
+ struc2 = (char *) space(n2+1);
+ /* char * struct_const; */
+
+ for (k=MIN2(i0,1); k<=i; k++) if (!st1[k-1]) st1[k-1] = '.';
+ /* for (k=j0; k<=j; k++) if (!st2[k-1]) st2[k-1] = struc_loop[k-1];*/ /* '.'; normal */
+ /* char * struct_const; */
+ /* struct_const = (char *) space(sizeof(char)*(n2+1)); */
+ for (k=1; k<=n2; k++) {
+ if (!st2[k-1]) st2[k-1] = struc_loop[k-1];/* '.'; */
+ struc2[k-1] = st2[k-1];/* '.'; */
+ /* if (k>=j0 && k<=j){ */
+ /* struct_const[k-1]='x'; */
+ /* } */
+ /* else{ */
+ /* if(k<j0) {struct_const[k-1]='<';} */
+ /* if(k>j) {struct_const[k-1]='>';} */
+ /* } */
+ }
+ char duplexseq_1[j];
+ char duplexseq_2[n2-j0+2];
+ if(j0<n2){
+ strncpy(duplexseq_1, snoseq, j-1);
+ strcpy(duplexseq_2, snoseq + j0);
+ duplexseq_1[j-1]='\0';
+ duplexseq_2[n2-j0+1]='\0';
+ duplexT temp;
+ temp=duplexfold(duplexseq_1, duplexseq_2);
+ *Loop_D = MIN2(0,-410 + (int) 100 * temp.energy);
+ if(*Loop_D){
+ int l1,ibegin, iend, jbegin, jend;
+ l1=strchr(temp.structure, '&')-temp.structure;
+ ibegin=temp.i-l1;
+ iend =temp.i-1;
+ jbegin=temp.j;
+ jend =temp.j+strlen(temp.structure)-l1-2-1;
+ for(k=ibegin+1; k<=iend+1; k++){
+ struc2[k-1]=temp.structure[k-ibegin-1];
+ }
+ for(k=jbegin+j0; k<=jend+j0; k++){
+ struc2[k-1]=temp.structure[l1+k-j0-jbegin+1];
+ }
+ }
+ free(temp.structure);
+ }
+ strcpy(struc, st1+MAX2(i0,1)); strcat(struc, "&");
+ /* strcat(struc, st2); */
+ strncat(struc, struc2+5, strlen(struc2)-10);
+ free(struc2);
+ free(struc_loop);
+ free(st1); free(st2);
+
+ for (i=0; i<=n1; i++) {free(r[i]);free(c[i]);}
+ free(c);
+ free(r);
+ free(S1);free(S2);free(SS1);free(SS2);
+ /* free_arrays(); */
+ return struc;
+}
+
+PRIVATE int covscore(const int *types, int n_seq) {
+ /* calculate co-variance bonus for a pair depending on */
+ /* compensatory/consistent mutations and incompatible seqs */
+ /* should be 0 for conserved pairs, >0 for good pairs */
+#define NONE -10000 /* score for forbidden pairs */
+ int k,l,s,score, pscore;
+ int dm[7][7]={{0,0,0,0,0,0,0}, /* hamming distance between pairs */
+ {0,0,2,2,1,2,2} /* CG */,
+ {0,2,0,1,2,2,2} /* GC */,
+ {0,2,1,0,2,1,2} /* GU */,
+ {0,1,2,2,0,2,1} /* UG */,
+ {0,2,2,1,2,0,2} /* AU */,
+ {0,2,2,2,1,2,0} /* UA */};
+
+ int pfreq[8]={0,0,0,0,0,0,0,0};
+ for (s=0; s<n_seq; s++)
+ pfreq[types[s]]++;
+
+ if (pfreq[0]*2>n_seq) return NONE;
+ for (k=1,score=0; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
+ for (l=k+1; l<=6; l++)
+ /* scores for replacements between pairtypes */
+ /* consistent or compensatory mutations score 1 or 2 */
+ score += pfreq[k]*pfreq[l]*dm[k][l];
+
+ /* counter examples score -1, gap-gap scores -0.25 */
+ pscore = cv_fact *
+ ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
+ return pscore;
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE short * aliencode_seq(const char *sequence) {
+ unsigned int i,l;
+ short *Stemp;
+ l = strlen(sequence);
+ Stemp = (short *) space(sizeof(short)*(l+2));
+ Stemp[0] = (short) l;
+
+ /* make numerical encoding of sequence */
+ for (i=1; i<=l; i++)
+ Stemp[i]= (short) encode_char(toupper(sequence[i-1]));
+
+ /* for circular folding add first base at position n+1 */
+ /* Stemp[l+1] = Stemp[1]; */
+
+ return Stemp;
+}
+
+PRIVATE short * encode_seq(const char *sequence) {
+ unsigned int i,l;
+ short *S;
+ l = strlen(sequence);
+extern double nc_fact;
+ S = (short *) space(sizeof(short)*(l+2));
+ S[0] = (short) l;
+
+ /* make numerical encoding of sequence */
+ for (i=1; i<=l; i++)
+ S[i]= (short) encode_char(toupper(sequence[i-1]));
+ /* for circular folding add first base at position n+1 */
+ S[l+1] = S[1];
+
+ return S;
+}
+
+PRIVATE void encode_seqs(const char *s1, const char *s2) {
+ unsigned int i,l;
+
+ l = strlen(s1);
+ S1 = encode_seq(s1);
+ SS1= (short *) space(sizeof(short)*(l+1));
+ /* SS1 exists only for the special X K and I bases and energy_set!=0 */
+
+ for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
+ SS1[i] = alias[S1[i]]; /* for mismatches of nostandard bases */
+ }
+
+ l = strlen(s2);
+ S2 = encode_seq(s2);
+ SS2= (short *) space(sizeof(short)*(l+1));
+ /* SS2 exists only for the special X K and I bases and energy_set!=0 */
+
+ for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
+ SS2[i] = alias[S2[i]]; /* for mismatches of nostandard bases */
+ }
+}
+
+PRIVATE int compare(const void *sub1, const void *sub2) {
+ int d;
+ if (((snoopT *) sub1)->energy > ((snoopT *) sub2)->energy)
+ return 1;
+ if (((snoopT *) sub1)->energy < ((snoopT *) sub2)->energy)
+ return -1;
+ d = ((snoopT *) sub1)->i - ((snoopT *) sub2)->i;
+ if (d!=0) return d;
+ return ((snoopT *) sub1)->j - ((snoopT *) sub2)->j;
+}
+
+
+
git://source.jalview.org / jabaws.git / blobdiff