binaries/src/ViennaRNA/Progs/RNA2Dfold.c

   1
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <math.h>
   5 #include <ctype.h>
   6 #include <unistd.h>
   7 #include <string.h>
   8 #include "fold.h"
   9 #include "part_func.h"
  10 #include "fold_vars.h"
  11 #include "PS_dot.h"
  12 #include "utils.h"
  13 #include "energy_par.h"
  14 #include "read_epars.h"
  15 #include "2Dfold.h"
  16 #include "2Dpfold.h"
  17 #include "RNA2Dfold_cmdl.h"
  18
  19 #ifdef _OPENMP
  20 #include <omp.h>
  21 #endif
  22
  23 typedef struct nbhoods{
  24   int k;
  25   int l;
  26   struct nbhoods *next;
  27 } nbhoods;
  28
  29 int main(int argc, char *argv[]){
  30   struct        RNA2Dfold_args_info args_info;
  31   unsigned int  input_type;
  32   char *string, *input_string, *orig_sequence;
  33   char *mfe_structure=NULL, *structure1=NULL, *structure2=NULL, *reference_struc1=NULL, *reference_struc2=NULL;
  34   char  *ParamFile=NULL;
  35   int   i, j, length, l;
  36   double min_en;
  37   double kT, sfact=1.07;
  38   int   pf=0,istty;
  39   int noconv=0;
  40   int circ=0;
  41   int maxDistance1 = -1;
  42   int maxDistance2 = -1;
  43   int do_backtrack = 1;
  44   int stBT = 0;
  45   int nstBT = 0;
  46   string=NULL;
  47   dangles = 2;
  48   struct nbhoods *neighborhoods = NULL;
  49   struct nbhoods *neighborhoods_cur = NULL;
  50
  51   string = input_string = orig_sequence = NULL;
  52   /*
  53   #############################################
  54   # check the command line prameters
  55   #############################################
  56   */
  57   if(RNA2Dfold_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
  58
  59   /* temperature */
  60   if(args_info.temp_given)
  61     temperature = args_info.temp_arg;
  62
  63   /* max distance to 1st reference structure */
  64   if(args_info.maxDist1_given)
  65     maxDistance1 = args_info.maxDist1_arg;
  66
  67   /* max distance to 2nd reference structure */
  68   if(args_info.maxDist2_given)
  69     maxDistance2 = args_info.maxDist2_arg;
  70
  71   /* compute partition function and boltzmann probabilities */
  72   if(args_info.partfunc_given)
  73     pf = 1;
  74
  75   /* do stachastic backtracking */
  76   if(args_info.stochBT_given){
  77     pf = 1;
  78     stBT = 1;
  79     nstBT = args_info.stochBT_arg;
  80   }
  81
  82   if(args_info.noTetra_given)
  83     tetra_loop=0;
  84
  85   /* assume RNA sequence to be circular */
  86   if(args_info.circ_given)
  87     circ=1;
  88
  89   /* dangle options */
  90   if(args_info.dangles_given){
  91     if((args_info.dangles_arg != 0) && (args_info.dangles_arg != 2))
  92       warn_user("required dangle model not implemented, falling back to default dangles=2");
  93     else
  94       dangles = args_info.dangles_arg;
  95   }
  96   /* set number of threads for parallel computation */
  97   if(args_info.numThreads_given)
  98 #ifdef _OPENMP
  99   omp_set_num_threads(args_info.numThreads_arg);
 100 #else
 101   nrerror("\'j\' option is available only if compiled with OpenMP support!");
 102 #endif
 103
 104   /* get energy parameter file name */
 105   if(args_info.parameterFile_given)
 106     ParamFile = strdup(args_info.parameterFile_arg);
 107
 108   /* do not allow GU pairs ? */
 109   if(args_info.noGU_given)
 110     noGU = 1;
 111
 112   /* do not allow GU pairs at the end of helices? */
 113   if(args_info.noClosingGU_given)
 114     no_closingGU = 1;
 115
 116   /* pf scaling factor */
 117   if(args_info.pfScale_given)
 118     sfact = args_info.pfScale_arg;
 119
 120   /* do not backtrack structures ? */
 121   if(args_info.noBT_given)
 122     do_backtrack = 0;
 123
 124   for (i = 0; i < args_info.neighborhood_given; i++){
 125     int kappa, lambda;
 126     kappa = lambda = 0;
 127     if(sscanf(args_info.neighborhood_arg[i], "%d:%d", &kappa, &lambda) == 2);
 128     if ((kappa>-2) && (lambda>-2)){
 129       if(neighborhoods_cur != NULL){
 130         neighborhoods_cur->next = (nbhoods *)space(sizeof(nbhoods));
 131         neighborhoods_cur = neighborhoods_cur->next;
 132       }
 133       else{
 134         neighborhoods = (nbhoods *)space(sizeof(nbhoods));
 135         neighborhoods_cur = neighborhoods;
 136       }
 137       neighborhoods_cur->k = kappa;
 138       neighborhoods_cur->l = lambda;
 139       neighborhoods_cur->next = NULL;
 140     }
 141   }
 142   /* free allocated memory of command line data structure */
 143   RNA2Dfold_cmdline_parser_free (&args_info);
 144
 145   /*
 146   #############################################
 147   # begin actual program code
 148   #############################################
 149   */
 150   if (ParamFile != NULL)
 151     read_parameter_file(ParamFile);
 152
 153   istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
 154
 155  /*
 156   #############################################
 157   # main loop, continue until end of file
 158   #############################################
 159   */
 160   do {
 161     if (istty)
 162       print_tty_input_seq_str("Input strings\n1st line: sequence (upper or lower case)\n2nd + 3rd line: reference structures (dot bracket notation)\n@ to quit\n");
 163
 164     while((input_type = get_input_line(&input_string, 0)) & VRNA_INPUT_FASTA_HEADER){
 165       printf(">%s\n", input_string); /* print fasta header if available */
 166       free(input_string);
 167     }
 168
 169     /* break on any error, EOF or quit request */
 170     if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)){ break;}
 171     /* else assume a proper sequence of letters of a certain alphabet (RNA, DNA, etc.) */
 172     else{
 173       length = (int)    strlen(input_string);
 174       string = strdup(input_string);
 175       free(input_string);
 176     }
 177
 178     mfe_structure = (char *) space((unsigned) length+1);
 179     structure1 = (char *) space((unsigned) length+1);
 180     structure2 = (char *) space((unsigned) length+1);
 181
 182     input_type = get_input_line(&input_string, VRNA_INPUT_NOSKIP_COMMENTS);
 183     if(input_type & VRNA_INPUT_QUIT){ break;}
 184     else if((input_type & VRNA_INPUT_MISC) && (strlen(input_string) > 0)){
 185       reference_struc1 = strdup(input_string);
 186       free(input_string);
 187       if(strlen(reference_struc1) != length)
 188         nrerror("sequence and 1st reference structure have unequal length");
 189     }
 190     else nrerror("1st reference structure missing\n");
 191     strncpy(structure1, reference_struc1, length);
 192
 193     input_type = get_input_line(&input_string, VRNA_INPUT_NOSKIP_COMMENTS);
 194     if(input_type & VRNA_INPUT_QUIT){ break;}
 195     else if((input_type & VRNA_INPUT_MISC) && (strlen(input_string) > 0)){
 196       reference_struc2 = strdup(input_string);
 197       free(input_string);
 198       if(strlen(reference_struc2) != length)
 199         nrerror("sequence and 2nd reference structure have unequal length");
 200     }
 201     else nrerror("2nd reference structure missing\n");
 202     strncpy(structure2, reference_struc2, length);
 203
 204     /* convert DNA alphabet to RNA if not explicitely switched off */
 205     if(!noconv) str_DNA2RNA(string);
 206     /* store case-unmodified sequence */
 207     orig_sequence = strdup(string);
 208     /* convert sequence to uppercase letters only */
 209     str_uppercase(string);
 210
 211     if (istty)  printf("length = %d\n", length);
 212
 213     min_en = (circ) ? circfold(string, mfe_structure) : fold(string, mfe_structure);
 214
 215     printf("%s\n%s", orig_sequence, mfe_structure);
 216
 217     if (istty)
 218       printf("\n minimum free energy = %6.2f kcal/mol\n", min_en);
 219     else
 220       printf(" (%6.2f)\n", min_en);
 221
 222     printf("%s (%6.2f) <ref 1>\n", structure1, (circ) ? energy_of_circ_structure(string, structure1, 0) : energy_of_structure(string,structure1, 0));
 223     printf("%s (%6.2f) <ref 2>\n", structure2, (circ) ? energy_of_circ_structure(string, structure2, 0) : energy_of_structure(string,structure2, 0));
 224
 225     /* get all variables need for the folding process (some memory will be preallocated here too) */
 226     TwoDfold_vars *mfe_vars = get_TwoDfold_variables(string, structure1, structure2, circ);
 227     mfe_vars->do_backtrack = do_backtrack;
 228     TwoDfold_solution *mfe_s = TwoDfoldList(mfe_vars, maxDistance1, maxDistance2);
 229
 230     if(!pf){
 231 #ifdef COUNT_STATES
 232       printf("k\tl\tn\tMFE\tMFE-structure\n");
 233       for(i = 0; mfe_s[i].k != INF; i++){
 234         printf("%d\t%d\t%lu\t%6.2f\t%s\n", mfe_s[i].k, mfe_s[i].l, mfe_vars->N_F5[length][mfe_s[i].k][mfe_s[i].l/2], mfe_s[i].en, mfe_s[i].s);
 235         if(mfe_s[i].s) free(mfe_s[i].s);
 236       }
 237       free(mfe_s);
 238 #else
 239       printf("k\tl\tMFE\tMFE-structure\n");
 240       for(i = 0; mfe_s[i].k != INF; i++){
 241         printf("%d\t%d\t%6.2f\t%s\n", mfe_s[i].k, mfe_s[i].l, mfe_s[i].en, mfe_s[i].s);
 242         if(mfe_s[i].s) free(mfe_s[i].s);
 243       }
 244       free(mfe_s);
 245 #endif
 246     }
 247
 248     if(pf){
 249       int maxD1 = (int) mfe_vars->maxD1;
 250       int maxD2 = (int) mfe_vars->maxD2;
 251       float mmfe = INF;
 252       double Q;
 253       for(i = 0; mfe_s[i].k != INF; i++){
 254         if(mmfe > mfe_s[i].en) mmfe = mfe_s[i].en;
 255       }
 256       kT = (temperature+K0)*GASCONST/1000.0; /* in Kcal */
 257       pf_scale = exp(-(sfact*mmfe)/kT/length);
 258       if (length>2000)
 259         fprintf(stdout, "scaling factor %f\n", pf_scale);
 260
 261       /* get all variables need for the folding process (some memory will be preallocated there too) */
 262       /* TwoDpfold_vars *q_vars = get_TwoDpfold_variables_from_MFE(mfe_vars); */
 263       /* we dont need the mfe vars and arrays anymore, so we can savely free their occupying memory */
 264       destroy_TwoDfold_variables(mfe_vars);
 265       TwoDpfold_vars *q_vars = get_TwoDpfold_variables(string, structure1, structure2, circ);
 266
 267       TwoDpfold_solution *pf_s = TwoDpfoldList(q_vars, maxD1, maxD2);
 268
 269       Q = 0.;
 270
 271       for(i = 0; pf_s[i].k != INF; i++){
 272         Q += pf_s[i].q;
 273       }
 274
 275       double fee = (-log(Q)-length*log(pf_scale))*kT;
 276
 277       if(!stBT){
 278         printf("free energy of ensemble = %6.2f kcal/mol\n",fee);
 279         printf("k\tl\tP(neighborhood)\tP(MFE in neighborhood)\tP(MFE in ensemble)\tMFE\tE_gibbs\tMFE-structure\n");
 280         for(i=0; pf_s[i].k != INF;i++){
 281           float free_energy = (-log((float)pf_s[i].q)-length*log(pf_scale))*kT;
 282           if((pf_s[i].k != mfe_s[i].k) || (pf_s[i].l != mfe_s[i].l))
 283             nrerror("This should never happen!");
 284           fprintf(stdout,
 285                   "%d\t%d\t%2.8f\t%2.8f\t%2.8f\t%6.2f\t%6.2f\t%s\n",
 286                   pf_s[i].k,
 287                   pf_s[i].l,
 288                   (float)(pf_s[i].q)/(float)Q,
 289                   exp((free_energy-mfe_s[i].en)/kT),
 290                   exp((fee-mfe_s[i].en)/kT),
 291                   mfe_s[i].en,
 292                   free_energy,
 293                   mfe_s[i].s);
 294         }
 295       }
 296       else{
 297         init_rand();
 298         printf("k\tl\ten\tstructure\n");
 299         if(neighborhoods != NULL){
 300           nbhoods *tmp, *tmp2;
 301           for(tmp = neighborhoods; tmp != NULL; tmp = tmp->next){
 302             int k,l;
 303             k = tmp->k;
 304             l = tmp->l;
 305             for(i = 0; i < nstBT; i++){
 306               char *s = TwoDpfold_pbacktrack(q_vars, k, l);
 307               printf("%d\t%d\t%6.2f\t%s\n", k, l, q_vars->circ ? energy_of_circ_structure(q_vars->sequence, s, 0) : energy_of_structure(q_vars->sequence, s, 0), s);
 308             }
 309           }
 310         }
 311         else{
 312           for(i=0; pf_s[i].k != INF;i++){
 313             for(l = 0; l < nstBT; l++){
 314               char *s = TwoDpfold_pbacktrack(q_vars, pf_s[i].k, pf_s[i].l);
 315               printf("%d\t%d\t%6.2f\t%s\n", pf_s[i].k, pf_s[i].l, q_vars->circ ? energy_of_circ_structure(q_vars->sequence, s, 0) : energy_of_structure(q_vars->sequence, s, 0), s);
 316             }
 317           }
 318         }
 319       }
 320       free_pf_arrays();
 321
 322       for(i=0; mfe_s[i].k != INF;i++){
 323         if(mfe_s[i].s) free(mfe_s[i].s);
 324       }
 325       free(pf_s);
 326       free(mfe_s);
 327       /* destroy the q_vars */
 328       destroy_TwoDpfold_variables(q_vars);
 329     }
 330     else
 331       destroy_TwoDfold_variables(mfe_vars);
 332
 333     free_arrays();
 334     free(string);
 335     free(orig_sequence);
 336     free(mfe_structure);
 337     free(structure1);
 338     free(structure2);
 339     free(reference_struc1);
 340     free(reference_struc2);
 341     string = orig_sequence = mfe_structure = NULL;
 342   } while (1);
 343   return 0;
 344 }