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