binaries/src/ViennaRNA/Progs/RNAheat.c

   1 /*
   2                       Heat Capacity of RNA molecule
   3
   4                     c Ivo Hofacker and Peter Stadler
   5                           Vienna RNA package
   6
   7
   8             calculates specific heat using C = - T d^2/dT^2 G(T)
   9 */
  10
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <ctype.h>
  15 #include <math.h>
  16 #include <unistd.h>
  17 #include "utils.h"
  18 #include "fold_vars.h"
  19 #include "fold.h"
  20 #include "params.h"
  21 #include "part_func.h"
  22 #include "read_epars.h"
  23 #include "RNAheat_cmdl.h"
  24
  25 /*@unused@*/
  26 static char rcsid[] = "$Id: RNAheat.c,v 1.12 2000/09/28 11:23:14 ivo Rel $";
  27
  28 #define MAXWIDTH  201
  29 #define GASCONST  1.98717  /* in [cal/K] */
  30 #define K0        273.15
  31
  32 PRIVATE float F[MAXWIDTH];
  33 PRIVATE float ddiff(float f[], float h, int m);
  34 PRIVATE void heat_capacity(char *string, float T_min, float T_max, float h, int m);
  35
  36 int main(int argc, char *argv[]){
  37   struct  RNAheat_args_info args_info;
  38   char                      *string, *input_string, *ns_bases, *c, *ParamFile, *rec_sequence, *rec_id, **rec_rest, *orig_sequence;
  39   int                       i, length, l, sym;
  40   float                     T_min, T_max, h;
  41   int                       mpoints, istty, noconv = 0;
  42   unsigned int              input_type;
  43   unsigned int              rec_type, read_opt;
  44
  45   string    = ParamFile = ns_bases = NULL;
  46   T_min     = 0.;
  47   T_max     = 100.;
  48   h         = 1;
  49   mpoints   = 2;
  50   dangles   = 2;   /* dangles can be 0 (no dangles) or 2, default is 2 */
  51   rec_type  = read_opt = 0;
  52   rec_id    = rec_sequence = orig_sequence = NULL;
  53   rec_rest  = NULL;
  54
  55   /*
  56   #############################################
  57   # check the command line parameters
  58   #############################################
  59   */
  60   if(RNAheat_cmdline_parser(argc, argv, &args_info) != 0) exit(1);
  61
  62   /* do not take special tetra loop energies into account */
  63   if(args_info.noTetra_given)     tetra_loop=0;
  64   /* set dangle model */
  65   if(args_info.dangles_given){
  66     if((args_info.dangles_arg != 0) && (args_info.dangles_arg != 2))
  67       warn_user("required dangle model not implemented, falling back to default dangles=2");
  68     else
  69       dangles = args_info.dangles_arg;
  70   }
  71   /* do not allow weak pairs */
  72   if(args_info.noLP_given)        noLonelyPairs = 1;
  73   /* do not allow wobble pairs (GU) */
  74   if(args_info.noGU_given)        noGU = 1;
  75   /* do not allow weak closing pairs (AU,GU) */
  76   if(args_info.noClosingGU_given) no_closingGU = 1;
  77   /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
  78   if(args_info.noconv_given)      noconv = 1;
  79   /* set energy model */
  80   if(args_info.energyModel_given) energy_set = args_info.energyModel_arg;
  81   /* take another energy parameter set */
  82   if(args_info.paramFile_given)   ParamFile = strdup(args_info.paramFile_arg);
  83   /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
  84   if(args_info.nsp_given)         ns_bases = strdup(args_info.nsp_arg);
  85   /* Tmin */
  86   if(args_info.Tmin_given)        T_min = args_info.Tmin_arg;
  87   /* Tmax */
  88   if(args_info.Tmax_given)        T_max = args_info.Tmax_arg;
  89   /* step size */
  90   if(args_info.stepsize_given)    h = args_info.stepsize_arg;
  91   /* ipoints */
  92   if(args_info.ipoints_given){
  93     mpoints = args_info.ipoints_arg;
  94     if (mpoints < 1)    mpoints = 1;
  95     if (mpoints > 100)  mpoints = 100;
  96   }
  97
  98   /* free allocated memory of command line data structure */
  99   RNAheat_cmdline_parser_free (&args_info);
 100
 101   /*
 102   #############################################
 103   # begin initializing
 104   #############################################
 105   */
 106   if (ParamFile!=NULL) read_parameter_file(ParamFile);
 107
 108   if (ns_bases != NULL) {
 109     nonstandards = space(33);
 110     c=ns_bases;
 111     i=sym=0;
 112     if (*c=='-') {
 113       sym=1; c++;
 114     }
 115     while (*c!='\0') {
 116       if (*c!=',') {
 117         nonstandards[i++]=*c++;
 118         nonstandards[i++]=*c;
 119         if ((sym)&&(*c!=*(c-1))) {
 120           nonstandards[i++]=*c;
 121           nonstandards[i++]=*(c-1);
 122         }
 123       }
 124       c++;
 125     }
 126   }
 127
 128   istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
 129
 130   read_opt |= VRNA_INPUT_NO_REST;
 131   if(istty){
 132     print_tty_input_seq();
 133     read_opt |= VRNA_INPUT_NOSKIP_BLANK_LINES;
 134   }
 135
 136   /*
 137   #############################################
 138   # main loop: continue until end of file
 139   #############################################
 140   */
 141   while(
 142     !((rec_type = read_record(&rec_id, &rec_sequence, &rec_rest, read_opt))
 143         & (VRNA_INPUT_ERROR | VRNA_INPUT_QUIT))){
 144     /*
 145     ########################################################
 146     # init everything according to the data we've read
 147     ########################################################
 148     */
 149     if(rec_id && !istty) printf("%s\n", rec_id);
 150
 151     length = (int)strlen(rec_sequence);
 152
 153     /* convert DNA alphabet to RNA if not explicitely switched off */
 154     if(!noconv) str_DNA2RNA(rec_sequence);
 155     /* store case-unmodified sequence */
 156     orig_sequence = strdup(rec_sequence);
 157     /* convert sequence to uppercase letters only */
 158     str_uppercase(rec_sequence);
 159
 160     if(istty) printf("length = %d\n", length);
 161     /*
 162     ########################################################
 163     # done with 'stdin' handling
 164     ########################################################
 165     */
 166
 167     heat_capacity(rec_sequence, T_min, T_max, h, mpoints);
 168     (void) fflush(stdout);
 169
 170     /* clean up */
 171     if(rec_id) free(rec_id);
 172     free(rec_sequence);
 173     free(orig_sequence);
 174     rec_id = rec_sequence = orig_sequence = NULL;
 175     rec_rest = NULL;
 176     /* print user help for the next round if we get input from tty */
 177
 178     if(istty) print_tty_input_seq();
 179   }
 180   return EXIT_SUCCESS;
 181 }
 182
 183 PRIVATE void heat_capacity(char *string, float T_min, float T_max,
 184                           float h, int m)
 185 {
 186    int length, i;
 187    char *structure;
 188    float hc, kT, min_en;
 189
 190    length = (int) strlen(string);
 191
 192    do_backtrack = 0;
 193
 194    temperature = T_min -m*h;
 195    /* initialize_fold(length); <- obsolete */
 196    structure = (char *) space((unsigned) length+1);
 197    min_en = fold(string, structure);
 198    free(structure); free_arrays();
 199    kT = (temperature+K0)*GASCONST/1000;    /* in kcal */
 200    pf_scale = exp(-(1.07*min_en)/kT/length );
 201    /* init_pf_fold(length); <- obsolete */
 202     pf_paramT       *pf_parameters = NULL;
 203     model_detailsT  md;
 204     set_model_details(&md);
 205
 206    for (i=0; i<2*m+1; i++) {
 207       if(pf_parameters) free(pf_parameters);
 208       pf_parameters = get_boltzmann_factors(temperature, 1.0, md, pf_scale);
 209       F[i] = pf_fold_par(string, NULL, pf_parameters, 0, 0, 0);   /* T_min -2h */
 210       temperature += h;
 211       kT = (temperature+K0)*GASCONST/1000;
 212       pf_scale=exp(-(F[i]/length +h*0.00727)/kT); /* try to extrapolate F */
 213       update_pf_params(length);
 214    }
 215    while (temperature <= (T_max+m*h+h)) {
 216
 217       hc = - ddiff(F,h,m)* (temperature +K0 - m*h -h);
 218       printf("%g   %g\n", (temperature-m*h-h), hc);
 219
 220       for (i=0; i<2*m; i++)
 221          F[i] = F[i+1];
 222
 223       F[2*m] = pf_fold_par(string, NULL, pf_parameters, 0, 0, 0);
 224 /*       printf("%g\n", F[2*m]);*/
 225       temperature += h;
 226       kT = (temperature+K0)*GASCONST/1000;
 227       pf_scale=exp(-(F[i]/length +h*0.00727)/kT);
 228       if(pf_parameters) free(pf_parameters);
 229       pf_parameters = get_boltzmann_factors(temperature, 1.0, md, pf_scale);
 230       update_pf_params(length);
 231    }
 232    free_pf_arrays();
 233 }
 234
 235 /* ------------------------------------------------------------------------- */
 236
 237 PRIVATE float ddiff(float f[], float h, int m)
 238 {
 239    float fp;
 240    int i;
 241    float A, B;
 242    A = (float)(m*(m+1)*(2*m+1)/3 );                            /* 2*sum(x^2) */
 243    B = (float)(m*(m+1)*(2*m+1) ) * (float)(3*m*m+3*m-1) /15.;  /* 2*sum(x^4) */
 244
 245    fp=0.;
 246    for (i=0; i<2*m+1; i++)
 247       fp += f[i]*( A - (float) ( (2*m+1)*(i-m)*(i-m)) );
 248
 249    fp /= ( ( A*A - B*( (float)(2*m+1) ) )*h*h/2. );
 250    return (float)fp;
 251
 252 }
 253