binaries/src/ViennaRNA/Progs/RNAdistance.c

   1 /*
   2                     Distances of Secondary Structures
   3            Walter Fontana, Ivo L Hofacker, Peter F Stadler
   4                           Vienna RNA Package
   5 */
   6 #include <stdio.h>
   7 #include <stdlib.h>
   8 #include <math.h>
   9 #include <ctype.h>
  10 #include <unistd.h>
  11 #include <string.h>
  12 #include "dist_vars.h"
  13 #include "RNAstruct.h"
  14 #include "treedist.h"
  15 #include "stringdist.h"
  16 #include "utils.h"
  17 #include "data_structures.h"
  18 #include "RNAdistance_cmdl.h"
  19
  20 #define MAXNUM      1000    /* max number of structs for distance matrix */
  21
  22 #define PUBLIC
  23 #define PRIVATE     static
  24 /*@unused@*/
  25 static char rcsid[] = "$Id: RNAdistance.c,v 1.8 2005/07/24 08:35:15 ivo Exp $";
  26 PRIVATE void command_line(int argc, char *argv[]);
  27 PRIVATE int parse_input(char *line);
  28 PRIVATE int check_tree(char *line, char alpha[]);
  29 PRIVATE int check_brackets(char *line);
  30 PRIVATE void print_aligned_lines(FILE *somewhere);
  31
  32 PRIVATE char  ruler[] ="....,....1....,....2....,....3....,....4"
  33 "....,....5....,....6....,....7....,....8";
  34 PRIVATE int types=1;
  35 PRIVATE int task;
  36 PRIVATE int taxa_list;
  37 PRIVATE char outfile[FILENAME_MAX_LENGTH], *list_title;
  38
  39 PRIVATE char ttype[10]="f";
  40 PRIVATE int n=0;
  41
  42 int main(int argc, char *argv[])
  43 {
  44   char     *line=NULL, *xstruc, *cc;
  45   Tree     *T[10][MAXNUM];
  46   int       tree_types = 0, ttree;
  47   swString *S[10][MAXNUM];
  48   char     *P[MAXNUM];   /* structures for base pair distances */
  49   int       string_types = 0, tstr;
  50   int       i,j, tt, istty, type;
  51   int       it, is;
  52   FILE     *somewhere=NULL;
  53
  54   command_line(argc, argv);
  55
  56   if((outfile[0]=='\0')&&(task==2)&&(edit_backtrack))
  57     strcpy(outfile,"backtrack.file");
  58   if(outfile[0]!='\0') somewhere = fopen(outfile,"w");
  59   if (somewhere==NULL) somewhere = stdout;
  60   istty = isatty(fileno(stdin))&&isatty(fileno(stdout));
  61
  62   do {
  63     if ((istty)&&(n==0)) {
  64       printf("\nInput structure;  @ to quit\n");
  65       printf("%s\n", ruler);
  66     }
  67     do {
  68       if (line!=NULL) free(line);
  69       line=get_line(stdin);
  70     } while ((type=parse_input(line))==0);
  71
  72     if (((type==999)||(type==888))&&(task==2)) {  /* do matrices */
  73       if (taxa_list)
  74         printf("* END of taxa list\n");
  75
  76       ttree = 0; tstr = 0;
  77       for (tt=0; tt< types; tt++) {
  78         printf("> %c   %d\n", ttype[tt], n);
  79         if(islower(ttype[tt])) {
  80           for (i=1; i<n; i++) {
  81             for (j=0; j<i; j++) {
  82               printf("%g ",tree_edit_distance(T[ttree][i], T[ttree][j]));
  83               if(edit_backtrack) {
  84                 fprintf(somewhere,"%d %d",i+1,j+1);
  85                 if (ttype[tt]=='f') unexpand_aligned_F(aligned_line);
  86                 print_aligned_lines(somewhere);
  87               }
  88             }
  89             printf("\n");
  90           }
  91           printf("\n");
  92           for (i=0; i<n; i++) free_tree(T[ttree][i]);
  93           ttree++;
  94         }
  95         if (ttype[tt]=='P') {
  96           for (i=1; i<n; i++) {
  97             for (j=0; j<i; j++)
  98               printf("%g ", (float) bp_distance(P[i], P[j]));
  99             printf("\n");
 100           }
 101           printf("\n");
 102           for (i=0; i<n; i++) free(P[i]);
 103         }
 104         else if(isupper(ttype[tt])) {
 105           for (i=1; i<n; i++) {
 106             for (j=0; j<i; j++) {
 107               printf("%g ",string_edit_distance(S[tstr][i], S[tstr][j]));
 108               if (edit_backtrack) {
 109                 fprintf(somewhere,"%d %d",i+1,j+1);
 110                 if (ttype[tt]=='F') unexpand_aligned_F(aligned_line);
 111                 print_aligned_lines(somewhere);
 112               }
 113             }
 114             printf("\n");
 115           }
 116           printf("\n");
 117           for (i=0; i<n; i++) free(S[tstr][i]);
 118           tstr++;
 119         }
 120       }
 121       fflush(stdout);
 122       if (type==888) {  /* do another distance matrix */
 123         n = 0;
 124         printf("%s\n", list_title);
 125         free(list_title);
 126         continue;
 127       }
 128     }                       /* END do matrices */
 129     if (type==999) {   /* finito */
 130       if (outfile[0]!='\0') fclose(somewhere);
 131       return 0;
 132     }
 133
 134     if (type<0) {
 135       xstruc = add_root(line);
 136       free(line);
 137       line = xstruc;
 138       type = -type;
 139     }
 140     if (type==2) {
 141       xstruc = unexpand_Full(line);
 142       free(line);
 143       line = xstruc;
 144       type=1;
 145     }
 146     tree_types   = 0;
 147     string_types = 0;
 148     for(tt=0; tt < types; tt++) {
 149       switch(ttype[tt]){
 150       case 'f' :
 151       case 'F' :
 152         if (type!=1) nrerror("Can't convert back to full structure");
 153         xstruc = expand_Full(line);
 154         if(islower(ttype[tt])) {  /* tree_edit */
 155           T[tree_types++][n] = make_tree(xstruc);
 156         }
 157         if(isupper(ttype[tt])) { /* string edit */
 158           S[string_types++][n] = Make_swString(xstruc);
 159         }
 160         free(xstruc);
 161         break;
 162       case 'P':
 163         if (type!=1) nrerror("Can't convert back to full structure");
 164         P[n] = strdup(line);
 165         break;
 166       case 'h' :
 167       case 'H' :
 168         switch (type) {
 169         case 1:
 170           xstruc = b2HIT(line);
 171           if(islower(ttype[tt])) {
 172             T[tree_types++][n] = make_tree(xstruc);
 173           }
 174           if(isupper(ttype[tt])) {
 175             S[string_types++][n] = Make_swString(xstruc);
 176           }
 177           free(xstruc);
 178           break;
 179         default:
 180           nrerror("Can't convert to HIT structure");
 181         }
 182         break;
 183       case 'c' :
 184       case 'C' :
 185         switch (type) {
 186         case 1:
 187           cc = b2C(line);
 188           xstruc = expand_Shapiro(cc);
 189           free(cc);
 190           break;
 191         case 4:
 192           cc = expand_Shapiro(line);
 193           xstruc = unweight(cc);
 194           free(cc);
 195           break;
 196         case 3:
 197           xstruc = unweight(line);
 198           break;
 199         default:
 200           nrerror("Unknown structure representation");
 201           exit(0);
 202         }
 203         if(islower(ttype[tt])) {
 204           T[tree_types++][n] = make_tree(xstruc);
 205         }
 206         if(isupper(ttype[tt])) {
 207           S[string_types++][n] = Make_swString(xstruc);
 208         }
 209         free(xstruc);
 210         break;
 211       case 'w' :
 212       case 'W' :
 213         if (type==1) {
 214           xstruc = b2Shapiro(line);
 215           if(islower(ttype[tt])) {
 216             T[tree_types++][n] = make_tree(xstruc);
 217           }
 218           if(isupper(ttype[tt])) {
 219             S[string_types++][n] = Make_swString(xstruc);
 220           }
 221           free(xstruc);
 222         }
 223         else {
 224           if(islower(ttype[tt])) {
 225             T[tree_types++][n] = make_tree(line);
 226           }
 227           if(isupper(ttype[tt])) {
 228             S[string_types++][n] = Make_swString(line);
 229           }
 230         }
 231         break;
 232       default:
 233         nrerror("Unknown distance type");
 234       }
 235     }
 236     n++;
 237     switch (task) {
 238       float     dist;
 239     case 1:
 240       if (n==2) {
 241         for (it=0, is=0, i=0; i<types; i++) {
 242           if(islower(ttype[i])) {
 243             dist = tree_edit_distance(T[it][0], T[it][1]);
 244             free_tree(T[it][0]); free_tree(T[it][1]);
 245             it++;
 246           }
 247           else if (ttype[i]=='P') {
 248             dist = (float) bp_distance(P[0], P[1]);
 249             free(P[0]); free(P[1]);
 250           }
 251           else /* isupper(ttype[i]) */ {
 252             dist = string_edit_distance(S[is][0], S[is][1]);
 253             free(S[is][0]); free(S[is][1]);
 254             is++;
 255           }
 256           printf("%c: %g  ", ttype[i], dist);
 257           if ((edit_backtrack)&&(ttype[i]!='P')) {
 258             if (ttype[i]=='f') unexpand_aligned_F(aligned_line);
 259             print_aligned_lines(somewhere);
 260           }
 261         }
 262         printf("\n");
 263         n=0;
 264       }
 265       break;
 266     case 3:
 267       if (n>1) {
 268         for (it=0, is=0, i=0; i<types; i++) {
 269           if(islower(ttype[i])) {
 270             dist = tree_edit_distance(T[it][1], T[it][0]);
 271             free_tree(T[it][1]);
 272             it++;
 273           }
 274           else if (ttype[i]=='P') {
 275             dist = (float) bp_distance(P[0], P[1]);
 276             free(P[1]);
 277           }
 278           else /* if(isupper(ttype[i])) */ {
 279             dist = string_edit_distance(S[is][0], S[is][1]);
 280             free(S[is][1]);
 281             is++;
 282           }
 283           printf("%c: %g  ", ttype[i], dist);
 284           if ((edit_backtrack)&&(ttype[i]!='P')) {
 285             if (ttype[i]=='f') unexpand_aligned_F(aligned_line);
 286             print_aligned_lines(somewhere);
 287           }
 288         }
 289         printf("\n");
 290         n=1;
 291       }
 292       break;
 293     case 4:
 294       if (n>1) {
 295         for (it=0, is=0, i=0; i<types; i++) {
 296           if(islower(ttype[i])) {
 297             dist = tree_edit_distance(T[it][1], T[it][0]);
 298             free_tree(T[it][0]);
 299             T[it][0] = T[it][1];
 300             it++;
 301           }
 302           else if (ttype[i]=='P') {
 303             dist = (float) bp_distance(P[0], P[1]);
 304             free(P[0]); P[0] = P[1];
 305           }
 306           else /* if(isupper(ttype[i])) */ {
 307             dist = string_edit_distance(S[is][0], S[is][1]);
 308             free(S[is][0]);
 309             S[is][0] = S[is][1];
 310             is++;
 311           }
 312           printf("%c: %g  ", ttype[i], dist);
 313           if ((edit_backtrack)&&(ttype[i]!='P')) {
 314             if (ttype[i]=='f') unexpand_aligned_F(aligned_line);
 315             print_aligned_lines(somewhere);
 316           }
 317         }
 318         printf("\n");
 319         n=1;
 320       }
 321       break;
 322     }
 323     fflush(stdout);
 324   } while(type!=999);
 325   return 0;
 326 }
 327
 328 /*--------------------------------------------------------------------------*/
 329
 330 PRIVATE int parse_input(char *line)
 331 {
 332   int type, rooted=0, i, xx;
 333   char *cp;
 334
 335   if (line==NULL) return 999;
 336   if (line[0]=='*') {
 337     if (taxa_list==0) {
 338       if (task==2) taxa_list=1;
 339       printf("%s\n", line);
 340       return 0;
 341     } else {
 342       list_title = strdup(line);
 343       return 888;
 344     }
 345   }
 346   if (line[0]=='>') {
 347     if (taxa_list)
 348       printf("%d :%s\n", n+1, line+1);
 349     else printf("%s\n", line);
 350     return 0;
 351   }
 352
 353   cp = strchr(line,' ');
 354   if (cp) *cp='\0';      /* get rid of junk at the end of line */
 355
 356   switch (line[0]) {
 357   case '.' :
 358     type = 1;
 359     break;
 360
 361   case '(' :            /* it's a tree */
 362     i=1;
 363     xx = 0;
 364     type = 4;           /* coarse */
 365     rooted = 0;
 366     while (line[i])  {
 367       if (line[i]=='.'){
 368         type = 1;     /* Full */
 369         break;
 370       }
 371       if ( (line[i]=='U')||(line[i]=='P') ) {
 372         type = 2;     /* FFull */
 373         xx = 1;
 374         break;
 375       }
 376       if (line[i]=='S') {
 377         type = 3;
 378         xx = 1;
 379         break;        /* Shapiro tree */
 380       }
 381       if ((line[i]!='(')&&(line[i]!=')')) xx = 1;
 382       i++;
 383     }
 384     if(!xx) type =1;
 385
 386     rooted = (line[strlen(line)-2]=='R');
 387     break;
 388   case '@' :
 389     return 999;
 390
 391   default:
 392     return 0;
 393   }
 394
 395   switch (type) {
 396   case 1:
 397     if(check_brackets(line))
 398       return 1;
 399     break;
 400   case 2:
 401     if(check_tree(line,"UP") ) {
 402       if(rooted) return 2;
 403       else return -2;
 404     }
 405     break;
 406   case 3:
 407     if(check_tree(line,"HBIMSE") ){
 408       if(rooted) return -3;
 409       else return -3;
 410     }
 411     break;
 412   case 4:
 413     if(check_tree(line,"HBIM") ){
 414       if(rooted) return 4;
 415       else return -4;
 416     }
 417     break;
 418   }
 419   return 0;
 420 }
 421
 422 /*--------------------------------------------------------------------------*/
 423
 424
 425 PRIVATE int check_tree(char *line, char alpha[])
 426 {
 427   int n, i, o;
 428   char *pos;
 429
 430   n = (int) strlen(line);
 431
 432   if( line[0] != '(' ) return 0;
 433   i=o=1;
 434   while( line[i] ){
 435     while( line[i] == '(' ){
 436       o++;
 437       i++;
 438     }
 439     pos=strchr(alpha, (int)line[i]);
 440     if (pos) {
 441       while (isdigit((int) line[++i]));
 442       if (line[i]!=')') return 0;
 443     }
 444     if (line[i]=='R') {
 445       i++;
 446       if ((i!=n-1)||(line[i]!=')')) return 0;
 447     }
 448     if (line[i]==')') {
 449       o--;
 450       if(o< 0) return 0;
 451       if( (i<n)&&(line[i+1]==')') ) return 0;
 452       i++;
 453     }
 454     else return 0;
 455   }
 456   if(o>0) return 0;
 457   return 1;
 458 }
 459
 460 /*--------------------------------------------------------------------------*/
 461
 462
 463 PRIVATE int check_brackets(char *line)
 464 {
 465   int i,o;
 466
 467   i=o=0;
 468   while( line[i] ){
 469     switch(line[i]) {
 470     case '(' :
 471       o++;
 472       i++;
 473       break;
 474     case '.' :
 475       i++;
 476       break;
 477     case ')' :
 478       i++;
 479       o--;
 480       if(o<0) return 0;
 481       break;
 482     default:
 483       return 0;
 484     }
 485   }
 486   if (o>0) return 0;
 487   return 1;
 488 }
 489
 490 /*--------------------------------------------------------------------------*/
 491
 492
 493 PRIVATE void command_line(int argc, char *argv[])
 494 {
 495   int i;
 496   struct        RNAdistance_args_info args_info;
 497
 498   /* default values */
 499   edit_backtrack = 0;
 500   types=1; ttype[0]='f'; task=1;
 501
 502   /*
 503   #############################################
 504   # check the command line parameters
 505   #############################################
 506   */
 507   if(RNAdistance_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
 508
 509   /* use specified distance representations */
 510   if(args_info.distance_given){
 511     strncpy(ttype, args_info.distance_arg,9);
 512     types = (int)strlen(ttype);
 513   }
 514
 515   if(args_info.compare_given){
 516     switch(args_info.compare_arg[0]){
 517       case 'p': task=1; break;
 518       case 'm': task=2; break;
 519       case 'f': task=3; break;
 520       case 'c': task=4; break;
 521       default:  RNAdistance_cmdline_parser_print_help();
 522                 exit(EXIT_FAILURE);
 523     }
 524   }
 525
 526   if(args_info.shapiro_given)
 527     cost_matrix = 1;
 528
 529   if(args_info.backtrack_given){
 530     if(strcmp(args_info.backtrack_arg, "none")){
 531       strncpy(outfile, args_info.backtrack_arg, FILENAME_MAX_LENGTH-1);
 532     }
 533     edit_backtrack = 1;
 534   }
 535
 536   /* free allocated memory of command line data structure */
 537   RNAdistance_cmdline_parser_free (&args_info);
 538 }
 539
 540 /*--------------------------------------------------------------------------*/
 541
 542 PRIVATE void print_aligned_lines(FILE *somewhere)
 543 {
 544   if (edit_backtrack) {
 545     fprintf(somewhere, "\n%s\n%s\n", aligned_line[0], aligned_line[1]);
 546     fflush(somewhere);
 547   }
 548 }