--- /dev/null
+/*
+ search for sequences that
+ fold into a given target structure
+
+ c Ivo Hofacker
+ Vienna RNA package
+*/
+/* Last changed Time-stamp: <2009-03-19 10:43:32 ivo> */
+
+#define TDIST 0 /* use tree distance */
+#define PF 1 /* include support for partiton function */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <math.h>
+#include <float.h>
+#if PF
+#include "part_func.h"
+#endif
+#include "fold.h"
+#if TDIST
+#include "dist_vars.h"
+#include "treedist.h"
+#include "RNAstruct.h"
+#endif
+#include "utils.h"
+#include "fold_vars.h"
+#include "pair_mat.h"
+/*@unused@*/
+static char rcsid[] = "$Id: inverse.c,v 1.13 2009/03/19 09:45:28 ivo Exp $";
+#define PUBLIC
+#define PRIVATE static
+PRIVATE double adaptive_walk(char *start, const char *target);
+PRIVATE void shuffle(int *list, int len);
+PRIVATE void make_start(char* start, const char *structure);
+PRIVATE void make_ptable(const char *structure, int *table);
+PRIVATE void make_pairset(void);
+PRIVATE double mfe_cost(const char *, char*, const char *);
+PRIVATE double pf_cost(const char *, char *, const char *);
+PRIVATE char *aux_struct(const char* structure );
+
+/* for backward compatibility, make sure symbolset can hold 20 characters */
+PRIVATE char default_alpha[21] = "AUGC";
+PUBLIC char *symbolset = default_alpha;
+PUBLIC int give_up = 0;
+PUBLIC float final_cost = 0; /* when to stop inverse_pf_fold */
+PUBLIC int inv_verbose=0; /* print out substructure on which inverse_fold() fails */
+
+PRIVATE char pairset[2*MAXALPHA+1];
+PRIVATE int base, npairs;
+PRIVATE int nc2;
+
+/*-------------------------------------------------------------------------*/
+PRIVATE int fold_type;
+#if TDIST
+PRIVATE Tree *T0;
+#endif
+PRIVATE double cost2;
+
+PRIVATE double adaptive_walk(char *start, const char *target)
+{
+#ifdef DUMMY
+ printf("%s\n%s %c\n", start, target, backtrack_type );
+ return 0.;
+#endif
+ int i,j,p,tt,w1,w2, n_pos, len, flag;
+ long walk_len;
+ char *string, *string2, *cstring, *structure, *struct2;
+ int *mut_pos_list, mut_sym_list[MAXALPHA+1], mut_pair_list[2*MAXALPHA+1];
+ int *w1_list, *w2_list, mut_position, symbol, bp;
+ int *target_table, *test_table;
+ char cont;
+ double cost, current_cost, ccost2;
+ double (*cost_function)(const char *, char *, const char *);
+
+ len = strlen(start);
+ if (strlen(target)!=len) {
+ fprintf(stderr, "%s\n%s\n", start, target);
+ nrerror("adaptive_walk: start and target have unequal length");
+ }
+ string = (char *) space(sizeof(char)*(len+1));
+ cstring = (char *) space(sizeof(char)*(len+1));
+ string2 = (char *) space(sizeof(char)*(len+1));
+ structure = (char *) space(sizeof(char)*(len+1));
+ struct2 = (char *) space(sizeof(char)*(len+1));
+ mut_pos_list = (int *) space(sizeof(int)*len);
+ w1_list = (int *) space(sizeof(int)*len);
+ w2_list = (int *) space(sizeof(int)*len);
+ target_table = (int *) space(sizeof(int)*len);
+ test_table = (int *) space(sizeof(int)*len);
+
+ make_ptable(target, target_table);
+
+ for (i=0; i<base; i++) mut_sym_list[i] = i;
+ for (i=0; i<npairs; i++) mut_pair_list[i] = i;
+
+ for (i=0; i<len; i++)
+ string[i] = (islower(start[i]))?toupper(start[i]):start[i];
+ walk_len = 0;
+
+ if (fold_type==0) cost_function = mfe_cost;
+ else cost_function = pf_cost;
+
+ cost = cost_function(string, structure, target);
+
+ if (fold_type==0) ccost2=cost2;
+ else { ccost2 = -1.; cost2=0; }
+
+ strcpy(cstring, string);
+ current_cost = cost;
+
+ if (cost>0) do {
+ cont=0;
+
+ if (fold_type==0) { /* min free energy fold */
+ make_ptable(structure, test_table);
+ for (j=w1=w2=flag=0; j<len; j++)
+ if ((tt=target_table[j])!=test_table[j]) {
+ if ((tt<j)&&(isupper(start[j]))) w1_list[w1++] = j; /* incorrectly paired */
+ if ((flag==0)&&(j>0))
+ if ((target_table[j-1]<j-1)&&isupper(start[j-1]))
+ w2_list[w2++] = j-1; /* adjacent to incorrect position */
+ if (w2>1) if (w2_list[w2-2]==w2_list[w2-1]) w2--;
+
+ flag = 1;
+ } else {
+ if (flag==1) if ((tt<j)&&isupper(start[j]))
+ w2_list[w2++] = j; /* adjacent to incorrect position */
+ flag = 0;
+ }
+ shuffle(w1_list, w1);
+ shuffle(w2_list, w2);
+ for (j=n_pos=0; j<w1; j++) mut_pos_list[n_pos++] = w1_list[j];
+ for (j=0; j<w2; j++) mut_pos_list[n_pos++] = w2_list[j];
+ } else { /* partition_function */
+ for (j=n_pos=0; j<len; j++) if (isupper(start[j]))
+ if (target_table[j]<=j) mut_pos_list[n_pos++] = j;
+ shuffle(mut_pos_list, n_pos);
+ }
+
+ string2[0]='\0';
+ for (mut_position=0; mut_position<n_pos; mut_position++){
+
+ strcpy(string, cstring);
+ shuffle(mut_sym_list, base);
+ shuffle(mut_pair_list, npairs);
+
+ i = mut_pos_list[mut_position];
+
+ if (target_table[i]<0) /* unpaired base */
+ for (symbol=0;symbol<base;symbol++) {
+
+ if(cstring[i]==
+ symbolset[mut_sym_list[symbol]]) continue;
+
+ string[i] = symbolset[mut_sym_list[symbol]];
+
+ cost = cost_function(string, structure, target);
+
+ if ( cost + DBL_EPSILON < current_cost ) break;
+ if (( cost == current_cost)&&(cost2<ccost2)){
+ strcpy(string2, string);
+ strcpy(struct2, structure);
+ ccost2 = cost2;
+ }
+ }
+ else /* paired base */
+ for (bp=0; bp<npairs; bp++) {
+ j = target_table[i];
+ p = mut_pair_list[bp]*2;
+ if ((cstring[i] == pairset[p]) &&
+ (cstring[j] == pairset[p+1]))
+ continue;
+ string[i] = pairset[p];
+ string[j] = pairset[p+1];
+
+ cost = cost_function(string, structure, target);
+
+ if ( cost < current_cost ) break;
+ if (( cost == current_cost)&&(cost2<ccost2)){
+ strcpy(string2, string);
+ strcpy(struct2, structure);
+ ccost2 = cost2;
+ }
+ }
+
+ if ( cost < current_cost ) {
+ strcpy(cstring, string);
+ current_cost = cost;
+ ccost2 = cost2;
+ walk_len++;
+ if (cost>0) cont=1;
+ break;
+ }
+ }
+ if ((current_cost>0)&&(cont==0)&&(string2[0])) {
+ /* no mutation that decreased cost was found,
+ but the the sequence in string2 decreases cost2 while keeping
+ cost constant */
+ strcpy(cstring, string2);
+ strcpy(structure, struct2);
+ nc2++; cont=1;
+ }
+ } while (cont);
+
+ for (i=0; i<len; i++) if (isupper(start[i])) start[i]=cstring[i];
+
+#if TDIST
+ if (fold_type==0) { free_tree(T0); T0=NULL; }
+#endif
+ free(test_table);
+ free(target_table);
+ free(mut_pos_list);
+ free(w2_list);
+ free(w1_list);
+ free(struct2);
+ free(structure);
+ free(string2);
+ free(cstring);
+ free(string);
+
+ return current_cost;
+}
+
+/*-------------------------------------------------------------------------*/
+
+/* shuffle produces a ronaom list by doing len exchanges */
+PRIVATE void shuffle(int *list, int len)
+{
+ int i, rn;
+
+ for (i=0;i<len;i++) {
+ int temp;
+ rn = i + (int) (urn()*(len-i)); /* [i..len-1] */
+ /* swap element i and rn */
+ temp = list[i];
+ list[i] = list[rn];
+ list[rn] = temp;
+ }
+}
+
+/*-------------------------------------------------------------------------*/
+
+PRIVATE void make_ptable(const char *structure, int *table)
+{
+ int i,j,hx;
+ int *stack;
+
+ hx=0;
+ stack = (int *) space(sizeof(int)*(strlen(structure)+1));
+
+ for (i=0; i<strlen(structure); i++) {
+ switch (structure[i]) {
+ case '.':
+ table[i]= -1;
+ break;
+ case '(':
+ stack[hx++]=i;
+ break;
+ case ')':
+ j = stack[--hx];
+ if (hx<0) {
+ fprintf(stderr, "%s\n", structure);
+ nrerror("unbalanced brackets in make_ptable");
+ }
+ table[i]=j;
+ table[j]=i;
+ break;
+ }
+ }
+ if (hx!=0) {
+ fprintf(stderr, "%s\n", structure);
+ nrerror("unbalanced brackets in make_ptable");
+ }
+ free(stack);
+}
+
+/*-------------------------------------------------------------------------*/
+
+#define WALK(i,j) \
+ strncpy(wstruct, structure+i, j-i+1); \
+ wstruct[j-i+1]='\0'; \
+ strncpy(wstring, string+i, j-i+1); \
+ wstring[j-i+1]='\0'; \
+ dist=adaptive_walk(wstring, wstruct); \
+ strncpy(string+i, wstring, j-i+1); \
+ if ((dist>0)&&(give_up)) goto adios
+
+PUBLIC float inverse_fold(char *start, char *structure)
+{
+ int i, j, jj, len, o;
+ int *pt;
+ char *string, *wstring, *wstruct, *aux;
+ double dist=0;
+
+ nc2 = j = o = fold_type = 0;
+
+ len = strlen(structure);
+ if (strlen(start)!=len) {
+ fprintf(stderr, "%s\n%s\n", start, structure);
+ nrerror("inverse_fold: start and structure have unequal length");
+ }
+ string = (char *) space(len+1);
+ wstring = (char *) space(len+1);
+ wstruct = (char *) space(len+1);
+ pt = (int *) space(sizeof(int)*(len+2));
+ pt[len] = len+1;
+
+ aux = aux_struct(structure);
+ strcpy(string, start);
+ make_pairset();
+ make_start(string, structure);
+
+ make_ptable(structure, pt);
+
+ while (j<len) {
+ while ((j<len)&&(structure[j]!=')')) {
+ if (aux[j]=='[') o++;
+ if (aux[j]==']') o--;
+ j++;
+ }
+ i=j;
+ while ((i>0) && structure[--i]!='(');
+ if (structure[i]=='.') { /* no pair found -> open chain */
+ WALK(0,len-1);
+ }
+
+ if (aux[i]!='[') { i--; j++;}
+ while (pt[j]==i) {
+ backtrack_type='C';
+ if (aux[i]!='[') {
+ while (aux[--i]!='[');
+ while (aux[++j]!=']');
+ /* WALK(i,j); */
+ }
+ WALK(i,j);
+ o--;
+ jj = j; i--;
+ while (aux[++j]=='.');
+ while ((i>=0)&&(aux[i]=='.')) i--;
+ if (pt[j]!=i) {
+ backtrack_type = (o==0)? 'F' : 'M';
+ if (j-jj>8) { WALK((i+1),(jj)); }
+ WALK((i+1), (j-1));
+ while ((i>=0) &&(aux[i]==']')) {
+ i=pt[i]-1;
+ while ((i>=0)&&(aux[i]=='.')) i--;
+ WALK((i+1), (j-1));
+ }
+ }
+ }
+ }
+ adios:
+ backtrack_type='F';
+ if ((dist>0)&&(inv_verbose)) printf("%s\n%s\n", wstring, wstruct);
+ /*if ((dist==0)||(give_up==0))*/ strcpy(start, string);
+ free(wstring); free(wstruct);
+ free(string); free(aux);
+ free(pt);
+/* if (dist>0) printf("%3d \n", nc2); */
+ return dist;
+}
+
+/*-------------------------------------------------------------------------*/
+
+PUBLIC float inverse_pf_fold(char *start, char *target)
+{
+ double dist;
+ int dang;
+
+ dang=dangles;
+ if (dangles!=0) dangles=2;
+
+ update_fold_params(); /* make sure there is a valid pair matrix */
+ make_pairset();
+ make_start(start, target);
+ fold_type=1;
+ do_backtrack = 0;
+ dist = adaptive_walk(start, target);
+ dangles=dang;
+ return (dist+final_cost);
+}
+
+/*-------------------------------------------------------------------------*/
+
+PRIVATE void make_start(char* start, const char *structure)
+{
+ int i,j,k,l,r,length;
+ int *table, *S, sym[MAXALPHA], ss;
+
+ length=strlen(start);
+ table = (int *) space(sizeof(int)*length);
+ S = (int *) space(sizeof(int)*length);
+
+ make_ptable(structure, table);
+ for (i=0; i<strlen(start); i++) S[i] = encode_char(toupper(start[i]));
+ for (i=0; i<strlen(symbolset); i++) sym[i] = i;
+
+ for (k=0; k<length; k++) {
+ if (table[k]<k) continue;
+ if (((urn()<0.5) && isupper(start[k])) ||
+ islower(start[table[k]])) {
+ i = table[k]; j = k;
+ } else {
+ i = k; j = table[k];
+ }
+
+ if (!pair[S[i]][S[j]]) { /* make a valid pair by mutating j */
+ shuffle(sym, (int) base);
+ for (l=0; l<base; l++) {
+ ss = encode_char(symbolset[sym[l]]);
+ if (pair[S[i]][ss]) break;
+ }
+ if (l==base) { /* nothing pairs start[i] */
+ r = 2*int_urn(0, npairs-1);
+ start[i] = pairset[r];
+ start[j] = pairset[r+1];
+ } else start[j] = symbolset[sym[l]];
+ }
+ }
+ free(table);
+ free(S);
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void make_pairset(void)
+{
+ int i,j;
+ int sym[MAXALPHA];
+
+ make_pair_matrix();
+ base = strlen(symbolset);
+
+ for (i=0; i< base; i++) sym[i] = encode_char(symbolset[i]);
+
+ for (i=npairs=0; i< base; i++)
+ for (j=0; j<base; j++)
+ if (pair[sym[i]][sym[j]]) {
+ pairset[npairs++] = symbolset[i];
+ pairset[npairs++] = symbolset[j];
+ }
+ npairs /= 2;
+ if (npairs==0) nrerror("No pairs in this alphabet!");
+}
+/*---------------------------------------------------------------------------*/
+
+PRIVATE double mfe_cost(const char *string, char *structure, const char *target)
+{
+#if TDIST
+ Tree *T1;
+ char *xstruc;
+#endif
+ double energy, distance;
+
+ if (strlen(string)!=strlen(target)) {
+ fprintf(stderr, "%s\n%s\n", string, target);
+ nrerror("unequal length in mfe_cost");
+ }
+ energy = fold(string, structure);
+#if TDIST
+ if (T0 == NULL) {
+ xstruc = expand_Full(target);
+ T0=make_tree(xstruc);
+ free(xstruc);
+ }
+
+ xstruc = expand_Full(structure);
+ T1=make_tree(xstruc);
+ distance = tree_edit_distance(T0,T1);
+ free(xstruc);
+ free_tree(T1);
+#else
+ distance = (double) bp_distance(target, structure);
+#endif
+ cost2 = energy_of_structure(string, target, 0) - energy;
+ return (double) distance;
+}
+/*---------------------------------------------------------------------------*/
+
+PRIVATE double pf_cost(const char *string, char *structure, const char *target)
+{
+#if PF
+ double f, e;
+
+ f = pf_fold(string, structure);
+ e = energy_of_structure(string, target, 0);
+ return (double) (e-f-final_cost);
+#else
+ nrerror("this version not linked with pf_fold");
+ return 0;
+#endif
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE char *aux_struct(const char* structure )
+{
+ int *match_paren;
+ int i, o, p;
+ char *string;
+
+ string = (char *) space(sizeof(char)*(strlen(structure)+1));
+ match_paren = (int *) space(sizeof(int)*(strlen(structure)/2+1));
+ strcpy(string, structure);
+
+ i = o = 0;
+ while (string[i]) {
+ switch (string[i]) {
+ case '.': break;
+ case '(':
+ match_paren[++o]=i;
+ break;
+ case ')':
+ p=i;
+ while ((string[p+1]==')')&&(match_paren[o-1]==match_paren[o]-1)) {
+ p++; o--;
+ }
+ string[p]=']';
+ i=p;
+ string[match_paren[o]]='[';
+ o--;
+ break;
+ default:
+ nrerror("Junk in structure at aux_structure\n");
+ }
+ i++;
+ }
+ free(match_paren);
+ return(string);
+}
git://source.jalview.org / jabaws.git / blobdiff