binaries/src/ViennaRNA/Progs/RNAinverse.c

   1 /*
   2             Interactive access to inverse folding routines
   3                     c Ivo Hofacker, Peter Stadler
   4                           Vienna RNA Package
   5 */
   6 /* Last changed Time-stamp: <2000-09-28 14:58:05 ivo> */
   7 #include <stdio.h>
   8 #include <stdlib.h>
   9 #include <math.h>
  10 #include <ctype.h>
  11 #include <unistd.h>
  12 #include <string.h>
  13 #include "inverse.h"
  14 #include "fold_vars.h"
  15 #include "fold.h"
  16 #include "part_func.h"
  17 #include "utils.h"
  18 #include "read_epars.h"
  19 #include "RNAinverse_cmdl.h"
  20
  21 #ifdef dmalloc
  22 #include  "/usr/local/include/dmalloc.h"
  23 #define space(X) calloc(1,(X))
  24 #endif
  25
  26 /*@unused@*/
  27 static char rcsid[] = "$Id: RNAinverse.c,v 1.11 2000/09/28 12:59:21 ivo Rel $";
  28
  29 #define  REPEAT_DEFAULT  100
  30
  31 PRIVATE void usage(void);
  32 extern int inv_verbose;
  33
  34 int main(int argc, char *argv[]){
  35   struct  RNAinverse_args_info args_info;
  36   int     input_type;
  37   char    *input_string, *start, *structure, *rstart, *str2, *line;
  38   char    *ParamFile=NULL, *c, *ns_bases;
  39   int     i,j, length, l, hd, sym;
  40   double  energy=0., kT;
  41   int     pf, mfe, istty;
  42   int     repeat, found;
  43
  44   dangles = 2;
  45   do_backtrack = 0;
  46   pf = 0;
  47   mfe = 1;
  48   repeat = 0;
  49   input_type = 0;
  50   input_string = ns_bases = NULL;
  51   init_rand();
  52
  53   /*
  54   #############################################
  55   # check the command line parameters
  56   #############################################
  57   */
  58   if(RNAinverse_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
  59   /* temperature */
  60   if(args_info.temp_given)        temperature = args_info.temp_arg;
  61   /* do not take special tetra loop energies into account */
  62   if(args_info.noTetra_given)     tetra_loop=0;
  63   /* set dangle model */
  64   if(args_info.dangles_given){
  65     if((args_info.dangles_arg < 0) || (args_info.dangles_arg > 3))
  66       warn_user("required dangle model not implemented, falling back to default dangles=2");
  67     else
  68       dangles = args_info.dangles_arg;
  69   }
  70   /* do not allow wobble pairs (GU) */
  71   if(args_info.noGU_given)        noGU = 1;
  72   /* do not allow weak closing pairs (AU,GU) */
  73   if(args_info.noClosingGU_given) no_closingGU = 1;
  74   /* set energy model */
  75   if(args_info.energyModel_given) energy_set = args_info.energyModel_arg;
  76   /* take another energy parameter set */
  77   if(args_info.paramFile_given)   ParamFile = strdup(args_info.paramFile_arg);
  78   /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
  79   if(args_info.nsp_given)         ns_bases = strdup(args_info.nsp_arg);
  80   /* alter the alphabet */
  81   if(args_info.alphabet_given){
  82     symbolset=args_info.alphabet_arg;
  83     /* symbolset should only have uppercase characters */
  84     for (l = 0; l < (int)strlen(symbolset); l++)
  85       symbolset[l] = toupper(symbolset[l]);
  86   }
  87   /* set function for optimization */
  88   if(args_info.function_given){
  89     if(strlen(args_info.function_arg) > 2){
  90       RNAinverse_cmdline_parser_print_help(); exit(EXIT_FAILURE);
  91     }
  92     else{
  93       if((*args_info.function_arg == 'm') || (*(args_info.function_arg+1) == 'm')) mfe = 1;
  94       if((*args_info.function_arg == 'p') || (*(args_info.function_arg+1) == 'p')) pf = 1;
  95     }
  96   }
  97   /* set repeat */
  98   if(args_info.repeat_given)      repeat = args_info.repeat_arg;
  99   /* set final cost */
 100   if(args_info.final_given)       final_cost = args_info.final_arg;
 101   /* do we wannabe verbose */
 102   if(args_info.verbose_given)     inv_verbose = 1;
 103
 104   /* free allocated memory of command line data structure */
 105   RNAinverse_cmdline_parser_free (&args_info);
 106
 107   kT = (temperature+273.15)*1.98717/1000.0;
 108
 109   istty = (isatty(fileno(stdout))&&isatty(fileno(stdin)));
 110
 111   if (ParamFile!=NULL)
 112     read_parameter_file(ParamFile);
 113
 114   give_up = (repeat<0);
 115
 116   do {
 117     /*
 118     ########################################################
 119     # handle user input from 'stdin'
 120     ########################################################
 121     */
 122     if(istty)
 123       print_tty_input_seq_str("Input structure & start string\n"
 124                               "(lower case letters for const positions) and 0 or empty line for random start string\n");
 125
 126     input_type = get_multi_input_line(&input_string, 0);
 127     /* we are waiting for a structure (i.e. something like a constraint) so we skip all sequences, fasta-headers and misc lines */
 128     while(input_type & (VRNA_INPUT_SEQUENCE | VRNA_INPUT_MISC | VRNA_INPUT_FASTA_HEADER)){
 129       if(!istty && (input_type & VRNA_INPUT_FASTA_HEADER)) printf(">%s\n", input_string);
 130       free(input_string); input_string = NULL;
 131       input_type = get_multi_input_line(&input_string, 0);
 132     }
 133     if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)) break;
 134
 135     if(input_type & (VRNA_INPUT_CONSTRAINT)){
 136       structure = (char *)space(sizeof(char) * (strlen(input_string) + 1));
 137       (void)sscanf(input_string, "%s", structure); /* scanf gets rid of trailing junk */
 138       length = (int)strlen(structure);
 139       free(input_string); input_string = NULL;
 140       input_type = get_multi_input_line(&input_string, VRNA_INPUT_NOSKIP_BLANK_LINES | VRNA_INPUT_NOSKIP_COMMENTS);
 141     }
 142     if(input_type & VRNA_INPUT_QUIT) break;
 143
 144     start = (char *)space(sizeof(char) * (length+1));
 145     /* now we assume to get a sequence (input_string may be empty as well) */
 146     if(input_type & VRNA_INPUT_SEQUENCE){
 147       (void)strncpy(start, input_string, length);
 148       start[length] = '\0';
 149       free(input_string); input_string = NULL;
 150     }
 151     /* fallback to empty start sequence */
 152     else start[0] = '\0';
 153
 154     /*
 155     ########################################################
 156     # done with 'stdin' handling
 157     ########################################################
 158     */
 159
 160     if (ns_bases != NULL) {
 161       nonstandards = space(33);
 162       c=ns_bases;
 163       i=sym=0;
 164       if (*c=='-') {
 165         sym=1; c++;
 166       }
 167       while (*c!='\0') {
 168         if (*c!=',') {
 169           nonstandards[i++]=*c++;
 170           nonstandards[i++]=*c;
 171           if ((sym)&&(*c!=*(c-1))) {
 172             nonstandards[i++]=*c;
 173             nonstandards[i++]=*(c-1);
 174           }
 175         }
 176         c++;
 177       }
 178     }
 179
 180     str2 = (char *) space((unsigned)length+1);
 181     if (istty) printf("length = %d\n", length);
 182
 183     if (repeat!=0) found = (repeat>0)? repeat : (-repeat);
 184     else found = 1;
 185
 186     /* initialize_fold(length); <- obsolete (hopefully commenting this out does not affect anything crucial ;) */
 187
 188     rstart = (char *) space((unsigned)length+1);
 189     while(found>0) {
 190       char *string;
 191       string = (char *) space((unsigned)length+1);
 192       strcpy(string, start);
 193       for (i=0; i<length; i++) {
 194         /* lower case characters are kept fixed, any other character
 195            not in symbolset is replaced by a random character */
 196         if (islower(string[i])) continue;
 197
 198         if (string[i]=='\0' || (strchr(symbolset,string[i])==NULL))
 199           string[i]=symbolset[int_urn(0,strlen(symbolset)-1)];
 200       }
 201       strcpy(rstart, string); /* remember start string */
 202
 203       if (mfe) {
 204         energy = inverse_fold(string, structure);
 205         if( (repeat>=0) || (energy<=0.0) ) {
 206           found--;
 207           hd = hamming(rstart, string);
 208           printf("%s  %3d", string, hd);
 209           if (energy>0) { /* no solution found */
 210             printf("   d= %g\n", energy);
 211             if(istty) {
 212               energy = fold(string,str2);
 213               printf("%s\n", str2);
 214             }
 215           } else printf("\n");
 216         }
 217       }
 218       if (pf) {
 219         if (!(mfe && give_up && (energy>0))) {
 220           /* unless we gave up in the mfe part */
 221           double prob, min_en, sfact=1.07;
 222
 223           /* get a reasonable pf_scale */
 224           min_en = fold(string,str2);
 225           pf_scale = exp(-(sfact*min_en)/kT/length);
 226           /* init_pf_fold(length); <- obsolete (hopefully commenting this out does not affect anything crucial ;) */
 227
 228           energy = inverse_pf_fold(string, structure);
 229           prob = exp(-energy/kT);
 230           hd = hamming(rstart, string);
 231           printf("%s  %3d  (%g)\n", string, hd, prob);
 232           free_pf_arrays();
 233         }
 234         if (!mfe) found--;
 235       }
 236       (void) fflush(stdout);
 237       free(string);
 238     }
 239     free(rstart);
 240     free_arrays();
 241
 242     free(structure);
 243     free(str2);
 244     free(start);
 245     (void) fflush(stdout);
 246   } while (1);
 247   return 0;
 248 }
 249
 250
 251 PRIVATE void usage(void)
 252 {
 253   nrerror("usage: RNAinverse [-F[mp]] [-a ALPHABET] [-R [repeats]] [-f final]\n"
 254           "                  [-T temp] [-4] [-d[2]] [-noGU] [-P paramfile] [-e e_set] [-v]");
 255 }