binaries/src/ViennaRNA/Progs/RNAplfold.c

   1 /* Last changed Time-stamp: <2008-07-02 17:10:04 berni> */
   2 /*
   3                   Ineractive Access to folding Routines
   4
   5                   c Ivo L Hofacker
   6                   Vienna RNA package
   7 */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <math.h>
  12 #include <ctype.h>
  13 #include <unistd.h>
  14 #include <string.h>
  15 #include "fold.h"
  16 #include "part_func.h"
  17 #include "fold_vars.h"
  18 #include "utils.h"
  19 #include "PS_dot.h"
  20 #include "energy_const.h"
  21 #include "read_epars.h"
  22 #include "LPfold.h"
  23 #include "params.h"
  24 #include "RNAplfold_cmdl.h"
  25
  26 /*@unused@*/
  27 static char rcsid[] = "$Id: RNAplfold.c,v 1.10 2008/07/02 15:34:24 ivo Exp $";
  28
  29 int unpaired;
  30 PRIVATE void putout_pup(double *pup,int length, int winsize, char *name);
  31 PRIVATE void putoutpU_G(double **pU,int length, int ulength, FILE *fp);
  32 PRIVATE void putoutphakim_u(double **pU,int length, int ulength, FILE *fp);
  33
  34 /*--------------------------------------------------------------------------*/
  35 int main(int argc, char *argv[]){
  36   struct        RNAplfold_args_info args_info;
  37   unsigned int  error = 0;
  38   char          fname[FILENAME_MAX_LENGTH], ffname[FILENAME_MAX_LENGTH], *c, *structure, *ParamFile, *ns_bases, *rec_sequence, *rec_id, **rec_rest, *orig_sequence;
  39   unsigned int  input_type;
  40   int           i, length, l, sym, r, istty, winsize, pairdist;
  41   float         cutoff;
  42   int           tempwin, temppair, tempunpaired;
  43   FILE          *pUfp = NULL, *spup = NULL;
  44   double        **pup = NULL; /*prob of being unpaired, lengthwise*/
  45   int           noconv, plexoutput, simply_putout, openenergies, binaries;
  46   plist         *pl, *dpp = NULL;
  47   unsigned int  rec_type, read_opt;
  48   double        betaScale;
  49   pf_paramT     *pf_parameters;
  50   model_detailsT  md;
  51
  52   dangles       = 2;
  53   cutoff        = 0.01;
  54   winsize       = 70;
  55   pairdist      = 0;
  56   unpaired      = 0;
  57   betaScale     = 1.;
  58   simply_putout = plexoutput = openenergies = noconv = 0;binaries=0;
  59   tempwin       = temppair = tempunpaired = 0;
  60   structure     = ParamFile = ns_bases = NULL;
  61   rec_type      = read_opt = 0;
  62   rec_id        = rec_sequence = orig_sequence = NULL;
  63   rec_rest      = NULL;
  64   pf_parameters = NULL;
  65
  66   set_model_details(&md);
  67
  68   /*
  69   #############################################
  70   # check the command line parameters
  71   #############################################
  72   */
  73   if(RNAplfold_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
  74   /* temperature */
  75   if(args_info.temp_given)              temperature = args_info.temp_arg;
  76   /* do not take special tetra loop energies into account */
  77   if(args_info.noTetra_given)           md.special_hp = tetra_loop=0;
  78   /* set dangle model */
  79   if(args_info.dangles_given){
  80     if((args_info.dangles_arg != 0) && (args_info.dangles_arg != 2))
  81       warn_user("required dangle model not implemented, falling back to default dangles=2");
  82     else
  83       md.dangles = dangles = args_info.dangles_arg;
  84   }
  85   /* do not allow weak pairs */
  86   if(args_info.noLP_given)              md.noLP = noLonelyPairs = 1;
  87   /* do not allow wobble pairs (GU) */
  88   if(args_info.noGU_given)              md.noGU = noGU = 1;
  89   /* do not allow weak closing pairs (AU,GU) */
  90   if(args_info.noClosingGU_given)       md.noGUclosure = no_closingGU = 1;
  91   /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
  92   if(args_info.noconv_given)            noconv = 1;
  93   /* set energy model */
  94   if(args_info.energyModel_given)       energy_set = args_info.energyModel_arg;
  95   /* take another energy parameter set */
  96   if(args_info.paramFile_given)         ParamFile = strdup(args_info.paramFile_arg);
  97   /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
  98   if(args_info.nsp_given)               ns_bases = strdup(args_info.nsp_arg);
  99   /* set the maximum base pair span */
 100   if(args_info.span_given)              pairdist = args_info.span_arg;
 101   /* set the pair probability cutoff */
 102   if(args_info.cutoff_given)            cutoff = args_info.cutoff_arg;
 103   /* set the windowsize */
 104   if(args_info.winsize_given)           winsize = args_info.winsize_arg;
 105   /* set the length of unstructured region */
 106   if(args_info.ulength_given)           unpaired = args_info.ulength_arg;
 107   /* compute opening energies */
 108   if(args_info.opening_energies_given)  openenergies = 1;
 109   /* print output on the fly */
 110   if(args_info.print_onthefly_given)    simply_putout = 1;
 111   /* turn on RNAplex output */
 112   if(args_info.plex_output_given)       plexoutput = 1;
 113   /* turn on binary output*/
 114   if(args_info.binaries_given)          binaries = 1;
 115
 116   if(args_info.betaScale_given)         betaScale = args_info.betaScale_arg;
 117
 118   /* check for errorneous parameter options */
 119   if((pairdist < 0) || (cutoff < 0.) || (unpaired < 0) || (winsize < 0)){
 120     RNAplfold_cmdline_parser_print_help();
 121     exit(EXIT_FAILURE);
 122   }
 123
 124   /* free allocated memory of command line data structure */
 125   RNAplfold_cmdline_parser_free(&args_info);
 126
 127   /*
 128   #############################################
 129   # begin initializing
 130   #############################################
 131   */
 132   if (ParamFile != NULL)
 133     read_parameter_file(ParamFile);
 134
 135   if (ns_bases != NULL) {
 136     nonstandards = space(33);
 137     c=ns_bases;
 138     i=sym=0;
 139     if (*c=='-') {
 140       sym=1; c++;
 141     }
 142     while (*c!='\0') {
 143       if (*c!=',') {
 144         nonstandards[i++]=*c++;
 145         nonstandards[i++]=*c;
 146         if ((sym)&&(*c!=*(c-1))) {
 147           nonstandards[i++]=*c;
 148           nonstandards[i++]=*(c-1);
 149         }
 150       }
 151       c++;
 152     }
 153   }
 154
 155   /* check parameter options again and reset to reasonable values if needed */
 156   if(openenergies && !unpaired) unpaired  = 31;
 157   if(pairdist == 0)             pairdist  = winsize;
 158   if(pairdist > winsize){
 159     fprintf(stderr, "pairdist (-L %d) should be <= winsize (-W %d);"
 160             "Setting pairdist=winsize\n",pairdist, winsize);
 161     pairdist = winsize;
 162   }
 163   if(dangles % 2){
 164     warn_user("using default dangles = 2");
 165     dangles = 2;
 166   }
 167
 168   istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
 169   read_opt |= VRNA_INPUT_NO_REST;
 170   if(istty){
 171     print_tty_input_seq();
 172     read_opt |= VRNA_INPUT_NOSKIP_BLANK_LINES;
 173   }
 174
 175   /*
 176   #############################################
 177   # main loop: continue until end of file
 178   #############################################
 179   */
 180   while(
 181     !((rec_type = read_record(&rec_id, &rec_sequence, &rec_rest, read_opt))
 182         & (VRNA_INPUT_ERROR | VRNA_INPUT_QUIT))){
 183
 184     /*
 185     ########################################################
 186     # init everything according to the data we've read
 187     ########################################################
 188     */
 189     if(rec_id){
 190       if(!istty) printf("%s\n", rec_id);
 191       (void) sscanf(rec_id, ">%" XSTR(FILENAME_ID_LENGTH) "s", fname);
 192     }
 193     else fname[0] = '\0';
 194
 195     length    = (int)strlen(rec_sequence);
 196     structure = (char *) space((unsigned) length+1);
 197
 198     /* convert DNA alphabet to RNA if not explicitely switched off */
 199     if(!noconv) str_DNA2RNA(rec_sequence);
 200     /* store case-unmodified sequence */
 201     orig_sequence = strdup(rec_sequence);
 202     /* convert sequence to uppercase letters only */
 203     str_uppercase(rec_sequence);
 204
 205     if(istty) printf("length = %d\n", length);
 206
 207     /*
 208     ########################################################
 209     # done with 'stdin' handling
 210     ########################################################
 211     */
 212
 213     if(length > 1000000){
 214       if(!simply_putout && !unpaired){
 215         printf("Switched to simple output mode!!!\n");
 216         simply_putout = 1;
 217       }
 218     }
 219     if(unpaired && simply_putout){
 220       printf("Output simplification not possible if unpaired is switched on\n");
 221       simply_putout = 0;
 222     }
 223
 224     /* restore winsize if altered before */
 225     if(tempwin != 0){
 226       winsize = tempwin;
 227       tempwin = 0;
 228     }
 229     /* restore pairdist if altered before */
 230     if(temppair != 0){
 231       pairdist = temppair;
 232       temppair = 0;
 233     }
 234     /* restore ulength if altered before */
 235     if(tempunpaired != 0){
 236       unpaired      = tempunpaired;
 237       tempunpaired  = 0;
 238     }
 239
 240     /* adjust winsize, pairdist and ulength if necessary */
 241     if(length < winsize){
 242       fprintf(stderr, "WARN: window size %d larger than sequence length %d\n", winsize, length);
 243       tempwin = winsize;
 244       winsize = length;
 245       if (pairdist>winsize) {
 246         temppair=pairdist;
 247         pairdist=winsize;
 248       }
 249       if (unpaired>winsize) {
 250         tempunpaired=unpaired;
 251         unpaired=winsize;
 252       }
 253     }
 254
 255     /*
 256     ########################################################
 257     # begin actual computations
 258     ########################################################
 259     */
 260
 261     if (length >= 5){
 262       /* construct output file names */
 263       char fname1[FILENAME_MAX_LENGTH], fname2[FILENAME_MAX_LENGTH], fname3[FILENAME_MAX_LENGTH], fname4[FILENAME_MAX_LENGTH], fname_t[FILENAME_MAX_LENGTH];
 264
 265       strcpy(fname_t, (fname[0] != '\0') ? fname : "plfold");
 266
 267       strcpy(fname1, fname_t);
 268       strcpy(fname2, fname_t);
 269       strcpy(fname3, fname_t);
 270       strcpy(fname4, fname_t);
 271       strcpy(ffname, fname_t);
 272
 273       strcat(fname1, "_lunp");
 274       strcat(fname2, "_basepairs");
 275       strcat(fname3, "_uplex");
 276       if(binaries){
 277         strcat(fname4, "_openen_bin");
 278       }
 279       else{
 280         strcat(fname4, "_openen");
 281       }
 282       strcat(ffname, "_dp.ps");
 283
 284       pf_parameters = get_boltzmann_factors(temperature, betaScale, md, -1);
 285
 286       if(unpaired > 0){
 287         pup       =(double **)  space((length+1)*sizeof(double *));
 288         pup[0]    =(double *)   space(sizeof(double)); /*I only need entry 0*/
 289         pup[0][0] = unpaired;
 290       }
 291
 292       pUfp = spup = NULL;
 293       if(simply_putout){
 294         spup = fopen(fname2, "w");
 295         pUfp = (unpaired > 0) ? fopen(fname1, "w") : NULL;
 296
 297         pl = pfl_fold_par(rec_sequence, winsize, pairdist, cutoff, pup, &dpp, pUfp, spup, pf_parameters);
 298
 299         if(pUfp != NULL)  fclose(pUfp);
 300         if(spup != NULL)  fclose(spup);
 301       }
 302       else{
 303         pl = pfl_fold_par(rec_sequence, winsize, pairdist, cutoff, pup, &dpp, pUfp, spup, pf_parameters);
 304         PS_dot_plot_turn(orig_sequence, pl, ffname, pairdist);
 305         if (unpaired > 0){
 306           if(plexoutput){
 307             pUfp = fopen(fname3, "w");
 308             putoutphakim_u(pup,length, unpaired, pUfp);
 309             fclose(pUfp);
 310           }
 311           pUfp = fopen(openenergies ? fname4 : fname1, "w");
 312           if(binaries){
 313             putoutpU_prob_bin_par(pup, length, unpaired, pUfp, openenergies, pf_parameters);
 314           }
 315           else{
 316             putoutpU_prob_par(pup, length, unpaired, pUfp, openenergies, pf_parameters);
 317           }
 318           fclose(pUfp);
 319         }
 320       }
 321       if(pl) free(pl);
 322       if(unpaired > 0){
 323         free(pup[0]);
 324         free(pup);
 325       }
 326
 327       free(pf_parameters);
 328     }
 329     (void) fflush(stdout);
 330
 331     /* clean up */
 332     if(rec_id) free(rec_id);
 333     free(rec_sequence);
 334     free(orig_sequence);
 335     free(structure);
 336     rec_id = rec_sequence = orig_sequence = NULL;
 337     rec_rest = NULL;
 338     /* print user help for the next round if we get input from tty */
 339
 340     if(istty) print_tty_input_seq();
 341   }
 342   return EXIT_SUCCESS;
 343 }
 344
 345 /* additional output functions */
 346
 347 PRIVATE void putout_pup(double *pup,int length, int winsize, char *name) {
 348   int i;
 349   float factor;
 350   float tfact;
 351   FILE *FP;
 352
 353   FP=fopen(name,"w");
 354   fprintf(FP,"&prob of being unpaired between i-%d and i\n",unpaired);
 355   fflush(NULL);
 356   for (i=unpaired; i<=length; i++) {
 357     factor=0.;
 358     if (i<winsize) {
 359       factor=1./(i-unpaired+1);
 360     }
 361     if (i>length-winsize+unpaired-1) {
 362       tfact=1./(length-i+1);
 363       if (tfact>factor) {
 364         factor=tfact;
 365       }
 366
 367     }
 368     else {
 369       tfact=1./(winsize-unpaired+1);
 370       if (tfact>factor) {
 371         factor=tfact;
 372       }
 373     }
 374     fprintf(FP,"%d %.6f\n",i,pup[i]*factor);
 375   }
 376   fclose(FP);
 377
 378 }
 379 PRIVATE void putoutpU_G(double **pU,int length, int ulength, FILE *fp) {
 380   /*put out unpaireds */
 381   int i,k;
 382   fprintf(fp,"#unpaired probabilities\n #i\tl=");
 383   for (i=1; i<=ulength; i++) {
 384     fprintf(fp,"%d\t", i);
 385   }
 386   fprintf(fp,"\n");
 387   for (k=1; k<=length; k++){
 388     fprintf(fp,"%d\t",k);
 389     for (i=1; i<=ulength; i++) {
 390       if (k+(i-1)<=length) fprintf(fp,"%.7g\t",pU[k+(i-1)][i]);
 391     }
 392     fprintf(fp,"\n");
 393     free(pU[k]);
 394   }
 395   free(pU[0]);
 396   free(pU);
 397   fflush(fp);
 398 }
 399 PRIVATE void putoutphakim_u(double **pU,int length, int ulength, FILE *fp) {
 400   /*put out Fopen in dekacalories per mol, and F(cond,open) also in dekacal*/
 401   int k;
 402
 403   float RT = (temperature+K0)*GASCONST;
 404   float p0;
 405   float pdep;
 406   int f0;
 407   int fdep;
 408
 409   fprintf(fp,"#energy necessary to unpair as well as to unpair if i-1 is unpaired also, if i+1 is unpaired also in dekacal/mol\n");
 410   for (k=1; k<=length; k++){
 411     fprintf(fp,"%d\t",k);
 412     p0=pU[k][1];
 413     f0=(int) -RT*log(p0)/10;
 414     fprintf(fp,"%d\t", f0);
 415     if (k>1) {
 416       pdep=pU[k][2]/pU[k-1][1];
 417       fdep=(int) -RT*log(pdep)/10;
 418       fprintf(fp,"%d\t",fdep);
 419
 420
 421     }
 422     else  fprintf(fp,"-0\t");
 423     if (k<length) {
 424       pdep=pU[k+1][2]/pU[k+1][1];
 425       fdep=(int) -RT*log(pdep)/10;
 426       fprintf(fp,"%d\t",fdep);
 427     }
 428     else  fprintf(fp,"-0\t");
 429     fprintf(fp,"\n");
 430   }
 431
 432   fflush(fp);
 433 }
 434