binaries/src/ViennaRNA/Progs/RNAsubopt.c

   1 /* Last changed Time-stamp: <2008-12-03 16:38:01 ivo> */
   2 /*
   3                 Ineractive access to suboptimal folding
   4
   5                            c Ivo L Hofacker
   6                           Vienna RNA package
   7 */
   8 #include <config.h>
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <unistd.h>
  12 #include <math.h>
  13 #include <ctype.h>
  14 #include <string.h>
  15 #include "part_func.h"
  16 #include "fold.h"
  17 #include "cofold.h"
  18 #include "fold_vars.h"
  19 #include "utils.h"
  20 #include "read_epars.h"
  21 #include "subopt.h"
  22 #include "params.h"
  23 #include "RNAsubopt_cmdl.h"
  24
  25 /*@unused@*/
  26 static char UNUSED rcsid[] = "$Id: RNAsubopt.c,v 1.20 2008/12/03 16:55:44 ivo Exp $";
  27
  28 PRIVATE char *tokenize(char *line);
  29 PRIVATE void putoutzuker(SOLUTION* zukersolution);
  30
  31 int main(int argc, char *argv[]){
  32   struct          RNAsubopt_args_info args_info;
  33   unsigned int    input_type;
  34   unsigned int    rec_type, read_opt;
  35   char            fname[FILENAME_MAX_LENGTH], *c, *input_string, *rec_sequence, *rec_id, **rec_rest, *orig_sequence;
  36   char            *cstruc, *structure, *ParamFile, *ns_bases;
  37   int             i, length, l, cl, sym, istty;
  38   double          deltaf, deltap, betaScale;
  39   int             delta, n_back, noconv, circular, dos, zuker;
  40   paramT          *P;
  41   pf_paramT       *pf_parameters;
  42   model_detailsT  md;
  43
  44   do_backtrack  = 1;
  45   dangles       = 2;
  46   betaScale     = 1.;
  47   delta         = 100;
  48   deltap = n_back = noconv = circular = dos = zuker = 0;
  49   rec_type      = read_opt = 0;
  50   rec_id        = rec_sequence = orig_sequence = NULL;
  51   rec_rest      = NULL;
  52   input_string  = c = cstruc = structure = ParamFile = ns_bases = NULL;
  53   pf_parameters = NULL;
  54   P             = NULL;
  55
  56   set_model_details(&md);
  57
  58   /*
  59   #############################################
  60   # check the command line parameters
  61   #############################################
  62   */
  63   if(RNAsubopt_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
  64   /* temperature */
  65   if(args_info.temp_given)        temperature = args_info.temp_arg;
  66   /* structure constraint */
  67   if(args_info.constraint_given)  fold_constrained=1;
  68   /* do not take special tetra loop energies into account */
  69   if(args_info.noTetra_given)     md.special_hp = tetra_loop=0;
  70   /* set dangle model */
  71   if(args_info.dangles_given){
  72     if((args_info.dangles_arg < 0) || (args_info.dangles_arg > 3))
  73       warn_user("required dangle model not implemented, falling back to default dangles=2");
  74     else
  75       md.dangles = dangles = args_info.dangles_arg;
  76   }
  77   /* do not allow weak pairs */
  78   if(args_info.noLP_given)        md.noLP = noLonelyPairs = 1;
  79   /* do not allow wobble pairs (GU) */
  80   if(args_info.noGU_given)        md.noGU = noGU = 1;
  81   /* do not allow weak closing pairs (AU,GU) */
  82   if(args_info.noClosingGU_given) md.noGUclosure = no_closingGU = 1;
  83   /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
  84   if(args_info.noconv_given)      noconv = 1;
  85   /* take another energy parameter set */
  86   if(args_info.paramFile_given)   ParamFile = strdup(args_info.paramFile_arg);
  87   /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
  88   if(args_info.nsp_given)         ns_bases = strdup(args_info.nsp_arg);
  89   /* energy range */
  90   if(args_info.deltaEnergy_given) delta = (int) (0.1+args_info.deltaEnergy_arg*100);
  91   /* energy range after post evaluation */
  92   if(args_info.deltaEnergyPost_given) deltap = args_info.deltaEnergyPost_arg;
  93   /* sorted output */
  94   if(args_info.sorted_given)      subopt_sorted = 1;
  95   /* assume RNA sequence to be circular */
  96   if(args_info.circ_given)        circular=1;
  97   /* stochastic backtracking */
  98   if(args_info.stochBT_given){
  99     n_back = args_info.stochBT_arg;
 100     init_rand();
 101   }
 102   if(args_info.betaScale_given)   betaScale = args_info.betaScale_arg;
 103   /* density of states */
 104   if(args_info.dos_given){
 105     dos = 1;
 106     print_energy = -999999;
 107   }
 108   /* logarithmic multiloop energies */
 109   if(args_info.logML_given) md.logML = logML = 1;
 110   /* zuker subopts */
 111   if(args_info.zuker_given) zuker = 1;
 112
 113   if(zuker){
 114     if(circular){
 115       warn_user("Sorry, zuker subopts not yet implemented for circfold");
 116       RNAsubopt_cmdline_parser_print_help();
 117       exit(1);
 118     }
 119     else if(n_back>0){
 120       warn_user("Can't do zuker subopts and stochastic subopts at the same time");
 121       RNAsubopt_cmdline_parser_print_help();
 122       exit(1);
 123     }
 124   }
 125
 126   /* free allocated memory of command line data structure */
 127   RNAsubopt_cmdline_parser_free(&args_info);
 128
 129   /*
 130   #############################################
 131   # begin initializing
 132   #############################################
 133   */
 134
 135   if (ParamFile != NULL) read_parameter_file(ParamFile);
 136
 137   if (ns_bases != NULL) {
 138     nonstandards = space(33);
 139     c=ns_bases;
 140     i=sym=0;
 141     if (*c=='-') {
 142       sym=1; c++;
 143     }
 144     while (*c!='\0') {
 145       if (*c!=',') {
 146         nonstandards[i++]=*c++;
 147         nonstandards[i++]=*c;
 148         if ((sym)&&(*c!=*(c-1))) {
 149           nonstandards[i++]=*c;
 150           nonstandards[i++]=*(c-1);
 151         }
 152       }
 153       c++;
 154     }
 155   }
 156
 157   istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
 158
 159   /* print user help if we get input from tty */
 160   if(istty){
 161     if(!zuker)
 162       printf("Use '&' to connect 2 sequences that shall form a complex.\n");
 163     if(fold_constrained){
 164       print_tty_constraint(VRNA_CONSTRAINT_DOT | VRNA_CONSTRAINT_X | VRNA_CONSTRAINT_ANG_BRACK | VRNA_CONSTRAINT_RND_BRACK);
 165       print_tty_input_seq_str("Input sequence (upper or lower case) followed by structure constraint\n");
 166     }
 167     else print_tty_input_seq();
 168   }
 169
 170   /* set options we wanna pass to read_record */
 171   if(istty)             read_opt |= VRNA_INPUT_NOSKIP_BLANK_LINES;
 172   if(!fold_constrained) read_opt |= VRNA_INPUT_NO_REST;
 173
 174   P = get_scaled_parameters(temperature, md);
 175
 176   /*
 177   #############################################
 178   # main loop: continue until end of file
 179   #############################################
 180   */
 181   while(
 182     !((rec_type = read_record(&rec_id, &rec_sequence, &rec_rest, read_opt))
 183         & (VRNA_INPUT_ERROR | VRNA_INPUT_QUIT))){
 184
 185     /*
 186     ########################################################
 187     # init everything according to the data we've read
 188     ########################################################
 189     */
 190     if(rec_id){
 191       /* if(!istty) printf("%s\n", rec_id); */
 192       (void) sscanf(rec_id, ">%" XSTR(FILENAME_ID_LENGTH) "s", fname);
 193     }
 194     else fname[0] = '\0';
 195
 196     cut_point = -1;
 197
 198     rec_sequence  = tokenize(rec_sequence); /* frees input_string and sets cut_point */
 199     length    = (int) strlen(rec_sequence);
 200     structure = (char *) space((unsigned) length+1);
 201
 202     /* parse the rest of the current dataset to obtain a structure constraint */
 203     if(fold_constrained){
 204       cstruc = NULL;
 205       int cp = cut_point;
 206       unsigned int coptions = (rec_id) ? VRNA_CONSTRAINT_MULTILINE : 0;
 207       coptions |= VRNA_CONSTRAINT_DOT | VRNA_CONSTRAINT_X | VRNA_CONSTRAINT_ANG_BRACK | VRNA_CONSTRAINT_RND_BRACK;
 208       getConstraint(&cstruc, (const char **)rec_rest, coptions);
 209       cstruc = tokenize(cstruc);
 210       if(cut_point != cp) nrerror("cut point in sequence and structure constraint differs");
 211       cl = (cstruc) ? (int)strlen(cstruc) : 0;
 212
 213       if(cl == 0)           warn_user("structure constraint is missing");
 214       else if(cl < length)  warn_user("structure constraint is shorter than sequence");
 215       else if(cl > length)  nrerror("structure constraint is too long");
 216
 217       if(cstruc) strncpy(structure, cstruc, sizeof(char)*(cl+1));
 218     }
 219
 220     /* convert DNA alphabet to RNA if not explicitely switched off */
 221     if(!noconv) str_DNA2RNA(rec_sequence);
 222     /* store case-unmodified sequence */
 223     orig_sequence = strdup(rec_sequence);
 224     /* convert sequence to uppercase letters only */
 225     str_uppercase(rec_sequence);
 226
 227     if(istty){
 228       if (cut_point == -1)
 229         printf("length = %d\n", length);
 230       else
 231         printf("length1 = %d\nlength2 = %d\n", cut_point-1, length-cut_point+1);
 232     }
 233
 234     /*
 235     ########################################################
 236     # begin actual computations
 237     ########################################################
 238     */
 239
 240     if((logML != 0 || dangles==1 || dangles==3) && dos == 0)
 241       if(deltap<=0) deltap = delta/100. + 0.001;
 242     if (deltap>0)
 243       print_energy = deltap;
 244
 245     /* stochastic backtracking */
 246     if(n_back>0){
 247       double mfe, kT;
 248       char *ss;
 249       if (fname[0] != '\0') printf(">%s\n", fname);
 250
 251       st_back=1;
 252       ss = (char *) space(strlen(rec_sequence)+1);
 253       strncpy(ss, structure, length);
 254       mfe = fold_par(rec_sequence, ss, P, fold_constrained, circular);
 255       kT = (betaScale*((temperature+K0)*GASCONST))/1000.; /* in Kcal */
 256       pf_scale = exp(-(1.03*mfe)/kT/length);
 257       strncpy(ss, structure, length);
 258
 259       pf_parameters = get_boltzmann_factors(temperature, betaScale, md, pf_scale);
 260       /* ignore return value, we are not interested in the free energy */
 261       (void) pf_fold_par(rec_sequence, ss, pf_parameters, do_backtrack, fold_constrained, circular);
 262       free(ss);
 263       for (i=0; i<n_back; i++) {
 264         char *s;
 265         s =(circular) ? pbacktrack_circ(rec_sequence) : pbacktrack(rec_sequence);
 266         printf("%s\n", s);
 267         free(s);
 268       }
 269       free_pf_arrays();
 270       free(pf_parameters);
 271     }
 272     /* normal subopt */
 273     else if(!zuker){
 274       /* first lines of output (suitable  for sort +1n) */
 275       if (fname[0] != '\0')
 276         printf("> %s [%d]\n", fname, delta);
 277
 278       subopt_par(rec_sequence, structure, P, delta, fold_constrained, circular, stdout);
 279       if (dos) {
 280         int i;
 281         for (i=0; i<= MAXDOS && i<=delta/10; i++) {
 282           printf("%4d %6d\n", i, density_of_states[i]);
 283         }
 284       }
 285     }
 286     /* Zuker suboptimals */
 287     else{
 288       SOLUTION *zr;
 289       int i;
 290       if (fname[0] != '\0') printf(">%s\n%s\n", fname, rec_sequence);
 291
 292       if (cut_point!=-1) {
 293         nrerror("Sorry, zuker subopts not yet implemented for cofold\n");
 294       }
 295       zr = zukersubopt(rec_sequence);
 296       putoutzuker(zr);
 297       (void)fflush(stdout);
 298       for (i=0; zr[i].structure; i++) {
 299         free(zr[i].structure);
 300       }
 301       free(zr);
 302     }
 303
 304     (void)fflush(stdout);
 305
 306     /* clean up */
 307     if(cstruc) free(cstruc);
 308     if(rec_id) free(rec_id);
 309     free(rec_sequence);
 310     free(orig_sequence);
 311     free(structure);
 312     /* free the rest of current dataset */
 313     if(rec_rest){
 314       for(i=0;rec_rest[i];i++) free(rec_rest[i]);
 315       free(rec_rest);
 316     }
 317     rec_id = rec_sequence = orig_sequence = structure = cstruc = NULL;
 318     rec_rest = NULL;
 319
 320     /* print user help for the next round if we get input from tty */
 321     if(istty){
 322       if(!zuker)
 323         printf("Use '&' to connect 2 sequences that shall form a complex.\n");
 324       if(fold_constrained){
 325         print_tty_constraint(VRNA_CONSTRAINT_DOT | VRNA_CONSTRAINT_X | VRNA_CONSTRAINT_ANG_BRACK | VRNA_CONSTRAINT_RND_BRACK);
 326         print_tty_input_seq_str("Input sequence (upper or lower case) followed by structure constraint\n");
 327       }
 328       else print_tty_input_seq();
 329     }
 330   }
 331   free(P);
 332   return 0;
 333 }
 334
 335 PRIVATE char *tokenize(char *line)
 336 {
 337   char *pos, *copy;
 338   int cut = -1;
 339
 340   copy = (char *) space(strlen(line)+1);
 341   (void) sscanf(line, "%s", copy);
 342   pos = strchr(copy, '&');
 343   if (pos) {
 344     cut = (int) (pos-copy)+1;
 345     if (cut >= strlen(copy)) cut = -1;
 346     if (strchr(pos+1, '&')) nrerror("more than one cut-point in input");
 347     for (;*pos;pos++) *pos = *(pos+1); /* splice out the & */
 348   }
 349   if (cut > -1) {
 350     if (cut_point==-1) cut_point = cut;
 351     else if (cut_point != cut) {
 352       fprintf(stderr,"cut_point = %d cut = %d\n", cut_point, cut);
 353       nrerror("Sequence and Structure have different cut points.");
 354     }
 355   }
 356
 357   free(line);
 358   return copy;
 359 }
 360 PRIVATE void putoutzuker(SOLUTION* zukersolution) {
 361   int i;
 362   printf("%s [%.2f]\n",zukersolution[0].structure,zukersolution[0].energy/100.);
 363   for(i=1; zukersolution[i].structure; i++) {
 364     printf("%s [%.2f]\n", zukersolution[i].structure, zukersolution[i].energy/100.);
 365   }
 366   return;
 367 }