/*
minimum free energy consensus
RNA secondary structure prediction
with maximum distance base pairs
c Ivo Hofacker, Stephan Bernhart
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/energy_par.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/pair_mat.h"
#include "ViennaRNA/params.h"
#include "ViennaRNA/ribo.h"
#include "ViennaRNA/alifold.h"
#include "ViennaRNA/fold.h"
#include "ViennaRNA/loop_energies.h"
#ifdef _OPENMP
#include <omp.h>
#endif
#define PAREN
#define STACK_BULGE1 1 /* stacking energies for bulges of size 1 */
#define NEW_NINIO 1 /* new asymetry penalty */
#define MAXSECTORS 500 /* dimension for a backtrack array */
#define LOCALITY 0. /* locality parameter for base-pairs */
#define UNIT 100
#define MINPSCORE -2 * UNIT
#define NONE -10000 /* score for forbidden pairs */
/*
#################################
# GLOBAL VARIABLES #
#################################
*/
/*
#################################
# PRIVATE VARIABLES #
#################################
*/
PRIVATE vrna_param_t *P = NULL;
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 *f3 = NULL; /* energy of 5' end */
PRIVATE int **fML = NULL; /* multi-loop auxiliary energy array */
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 int **pscore = NULL; /* precomputed array of pair types */
PRIVATE unsigned int length;
PRIVATE short **S = NULL;
PRIVATE short **S5 = NULL; /*S5[s][i] holds next base 5' of i in sequence s*/
PRIVATE short **S3 = NULL; /*Sl[s][i] holds next base 3' of i in sequence s*/
PRIVATE char **Ss = NULL;
PRIVATE unsigned short **a2s = NULL;
PRIVATE float **dm = NULL;
PRIVATE int olddm[7][7]= {{0,0,0,0,0,0,0}, /* hamming distance between pairs PRIVATE needed??*/
{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 */};
PRIVATE int energyout;
PRIVATE int energyprev;
#ifdef _OPENMP
/* NOTE: all variables are assumed to be uninitialized if they are declared as threadprivate
*/
#pragma omp threadprivate(P, c, cc, cc1, f3, fML, Fmi, DMLi, DMLi1, DMLi2, pscore, length, S, dm, S5, S3, Ss, a2s, energyout, energyprev)
#endif
/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
PRIVATE void initialize_aliLfold(int length, int maxdist);
PRIVATE void free_aliL_arrays(int maxdist);
PRIVATE void get_arrays(unsigned int size, int maxdist);
PRIVATE short *encode_seq(const char *sequence, short *s5, short *s3, char *ss, unsigned short *as);
PRIVATE void make_pscores(const char ** AS, const char *structure,int maxdist, int start);
PRIVATE int fill_arrays(const char **strings, int maxdist, char *structure);
PRIVATE char *backtrack(const char **strings, int start, int maxdist);
/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
PRIVATE void initialize_aliLfold(int length, int maxdist){
if (length<1) vrna_message_error("initialize_fold: argument must be greater 0");
get_arrays((unsigned) length, maxdist);
make_pair_matrix();
if(P) free(P);
vrna_md_t md;
set_model_details(&md);
P = vrna_params(&md);
}
/*--------------------------------------------------------------------------*/
PRIVATE void get_arrays(unsigned int size, int maxdist)
{
int i;
c = (int **)vrna_alloc(sizeof(int *)*(size+1));
fML = (int **)vrna_alloc(sizeof(int *)*(size+1));
pscore = (int **)vrna_alloc(sizeof(int *)*(size+1));
f3 = (int *) vrna_alloc(sizeof(int)*(size+2)); /* has to be one longer */
cc = (int *) vrna_alloc(sizeof(int)*(maxdist+5));
cc1 = (int *) vrna_alloc(sizeof(int)*(maxdist+5));
Fmi = (int *) vrna_alloc(sizeof(int)*(maxdist+5));
DMLi = (int *) vrna_alloc(sizeof(int)*(maxdist+5));
DMLi1 = (int *) vrna_alloc(sizeof(int)*(maxdist+5));
DMLi2 = (int *) vrna_alloc(sizeof(int)*(maxdist+5));
for (i=size; i>(int)size-maxdist-5 && i>=0; i--) {
c[i] = (int *) vrna_alloc(sizeof(int) *(maxdist+5));
fML[i] = (int *) vrna_alloc(sizeof(int) *(maxdist+5));
pscore[i] = (int *) vrna_alloc(sizeof(int )*(maxdist+5));
}
}
/*--------------------------------------------------------------------------*/
PRIVATE void free_aliL_arrays(int maxdist) {
int i;
for(i=0; i<maxdist+5 && i<=length; i++){
free(c[i]);
free(fML[i]);
free(pscore[i]);
}
free(c);
free(fML);
free(f3);
free(cc);
free(cc1);
free(pscore);
free(Fmi);
free(DMLi);
free(DMLi1);
free(DMLi2);
}
/*--------------------------------------------------------------------------*/
PUBLIC float aliLfold(const char **strings, char *structure, int maxdist) {
int length, energy, s, n_seq, i, j;
length = (int) strlen(strings[0]);
if (maxdist>length) maxdist = length;
initialize_aliLfold(length, maxdist);
for (s=0; strings[s]!=NULL; s++);
n_seq = s;
S = (short **) vrna_alloc(n_seq*sizeof(short *));
S5 = (short **) vrna_alloc(n_seq*sizeof(short *));
S3 = (short **) vrna_alloc(n_seq*sizeof(short *));
a2s = (unsigned short **) vrna_alloc(n_seq*sizeof(unsigned short *));
Ss = (char **) vrna_alloc(n_seq*sizeof(char *));
for (s=0; s<n_seq; s++) {
if (strlen(strings[s]) != length) vrna_message_error("uneqal seqence lengths");
S5[s] = (short *) vrna_alloc((length+2)*sizeof(short));
S3[s] = (short *) vrna_alloc((length+2)*sizeof(short));
a2s[s] = (unsigned short *) vrna_alloc((length+2)*sizeof(unsigned short));
Ss[s] = (char *) vrna_alloc((length+2)*sizeof(char));
S[s] = encode_seq(strings[s], S5[s],S3[s],Ss[s],a2s[s]);
}
if (ribo) {
if (RibosumFile !=NULL) dm=readribosum(RibosumFile);
else dm=get_ribosum(strings, n_seq, S[0][0]);
}
else { /*use usual matrix*/
dm=(float **)vrna_alloc(7*sizeof(float*));
for (i=0; i<7;i++) {
dm[i]=(float *)vrna_alloc(7*sizeof(float));
for (j=0; j<7; j++)
dm[i][j] = (float) olddm[i][j];
}
}
for (i=length; i>=(int)length-(int)maxdist-4 && i>0; i--)
make_pscores((const char **) strings,structure,maxdist,i);
energy = fill_arrays(strings, maxdist, structure);
free_aliL_arrays(maxdist);
return (float) energy/100.;
}
PRIVATE int fill_arrays(const char **strings, int maxdist, char *structure) {
/* fill "c", "fML" and "f3" arrays and return optimal energy */
int i, j, k, length, energy;
int decomp, new_fML,MLenergy ;
int *type, type_2, tt, s, n_seq, lastf, lastf2, thisj, lastj;
lastf = lastf2 = INF;
/* int bonus=0;*/
length = (int) strlen(strings[0]);
for (s=0; strings[s]!=NULL; s++);
n_seq = s;
type = (int *) vrna_alloc(n_seq*sizeof(int));
for (j=0; j<maxdist+5; j++)
Fmi[j]=DMLi[j]=DMLi1[j]=DMLi2[j]=INF;
for (j=length; j>length-maxdist-3; j--) {
for (i=(length-maxdist-2>0)?length-maxdist-2:1 ; i<j; i++)
c[i][j-i] = fML[i][j-i] = INF;
}
for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
for (j = i+1; j<=length && j<=i+TURN; j++) {
c[i][j-i]=fML[i][j-i]=INF;
}
for (j = i+TURN+1; j <= length && j <= i+maxdist; j++) {
int p, q, psc;
/* bonus = 0;*/
for (s=0; s<n_seq; s++) {
type[s] = pair[S[s][i]][S[s][j]];
if (type[s]==0) type[s]=7;
}
psc = pscore[i][j-i];
if (psc>=cv_fact*MINPSCORE) { /* we have a pair 2 consider */
int new_c=0, stackEnergy=INF;
/* hairpin ----------------------------------------------*/
for (new_c=s=0; s<n_seq; s++){
if((a2s[s][j-1] - a2s[s][i]) < 3) new_c += 600;
else new_c += E_Hairpin(a2s[s][j-1]-a2s[s][i],type[s],S3[s][i],S5[s][j],Ss[s]+(a2s[s][i-1]), P);
}
/*--------------------------------------------------------
check for elementary structures involving more than one
closing pair.
--------------------------------------------------------*/
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[p][q-p]<MINPSCORE) continue;
for (energy = s=0; s<n_seq; s++) {
type_2 = pair[S[s][q]][S[s][p]]; /* q,p not p,q! */
if (type_2 == 0) type_2 = 7;
energy += E_IntLoop(a2s[s][p-1]-a2s[s][i],
a2s[s][j-1]-a2s[s][q],
type[s],
type_2,
S3[s][i],
S5[s][j],
S5[s][p],
S3[s][q],
P);
}
new_c = MIN2(energy+c[p][q-p], new_c);
if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
} /* end q-loop */
} /* end p-loop */
/* multi-loop decomposition ------------------------*/
decomp = DMLi1[j-1-(i+1)];
if (dangles) {
for (s=0; s<n_seq; s++) {
tt = rtype[type[s]];
decomp += E_MLstem(tt, S5[s][j], S3[s][i], P);
}
}
else{
for(s=0; s<n_seq; s++){
tt = rtype[type[s]];
decomp += E_MLstem(tt, -1, -1, P);
}
}
MLenergy = decomp + n_seq*P->MLclosing;
new_c = MIN2(new_c, MLenergy);
new_c = MIN2(new_c, cc1[j-1-(i+1)]+stackEnergy);
cc[j-i] = new_c - psc; /* add covariance bonnus/penalty */
if (noLonelyPairs)
c[i][j-i] = cc1[j-1-(i+1)]+stackEnergy-psc;
else
c[i][j-i] = cc[j-i];
} /* end >> if (pair) << */
else c[i][j-i] = INF;
/* done with c[i,j], now compute fML[i,j] */
/* free ends ? -----------------------------------------*/
new_fML = fML[i+1][j-i-1]+n_seq*P->MLbase;
new_fML = MIN2(fML[i][j-1-i]+n_seq*P->MLbase, new_fML);
energy = c[i][j-i]/*+P->MLintern[type]*/;
if(dangles){
for (s=0; s<n_seq; s++) {
energy += E_MLstem(type[s], (i > 1) ? S5[s][i] : -1, (j < length) ? S3[s][j] : -1, P);
}
}
else{
for (s=0; s<n_seq; s++) {
energy += E_MLstem(type[s], -1, -1, P);
}
}
new_fML = MIN2(energy, new_fML);
/* modular decomposition -------------------------------*/
for (decomp = INF, k = i+1+TURN; k <= j-2-TURN; k++)
decomp = MIN2(decomp, Fmi[k-i]+fML[k+1][j-k-1]);
DMLi[j-i] = decomp; /* store for use in ML decompositon */
new_fML = MIN2(new_fML,decomp);
fML[i][j-i] = Fmi[j-i] = new_fML; /* substring energy */
} /* for (j...) */
/* calculate energies of 5' and 3' fragments */
{
static int do_backtrack = 0, prev_i=0;
static char * prev=NULL;
char *ss;
int thisf=0;
f3[i] = f3[i+1];
for (j=i+TURN+1; j<length && j<=i+maxdist; j++) {
if(c[i][j-i]<INF) {
/* if (c[j+1]<INF) {*/
energy = c[i][j-i];
if(dangles){
for(s = 0; s < n_seq; s++){
tt = pair[S[s][i]][S[s][j]];
if(tt==0) tt=7;
energy += E_ExtLoop(tt, (i>1) ? S5[s][i] : -1, S3[s][j], P);
}
}
else{
for(s = 0; s < n_seq; s++){
tt = pair[S[s][i]][S[s][j]];
if(tt==0) tt=7;
energy += E_ExtLoop(tt, -1, -1, P);
}
}
if (energy/(j-i+1) < thisf){
thisf = energy/(j-i+1);
thisj = j;
}
energy += f3[j+1];
if(f3[i] > energy){
f3[i] = energy;
}
}
}
if(length <= i+maxdist){
j = length;
if(c[i][j-i]<INF) {
energy = c[i][j-i];
if(dangles){
for (s=0; s<n_seq; s++) {
tt = pair[S[s][i]][S[s][j]];
if(tt==0) tt=7;
energy += E_ExtLoop(tt, (i>1) ? S5[s][i] : -1, -1, P);
}
}
else{
for (s=0; s<n_seq; s++) {
tt = pair[S[s][i]][S[s][j]];
if(tt==0) tt=7;
energy += E_ExtLoop(tt, -1, -1, P);
}
}
/* thisf=MIN2(energy/(j-i+1),thisf); ???*/
if (energy/(j-i+1) < thisf){
thisf = energy/(j-i+1);
thisj = j;
}
f3[i] = MIN2(f3[i], energy);
}
}
/* backtrack partial structure */
/* if (i+maxdist<length) {*/
if (i<length-1){
if (f3[i] != f3[i+1]) {
do_backtrack = 1;
backtrack_type = 'F';
if (prev_i==0) {
prev = backtrack(strings, i , MIN2(maxdist,length-i));
prev_i = i;
do_backtrack = 0;
lastf2 = lastf;
energyprev = f3[i];
}
}
else if((thisf < lastf) && (thisf < lastf2) && ((thisf/(n_seq*100)) < -0.01)){ /*?????????*/
do_backtrack = 2;
backtrack_type = 'C';
}
else if (do_backtrack){
if(do_backtrack == 1){
ss = backtrack(strings, i+1 , MIN2(maxdist,length-i)/*+1*/);
energyout = f3[i] - f3[i+strlen(ss)-1];/*??*/
}
else {
ss = backtrack(strings, i+1 , lastj-i-2);
energyout=c[i+1][lastj-(i+1)];
if(dangles){
for (s=0; s<n_seq; s++) {
int type;
type = pair[S[s][i+1]][S[s][lastj-i]]; if (type==0) type=7;
energyout += E_ExtLoop(type, (i>1) ? S5[s][i+1] : -1, S3[s][lastj-i], P);
}
}
else{
for (s=0; s<n_seq; s++) {
int type;
type = pair[S[s][i+1]][S[s][lastj-i]]; if (type==0) type=7;
energyout += E_ExtLoop(type, -1, -1, P);
}
}
}
if((prev_i + strlen(prev) > i+1+strlen(ss)) || (do_backtrack==2)){
char *outstr = (char *)vrna_alloc(sizeof(char) * (strlen(prev)+1));
strncpy(outstr, strings[0]+prev_i-1, strlen(prev));
outstr[strlen(prev)] = '\0';
if (csv==1) printf("%s , %6.2f, %4d, %4d\n",prev, energyprev/(100.*n_seq), prev_i,prev_i + (int)strlen(prev)-1);
/* if(do_backtrack==1)*/
else {
printf("%s (%6.2f) %4d - %4d\n",prev, energyprev/(100.*n_seq), prev_i,prev_i + (int)strlen(prev)-1);
}
free(outstr);
}
free(prev);
prev = ss;
energyprev = energyout;
prev_i = i+1;
do_backtrack = 0;
backtrack_type='F';
}
}
lastf2 = lastf;
lastf = thisf;
lastj = thisj;
if (i==1) {
char *outstr = NULL;
if (prev) {
outstr = (char *)vrna_alloc(sizeof(char) *(strlen(prev) + 1));
strncpy(outstr, strings[0]+prev_i-1, strlen(prev));
outstr[strlen(prev)] = '\0';
if(csv==1)
printf("%s ,%6.2f, %4d, %4d\n", prev, (energyprev)/(100.*n_seq), prev_i,prev_i + (int)strlen(prev)-1);
else{
printf("%s (%6.2f) %4d - %4d\n", prev, (energyprev)/(100.*n_seq), prev_i,prev_i + (int)strlen(prev)-1);
}
}
if ((f3[prev_i] != f3[1]) || !prev){
ss = backtrack(strings, i , maxdist);
if(outstr) free(outstr);
outstr = (char *)vrna_alloc(sizeof(char) * (strlen(ss) + 1));
strncpy(outstr, strings[0], strlen(ss));
outstr[strlen(ss)] = '\0';
printf("%s \n", outstr);
if(csv==1)
printf("%s ,%6.2f ,%4d ,%4d\n", ss, (f3[1]-f3[1 + strlen(ss)-1])/(100.*n_seq), 1, (int)strlen(ss)-1);
else{
printf("%s (%6.2f) %4d - %4d\n", ss, (f3[1]-f3[1 + strlen(ss)-1])/(100.*n_seq), 1, (int)strlen(ss)-1);
}
free(ss);
}
if(prev) free(prev);
if(outstr) free(outstr);
}
}
{
int ii, *FF; /* rotate the auxilliary arrays */
FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
FF = cc1; cc1=cc; cc=FF;
for (j=0; j< maxdist+5; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
if (i<=length-maxdist-4) {
c[i-1] = c[i+maxdist+4]; c[i+maxdist+4] = NULL;
fML[i-1] = fML[i+maxdist+4]; fML[i+maxdist+4]=NULL;
pscore[i-1] = pscore[i+maxdist+4]; pscore[i+maxdist+4] = NULL;
if(i > 1)
make_pscores((const char**) strings, structure, maxdist, i-1);
for(ii=0; ii<maxdist+5; ii++) {
c[i-1][ii] = fML[i-1][ii] = INF;
}
}
}
}
return f3[1];
}
PRIVATE char * backtrack(const char **strings, int start, int maxdist) {
/*------------------------------------------------------------------
trace back through the "c", "f3" 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.
------------------------------------------------------------------*/
sect sector[MAXSECTORS]; /* backtracking sectors */
int i, j, k, energy;
int *type, type_2, tt, n_seq;
/*int bonus;*/
int s=0, ss;
char *structure;
for (s=0; strings[s]!=NULL; s++);
n_seq = s;
type = (int *) vrna_alloc(n_seq*sizeof(int));
s=0;
length = strlen(strings[0]);
sector[++s].i = start;
sector[s].j = MIN2(length, start+maxdist+1);
sector[s].ml = (backtrack_type=='M') ? 1 : ((backtrack_type=='C')?2:0);
structure = (char *) vrna_alloc((MIN2(length-start, maxdist)+3)*sizeof(char));
for (i=0; i<=MIN2(length-start, maxdist); i++) structure[i] = '.';
while (s>0) {
int ml, fij, cij, traced, i1, j1, mm, p, q, jj=0;
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) {
structure[i-start] = '(';
structure[j-start] = ')';
goto repeat1;
}
if (j < i+TURN+1) continue; /* no more pairs in this interval */
fij = (ml)? fML[i][j-i] : f3[i];
if (ml == 0) { /* backtrack in f3 */
if (fij == f3[i+1]) {
sector[++s].i = i+1;
sector[s].j = j;
sector[s].ml = ml;
continue;
}
/* i is paired. Find pairing partner */
for (k=i+TURN+1,traced=0; k<=j; k++) {
int cc;
jj = k+1;
cc = c[i][k-(i)];
if (cc<INF) {
if(dangles){
for (ss=0; ss<n_seq; ss++) {
type[ss] = pair[S[ss][i]][S[ss][k]];
if (type[ss]==0) type[ss]=7;
cc += E_ExtLoop(type[ss], (i>1) ? S5[ss][i] : -1, (k<length) ? S3[ss][k] : -1, P);
}
}
else{
for (ss=0; ss<n_seq; ss++) {
type[ss] = pair[S[ss][i]][S[ss][k]];
if (type[ss]==0) type[ss]=7;
cc += E_ExtLoop(type[ss], -1, -1, P);
}
}
if (fij == cc + f3[k+1]) traced=i;
}
if (traced) break;
}
if (!traced) vrna_message_error("backtrack failed in f3");
if (j==length) { /* backtrack only one component, unless j==length */
sector[++s].i = jj;
sector[s].j = j;
sector[s].ml = ml;
}
i=traced; j=k;
structure[i-start] = '('; structure[j-start] = ')';
goto repeat1;
}
else { /* trace back in fML array */
if (fML[i][j-1-i]+n_seq*P->MLbase == fij) { /* 3' end is unpaired */
sector[++s].i = i;
sector[s].j = j-1;
sector[s].ml = ml;
continue;
}
if (fML[i+1][j-(i+1)]+n_seq*P->MLbase == fij) { /* 5' end is unpaired */
sector[++s].i = i+1;
sector[s].j = j;
sector[s].ml = ml;
continue;
}
cij = c[i][j-i] ;
if(dangles){
for (ss=0; ss<n_seq; ss++) {
tt = pair[S[ss][i]][S[ss][j]];
if (tt==0) tt=7;
cij += E_MLstem(tt, (i>1) ? S5[ss][i] : -1, (j<length) ? S3[ss][j] : -1, P);
}
}
else{
for (ss=0; ss<n_seq; ss++) {
tt = pair[S[ss][i]][S[ss][j]];
if (tt==0) tt=7;
cij += E_MLstem(tt, -1, -1, P);
}
}
if(fij==cij){
/* found a pair */
structure[i-start] = '('; structure[j-start] = ')';
goto repeat1;
}
for (k = i+1+TURN; k <= j-2-TURN; k++)
if (fij == (fML[i][k-i]+fML[k+1][j-(k+1)]))
break;
sector[++s].i = i;
sector[s].j = k;
sector[s].ml = ml;
sector[++s].i = k+1;
sector[s].j = j;
sector[s].ml = ml;
if (k>j-2-TURN) vrna_message_error("backtrack failed in fML");
continue;
}
repeat1:
/*----- begin of "repeat:" -----*/
if (canonical) cij = c[i][j-i];
for (ss=0; ss<n_seq; ss++) {
type[ss] = pair[S[ss][i]][S[ss][j]];
if (type[ss]==0) type[ss] = 7;
}
/* bonus = 0;*/
if (noLonelyPairs)
if (cij == c[i][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[S[ss][j-1]][S[ss][i+1]]; /* j,i not i,j */
if (type_2==0) type_2 = 7;
cij -= P->stack[type[ss]][type_2];
}
cij += pscore[i][j-i];
structure[i+1-start] = '('; structure[j-1-start] = ')';
i++; j--;
canonical=0;
goto repeat1;
}
canonical = 1;
cij += pscore[i][j-i];
{
int cc=0;
for (ss=0; ss<n_seq; ss++){
if((a2s[ss][j-1] - a2s[ss][i]) < 3) cc += 600;
else cc += E_Hairpin(a2s[ss][j-1] - a2s[ss][i], type[ss], S3[ss][i], S5[ss][j], Ss[ss] + a2s[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--) {
if (c[p][q-p]>=INF) continue;
for (ss=energy=0; ss<n_seq; ss++) {
type_2 = pair[S[ss][q]][S[ss][p]]; /* q,p not p,q */
if (type_2==0) type_2 = 7;
energy += E_IntLoop(a2s[ss][p-1] - a2s[ss][i],
a2s[ss][j-1] - a2s[ss][q],
type[ss],
type_2,
S3[ss][i],
S5[ss][j],
S5[ss][p],
S3[ss][q],
P);
}
traced = (cij == energy+c[p][q-p]);
if (traced) {
structure[p-start] = '(';
structure[q-start] = ')';
i = p, j = q;
goto repeat1;
}
}
}
/* end of repeat: --------------------------------------------------*/
/* (i.j) must close a multi-loop */
mm = n_seq*P->MLclosing;
if(dangles){
for (ss=0; ss<n_seq; ss++) {
tt = rtype[type[ss]];
mm += E_MLstem(tt, S5[ss][j],S3[ss][i], P);
}
}
else{
for (ss=0; ss<n_seq; ss++) {
tt = rtype[type[ss]];
mm += E_MLstem(tt, -1, -1, P);
}
}
i1 = i+1; j1 = j-1;
sector[s+1].ml = sector[s+2].ml = 1;
for (k = i+TURN+2; k < j-TURN-2; k++){
if(cij == fML[i+1][k-(i+1)] + fML[k+1][j-1-(k+1)] + mm) break;
}
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");
}
}
if (start+maxdist<length) {
for (i=strlen(structure); i>0 && structure[i-1] == '.'; i--)
structure[i] = '\0';
}
return structure;
}
/*---------------------------------------------------------------------------*/
PRIVATE short *encode_seq(const char *sequence, short *s5, short *s3, char *ss, unsigned short *as){
unsigned int i,l;
short *S;
unsigned short p;
l = strlen(sequence);
S = (short *) vrna_alloc(sizeof(short)*(l+2));
S[0] = (short) l;
s5[0]=s5[1]=0;
/* make numerical encoding of sequence */
for (i=1; i<=l; i++) {
short ctemp = (short)encode_char(toupper(sequence[i-1]));
S[i] = ctemp ;
}
if (oldAliEn) {
/*use alignment sequences in all energy evaluations*/
ss[0]=sequence[0];
for (i=1; i<l; i++) {
s5[i]=S[i-1];
s3[i]=S[i+1];
ss[i]= sequence[i];
as[i]=i;
}
ss[l] = sequence[l];
as[l]=l;
s5[l]=S[l-1];
s3[l]=0;
S[l+1] = S[1];
s5[1]=0;
if (1) {
s5[1]=S[l];
s3[l]=S[1];
ss[l+1]=S[1];
}
return S;
}
else {
if (1) {
for (i=l; i>0; i--) {
char c5;
c5=sequence[i-1];
if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.')) continue;
s5[1] = S[i];
break;
}
for (i=1; i<=l; i++) {
char c3;
c3 = sequence[i-1];
if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.')) continue;
s3[l] = S[i];
break;
}
} else s5[1]=s3[l]=0;
for (i=1,p=0; i<=l; i++) {
char c5;
c5=sequence[i-1];
if ((c5=='-')||(c5=='_')||(c5=='~')||(c5=='.'))
s5[i+1]=s5[i];
else { /* no gap */
ss[p++]=sequence[i-1]; /*start at 0!!*/
s5[i+1]=S[i];
}
as[i]=p;
}
for (i=l; i>=1; i--) {
char c3;
c3=sequence[i-1];
if ((c3=='-')||(c3=='_')||(c3=='~')||(c3=='.'))
s3[i-1]=s3[i];
else
s3[i-1]=S[i];
}
}
return S;
}
PRIVATE double cov_score(const char **AS, int i, int j) {
int n_seq,k,l,s;
double score;
int pfreq[8]={0,0,0,0,0,0,0,0};
for (n_seq=0; AS[n_seq]!=NULL; n_seq++);
for (s=0; s<n_seq; s++) {
int type;
if (S[s][i]==0 && S[s][j]==0) type = 7; /* gap-gap */
else {
if ((AS[s][i] == '~')||(AS[s][j] == '~')) type = 7;
else type = pair[S[s][i]][S[s][j]];
}
pfreq[type]++;
}
if (pfreq[0]*2+pfreq[7]>n_seq)
return NONE;
else
for (k=1,score=0.; k<=6; k++) /* ignore pairtype 7 (gap-gap) */
for (l=k; 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 */
return cv_fact * ((UNIT*score)/n_seq - nc_fact*UNIT*(pfreq[0] + pfreq[7]*0.25));
}
PRIVATE void make_pscores(const char ** AS,
const char *structure, int maxd, int i) {
/* 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 */
int n,j,l;
n=S[0][0]; /* length of seqs */
/*first allocate space:*/
pscore[i]=(int *)vrna_alloc((maxd+5)*sizeof(int));
/* pscore[start]-=start;*/
/*fill pscore[start], too close*/
for (j=i+1; (j<i+TURN+1) && (j<=n); j++) {
pscore[i][j-i] = NONE;
}
for (j=i+TURN+1; ((j<=n) && (j<=i+maxd)); j++) {
pscore[i][j-i] = cov_score(AS, i, j);
}
if (noLonelyPairs) { /* remove unwanted lonely pairs */
int otype=0, ntype=0;
for (j=i+TURN; ((j<n)&&(j<i+maxd)); j++) {
if ((i>1) && (j<n)) otype = cov_score(AS, i-1, j+1);
if (i<n) ntype=pscore[i+1][j-1-(i+1)];
else ntype=NONE;
if ((otype<-4*UNIT)&&(ntype<-4*UNIT)) /* worse than 2 counterex */
pscore[i][j-i] = NONE; /* i.j can only form isolated pairs */
}
}
if (fold_constrained&&(structure!=NULL)) {
int psij, hx, *stack;
stack = (int *) vrna_alloc(sizeof(int)*(n+1));
hx=psij=0;
/* for(hx=0, j=i+TURN; ((j<=i+maxd)&&(j<=n)); j++) {*/
switch (structure[i-1]) {
case 'x': /* can't pair */
for (l=i+TURN+1; l<=i+maxd; l++) pscore[i][l-i] = NONE;
break;
case '(':
hx=1;
psij=1;
for (l=i+1; l<=i+maxd; l++) {
switch (structure[l-1]) {
case '(':
hx++;
pscore[i][l-i] = NONE;
break;
case ')':
hx--;
if (hx!=0) pscore[i][l-i] = NONE;
break;
default:
pscore[i][l-i] = NONE;
}
/* fallthrough */
}
case ')':
for (l=i+TURN+1; l<=i+maxd; l++) pscore[i][l-i] = NONE;
break;
case '>':
for (l=i+TURN+1; l<=i+maxd; l++) pscore[i][l-i] = NONE;
break;
}
if (!psij) for (l=i+1; l<=i+maxd; l++) { /*no '(' constraint on i*/
switch (structure[l-1]) {
case '(':
pscore[i][l-i] = NONE;
break;
case '<':
pscore[i][l-i] = NONE;
break;
case 'x':
pscore[i][l-i] = NONE;
break;
case ')':
pscore[i][l-i] = NONE;
break;
}
}
if (hx!=0) {
vrna_message_error("%s\nunbalanced brackets in constraint string", structure);
}
free(stack);
}
}