binaries/src/ViennaRNA/Progs/RNAeval.c

   1 /* Last changed Time-stamp: <2008-09-02 10:47:24 ivo> */
   2 /*
   3
   4           Calculate Energy of given Sequences and Structures
   5                            c Ivo L Hofacker
   6                           Vienna RNA Pckage
   7 */
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <math.h>
  11 #include <ctype.h>
  12 #include <unistd.h>
  13 #include <string.h>
  14 #include <sys/types.h>
  15 #include "fold_vars.h"
  16 #include "fold.h"
  17 #include "utils.h"
  18 #include "read_epars.h"
  19 #include "RNAeval_cmdl.h"
  20
  21 #ifdef __GNUC__
  22 #define UNUSED __attribute__ ((unused))
  23 #else
  24 #define UNUSED
  25 #endif
  26
  27 /*@unused@*/
  28 static char UNUSED rcsid[]="$Id: RNAeval.c,v 1.10 2008/10/09 07:08:40 ivo Exp $";
  29
  30 PRIVATE char *costring(char *string);
  31 PRIVATE char *tokenize(char *line);
  32
  33 int main(int argc, char *argv[]){
  34   struct RNAeval_args_info  args_info;
  35   char                      *string, *structure, *orig_sequence, *tmp;
  36   char                      *rec_sequence, *rec_id, **rec_rest;
  37   char                      fname[FILENAME_MAX_LENGTH];
  38   char                      *ParamFile;
  39   int                       i, length1, length2;
  40   float                     energy;
  41   int                       istty;
  42   int                       circular=0;
  43   int                       noconv=0;
  44   int                       verbose = 0;
  45   unsigned int              rec_type, read_opt;
  46
  47   string  = orig_sequence = ParamFile = NULL;
  48   gquad   = 0;
  49   dangles = 2;
  50
  51   /*
  52   #############################################
  53   # check the command line parameters
  54   #############################################
  55   */
  56   if(RNAeval_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
  57   /* temperature */
  58   if(args_info.temp_given)        temperature = args_info.temp_arg;
  59   /* do not take special tetra loop energies into account */
  60   if(args_info.noTetra_given)     tetra_loop=0;
  61   /* set dangle model */
  62   if(args_info.dangles_given){
  63     if((args_info.dangles_arg < 0) || (args_info.dangles_arg > 3))
  64       warn_user("required dangle model not implemented, falling back to default dangles=2");
  65     else
  66       dangles = args_info.dangles_arg;
  67   }
  68   /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
  69   if(args_info.noconv_given)      noconv = 1;
  70   /* set energy model */
  71   if(args_info.energyModel_given) energy_set = args_info.energyModel_arg;
  72   /* take another energy parameter set */
  73   if(args_info.paramFile_given)   ParamFile = strdup(args_info.paramFile_arg);
  74   /* assume RNA sequence to be circular */
  75   if(args_info.circ_given)        circular=1;
  76   /* logarithmic multiloop energies */
  77   if(args_info.logML_given)       logML = 1;
  78   /* be verbose */
  79   if(args_info.verbose_given)     verbose = 1;
  80   /* gquadruplex support */
  81   if(args_info.gquad_given)       gquad = 1;
  82
  83   /* free allocated memory of command line data structure */
  84   RNAeval_cmdline_parser_free (&args_info);
  85
  86   /*
  87   #############################################
  88   # begin initializing
  89   #############################################
  90   */
  91
  92   if (ParamFile!=NULL) read_parameter_file(ParamFile);
  93
  94   rec_type      = read_opt = 0;
  95   rec_id        = rec_sequence = NULL;
  96   rec_rest      = NULL;
  97   istty         = isatty(fileno(stdout)) && isatty(fileno(stdin));
  98
  99   if(circular && gquad){
 100     nrerror("G-Quadruplex support is currently not available for circular RNA structures");
 101   }
 102
 103   /* set options we wanna pass to read_record */
 104   if(istty){
 105     read_opt |= VRNA_INPUT_NOSKIP_BLANK_LINES;
 106     print_tty_input_seq_str("Use '&' to connect 2 sequences that shall form a complex.\n"
 107                             "Input sequence (upper or lower case) followed by structure");
 108   }
 109
 110   /*
 111   #############################################
 112   # main loop: continue until end of file
 113   #############################################
 114   */
 115   while(
 116     !((rec_type = read_record(&rec_id, &rec_sequence, &rec_rest, read_opt))
 117         & (VRNA_INPUT_ERROR | VRNA_INPUT_QUIT))){
 118
 119     if(rec_id){
 120       (void) sscanf(rec_id, ">%" XSTR(FILENAME_ID_LENGTH) "s", fname);
 121     }
 122     else fname[0] = '\0';
 123
 124     cut_point = -1;
 125
 126     string    = tokenize(rec_sequence);
 127     length2   = (int) strlen(string);
 128     tmp       = extract_record_rest_structure((const char **)rec_rest, 0, (rec_id) ? VRNA_OPTION_MULTILINE : 0);
 129
 130     if(!tmp)
 131       nrerror("structure missing");
 132
 133     structure = tokenize(tmp);
 134     length1   = (int) strlen(structure);
 135     if(length1 != length2)
 136       nrerror("structure and sequence differ in length!");
 137
 138     free(tmp);
 139
 140     /* convert DNA alphabet to RNA if not explicitely switched off */
 141     if(!noconv) str_DNA2RNA(string);
 142     /* store case-unmodified sequence */
 143     orig_sequence = strdup(string);
 144     /* convert sequence to uppercase letters only */
 145     str_uppercase(string);
 146
 147     if(istty){
 148       if (cut_point == -1)
 149         printf("length = %d\n", length1);
 150       else
 151         printf("length1 = %d\nlength2 = %d\n", cut_point-1, length1-cut_point+1);
 152     }
 153
 154     if(gquad)
 155       energy = energy_of_gquad_structure(string, structure, verbose);
 156     else
 157       energy = (circular) ? energy_of_circ_structure(string, structure, verbose) : energy_of_structure(string, structure, verbose);
 158
 159     if (cut_point == -1)
 160       printf("%s\n%s", orig_sequence, structure);
 161     else {
 162       char *pstring, *pstruct;
 163       pstring = costring(orig_sequence);
 164       pstruct = costring(structure);
 165       printf("%s\n%s", pstring,  pstruct);
 166       free(pstring);
 167       free(pstruct);
 168     }
 169     if (istty)
 170       printf("\n energy = %6.2f\n", energy);
 171     else
 172       printf(" (%6.2f)\n", energy);
 173
 174     /* clean up */
 175     (void) fflush(stdout);
 176     if(rec_id) free(rec_id);
 177     free(rec_sequence);
 178     free(structure);
 179     /* free the rest of current dataset */
 180     if(rec_rest){
 181       for(i=0;rec_rest[i];i++) free(rec_rest[i]);
 182       free(rec_rest);
 183     }
 184     rec_id = rec_sequence = structure = NULL;
 185     rec_rest = NULL;
 186
 187     free(string);
 188     free(orig_sequence);
 189     string = orig_sequence = NULL;
 190
 191     /* print user help for the next round if we get input from tty */
 192     if(istty){
 193       print_tty_input_seq_str("Use '&' to connect 2 sequences that shall form a complex.\n"
 194                               "Input sequence (upper or lower case) followed by structure");
 195     }
 196   }
 197   return EXIT_SUCCESS;
 198 }
 199
 200 PRIVATE char *tokenize(char *line)
 201 {
 202   char *token, *copy, *ctmp;
 203   int cut = -1;
 204
 205   copy = (char *) space(strlen(line)+1);
 206   ctmp = (char *) space(strlen(line)+1);
 207   (void) sscanf(line, "%s", copy);
 208   ctmp[0] = '\0';
 209   token = strtok(copy, "&");
 210   cut = strlen(token)+1;
 211   while (token) {
 212     strcat(ctmp, token);
 213     token = strtok(NULL, "&");
 214   }
 215   if (cut > strlen(ctmp)) cut = -1;
 216   if (cut > -1) {
 217     if (cut_point==-1) cut_point = cut;
 218     else if (cut_point != cut) {
 219       fprintf(stderr,"cut_point = %d cut = %d\n", cut_point, cut);
 220       nrerror("Sequence and Structure have different cut points.");
 221     }
 222   }
 223   free(copy);
 224
 225   return ctmp;
 226 }
 227
 228 PRIVATE char *costring(char *string)
 229 {
 230   char *ctmp;
 231   int len;
 232
 233   len = strlen(string);
 234   ctmp = (char *)space((len+2) * sizeof(char));
 235   /* first sequence */
 236   (void) strncpy(ctmp, string, cut_point-1);
 237   /* spacer */
 238   ctmp[cut_point-1] = '&';
 239   /* second sequence */
 240   (void) strcat(ctmp, string+cut_point-1);
 241
 242   return ctmp;
 243 }