binaries/src/ViennaRNA/Progs/RNAduplex.c

   1 /* Last changed Time-stamp: <2005-07-23 16:50:24 ivo> */
   2 /*
   3              Compute duplex structure of two RNA strands
   4
   5                            c Ivo L Hofacker
   6                           Vienna RNA package
   7 */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <ctype.h>
  12 #include <unistd.h>
  13 #include <string.h>
  14 #include "fold.h"
  15 #include "fold_vars.h"
  16 #include "utils.h"
  17 #include "read_epars.h"
  18 #include "subopt.h"
  19 #include "duplex.h"
  20 #include "RNAduplex_cmdl.h"
  21
  22 PRIVATE void  print_struc(duplexT const *dup);
  23 /*@unused@*/
  24 static char rcsid[] = "$Id: RNAduplex.c,v 1.5 2007/08/26 09:41:12 ivo Exp $";
  25
  26 /*--------------------------------------------------------------------------*/
  27
  28 int main(int argc, char *argv[]){
  29   struct        RNAduplex_args_info args_info;
  30   char          fname[FILENAME_MAX_LENGTH], *input_string, *s1, *s2, *orig_s1, *orig_s2, *c, *ParamFile=NULL, *ns_bases=NULL;
  31   unsigned int  input_type;
  32   int           i, l, sym, r;
  33   int           istty, delta=-1;
  34   int           noconv=0;
  35
  36   s1 = s2 = orig_s1 = orig_s2 = NULL;
  37
  38   dangles = 2;
  39
  40   /*
  41   #############################################
  42   # check the command line parameters
  43   #############################################
  44   */
  45   if(RNAduplex_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
  46   /* temperature */
  47   if(args_info.temp_given)        temperature = args_info.temp_arg;
  48   /* do not take special tetra loop energies into account */
  49   if(args_info.noTetra_given)     tetra_loop=0;
  50   /* set dangle model */
  51   if(args_info.dangles_given){
  52     if((args_info.dangles_arg < 0) || (args_info.dangles_arg > 3))
  53       warn_user("required dangle model not implemented, falling back to default dangles=2");
  54     else
  55       dangles = args_info.dangles_arg;
  56   }
  57   /* do not allow weak pairs */
  58   if(args_info.noLP_given)        noLonelyPairs = 1;
  59   /* do not allow wobble pairs (GU) */
  60   if(args_info.noGU_given)        noGU = 1;
  61   /* do not allow weak closing pairs (AU,GU) */
  62   if(args_info.noClosingGU_given) no_closingGU = 1;
  63   /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
  64   if(args_info.noconv_given)      noconv = 1;
  65   /* take another energy parameter set */
  66   if(args_info.paramFile_given)   ParamFile = strdup(args_info.paramFile_arg);
  67   /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
  68   if(args_info.nsp_given)         ns_bases = strdup(args_info.nsp_arg);
  69   /*energy range */
  70   if(args_info.deltaEnergy_given) delta = (int) (0.1+args_info.deltaEnergy_arg*100);
  71   /* sorted output */
  72   if(args_info.sorted_given)      subopt_sorted = 1;
  73
  74   /* free allocated memory of command line data structure */
  75   RNAduplex_cmdline_parser_free (&args_info);
  76
  77   /*
  78   #############################################
  79   # begin initializing
  80   #############################################
  81   */
  82   if (ParamFile != NULL)
  83     read_parameter_file(ParamFile);
  84
  85   if (ns_bases != NULL) {
  86     nonstandards = space(33);
  87     c=ns_bases;
  88     i=sym=0;
  89     if (*c=='-') {
  90       sym=1; c++;
  91     }
  92     while (*c!='\0') {
  93       if (*c!=',') {
  94         nonstandards[i++]=*c++;
  95         nonstandards[i++]=*c;
  96         if ((sym)&&(*c!=*(c-1))) {
  97           nonstandards[i++]=*c;
  98           nonstandards[i++]=*(c-1);
  99         }
 100       }
 101       c++;
 102     }
 103   }
 104   istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
 105
 106   /*
 107   #############################################
 108   # main loop: continue until end of file
 109   #############################################
 110   */
 111   do{
 112     duplexT mfe, *subopt;
 113     /*
 114     ########################################################
 115     # handle user input from 'stdin'
 116     ########################################################
 117     */
 118     if(istty) print_tty_input_seq_str("Input two sequences (one line each)");
 119
 120     /* extract filename from fasta header if available */
 121     fname[0] = '\0';
 122     while((input_type = get_input_line(&input_string, 0)) == VRNA_INPUT_FASTA_HEADER){
 123       printf(">%s\n", input_string);
 124       (void) sscanf(input_string, "%" XSTR(FILENAME_ID_LENGTH) "s", fname);
 125       free(input_string);
 126     }
 127
 128     /* break on any error, EOF or quit request */
 129     if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)){ break;}
 130     /* else assume a proper sequence of letters of a certain alphabet (RNA, DNA, etc.) */
 131     else{
 132       s1 = strdup(input_string);
 133       free(input_string);
 134     }
 135
 136     /* get second sequence */
 137     while((input_type = get_input_line(&input_string, 0)) == VRNA_INPUT_FASTA_HEADER){
 138       printf(">%s\n", input_string);
 139       free(input_string);
 140     }
 141     /* break on any error, EOF or quit request */
 142     if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)){ break;}
 143     /* else assume a proper sequence of letters of a certain alphabet (RNA, DNA, etc.) */
 144     else{
 145       s2 = strdup(input_string);
 146       free(input_string);
 147     }
 148
 149     /* convert DNA alphabet to RNA if not explicitely switched off */
 150     if(!noconv){
 151       str_DNA2RNA(s1);
 152       str_DNA2RNA(s2);
 153     }
 154     /* store case-unmodified sequence */
 155     orig_s1 = strdup(s1);
 156     orig_s2 = strdup(s2);
 157     /* convert sequence to uppercase letters only */
 158     str_uppercase(s1);
 159     str_uppercase(s2);
 160
 161     if (istty) printf("lengths = %d,%d\n", strlen(s1), strlen(s2));
 162
 163     /*
 164     ########################################################
 165     # done with 'stdin' handling, now init everything properly
 166     ########################################################
 167     */
 168     update_fold_params();
 169
 170      /*
 171     ########################################################
 172     # begin actual computations
 173     ########################################################
 174     */
 175     if (delta>=0) {
 176       duplexT *sub;
 177       subopt = duplex_subopt(s1, s2, delta, 5);
 178       for (sub=subopt; sub->i >0; sub++) {
 179         print_struc(sub);
 180         free(sub->structure);
 181       }
 182       free(subopt);
 183     }
 184     else {
 185       mfe = duplexfold(s1, s2);
 186       print_struc(&mfe);
 187       free(mfe.structure);
 188     }
 189     (void) fflush(stdout);
 190     free(s1);
 191     free(s2);
 192     free(orig_s1);
 193     free(orig_s2);
 194     s1 = s2 = orig_s1 = orig_s2 = NULL;
 195   } while (1);
 196   return 0;
 197 }
 198
 199 PRIVATE void print_struc(duplexT const *dup) {
 200   int l1;
 201   l1 = strchr(dup->structure, '&')-dup->structure;
 202   printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n", dup->structure, dup->i+1-l1,
 203          dup->i, dup->j, dup->j+strlen(dup->structure)-l1-2, dup->energy);
 204 }