/*
minimum free energy
RNA secondary structure prediction
c Ivo Hofacker, Chrisoph Flamm
original implementation by
Walter Fontana
Vienna RNA package
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include "ViennaRNA/utils.h"
#include "ViennaRNA/structure_utils.h"
#include "ViennaRNA/energy_par.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/pair_mat.h"
#include "ViennaRNA/params.h"
#include "ViennaRNA/snofold.h"
#include "ViennaRNA/loop_energies.h"
#ifdef __GNUC__
#define INLINE inline
#else
#define INLINE
#endif
#define PAREN
#define PUBLIC
#define PRIVATE static
#define STACK_BULGE1 1 /* stacking energies for bulges of size 1 */
#define NEW_NINIO 1 /* new asymetry penalty */
/*@unused@*/
PRIVATE void get_arrays(unsigned int size);
/* PRIVATE int stack_energy(int i, const char *string); */
PRIVATE void make_ptypes(const short *S, const char *structure);
PRIVATE void encode_seq(const char *sequence);
PRIVATE void backtrack(const char *sequence, int s);
PRIVATE int fill_arrays(const char *sequence, const int max_asymm, const int threshloop,
const int min_s2, const int max_s2, const int half_stem, const int max_half_stem);
/*@unused@*/
/* alifold */
PRIVATE void alisnoinitialize_fold(const int length);
PRIVATE void make_pscores(const short *const* S, const char *const* AS,int n_seq, const char *structure);
PRIVATE int *pscore; /* precomputed array of pair types */
PRIVATE short **Sali;
PRIVATE int alifill_arrays(const char **string, const int max_asymm, const int threshloop,
const int min_s2, const int max_s2, const int half_stem,
const int max_half_stem);
PRIVATE void aliget_arrays(unsigned int size);
PRIVATE short * aliencode_seq(const char *sequence);
PRIVATE int alibacktrack(const char **strings, int s);
#define UNIT 100
#define MINPSCORE -2 * UNIT
/* end alifold */
#define MAXSECTORS 500 /* dimension for a backtrack array */
#define LOCALITY 0. /* locality parameter for base-pairs */
#define MIN2(A, B) ((A) < (B) ? (A) : (B))
#define MAX2(A, B) ((A) > (B) ? (A) : (B))
#define SAME_STRAND(I,J) (((I)>=cut_point)||((J)<cut_point))
PRIVATE vrna_param_t *P = NULL;
PRIVATE int *indx = NULL; /* index for moving in the triangle matrices c[] and fMl[]*/
PRIVATE int *c = NULL; /* energy array, given that i-j pair */
PRIVATE int *cc = NULL; /* linear array for calculating canonical structures */
PRIVATE int *cc1 = NULL; /* " " */
PRIVATE int *Fmi = NULL; /* holds row i of fML (avoids jumps in memory) */
PRIVATE int *DMLi = NULL; /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j]) */
PRIVATE int *DMLi1 = NULL; /* MIN(fML[i+1,k]+fML[k+1,j]) */
PRIVATE int *DMLi2 = NULL; /* MIN(fML[i+2,k]+fML[k+1,j]) */
PRIVATE char *ptype = NULL; /* precomputed array of pair types */
PRIVATE short *S = NULL, *S1 = NULL;
PRIVATE int init_length=-1;
PRIVATE int *mLoop = NULL; /*contains the minimum of c for a xy range*/
PRIVATE folden **foldlist = NULL;
PRIVATE folden **foldlist_XS = NULL;
PRIVATE int *BP = NULL; /* contains the structure constrainsts: BP[i]
-1: | = base must be paired
-2: < = base must be paired with j<i
-3: > = base must be paired with j>i
-4: x = base must not pair
positive int: base is paired with int */
static sect sector[MAXSECTORS]; /* stack of partial structures for backtracking */
PRIVATE char alpha[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
/* needed by cofold/eval */
/* PRIVATE int cut_in_loop(int i); */
/* PRIVATE int min_hairpin = TURN; */
/* some definitions to take circfold into account... */
/* PRIVATE int *fM2 = NULL;*/ /* fM2 = multiloop region with exactly two stems, extending to 3' end */
PUBLIC int Fc, FcH, FcI, FcM; /* parts of the exterior loop energies */
/*--------------------------------------------------------------------------*/
void snoinitialize_fold(const int length)
{
unsigned int n;
if (length<1) vrna_message_error("snoinitialize_fold: argument must be greater 0");
if (init_length>0) snofree_arrays(length);
get_arrays((unsigned) length);
init_length=length;
for (n = 1; n <= (unsigned) length; n++)
indx[n] = (n*(n-1)) >> 1; /* n(n-1)/2 */
snoupdate_fold_params();
}
PRIVATE void alisnoinitialize_fold(const int length)
{
unsigned int n;
if (length<1) vrna_message_error("snoinitialize_fold: argument must be greater 0");
if (init_length>0) snofree_arrays(length);
aliget_arrays((unsigned) length);
make_pair_matrix();
init_length=length;
for (n = 1; n <= (unsigned) length; n++)
indx[n] = (n*(n-1)) >> 1; /* n(n-1)/2 */
snoupdate_fold_params();
}
/*--------------------------------------------------------------------------*/
PRIVATE void get_arrays(unsigned int size)
{
indx = (int *) vrna_alloc(sizeof(int)*(size+1));
c = (int *) vrna_alloc(sizeof(int)*((size*(size+1))/2+2));
mLoop = (int *) vrna_alloc(sizeof(int)*((size*(size+1))/2+2));
ptype = (char *) vrna_alloc(sizeof(char)*((size*(size+1))/2+2));
cc = (int *) vrna_alloc(sizeof(int)*(size+2));
cc1 = (int *) vrna_alloc(sizeof(int)*(size+2));
Fmi = (int *) vrna_alloc(sizeof(int)*(size+1));
DMLi = (int *) vrna_alloc(sizeof(int)*(size+1));
DMLi1 = (int *) vrna_alloc(sizeof(int)*(size+1));
DMLi2 = (int *) vrna_alloc(sizeof(int)*(size+1));
if (base_pair) free(base_pair);
base_pair = (vrna_bp_stack_t *) vrna_alloc(sizeof(vrna_bp_stack_t)*(1+size/2));
/* extra array(s) for circfold() */
}
PRIVATE void aliget_arrays(unsigned int size)
{
indx = (int *) vrna_alloc(sizeof(int)*(size+1));
c = (int *) vrna_alloc(sizeof(int)*((size*(size+1))/2+2));
mLoop = (int *) vrna_alloc(sizeof(int)*((size*(size+1))/2+2));
pscore = (int *) vrna_alloc(sizeof(int)*((size*(size+1))/2+2));
ptype = (char *) vrna_alloc(sizeof(char)*((size*(size+1))/2+2));
cc = (int *) vrna_alloc(sizeof(int)*(size+2));
cc1 = (int *) vrna_alloc(sizeof(int)*(size+2));
Fmi = (int *) vrna_alloc(sizeof(int)*(size+1));
DMLi = (int *) vrna_alloc(sizeof(int)*(size+1));
DMLi1 = (int *) vrna_alloc(sizeof(int)*(size+1));
DMLi2 = (int *) vrna_alloc(sizeof(int)*(size+1));
if (base_pair) free(base_pair);
base_pair = (vrna_bp_stack_t *) vrna_alloc(sizeof(vrna_bp_stack_t)*(1+size/2));
/* extra array(s) for circfold() */
}
/*--------------------------------------------------------------------------*/
void snofree_arrays(const int length)
{
free(indx); free(c);free(cc); free(cc1);
free(ptype);free(mLoop);
int i;
for(i=length;i>-1;i--){
while(foldlist[i]!=NULL){
folden *n = foldlist[i];
foldlist[i] = foldlist[i]->next;
free(n);
}
free(foldlist[i]);
}
free(foldlist);
for(i=length;i>-1;i--){
while(foldlist_XS[i]!=NULL){
folden *n = foldlist_XS[i];
foldlist_XS[i] = foldlist_XS[i]->next;
free(n);
}
free(foldlist_XS[i]);
}
free(foldlist_XS);
free(base_pair); base_pair=NULL; free(Fmi);
free(DMLi); free(DMLi1);free(DMLi2);
free(BP);
init_length=0;
}
void alisnofree_arrays(const int length)
{
free(indx); free(c);free(cc); free(cc1);
free(ptype);free(mLoop);free(pscore);
int i;
for(i=length-1;i>-1;i--){
while(foldlist[i]!=NULL){
folden *n = foldlist[i];
foldlist[i] = foldlist[i]->next;
free(n);
}
free(foldlist[i]);
}
free(foldlist);
free(base_pair); base_pair=NULL; free(Fmi);
free(DMLi); free(DMLi1);free(DMLi2);
free(BP);
init_length=0;
}
/*--------------------------------------------------------------------------*/
void snoexport_fold_arrays(int **indx_p, int **mLoop_p, int **cLoop, folden ***fold_p, folden ***fold_p_XS) {
/* make the DP arrays available to routines such as subopt() */
*indx_p = indx; *mLoop_p = mLoop;
*cLoop = c; *fold_p = foldlist;*fold_p_XS=foldlist_XS;
}
/* void alisnoexport_fold_arrays(int **indx_p, int **mLoop_p, int **cLoop, folden ***fold_p, int **pscores) { */
/* /\* make the DP arrays available to routines such as subopt() *\/ */
/* *indx_p = indx; *mLoop_p = mLoop; */
/* *cLoop = c; *fold_p = foldlist; */
/* *pscores=pscore; */
/* } */
/*--------------------------------------------------------------------------*/
int snofold(const char *string, char *structure, const int max_assym, const int threshloop,
const int min_s2, const int max_s2, const int half_stem, const int max_half_stem) {
int length, energy, bonus, bonus_cnt, s;
/* Variable initialization */
bonus = 0;
bonus_cnt = 0;
s = 0;
length = (int) strlen(string);
S = encode_sequence(string, 0);
S1 = encode_sequence(string, 1);
/* structure = (char *) vrna_alloc((unsigned) length+1); */
if (length>init_length) snoinitialize_fold(length);
else if (fabs(P->temperature - temperature)>1e-6) snoupdate_fold_params();
/* encode_seq(string); */
BP = (int *)vrna_alloc(sizeof(int)*(length+2));
make_ptypes(S, structure);
energy=fill_arrays(string, max_assym, threshloop, min_s2, max_s2, half_stem, max_half_stem);
backtrack(string, s);
free(structure);
free(S); free(S1); /* free(BP); */
return energy;
}
PRIVATE void make_pscores(const short *const* S, const char *const* AS,
int n_seq, const char *structure) {
/* calculate co-variance bonus for each 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 n,i,j,k,l,s,score;
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 */};
n=Sali[0][0]; /* length of seqs */
for (i=1; i<n; i++) {
for (j=i+1; (j<i+TURN+1) && (j<=n); j++)
pscore[indx[j]+i] = NONE;
for (j=i+TURN+1; j<=n; j++) {
int pfreq[8]={0,0,0,0,0,0,0,0};
for (s=0; s<n_seq; s++) {
int type;
if (Sali[s][i]==0 && Sali[s][j]==0) type = 7; /* gap-gap */
else {
if ((AS[s][i] == '~')||(AS[s][j] == '~')) type = 7;
else type = pair[Sali[s][i]][Sali[s][j]];
}
pfreq[type]++;
}
if (pfreq[0]*2>n_seq) { pscore[indx[j]+i] = NONE; continue;}
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[indx[j]+i] = cv_fact *
((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
}
}
if (noLonelyPairs) /* remove unwanted pairs */
for (k=1; k<n-TURN-1; k++)
for (l=1; l<=2; l++) {
int type,ntype=0,otype=0;
i=k; j = i+TURN+l;
type = pscore[indx[j]+i];
while ((i>=1)&&(j<=n)) {
if ((i>1)&&(j<n)) ntype = pscore[indx[j+1]+i-1];
if ((otype<-4*UNIT)&&(ntype<-4*UNIT)) /* worse than 2 counterex */
pscore[indx[j]+i] = NONE; /* i.j can only form isolated pairs */
otype = type;
type = ntype;
i--; j++;
}
}
if (fold_constrained&&(structure!=NULL)) {
int psij, hx, hx2, *stack, *stack2;
stack = (int *) vrna_alloc(sizeof(int)*(n+1));
stack2 = (int *) vrna_alloc(sizeof(int)*(n+1));
for(hx=hx2=0, j=1; j<=n; j++) {
switch (structure[j-1]) {
case 'x': /* can't pair */
for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
break;
case '(':
stack[hx++]=j;
/* fallthrough */
case '[':
stack2[hx2++]=j;
/* fallthrough */
case '<': /* pairs upstream */
for (l=1; l<j-TURN; l++) pscore[indx[j]+l] = NONE;
break;
case ']':
if (hx2<=0) {
vrna_message_error("unbalanced brackets in constraints\n%s", structure);
}
i = stack2[--hx2];
pscore[indx[j]+i]=NONE;
break;
case ')':
if (hx<=0) {
vrna_message_error("unbalanced brackets in constraints\n%s", structure);
}
i = stack[--hx];
psij = pscore[indx[j]+i]; /* store for later */
for (k=j; k<=n; k++)
for (l=i; l<=j; l++)
pscore[indx[k]+l] = NONE;
for (l=i; l<=j; l++)
for (k=1; k<=i; k++)
pscore[indx[l]+k] = NONE;
for (k=i+1; k<j; k++)
pscore[indx[k]+i] = pscore[indx[j]+k] = NONE;
pscore[indx[j]+i] = (psij>0) ? psij : 0;
/* fallthrough */
case '>': /* pairs downstream */
for (l=j+TURN+1; l<=n; l++) pscore[indx[l]+j] = NONE;
break;
}
}
if (hx!=0) {
vrna_message_error("unbalanced brackets in constraint string\n%s", structure);
}
free(stack); free(stack2);
}
}
float alisnofold(const char **strings, const int max_assym, const int threshloop,
const int min_s2, const int max_s2, const int half_stem, const int max_half_stem) {
int s,n_seq, length, energy;
char * structure;
length = (int) strlen(strings[0]);
/* structure = (char *) vrna_alloc((unsigned) length+1); */
structure = NULL;
if (length>init_length) alisnoinitialize_fold(length);
if (fabs(P->temperature - temperature)>1e-6) snoupdate_fold_params();
for (s=0; strings[s]!=NULL; s++);
n_seq = s;
Sali = (short **) vrna_alloc(n_seq*sizeof(short *));
for (s=0; s<n_seq; s++) {
if (strlen(strings[s]) != length) vrna_message_error("uneqal seqence lengths");
Sali[s] = aliencode_seq(strings[s]);
}
make_pscores((const short **) Sali, (const char *const *) strings, n_seq, structure);
energy=alifill_arrays(strings, max_assym, threshloop, min_s2, max_s2, half_stem, max_half_stem);
alibacktrack((const char **)strings, 0);
for (s=0; s<n_seq; s++) free(Sali[s]);
free(Sali);
/* free(structure); */
/* free(S)*/; free(S1); /* free(BP); */
return (float) energy/100.;
}
PRIVATE int alifill_arrays(const char **strings, const int max_asymm, const int threshloop,
const int min_s2, const int max_s2, const int half_stem,
const int max_half_stem) {
int i, j, length, energy;
/* int decomp, new_fML; */
int *type, type_2;
int bonus,n_seq,s;
for (n_seq=0; strings[n_seq]!=NULL; n_seq++);
type = (int *) vrna_alloc(n_seq*sizeof(int));
length = strlen(strings[0]);
bonus=0;
/* max_separation = (int) ((1.-LOCALITY)*(double)(length-2));*/ /* not in use */
/* for (i=(j>TURN?(j-TURN):1); i<j; i++) { */
/* } */
for (i = (length)-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
for (j = i+TURN+1; j <= length; j++) {
int p, q, ij,psc;
ij = indx[j]+i;
for (s=0; s<n_seq; s++) {
type[s] = pair[Sali[s][i]][Sali[s][j]];
if (type[s]==0) type[s]=7;
}
psc = pscore[indx[j]+i];
if (psc>=MINPSCORE) { /* we have a pair */
int new_c=0, stackEnergy=INF; /* seems that new_c immer den minimum von cij enthaelt */
/* hairpin ----------------------------------------------*/
for (new_c=s=0; s<n_seq; s++)
new_c += E_Hairpin(j-i-1,type[s],Sali[s][i+1],Sali[s][j-1],strings[s]+i-1,P);
/*--------------------------------------------------------
check for elementary structures involving more than one
closing pair (interior loop).
--------------------------------------------------------*/
for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++) {
int minq = j-i+p-MAXLOOP-2;
if (minq<p+1+TURN) minq = p+1+TURN;
for (q = minq; q < j; q++) {
if (pscore[indx[q]+p]<MINPSCORE) continue;
if(abs((p-i) - (j-q)) > max_asymm) continue;
for (energy = s=0; s<n_seq; s++) {
type_2 = pair[Sali[s][q]][Sali[s][p]]; /* q,p not p,q! */
if (type_2 == 0) type_2 = 7;
energy += E_IntLoop(p-i-1, j-q-1, type[s], type_2,
Sali[s][i+1], Sali[s][j-1],
Sali[s][p-1], Sali[s][q+1],P);
}
new_c = MIN2(energy+c[indx[q]+p], new_c);
if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
} /* end q-loop */
} /* end p-loop */
/* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */
new_c = MIN2(new_c, cc1[j-1]+stackEnergy);
cc[j] = new_c - psc; /* add covariance bonnus/penalty */
c[ij]=cc[j];
} /* end >> if (pair) << */
else c[ij] = INF;
/* done with c[i,j], now compute fML[i,j] */
/* free ends ? -----------------------------------------*/
}
{
int *FF; /* rotate the auxilliary arrays */
FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
FF = cc1; cc1=cc; cc=FF;
for (j=1; j<=length; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
}
}
foldlist = (folden**) vrna_alloc((length)*sizeof(folden*));
for(i=0; i< length; i++){
foldlist[i]=(folden*) vrna_alloc(sizeof(folden));
foldlist[i]->next=NULL;
foldlist[i]->k=INF+1;
foldlist[i]->energy=INF;
}
folden* head; /* we save the stem loop information in a list like structure */
for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
int max_k, min_k;
max_k = MIN2(length-min_s2,i+max_half_stem+1);
min_k = MAX2(i+half_stem+1, length-max_s2);
for (j = i+TURN+1; j <= length; j++) {
int ij,a,b;
ij = indx[j]+i;
for(a=0; a< MISMATCH ;a++){
for(b=0; b< MISMATCH ; b++){
mLoop[ij]=MIN2(mLoop[ij], c[indx[j-a]+i+b]);
}
}
if(mLoop[ij]>=n_seq*threshloop){
mLoop[ij]=INF;
}
else{
if(j>=min_k-1 && j < max_k){ /* comment if out to recover the known behaviour */
head = (folden*) vrna_alloc(sizeof(folden));
head->k=j;
head->energy=mLoop[ij];
head->next=foldlist[i];
foldlist[i] = head;
}
}
}
}
free(type);
return mLoop[indx[length]+1];/* mLoop; */
}
PRIVATE int alibacktrack(const char **strings, int s) {
/*------------------------------------------------------------------
trace back through the "c", "f5" and "fML" arrays to get the
base pairing list. No search for equivalent structures is done.
This is fast, since only few structure elements are recalculated.
------------------------------------------------------------------*/
/* normally s=0.
If s>0 then s items have been already pushed onto the sector stack */
int i, j, length, energy;/* , new; */
int type_2;
int bonus,n_seq,*type; int b=0,cov_en = 0;
length = strlen(strings[0]);
for (n_seq=0; strings[n_seq]!=NULL; n_seq++);
type = (int *) vrna_alloc(n_seq*sizeof(int));
if (s==0) {
sector[++s].i = 1;
sector[s].j = length;
sector[s].ml = 2 ;
}
while (s>0) {
int ml, ss, cij, traced, i1, j1, p, q;
int canonical = 1; /* (i,j) closes a canonical structure */
i = sector[s].i;
j = sector[s].j;
ml = sector[s--].ml; /* ml is a flag indicating if backtracking is to
occur in the fML- (1) or in the f-array (0) */
if (ml==2) {
base_pair[++b].i = i;
base_pair[b].j = j;
goto repeat1;
}
if (j < i+TURN+1) continue; /* no more pairs in this interval */
repeat1:
/*----- begin of "repeat:" -----*/
if (canonical) cij = c[indx[j]+i];
for (ss=0; ss<n_seq; ss++) {
type[ss] = pair[Sali[ss][i]][Sali[ss][j]];
if (type[ss]==0) type[ss] = 7;
}
bonus = 0;
if (noLonelyPairs)
if (cij == c[indx[j]+i]) {
/* (i.j) closes canonical structures, thus
(i+1.j-1) must be a pair */
for (ss=0; ss<n_seq; ss++) {
type_2 = pair[Sali[ss][j-1]][Sali[ss][i+1]]; /* j,i not i,j */
if (type_2==0) type_2 = 7;
cij -= P->stack[type[ss]][type_2];
}
cij += pscore[indx[j]+i];
base_pair[++b].i = i+1;
base_pair[b].j = j-1;
cov_en += pscore[indx[j-1]+i+1];
i++; j--;
canonical=0;
goto repeat1;
}
canonical = 1;
cij += pscore[indx[j]+i];
{int cc=0;
for (ss=0; ss<n_seq; ss++)
cc += E_Hairpin(j-i-1, type[ss], Sali[ss][i+1], Sali[ss][j-1], strings[ss]+i-1,P);
if (cij == cc) /* found hairpin */
continue;
}
for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
int minq;
minq = j-i+p-MAXLOOP-2;
if (minq<p+1+TURN) minq = p+1+TURN;
for (q = j-1; q >= minq; q--) {
for (ss=energy=0; ss<n_seq; ss++) {
type_2 = pair[Sali[ss][q]][Sali[ss][p]]; /* q,p not p,q */
if (type_2==0) type_2 = 7;
energy += E_IntLoop(p-i-1, j-q-1, type[ss], type_2,
Sali[ss][i+1], Sali[ss][j-1],
Sali[ss][p-1], Sali[ss][q+1],P);
}
traced = (cij == energy+c[indx[q]+p]);
if (traced) {
base_pair[++b].i = p;
base_pair[b].j = q;
cov_en += pscore[indx[q]+p];
i = p, j = q;
goto repeat1;
}
}
}
/* end of repeat: --------------------------------------------------*/
/* (i.j) must close a multi-loop */
/* tt = rtype[type]; */
/* mm = bonus+P->MLclosing+P->MLintern[tt]; */
/* d5 = P->dangle5[tt][S1[j-1]]; */
/* d3 = P->dangle3[tt][S1[i+1]]; */
i1 = i+1; j1 = j-1;
sector[s+1].ml = sector[s+2].ml = 1;
/* if (k<=j-3-TURN) { /\\* found the decomposition *\\/ *\/ */
/* sector[++s].i = i1; */
/* sector[s].j = k; */
/* sector[++s].i = k+1; */
/* sector[s].j = j1; */
/* } /\* else { *\/ */
/* vrna_message_error("backtracking failed in repeat"); */
/* } */
}
base_pair[0].i = b; /* save the total number of base pairs */
free(type);
return cov_en;
}
PRIVATE int fill_arrays(const char *string, const int max_asymm, const int threshloop,
const int min_s2, const int max_s2, const int half_stem, const int max_half_stem) {
int i, j, length, energy;
/* int decomp;*/ /*, new_fML; */
int no_close, type, type_2;
int bonus;
int min_c;
min_c=INF;
length = (int) strlen(string);
bonus=0;
/* max_separation = (int) ((1.-LOCALITY)*(double)(length-2)); */ /* not in use */
for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
/* printf("i=%d\t",i); */
for (j = i+TURN+1; j <= length; j++) {
/* printf("j=%d,",j); */
int p, q, ij;
ij = indx[j]+i;
type = ptype[ij];
bonus = 0;
energy = INF;
if ((BP[i]==j)||(BP[i]==-1)||(BP[i]==-2)) bonus -= BONUS;
if ((BP[j]==-1)||(BP[j]==-3)) bonus -= BONUS;
if ((BP[i]==-4)||(BP[j]==-4)) type=0;
no_close = (((type==3)||(type==4))&&no_closingGU);
/* if (j-i-1 > max_separation) type = 0; */ /* forces locality degree */
if (type) { /* we have a pair */
int new_c=0, stackEnergy=INF; /* seems that new_c immer den minimum von cij enthaelt */
/* hairpin ----------------------------------------------*/
if (no_close) new_c = FORBIDDEN;
else
new_c = E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], string+i-1,P); /* computes hair pin structure for subsequence i...j */
/*--------------------------------------------------------
check for elementary structures involving more than one
closing pair (interior loop).
--------------------------------------------------------*/
for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++) {
int minq = j-i+p-MAXLOOP-2;
if (minq<p+1+TURN) minq = p+1+TURN;
for (q = minq; q < j; q++) {
if(abs((p-i) - (j-q)) > max_asymm) continue;
type_2 = ptype[indx[q]+p];
if (type_2==0) continue;
type_2 = rtype[type_2];
if (no_closingGU)
if (no_close||(type_2==3)||(type_2==4))
if ((p>i+1)||(q<j-1)) continue; /* continue unless stack */
energy = E_IntLoop(p-i-1, j-q-1, type, type_2,
S1[i+1], S1[j-1], S1[p-1], S1[q+1],P);
new_c = MIN2(energy+c[indx[q]+p], new_c);
if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
} /* end q-loop */
} /* end p-loop */
/* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */
new_c = MIN2(new_c, cc1[j-1]+stackEnergy);
cc[j] = new_c;
c[ij] = new_c;
/* min_c=MIN2(min_c, c[ij]); */
} /* end >> if (pair) << */
else c[ij] = INF;
/* done with c[i,j], now compute fML[i,j] */
/* free ends ? -----------------------------------------*/
}
{
int *FF; /* rotate the auxilliary arrays */
FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
FF = cc1; cc1=cc; cc=FF;
for (j=1; j<=length; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
}
}
foldlist = (folden**) vrna_alloc((length+1)*sizeof(folden*));
foldlist_XS = (folden**) vrna_alloc((length+1)*sizeof(folden*));
/* linked list initialization*/
for(i=0; i<=length; i++){
foldlist[i]=(folden*) vrna_alloc(sizeof(folden));
foldlist[i]->next=NULL;
foldlist[i]->k=INF+1;
foldlist[i]->energy=INF;
foldlist_XS[i]=(folden*) vrna_alloc(sizeof(folden));
foldlist_XS[i]->next=NULL;
foldlist_XS[i]->k=INF+1;
foldlist_XS[i]->energy=INF;
}
folden* head; /* we save the stem loop information in a list like structure */
folden* head_XS;
for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
int max_k, min_k;
max_k = MIN2(length-min_s2,i+max_half_stem+1);
min_k = MAX2(i+half_stem+1, length-max_s2);
for (j = i+TURN+1; j <= length; j++) {
int ij,a,b;
ij = indx[j]+i;
for(a=0; a< MISMATCH ;a++){
for(b=0; b< MISMATCH ; b++){
mLoop[ij]=MIN2(mLoop[ij], c[indx[j-a]+i+b]);
/* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-2]+i]); */
/* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j]+i+1]); */
/* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-1]+i+1]); */
/* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-2]+i+1]); */
/* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j]+i+2]); */
/* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-1]+i+2]); */
/* #mLoop[ij]=MIN2(mLoop[ij], c[indx[j-2]+i+2]); */
}
}
min_c = MIN2(mLoop[ij] ,min_c);
if(mLoop[ij]>=threshloop){
mLoop[ij]=INF;
}
else{
if(j>=min_k-1 && j <= max_k){ /* comment if out to recover the known behaviour */
head = (folden*) vrna_alloc(sizeof(folden));
head->k=j;
head->energy=mLoop[ij];
head->next=foldlist[i];
foldlist[i] = head;
head_XS = (folden*) vrna_alloc(sizeof(folden));
head_XS->k=i;
head_XS->energy=mLoop[ij];
head_XS->next=foldlist_XS[j];
foldlist_XS[j] = head_XS;
}
}
}
}
/* int count=0; */
/* for(i=0; i< length; i++){ */
/* folden *temp; */
/* temp = foldlist[i]; */
/* while(temp->next){ */
/* count++; */
/* printf("count %d: i%d j%d energy %d \n", count, i, temp->k, temp->energy); */
/* temp=temp->next; */
/* } */
/* } */
/* printf("Count %d \n", count); */
/* count=0; */
/* for(i=length-1; i>=0; i--){ */
/* folden *temp; */
/* temp = foldlist_XS[i]; */
/* while(temp->next){ */
/* count++; */
/* printf("count %d: i%d j%d energy %d \n", count, temp->k,i, temp->energy); */
/* temp=temp->next; */
/* } */
/* } */
/* printf("Count %d \n", count); */
/* return mLoop[indx[length]+1]; */ /* mLoop; */
return min_c;
/* printf("\nmin_array = %d\n", min_c); */
/* return f5[length]; */
}
PRIVATE void backtrack(const char *string, int s) {
/*------------------------------------------------------------------
trace back through the "c", "f5" and "fML" arrays to get the
base pairing list. No search for equivalent structures is done.
This is fast, since only few structure elements are recalculated.
------------------------------------------------------------------*/
/* normally s=0.
If s>0 then s items have been already pushed onto the sector stack */
int i, j, /*k,*/ length, energy, new;
int no_close, type, type_2;/* , tt; */
int bonus;
int b=0;
length = strlen(string);
if (s==0) {
sector[++s].i = 1;
sector[s].j = length;
sector[s].ml = 2 ;
}
while (s>0) {
int ml, cij, traced, i1, j1, /*d3, d5, mm,*/ p, q;
int canonical = 1; /* (i,j) closes a canonical structure */
i = sector[s].i;
j = sector[s].j;
ml = sector[s--].ml; /* ml is a flag indicating if backtracking is to
occur in the fML- (1) or in the f-array (0) */
if (ml==2) {
base_pair[++b].i = i;
base_pair[b].j = j;
goto repeat1;
}
if (j < i+TURN+1) continue; /* no more pairs in this interval */
repeat1:
/*----- begin of "repeat:" -----*/
if (canonical) cij = c[indx[j]+i];
type = ptype[indx[j]+i];
bonus = 0;
if (fold_constrained) {
if ((BP[i]==j)||(BP[i]==-1)||(BP[i]==-2)) bonus -= BONUS;
if ((BP[j]==-1)||(BP[j]==-3)) bonus -= BONUS;
}
if (noLonelyPairs)
if (cij == c[indx[j]+i]) {
/* (i.j) closes canonical structures, thus
(i+1.j-1) must be a pair */
type_2 = ptype[indx[j-1]+i+1]; type_2 = rtype[type_2];
cij -= P->stack[type][type_2] + bonus;
base_pair[++b].i = i+1;
base_pair[b].j = j-1;
i++; j--;
canonical=0;
goto repeat1;
}
canonical = 1;
no_close = (((type==3)||(type==4))&&no_closingGU&&(bonus==0));
if (no_close) {
if (cij == FORBIDDEN) continue;
} else
if (cij == E_Hairpin(j-i-1, type, S1[i+1], S1[j-1],string+i-1,P)+bonus)
continue;
for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
int minq;
minq = j-i+p-MAXLOOP-2;
if (minq<p+1+TURN) minq = p+1+TURN;
for (q = j-1; q >= minq; q--) {
type_2 = ptype[indx[q]+p];
if (type_2==0) continue;
type_2 = rtype[type_2];
if (no_closingGU)
if (no_close||(type_2==3)||(type_2==4))
if ((p>i+1)||(q<j-1)) continue; /* continue unless stack */
energy = E_IntLoop(p-i-1, j-q-1, type, type_2,
S1[i+1], S1[j-1], S1[p-1], S1[q+1],P);
new = energy+c[indx[q]+p]+bonus;
traced = (cij == new);
if (traced) {
base_pair[++b].i = p;
base_pair[b].j = q;
i = p, j = q;
goto repeat1;
}
}
}
/* end of repeat: --------------------------------------------------*/
/* (i.j) must close a multi-loop */
/* tt = rtype[type]; */
/* mm = bonus+P->MLclosing+P->MLintern[tt]; */
/* d5 = P->dangle5[tt][S1[j-1]]; */
/* d3 = P->dangle3[tt][S1[i+1]]; */
i1 = i+1; j1 = j-1;
sector[s+1].ml = sector[s+2].ml = 1;
/* if (k<=j-3-TURN) { */ /* found the decomposition */
/* sector[++s].i = i1; */
/* sector[s].j = k; */
/* sector[++s].i = k+1; */
/* sector[s].j = j1; */
/* } else { */
/* vrna_message_error("backtracking failed in repeat"); */
/* } */
/* */
}
base_pair[0].i = b; /* save the total number of base pairs */
}
char *snobacktrack_fold_from_pair(const char *sequence, int i, int j) {
char *structure;
sector[1].i = i;
sector[1].j = j;
sector[1].ml = 2;
base_pair[0].i=0;
encode_seq(sequence);
backtrack(sequence, 1);
structure = vrna_db_from_bp_stack(base_pair, strlen(sequence));
free(S);free(S1);
return structure;
}
char *alisnobacktrack_fold_from_pair(const char **strings, int i, int j, int *cov) {
char *structure;
int n_seq, s, length;
length = (int) strlen(strings[0]);
for (s=0; strings[s]!=NULL; s++);
n_seq = s;
sector[1].i = i;
sector[1].j = j;
sector[1].ml = 2;
base_pair[0].i=0;
/* encode_seq(sequence); */
Sali = (short **) vrna_alloc(n_seq*sizeof(short *));
for (s=0; s<n_seq; s++) {
if (strlen(strings[s]) != length) vrna_message_error("uneqal seqence lengths");
Sali[s] = aliencode_seq(strings[s]);
}
*cov=alibacktrack(strings, 1);
structure = vrna_db_from_bp_stack(base_pair, length);
free(S);free(S1);
for (s=0; s<n_seq; s++) {
free(Sali[s]);
}
free(Sali);
return structure;
}
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
PRIVATE void encode_seq(const char *sequence) {
unsigned int i,l;
l = strlen(sequence);
S = (short *) vrna_alloc(sizeof(short)*(l+2));
S1= (short *) vrna_alloc(sizeof(short)*(l+2));
/* S1 exists only for the special X K and I bases and energy_set!=0 */
S[0] = (short) l;
for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
S[i]= (short) encode_char(toupper(sequence[i-1]));
S1[i] = alias[S[i]]; /* for mismatches of nostandard bases */
}
/* for circular folding add first base at position n+1 and last base at
position 0 in S1 */
S[l+1] = S[1]; S1[l+1]=S1[1]; S1[0] = S1[l];
}
PRIVATE short * aliencode_seq(const char *sequence) {
unsigned int i,l;
short *Stemp;
l = strlen(sequence);
Stemp = (short *) vrna_alloc(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;
}
/*---------------------------------------------------------------------------*/
PUBLIC void snoupdate_fold_params(void)
{
vrna_md_t md;
if(P)
free(P);
set_model_details(&md);
P = vrna_params(&md);
make_pair_matrix();
if (init_length < 0) init_length=0;
}
/*---------------------------------------------------------------------------*/
PRIVATE void make_ptypes(const short *S, const char *structure) {
int n,i,j,k,l;
n=S[0];
for (k=1; k<n-TURN; k++)
for (l=1; l<=2; l++) {
int type,ntype=0,otype=0;
i=k; j = i+TURN+l; if (j>n) continue;
type = pair[S[i]][S[j]];
while ((i>=1)&&(j<=n)) {
if ((i>1)&&(j<n)) ntype = pair[S[i-1]][S[j+1]];
if (noLonelyPairs && (!otype) && (!ntype))
type = 0; /* i.j can only form isolated pairs */
ptype[indx[j]+i] = (char) type;
otype = type;
type = ntype;
i--; j++;
}
}
if (fold_constrained&&(structure!=NULL)) {
constrain_ptypes(structure, (unsigned int)n, ptype, BP, TURN, 0);
}
}