binaries/src/ViennaRNA/Progs/RNApdist.c

   1 /*
   2                 Distances of Secondary Structure Ensembles
   3           Peter F Stadler, Ivo L Hofacker, Sebastian Bonhoeffer
   4                         Vienna RNA Package
   5 */
   6 #include <stdio.h>
   7 #include <stdlib.h>
   8 #include <math.h>
   9 #include <ctype.h>
  10 #include <errno.h>
  11 #include <string.h>
  12 #include <unistd.h>
  13 #include "part_func.h"
  14 #include "fold_vars.h"
  15 #include "dist_vars.h"
  16 #include "utils.h"
  17 #include "PS_dot.h"
  18 #include "read_epars.h"
  19 #include "profiledist.h"
  20 #include "RNApdist_cmdl.h"
  21
  22
  23 #define MAXLENGTH  10000
  24 #define MAXSEQ      1000
  25 /*@unused@*/
  26 static char rcsid[] = "$Id: RNApdist.c,v 1.8 2002/11/07 12:19:41 ivo Exp $";
  27
  28 PRIVATE void command_line(int argc, char *argv[]);
  29 PRIVATE void usage(void);
  30 PRIVATE void print_aligned_lines(FILE *somewhere);
  31
  32 PRIVATE char task;
  33 PRIVATE char outfile[FILENAME_MAX_LENGTH];
  34 PRIVATE char  ruler[] ="....,....1....,....2....,....3....,....4"
  35                        "....,....5....,....6....,....7....,....8";
  36 static int noconv = 0;
  37
  38 int main(int argc, char *argv[])
  39
  40 {
  41   float     *T[MAXSEQ];
  42   int        i,j, istty, n=0;
  43   int        type, length, taxa_list=0;
  44   float      dist;
  45   FILE      *somewhere=NULL;
  46   char      *structure;
  47   char      *line=NULL, fname[FILENAME_MAX_LENGTH], *list_title=NULL;
  48   plist     *pr_pl, *mfe_pl;
  49
  50   pr_pl = mfe_pl = NULL;
  51
  52   command_line(argc, argv);
  53
  54   if((outfile[0]=='\0')&&(task=='m')&&(edit_backtrack))
  55     strcpy(outfile,"backtrack.file");
  56   if (outfile[0]!='\0') somewhere = fopen(outfile,"w");
  57   if (somewhere==NULL) somewhere = stdout;
  58   istty   = (isatty(fileno(stdout))&&isatty(fileno(stdin)));
  59
  60   while (1) {
  61     if ((istty)&&(n==0)) {
  62       printf("\nInput sequence;  @ to quit\n");
  63       printf("%s\n", ruler);
  64     }
  65
  66     type = 0;
  67     do {  /* get sequence to fold */
  68       if (line!=NULL) free(line);
  69       *fname='\0';
  70       if ((line=get_line(stdin))==NULL) {type = 999; break;}
  71       if (line[0]=='@') type = 999;
  72       if (line[0]=='*') {
  73         if (taxa_list==0) {
  74           if (task=='m') taxa_list=1;
  75           printf("%s\n", line);
  76           type = 0;
  77         } else {
  78           list_title = strdup(line);
  79           type = 888;
  80         }
  81       }
  82       if (line[0]=='>') {
  83         if (sscanf(line,">%" XSTR(FILENAME_ID_LENGTH) "s", fname)!=0)
  84           strcat(fname, "_dp.ps");
  85         if (taxa_list)
  86           printf("%d : %s\n", n+1, line+1);
  87         else printf("%s\n",line);
  88         type = 0;
  89       }
  90       if (isalpha(line[0]))  {
  91         char *cp;
  92         cp =strchr(line,' ');
  93         if (cp) *cp='\0';
  94         type = 1;
  95       }
  96     } while(type==0);
  97
  98     if( (task == 'm')&&(type>800) ) {
  99       if (taxa_list)
 100         printf("* END of taxa list\n");
 101       printf("> p %d (pdist)\n",n);
 102       for (i=1; i<n; i++) {
 103         for (j=0; j<i; j++) {
 104           printf("%g ",profile_edit_distance(T[i], T[j]));
 105           if(edit_backtrack) fprintf(somewhere,"> %d %d\n",i+1,j+1);
 106           print_aligned_lines(somewhere);
 107         }
 108         printf("\n");
 109       }
 110       if (type==888) {  /* do another distance matrix */
 111         n = 0;
 112         printf("%s\n", list_title);
 113         free(list_title);
 114       }
 115     }
 116
 117     if(type>800) {
 118       for (i=0; i<n; i++)
 119         free_profile(T[i]);
 120       if (type == 888) continue;
 121       if (outfile[0]!='\0') (void) fclose(somewhere);
 122       if (line!= NULL) free(line);
 123       return 0; /* finito */
 124     }
 125
 126     length = (int) strlen(line);
 127     for (i=0; i<length; i++) {
 128       line[i]=toupper(line[i]);
 129       if (!noconv && line[i] == 'T') line[i] = 'U';
 130     }
 131
 132     /* init_pf_fold(length); <- obsolete */
 133     structure = (char *) space((length+1)*sizeof(char));
 134     (void) pf_fold(line,structure);
 135
 136     if (*fname=='\0')
 137       sprintf(fname, "%d_dp.ps", n+1);
 138
 139     /* PS_dot_plot(line, fname); <- NOT THREADSAFE and obsolete function! */
 140
 141     /* get pairlist of probability matrix */
 142     assign_plist_from_pr(&pr_pl, pr, length, 1e-5);
 143     /* no previous mfe call thus no mfe structure information known */
 144     mfe_pl = (plist *)space(sizeof(plist));
 145     mfe_pl[0].i = mfe_pl[0].j = 0;
 146
 147     /* call threadsafe dot plot printing function */
 148     PS_dot_plot_list(line, fname, pr_pl, mfe_pl, "");
 149
 150     T[n] = Make_bp_profile_bppm(pr, length);
 151     if((istty)&&(task=='m')) printf("%s\n",structure);
 152     free(structure);
 153     free(mfe_pl);
 154     free(pr_pl);
 155     free_pf_arrays();
 156
 157     n++;
 158     switch (task) {
 159     case 'p' :
 160       if (n==2) {
 161         dist = profile_edit_distance(T[0],T[1]);
 162         printf("%g\n",dist);
 163         print_aligned_lines(somewhere);
 164         free_profile(T[0]);
 165         free_profile(T[1]);
 166         n=0;
 167       }
 168       break;
 169     case 'f' :
 170       if (n>1) {
 171         dist = profile_edit_distance(T[1], T[0]);
 172         printf("%g\n",dist);
 173         print_aligned_lines(somewhere);
 174         free_profile(T[1]);
 175         n=1;
 176       }
 177       break;
 178     case 'c' :
 179       if (n>1) {
 180         dist = profile_edit_distance(T[1], T[0]);
 181         printf("%g\n",dist);
 182         print_aligned_lines(somewhere);
 183         free_profile(T[0]);
 184         T[0] = T[1];
 185         n=1;
 186       }
 187       break;
 188
 189     case 'm' :
 190       break;
 191
 192     default :
 193       nrerror("This can't happen.");
 194     }    /* END switch task */
 195     (void) fflush(stdout);
 196   }    /* END while */
 197   if (line !=NULL) free(line);
 198   return 0;
 199 }
 200
 201 /* ----------------------------------------------------------------- */
 202
 203 PRIVATE void command_line(int argc, char *argv[])
 204 {
 205
 206   int i, sym;
 207   char *ns_bases=NULL, *c;
 208   char *ParamFile=NULL;
 209   struct  RNApdist_args_info args_info;
 210
 211   task = 'p';
 212
 213   /*
 214   #############################################
 215   # check the command line parameters
 216   #############################################
 217   */
 218   if(RNApdist_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
 219
 220   /* temperature */
 221   if(args_info.temp_given)
 222     temperature = args_info.temp_arg;
 223
 224   /* do not take special tetra loop energies into account */
 225   if(args_info.noTetra_given)
 226     tetra_loop=0;
 227
 228   /* set dangle model */
 229   if(args_info.dangles_given){
 230     dangles = args_info.dangles_arg;
 231     if(dangles) dangles = 2;
 232   }
 233
 234   /* set energy model */
 235   if(args_info.energyModel_given)
 236     energy_set = args_info.energyModel_arg;
 237
 238   /* do not allow weak pairs */
 239   if(args_info.noLP_given)
 240     noLonelyPairs = 1;
 241
 242   /* do not allow wobble pairs (GU) */
 243   if(args_info.noGU_given)
 244     noGU = 1;
 245
 246   /* do not allow weak closing pairs (AU,GU) */
 247   if(args_info.noClosingGU_given)
 248     no_closingGU = 1;
 249
 250   /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
 251   if(args_info.noconv_given)
 252     noconv = 1;
 253
 254   /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
 255   if(args_info.nsp_given)
 256     ns_bases = strdup(args_info.nsp_arg);
 257
 258   /* take another energy parameter set */
 259   if(args_info.paramFile_given)
 260     ParamFile = strdup(args_info.paramFile_arg);
 261
 262   if(args_info.compare_given){
 263     switch(args_info.compare_arg[0]){
 264       case 'p':
 265       case 'm':
 266       case 'f':
 267       case 'c': task=args_info.compare_arg[0];
 268                 break;
 269       default:  RNApdist_cmdline_parser_print_help();
 270                 exit(EXIT_FAILURE);
 271     }
 272   }
 273
 274   if(args_info.backtrack_given){
 275     if(strcmp(args_info.backtrack_arg, "none")){
 276       strncpy(outfile, args_info.backtrack_arg, FILENAME_MAX_LENGTH-1);
 277     }
 278     edit_backtrack = 1;
 279   }
 280
 281   /* free allocated memory of command line data structure */
 282   RNApdist_cmdline_parser_free (&args_info);
 283
 284   /* do some preparations */
 285   if (ParamFile!=NULL)
 286     read_parameter_file(ParamFile);
 287
 288   if (ns_bases!=NULL) {
 289     nonstandards = space(33);
 290     c=ns_bases;
 291     i=sym=0;
 292     if (*c=='-') {
 293       sym=1; c++;
 294     }
 295     while (*c!='\0') {
 296       if (*c!=',') {
 297         nonstandards[i++]=*c++;
 298         nonstandards[i++]=*c;
 299         if ((sym)&&(*c!=*(c-1))) {
 300           nonstandards[i++]=*c;
 301           nonstandards[i++]=*(c-1);
 302         }
 303       }
 304       c++;
 305     }
 306   }
 307
 308 }
 309
 310 /* ------------------------------------------------------------------------- */
 311
 312 PRIVATE void usage(void)
 313 {
 314   nrerror("usage: RNApdist [-Xpmfc] [-B [file]] [-T temp] [-4] [-d] [-noGU]\n"
 315           "                [-noCloseGU] [-noLP] [-e e_set] [-P paramfile] [-nsp pairs]");
 316 }
 317
 318 /*--------------------------------------------------------------------------*/
 319
 320 PRIVATE void print_aligned_lines(FILE *somewhere)
 321 {
 322   if (edit_backtrack)
 323     fprintf(somewhere, "%s\n%s\n", aligned_line[0], aligned_line[1]);
 324 }
 325
 326 /*--------------------------------------------------------------------------*/