binaries/src/ViennaRNA/Progs/RNAaliduplex.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 "aln_util.h"
  18 #include "read_epars.h"
  19 #include "subopt.h"
  20 #include "duplex.h"
  21 #include "RNAaliduplex_cmdl.h"
  22
  23 /*@unused@*/
  24 static char rcsid[] = "$Id: RNAaliduplex.c,v 1.1 2007/08/26 10:08:44 ivo Exp $";
  25
  26 PRIVATE void  print_struc(duplexT const *dup);
  27
  28 #define MAX_NUM_NAMES    500
  29
  30 int main(int argc, char *argv[]){
  31   struct        RNAaliduplex_args_info args_info;
  32   char          *AS1[MAX_NUM_NAMES], *AS2[MAX_NUM_NAMES], *names1[MAX_NUM_NAMES], *names2[MAX_NUM_NAMES];
  33   char          *ParamFile=NULL, *c, *ns_bases=NULL;
  34   int           i, r, n_seq, n_seq2;
  35   int           istty, delta=-1, sym;
  36   int           noconv=0;
  37   duplexT       mfe, *subopt;
  38   FILE          *file1=NULL, *file2=NULL; /* input alignments */
  39
  40   dangles = 2;
  41
  42   /*
  43   #############################################
  44   # check the command line parameters
  45   #############################################
  46   */
  47   if(RNAaliduplex_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
  48   /* temperature */
  49   if(args_info.temp_given)        temperature = args_info.temp_arg;
  50   /* do not take special tetra loop energies into account */
  51   if(args_info.noTetra_given)     tetra_loop=0;
  52   /* set dangle model */
  53   if(args_info.dangles_given){
  54     if((args_info.dangles_arg < 0) || (args_info.dangles_arg > 3))
  55       warn_user("required dangle model not implemented, falling back to default dangles=2");
  56     else
  57      dangles = args_info.dangles_arg;
  58   }
  59   /* do not allow weak pairs */
  60   if(args_info.noLP_given)        noLonelyPairs = 1;
  61   /* do not allow wobble pairs (GU) */
  62   if(args_info.noGU_given)        noGU = 1;
  63   /* do not allow weak closing pairs (AU,GU) */
  64   if(args_info.noClosingGU_given) no_closingGU = 1;
  65   /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
  66   if(args_info.noconv_given)      noconv = 1;
  67   /* take another energy parameter set */
  68   if(args_info.paramFile_given)   ParamFile = strdup(args_info.paramFile_arg);
  69   /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
  70   if(args_info.nsp_given)         ns_bases = strdup(args_info.nsp_arg);
  71   /*energy range */
  72   if(args_info.deltaEnergy_given) delta = (int) (0.1+args_info.deltaEnergy_arg*100);
  73   /* sorted output */
  74   if(args_info.sorted_given)      subopt_sorted = 1;
  75
  76   /* check unnamed options a.k.a. filenames of input alignments */
  77   if(args_info.inputs_num == 2){
  78     file1 = fopen(args_info.inputs[0], "r");
  79     if(file1 == NULL){
  80       fprintf(stderr, "can't open %s\n", args_info.inputs[0]);
  81     }
  82     file2 = fopen(args_info.inputs[1], "r");
  83     if(file2 == NULL) {
  84       fprintf(stderr, "can't open %s\n", args_info.inputs[1]);
  85     }
  86   }
  87   else{
  88     RNAaliduplex_cmdline_parser_print_help();
  89     exit(1);
  90   }
  91
  92   if (!(file1 && file2)){
  93     fprintf(stderr, "No input files");
  94     RNAaliduplex_cmdline_parser_print_help();
  95     exit(1);
  96   }
  97   else{
  98     n_seq   = read_clustal(file1, AS1, names1);
  99     n_seq2  = read_clustal(file2, AS2, names2);
 100     fclose(file1);
 101     fclose(file2);
 102     if(n_seq != n_seq2) nrerror("unequal number of seqs in alignments");
 103   }
 104
 105   /* free allocated memory of command line data structure */
 106   RNAaliduplex_cmdline_parser_free (&args_info);
 107
 108   /*
 109   #############################################
 110   # begin initializing
 111   #############################################
 112   */
 113
 114
 115   if (ParamFile != NULL)
 116     read_parameter_file(ParamFile);
 117
 118   if (ns_bases != NULL) {
 119     nonstandards = space(33);
 120     c=ns_bases;
 121     i=sym=0;
 122     if (*c=='-') {
 123       sym=1; c++;
 124     }
 125     while (*c!='\0') {
 126       if (*c!=',') {
 127         nonstandards[i++]=*c++;
 128         nonstandards[i++]=*c;
 129         if ((sym)&&(*c!=*(c-1))) {
 130           nonstandards[i++]=*c;
 131           nonstandards[i++]=*(c-1);
 132         }
 133       }
 134       c++;
 135     }
 136   }
 137
 138   istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
 139
 140   /*
 141   #############################################
 142   # begin calculations
 143   #############################################
 144   */
 145   update_fold_params();
 146
 147   if (delta>=0) {
 148     duplexT *sub;
 149     subopt = aliduplex_subopt((const char **)AS1, (const char **)AS2, delta, 5);
 150     for (sub=subopt; sub->i >0; sub++) {
 151       print_struc(sub);
 152       free(sub->structure);
 153     }
 154     free(subopt);
 155   } else {
 156     mfe = aliduplexfold((const char **)AS1, (const char **)AS2);
 157     print_struc(&mfe);
 158     free(mfe.structure);
 159   }
 160   for (i=0; AS1[i]; i++) {
 161     free(AS1[i]); free(names1[i]);
 162     free(AS2[i]); free(names2[i]);
 163   }
 164   return 0;
 165 }
 166
 167 PRIVATE void print_struc(duplexT const *dup) {
 168   int l1;
 169   l1 = strchr(dup->structure, '&')-dup->structure;
 170   printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n", dup->structure, dup->i+1-l1,
 171          dup->i, dup->j, dup->j+strlen(dup->structure)-l1-2, dup->energy);
 172 }