binaries/src/ViennaRNA/Progs/RNAsnoop.c

   1 /* Last changed Time-stamp: <2007-12-05 13:55:42 ronny> */
   2 /*
   3                   Ineractive Access to folding Routines
   4
   5                   c Ivo L Hofacker
   6                   Vienna RNA package
   7 */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <math.h>
  12 #include <ctype.h>
  13 #include <unistd.h>
  14 #include <string.h>
  15 #include <time.h>
  16 #include "snofold.h"
  17 #include "fold.h"
  18 #include "energy_par.h"
  19 #include "snoop.h"
  20 #include "part_func.h"
  21 #include "fold_vars.h"
  22 #include "PS_dot.h"
  23 #include "utils.h"
  24 #include "aln_util.h"
  25 #include "RNAsnoop_cmdl.h"
  26 static void  aliprint_struc(snoopT *dup, const char **s1, const char **s2, char**,char**,int count, int);
  27 static void  print_struc(snoopT *dup, const char *s1, const char *s2, char*, char*, int,int);
  28 /* static void  print_struc_L(snoopT const *dup, const char *s1, const char *s2, char*, char*); */
  29
  30
  31 static void redraw_output(char *fname, char *output, int plfold_up_flag, char *suffix, char *access, char *sname, char *tname);
  32 static int ** read_plfold_i(char *fname, const int beg, const int end); /* read plfold profiles */
  33 static int ** read_rnaup(char *fname, const int beg, const int end);
  34 static int get_max_u(const char *s, char delim);
  35 extern int cut_point;
  36 /*@unused@*/
  37 static char UNUSED rcsid[] = "$Id: RNAfold.c,v 1.23 2007/12/05 13:04:08 ivo Exp $";
  38
  39 #define PRIVATE static
  40 #define MAX_NUM_NAMES    500
  41
  42 /*--------------------------------------------------------------------------*/
  43
  44 int main(int argc, char *argv[])
  45 {
  46   struct        RNAsnoop_args_info args_info;
  47   char *string_s, *line_s, *name_s, *temp_s; /*string for snoRNA*/
  48   char *string_t, *line_t, *temp_t; /*string for target RNA*/
  49   char *structure=NULL, *cstruc=NULL;
  50   char  *sname=NULL,*tname=NULL/*name of the sno RNa file, mRNA file respectively*/;
  51   char *name=NULL;
  52   char *access;char *result_file;
  53   char *output_directory;
  54   int fullStemEnergy;
  55   output_directory=NULL;
  56   result_file=NULL;
  57   access=NULL;
  58   char *suffix;
  59   suffix=NULL;
  60   FILE *sno, *mrna;
  61   int fast;fast=0;
  62   /* long int elapTicks; */
  63   /* clock_t Begin, End; */
  64   int plfold_up_flag=0;
  65   int   nice,i,r,l,length_s,/* , length_t,  */
  66    penalty,         /*extension penalty*/
  67    threshloop,         /*energy threshold on loop*/
  68    threshLE,         /*energy threshold on the S2 part*/
  69    threshRE,        /*energy threshold on the S1 part*/
  70    threshDE,         /*Total duplex energy threshold*/
  71    threshTE,        /*Total duplex + Loop energy*/
  72    threshSE,        /*Sum of all energies*/
  73    threshD,            /*lower stem energy*/
  74    half_stem,         /*minimal length of stem*/
  75    max_half_stem,        /*maximal length of stem*/
  76    max_s2                /*maximale position wo stem anfangen darf in 3->5*/,
  77    min_s2                 /*minimale position wo stem anfangen darf in 3->5*/,
  78    max_s1                /*minimal position wo stem enden darf in 3->5 */,
  79    min_s1                /*maximale position wo stem enden darf in 3->5*/,
  80    distance        /*distance between two subopts*/,
  81    min_d1                /*minimal distance between 5' sno and first duplex*/,
  82    min_d2                /*minimal distance between 3' sno and second duplex*/,
  83    max_asymm        /*maximal asymmetry in the stem interior loop*/,
  84    alignment_length/*maximal target RNA alignment length*/,
  85    alignment       /*flag to use ali version*/,
  86     delta,
  87    redraw          /*if used (option I) allow to redraw command line output into ps files */;
  88
  89   int  noconv=0;
  90   string_s=NULL;
  91   string_t=NULL;
  92   plfold_up_flag=0;
  93   alignment=0;
  94   redraw=0;
  95   nice=0;
  96   fast=0;
  97   snoop_subopt_sorted=0;
  98   /*
  99   #############################################
 100   # check the command line parameters
 101   #############################################
 102   */
 103   if(RNAsnoop_cmdline_parser(argc,argv, &args_info)!=0) exit(1);
 104   /* Redraw results from RNAsnoop */
 105   if(args_info.constraint_given)  fold_constrained=1;
 106   /* structure constraint */
 107   if(args_info.produce_ps_given)  redraw=1;
 108   /*output directory for RNAsnoop results */
 109   output_directory=strdup(args_info.output_directory_arg);
 110   /* Draw directly nice pictures from RNAsnoop*/
 111   if(args_info.direct_redraw_given)  nice=1;
 112   /*Accessibility files are coming from RNAup*/
 113   if(args_info.from_RNAup_given) {plfold_up_flag=2;access=strdup(args_info.from_RNAup_arg);}
 114   /*Suffix for the name of the accessibility files*/
 115   suffix=strdup(args_info.suffix_arg);
 116   /*Accessibility files are coming from RNAup*/
 117   if(args_info.from_RNAplfold_given) {plfold_up_flag=1;access=strdup(args_info.from_RNAplfold_arg);}
 118   /*Check if we are in alignment mode*/
 119   if(args_info.alignment_mode_given) alignment=1;
 120   /*Name of the stems file*/
 121   if(args_info.query_given) sname=strdup(args_info.query_arg);
 122   /*Name of the target file*/
 123   if(args_info.target_given) tname=strdup(args_info.target_arg);
 124   /*Delta for suboptimals*/
 125   delta=(int) (100*(args_info.energy_threshold_arg+0.1));
 126   /*fast folding and target search*/
 127   fast=args_info.fast_folding_arg;
 128   /*Penalty for duplex extension*/
 129   penalty=args_info.extension_cost_arg;
 130   /*threshold on loop energy*/
 131   threshloop=args_info.minimal_loop_energy_arg;
 132   /*threshold on right duplex*/
 133   threshRE=args_info.minimal_right_duplex_arg;
 134   /*threshold on left duplex*/
 135   threshLE=args_info.minimal_left_duplex_arg;
 136   /*threshold on minimal duplex*/
 137   threshDE=args_info.minimal_duplex_arg;
 138   /*threshold on duplex distance*/
 139   distance=args_info.duplex_distance_arg;
 140   /*threshold on minimal stem length*/
 141   half_stem=args_info.minimal_stem_length_arg;
 142   /*threshold on maximal stem length*/
 143   max_half_stem=args_info.maximal_stem_length_arg;
 144   /*threshold on minimal duplex box length*/
 145   min_s2=args_info.minimal_duplex_box_length_arg;
 146   /*threshold on maximal duplex box length*/
 147   max_s2=args_info.maximal_duplex_box_length_arg;
 148   /*threshold on minimal snoRNA stem loop length*/
 149   min_s1=args_info.minimal_snoRNA_stem_loop_length_arg;
 150   /*thrshold on maximal snoRNA stem loop length*/
 151   max_s1=args_info.maximal_snoRNA_stem_loop_length_arg;
 152   /*threshold on minimal snoRNA duplex length*/
 153   min_d1=args_info.minimal_snoRNA_duplex_length_arg;
 154   /*threshold on minimal snoRNA duplex length*/
 155   min_d2=args_info.minimal_snoRNA_duplex_length_arg;
 156   /*threshold on minimal duplex stem energy*/
 157   threshTE=args_info.minimal_duplex_stem_energy_arg;
 158   /*threshold on minimal total energy*/
 159   threshSE=args_info.minimal_total_energy_arg;
 160   /*threshold on maximal stem asymmetry*/
 161   max_asymm=args_info.maximal_stem_asymmetry_arg;
 162   /*threshold on minimal lower stem energy*/
 163   threshD=args_info.minimal_lower_stem_energy_arg;
 164   /*threshold on minimal lower stem energy*/
 165   alignment_length=args_info.alignmentLength_arg;
 166
 167   threshloop=MIN2(threshloop,0);
 168
 169 /*   if(plfold_up_flag && !fast){ */
 170 /*     printf("Sorry, no accessibility implementation with the non fast implementation\n"); */
 171 /*     printf("If you want to include accessibility information please run RNAsnoop with the -f option\n"); */
 172 /*     printf("If you want to run RNAsnoop with the slow algorithm, please remove run RNAsnoop without -P\n"); */
 173 /*     exit(0); */
 174 /*   } */
 175
 176
 177   if(plfold_up_flag==2 && suffix==NULL){
 178     printf("RNAsnoop needs a suffix (-S u1-to-30) for the RNAup accessibility file in order to localize them\n");
 179     exit(0);
 180   }
 181   if(redraw){
 182     if(!alignment){
 183       redraw_output(result_file, output_directory, plfold_up_flag, suffix, access,NULL,NULL);
 184     }
 185     else{
 186       redraw_output(result_file, output_directory, plfold_up_flag, suffix, access,sname,tname);
 187     }
 188
 189     /* readfile */
 190     /* parse lines */
 191     /* use current ploting function to output the results */
 192     exit(0);
 193   }
 194   /* istty = isatty(fileno(stdout))&&isatty(fileno(stdin)); */
 195
 196   min_s2+=5;
 197   max_s2+=5;
 198   min_s1+=5;
 199   max_s1+=5;
 200   min_d1+=5;
 201   min_d2+=5;
 202
 203
 204   if(!alignment) {
 205     if(tname==NULL || sname==NULL) RNAsnoop_cmdline_parser_print_help();
 206     sno=fopen(sname, "r");
 207     if(sno==NULL){printf("%s: Wrong snoRNA file name\n", sname);return 0;}
 208     mrna=fopen(tname, "r");
 209     if(mrna==NULL){printf("%s: Wrong target file name\n", tname);return 0;}
 210     do {                                /* main loop: continue until end of file */
 211       if ((line_s = get_line(sno))==NULL) {
 212         free(line_s);
 213         break;
 214       }
 215
 216       /*   skip comment lines and get filenames */
 217       while ((*line_s=='*')||(*line_s=='\0')||(*line_s=='>')) {
 218         if (*line_s=='>')
 219           name_s = (char *) space(strlen(line_s)+1);
 220         (void) sscanf(line_s,"%s",name_s);
 221         free(line_s);
 222         if ((line_s = get_line(sno))==NULL) {
 223           free(line_s);
 224           break;
 225         }
 226       }
 227       if(name_s == NULL) {printf("Your snoRNA sequence: \n%s\nhas no header. Please update your fasta file\n", line_s); exit(0);}
 228       /*   if ((line ==NULL) || (strcmp(line, "@") == 0)) break; */
 229       temp_s = (char *) space(strlen(line_s)+1);
 230       (void) sscanf(line_s,"%s",temp_s);
 231       free(line_s);
 232       length_s = (int) strlen(temp_s);
 233       for (l = 0; l < length_s; l++) {
 234         temp_s[l] = toupper(temp_s[l]);
 235         if (!noconv && temp_s[l] == 'T') temp_s[l] = 'U';
 236       }
 237       string_s = (char*) space(length_s + 11);
 238       strcpy(string_s, "NNNNN");
 239       strcat(string_s+5, temp_s);
 240       strcat(string_s+5+length_s, "NNNNN\0");
 241       free(temp_s);
 242       /* We declare the structure variable here as it will also contains the final stem structure */
 243       structure=(char *)space((unsigned) length_s+11);
 244       if(fold_constrained){
 245         cstruc = get_line(sno);
 246         if(cstruc!=NULL){
 247           int dn3=strlen(cstruc)-(length_s-10);
 248           strcpy(structure,".....");
 249           strcat(structure, cstruc);
 250           if(dn3>=0){
 251             strcat(structure,".....\0");
 252           }
 253           else{
 254             while(dn3++){
 255               strcat(structure,".");
 256             }
 257             strcat(structure,"\0");
 258           }
 259           /* Now we fold with constraints the  */
 260         }
 261       }
 262       fullStemEnergy = snofold(string_s, structure, max_asymm, threshloop, min_s2, max_s2, half_stem, max_half_stem);
 263       do{                                /* main loop for target continue until end of file */
 264         snoopT mfe;
 265         if ((line_t = get_line(mrna))==NULL) {
 266           /* free(line_t); */
 267           break;
 268         }
 269         char *name_t=NULL;
 270
 271         /* skip comment lines and get filenames */
 272         while ((*line_t=='*')||(*line_t=='\0')||(*line_t=='>')) {
 273           if (*line_t=='>'){
 274             printf("%s\n", name_s);
 275             name_t = (char*) space(strlen(line_t)+1);
 276             (void) sscanf(line_t,"%s",name_t);
 277
 278             printf("%s\n", name_t);
 279             free(line_t);
 280             /*             free(string_t); */
 281               }
 282           if ((line_t = get_line(mrna))==NULL) {
 283             free(line_t);
 284             break;
 285           }
 286         }
 287         if(name_t == NULL) {printf("Your target sequence: \n%s\nhas no header. Please update your fasta file\n", line_t); exit(0);}
 288         /*   if ((line ==NULL) || (strcmp(line, "@") == 0)) break; */
 289         temp_t = (char *) space(strlen(line_t)+1);
 290         (void) sscanf(line_t,"%s",temp_t);
 291         int length_t;
 292         length_t = (int) strlen(temp_t);
 293         for (l = 0; l < length_t; l++) {
 294           temp_t[l] = toupper(temp_t[l]);
 295           if (!noconv && temp_t[l] == 'T') temp_t[l] = 'U';
 296         }
 297         string_t =(char *) space(length_t +11);
 298         strcpy(string_t, "NNNNN");
 299         strcat(string_t+5, temp_t);
 300         strcat(string_t+5+length_t, "NNNNN");
 301         free(temp_t);
 302         char *name_output;
 303         name_output=NULL;
 304         if(nice){
 305           name_output = (char*) space(sizeof(char)*(length_t+length_s+2));
 306           strcpy(name_output, name_t+1);
 307           strcat(name_output, "_");
 308           strcat(name_output, name_s+1);
 309           name_output[length_t + length_s +1]='\0';
 310         }
 311         if(delta>=0){
 312           snoopT* subopt;
 313           snoopT* sub;
 314           if(!fast && !plfold_up_flag){
 315             subopt = snoop_subopt(string_t, string_s, delta, 5, penalty, threshloop,
 316                                   threshLE,threshRE, threshDE,threshTE,threshSE,threshD,distance,
 317                                   half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2,fullStemEnergy);
 318             if(subopt==NULL){printf("no target found under the given constraints\n");free(subopt); free(string_t);free(line_t);continue;}
 319             int count=0;
 320             for (sub=subopt; sub->structure !=NULL; sub++) {
 321               print_struc(sub, string_t, string_s, name_s, name_t, count++,nice);
 322               free(sub->structure);
 323             }
 324             free(subopt);
 325             free(string_t);
 326             free(line_t);
 327           }
 328           else if(!plfold_up_flag){
 329             Lsnoop_subopt_list(string_t, string_s, delta, 5, penalty, threshloop,
 330                                threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance,
 331                                half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length,name_output,fullStemEnergy);
 332             free(string_t);
 333             free(line_t);
 334           }
 335           else {
 336             if(plfold_up_flag==1){
 337               int **access_s1;
 338               char *file_s1;
 339               int s1_len;/* k;*/ /*,j; */
 340               s1_len=strlen(string_t);
 341               file_s1 = (char *) space(sizeof(char) * (strlen(name_t+1)+strlen(access)+9));
 342               strcpy(file_s1, access);
 343               strcat(file_s1, "/");
 344               strcat(file_s1, name_t+1);
 345               strcat(file_s1, "_openen");
 346               access_s1 = read_plfold_i(file_s1, 1, s1_len);
 347               if(fast) {
 348                 Lsnoop_subopt_list_XS(string_t, string_s, (const int**) access_s1, delta, 5, penalty, threshloop,
 349                                       threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance,
 350                                       half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length,name_output, fullStemEnergy);
 351               }
 352               else{
 353                 snoop_subopt_XS(string_t, string_s, (const int **) access_s1, delta, 5, penalty, threshloop,
 354                                 threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance,
 355                                 half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length,name_output,fullStemEnergy);
 356               }
 357 /*                for(j=0; j<s1_len; j++){  */
 358 /*                  for(k=0; k<access_s1[0][0];k++){  */
 359 /*                    printf("%d \t",access_s1[k][j]);  */
 360 /*                  }  */
 361 /*                  printf("\n");  */
 362 /*                }  */
 363               free(file_s1);free(string_t);free(line_t);
 364               int i = access_s1[0][0];
 365               while(--i>-1){
 366                 free(access_s1[i]);
 367               }
 368               free(access_s1);
 369             }
 370             else if(plfold_up_flag==2){
 371               int **access_s1;
 372               char *file_s1;
 373               int s1_len,k;/* ,j; */
 374               s1_len=strlen(string_t);
 375               file_s1 = (char *) space(sizeof(char) * (strlen(name_t+1)+strlen(suffix) + strlen(access)+3));
 376               strcpy(file_s1, access);
 377               strcat(file_s1, "/");
 378               strcat(file_s1, name_t+1);
 379               strcat(file_s1, "_");
 380               strcat(file_s1, suffix);
 381               access_s1 = read_rnaup(file_s1, 1, s1_len);
 382               if(fast){
 383                 Lsnoop_subopt_list_XS(string_t, string_s, (const int**) access_s1, delta, 5, penalty, threshloop,
 384                                       threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance,
 385                                       half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length, name_output, fullStemEnergy);
 386               }
 387               else{
 388                 snoop_subopt_XS(string_t, string_s, (const int**) access_s1, delta, 5, penalty, threshloop,
 389                                 threshLE,threshRE,threshDE,threshTE,threshSE,threshD,distance,
 390                                 half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, alignment_length,name_output, fullStemEnergy);
 391               }
 392
 393
 394 /*                for(j=0; j<s1_len; j++){  */
 395 /*                  for(k=0; k<access_s1[0][0];k++){  */
 396 /*                    printf("%d \t",access_s1[k][j]);  */
 397 /*                  }  */
 398 /*                  printf("\n");  */
 399 /*                }  */
 400               free(file_s1);free(string_t);free(line_t);
 401               k = access_s1[0][0];
 402               while(--k>-1){
 403                 free(access_s1[k]);
 404               }
 405               free(access_s1);
 406             }
 407           }
 408         }
 409         else{
 410           mfe=snoopfold(string_t, string_s, penalty, threshloop, threshLE, threshRE,threshDE,  threshD,
 411                         half_stem, max_half_stem, min_s2, max_s2, min_s1, max_s1, min_d1, min_d2, fullStemEnergy);
 412           if(mfe.energy < INF){
 413             print_struc(&mfe, string_t, string_s, name_s, name_t,0,1);
 414             free(mfe.structure);
 415           }
 416           free(line_t);
 417           length_t= (int) strlen(string_t);
 418           free(string_t);
 419         }
 420         free(name_t);
 421         if(nice){
 422           free(name_output);
 423         }
 424       } while(1);
 425       rewind(mrna);
 426       snofree_arrays(strlen(string_s));  /* free's base_pair */
 427       free(string_s);string_s=NULL;
 428       free(name_s);name_s=NULL;
 429     } while (1);
 430   }
 431   else{
 432     if(tname==NULL || sname==NULL) RNAsnoop_cmdline_parser_print_help();
 433     sno=fopen(sname, "r");
 434     if(sno==NULL){printf("%s: Wrong snoRNA file name\n", sname);return 0;}
 435     mrna=fopen(tname, "r");
 436     if(mrna==NULL){printf("%s: Wrong target file name\n", tname);return 0;}
 437     char *temp1[MAX_NUM_NAMES], *temp2[MAX_NUM_NAMES],*AS1[MAX_NUM_NAMES], *AS2[MAX_NUM_NAMES], *names1[MAX_NUM_NAMES], *names2[MAX_NUM_NAMES];
 438     int n_seq, n_seq2;
 439     n_seq =  read_clustal(mrna, temp1, names1);
 440     n_seq2 = read_clustal(sno, temp2, names2);
 441     if (n_seq != n_seq2){
 442       for (i=0; temp1[i]; i++) {
 443         free(temp1[i]);
 444         free(temp2[i]);
 445       }
 446       nrerror("unequal number of seqs in alignments");
 447     }
 448     for(i=0;temp1[i];i++){
 449       AS1[i] = (char*) space((strlen(temp1[i])+11)*sizeof(char));
 450       AS2[i] = (char*) space((strlen(temp2[i])+11)*sizeof(char));
 451       strcpy(AS1[i],"NNNNN");
 452       strcpy(AS2[i],"NNNNN");
 453       strcat(AS1[i],temp1[i]);
 454       strcat(AS2[i],temp2[i]);
 455       strcat(AS1[i],"NNNNN\0");
 456       strcat(AS2[i],"NNNNN\0");
 457     }
 458     for (i=0; temp1[i]; i++) {
 459       free(temp1[i]);
 460       free(temp2[i]);
 461     }
 462     AS1[n_seq]=NULL;
 463     AS2[n_seq]=NULL;
 464     update_fold_params();
 465     alisnofold((const char **)AS2, max_asymm, threshloop, min_s2, max_s2, half_stem, max_half_stem);
 466      snoopT struc;
 467      struc = alisnoopfold((const char**) AS1, (const char**) AS2,
 468                   penalty, threshloop,
 469                           threshLE,threshRE,threshDE,threshD,
 470                   half_stem, max_half_stem,
 471                   min_s2, max_s2, min_s1,
 472                           max_s1, min_d1, min_d2);
 473      if(!(struc.structure==NULL)){
 474        aliprint_struc(&struc, (const char**) AS1,(const char**) AS2, names1, names2,0,nice);
 475        free(struc.structure);
 476      }
 477
 478
 479
 480     snoopT* subopt = alisnoop_subopt((const char**) AS1, (const char**) AS2, delta, 5, penalty, threshloop,
 481                                       threshLE,threshRE, threshDE, threshD,
 482                                       threshTE,threshSE,distance,
 483                                       half_stem, max_half_stem, min_s2, max_s2,
 484                                       min_s1, max_s1, min_d1, min_d2);
 485
 486
 487     snoopT *sub;
 488     if(subopt==NULL){
 489       printf("no target found under the given constraints\n");
 490     }else{
 491       int count=1;
 492       for (sub=subopt; !(sub->structure ==NULL); sub++) {
 493         aliprint_struc(sub, (const char**) AS1,(const char**) AS2,names1,names2,count,nice);
 494         free(sub->structure);
 495         count++;
 496       }
 497     }
 498
 499      alisnofree_arrays(strlen(AS2[0]));
 500      free(subopt);
 501     for (i=0; temp1[i]; i++) {
 502       free(AS1[i]);
 503       free(AS2[i]);
 504       free(names1[i]);
 505       free(names2[i]);
 506     }
 507
 508   }
 509   fclose(sno);
 510   fclose(mrna);
 511   return 0;
 512 }
 513
 514 static void print_struc(snoopT  *dup, const char *s1, const char *s2, char *name_s, char *name_t, int count, int nice) {
 515   int l1;
 516   l1 = strchr(dup->structure, '&')-dup->structure;
 517   char *target_struct;
 518   int shift=0, n2;
 519   char* psoutput;
 520   psoutput = (char*) space(100*sizeof(char));
 521 /*   if(dup->i > strlen(s1)-10){ */
 522 /*         dup->i--; */
 523 /*         l1--; */
 524 /*   } */
 525 /*   if(dup->i-l1< 0){ */
 526 /*         l1--; */
 527 /*         shift++; */
 528 /*   } */
 529   target_struct = (char*) space(sizeof(char) * (strlen(dup->structure)+1));
 530   strncpy(target_struct, dup->structure+shift, l1);
 531   strncat(target_struct, dup->structure + (strchr(dup->structure, '&')-dup->structure), strlen(dup->structure) - (strchr(dup->structure, '&')-dup->structure));
 532   strcat(target_struct,"\0");
 533   char *target;
 534   target = (char *) space(l1+1);
 535   strncpy(target, (s1+dup->i-l1+5), l1);
 536   target[l1]='\0';
 537   char *s4;
 538   n2 = strlen(s2);
 539   s4 = (char*) space(sizeof(char) *(n2-9));
 540   strncpy(s4, s2+5, n2-10);
 541   s4[n2-10]='\0';
 542   printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + 4.1 ) (%5.2f) \n%s&%s\n",
 543         target_struct, dup->i+1-l1,
 544          dup->i, dup->u, dup->j+1, dup->j + (strrchr(dup->structure,'>') - strchr(dup->structure, '>'))+1,
 545          (dup->Loop_D + dup->Duplex_El + dup->Duplex_Er + dup->Loop_E) + 4.10,
 546          dup->Duplex_El, dup->Duplex_Er, dup->Loop_E, dup->Loop_D,dup->fullStemEnergy,target,s4);
 547   if(nice){
 548     char *temp_seq;
 549     char *temp_struc;
 550     temp_seq = (char*) space(sizeof(char)*(l1+n2-9));
 551     temp_struc = (char*) space(sizeof(char)*(l1+n2-9));
 552     strcpy(temp_seq, target);
 553     strcat(temp_seq, s4);
 554     strncpy(temp_struc, target_struct, l1);
 555     strcat(temp_struc, target_struct+l1+1);
 556     temp_seq[n2+l1-10]='\0';
 557     temp_struc[n2+l1-10]='\0';
 558     cut_point = l1+1;
 559     char str[16];char upos[16];
 560     char *temp;
 561     int length_name = strlen(name_t)+strlen(name_s);
 562     temp = (char *) space(sizeof(char)*(length_name+1));
 563     strcpy(temp,name_t+1);
 564     strcat(temp,"_");
 565     strcat(temp,name_s + 1);
 566     temp[length_name] = '\0';
 567     strcpy(psoutput,"sno_");
 568     sprintf(str,"%d",count);
 569     strcat(psoutput,str);
 570     sprintf(upos,"%d",dup->u);
 571     strcat(psoutput,"_u_");
 572     strcat(psoutput,upos);
 573     strcat(psoutput,"_");
 574     strcat(psoutput,temp);
 575     strcat(psoutput,".ps\0");
 576     PS_rna_plot_snoop_a(temp_seq, temp_struc, psoutput, NULL, NULL);
 577     cut_point = -1;
 578     free(temp_seq);free(temp_struc);free(psoutput);free(temp);
 579   }
 580 free(s4);
 581 free(target_struct);
 582 free(target);
 583 }
 584
 585
 586 static void aliprint_struc(snoopT  *dup, const char **s1, const char **s2, char **name_t, char **name_s,int count,int nice) {
 587   int l1;
 588   l1 = strchr(dup->structure, '&')-dup->structure;
 589   char *target_struct;
 590   int shift=0;
 591 /*   if(dup->i > strlen(s1[0])-10){ */
 592 /*     dup->i--; */
 593 /*     l1--; */
 594 /*   } */
 595 /*   if(dup->i-l1< 0){ */
 596 /*     l1--; */
 597 /*     shift++; */
 598 /*   } */
 599   int length_struct = strlen(dup->structure);
 600   target_struct = (char*) space(sizeof(char) * (length_struct+1));
 601   strncpy(target_struct, dup->structure+shift, l1);
 602   strncat(target_struct, dup->structure + (strchr(dup->structure, '&')-dup->structure), length_struct - (strchr(dup->structure, '&')-dup->structure));
 603   strcat(target_struct,"\0");
 604   /* get the corresponding alignment slice */
 605   int n_seq,s;
 606   for(s=0; s1[s]!=NULL; s++);
 607   n_seq=s;
 608   int n1,n2;
 609   n1=strlen(s1[0]);n2=strlen(s2[0]);
 610   char **target;
 611   target = (char**) space((n_seq+1)*sizeof(char*));
 612   for(s=0;s<n_seq;s++){
 613     target[s] = (char *) space((l1 + n2-8)*sizeof(char));
 614     strncpy(target[s], (s1[s]+dup->i-l1+5), l1);
 615     strcat(target[s],"&");
 616     strncat(target[s], (s2[s]+5), n2-10);
 617     target[s][l1+n2-9]='\0';
 618   }
 619   char * consens;
 620   consens =  consens_mis((const char**)target);
 621   consens[l1]='&';
 622
 623   printf("%s %3d,%-3d;%3d : %3d,%-3d (%5.2f = %5.2f + %5.2f + %5.2f + %5.2f + 4.1; duplex cov = %5.2f; stem cov = %5.2f )\n%s\n",
 624           dup->structure, dup->i+1-l1,
 625           dup->i, dup->u, dup->j+1, dup->j + (strrchr(dup->structure,'>') - strchr(dup->structure, '>'))+1,
 626           (dup->Loop_D + dup->Duplex_El + dup->Duplex_Er + dup->Loop_E)/n_seq + 4.10,
 627           dup->Duplex_El/n_seq, dup->Duplex_Er/n_seq, dup->Loop_E/n_seq, dup->Loop_D/n_seq,dup->pscd/n_seq, dup->psct/n_seq,consens );
 628   if(nice){
 629     char* psoutput;
 630     psoutput = (char*) space(100*sizeof(char));
 631
 632     char *temp_seq, *temp_struct, **temp_target;
 633     temp_seq    = (char*) space((l1 + n2 -9)*sizeof(char));
 634     temp_struct = (char*) space((l1 + n2 -9)*sizeof(char));
 635     temp_target = (char**) space((n_seq+1)*sizeof(char*));
 636     for(s=0;s<n_seq;s++){
 637       temp_target[s] = (char *) space((l1 + n2-9)*sizeof(char));
 638       strncpy(temp_target[s], (s1[s]+dup->i-l1+5), l1);
 639       strncat(temp_target[s], (s2[s]+5), n2-10);
 640       temp_target[s][l1+n2-10]='\0';
 641     }
 642     strncpy(temp_seq, consens, l1);
 643     strncpy(temp_struct, target_struct, l1);
 644     strcat(temp_seq, consens+l1+1);
 645     strcat(temp_struct, target_struct+l1+1);
 646     temp_seq[n2-10+l1]='\0';
 647     temp_struct[n2-10+l1]='\0';
 648     char str[16];
 649     char upos[16];
 650     char *temp;
 651     int length_name = strlen(name_t[0]) + strlen(name_s[0])+1;
 652     temp = (char *) space(sizeof(char)*(length_name+1));
 653     strcpy(temp,name_t[0]);
 654     strcat(temp,"_");
 655     strcat(temp, name_s[0]);
 656     temp[length_name] = '\0';
 657     strcpy(psoutput,"snoaln_");
 658     sprintf(str,"%d",count);
 659     strcat(psoutput,str);
 660     sprintf(upos,"%d",dup->u);
 661     strcat(psoutput,"_u_");
 662     strcat(psoutput,upos);
 663     strcat(psoutput,"_");
 664     for(s=0; s<length_name; s++){
 665       if(temp[s]=='/') temp[s]='-';
 666     }
 667     strcat(psoutput,temp);
 668     strcat(psoutput,".ps\0");
 669     /*   psoutput[strlen(temp)+4+strlen(str)+39]='\0'; */
 670     cut_point = l1+1;
 671     aliPS_color_aln(dup->structure, psoutput, (const char**) target, (const char**) name_t);
 672     psoutput[1]='t';
 673     psoutput[2]='r';
 674     PS_rna_plot_snoop_a(temp_seq, temp_struct, psoutput, NULL, (const char**) temp_target);
 675     cut_point = -1;
 676     free(psoutput);
 677     for(s=0; s<n_seq; s++){
 678       free(temp_target[s]);
 679     }
 680     free(temp);
 681     free(temp_seq);
 682     free(temp_struct);
 683     free(temp_target);
 684   }
 685   for(s=0; s<n_seq; s++){
 686     free(target[s]);
 687   }
 688   free(consens);
 689   free(target_struct);
 690   free(target);
 691 }
 692
 693 static int get_max_u(const char *s, char delim){
 694   char * pch;
 695   int total;
 696   total=0;
 697   pch = strchr(s, delim);
 698   total++;
 699   while(pch!=NULL){
 700     pch=strchr(pch+1, delim);
 701     total++;
 702   }
 703   return total-2;
 704 }
 705
 706
 707
 708 static int ** read_rnaup(char *fname, const int beg, const int end)
 709 {
 710   FILE *in = fopen(fname,"r"); /*  open RNA_up file */
 711
 712   int i,j;
 713   int **access;
 714   int offset, temp;
 715   temp=0; offset=0;
 716   int seq_pos;
 717   int beg_r, end_r;
 718   beg_r=beg;
 719   end_r=end;
 720
 721   if(in==NULL){
 722     printf("%s", fname);
 723     perror("RNAup File open error here, Computing next target");
 724
 725     exit(EXIT_FAILURE);
 726   }
 727   char tmp[2048]={0x0};
 728   /* char *ptr; */
 729
 730
 731   if(fgets(tmp,sizeof(tmp),in)==0){
 732     perror("Empty File");
 733   }
 734   if(strchr(tmp,'>')){
 735     fprintf(stderr,"file %s is not in RNAup format\n",fname);
 736     exit(EXIT_FAILURE);
 737   }
 738
 739   while(!strstr(fgets(tmp,sizeof(tmp),in), "pos")){};
 740  /*  (void) fgets(tmp,sizeof(tmp),in); //Datum  */
 741 /*   (void) fgets(tmp,sizeof(tmp),in); //white line */
 742 /*   (void) fgets(tmp,sizeof(tmp),in); //RNAup  */
 743 /*   (void) fgets(tmp,sizeof(tmp),in); //sequence length */
 744 /*   (void) fgets(tmp,sizeof(tmp),in); //sequence */
 745
 746   int dim_x;
 747   dim_x = get_max_u(tmp,'S'); /*  get unpaired regions by conting tabs in second line */
 748   access = (int**) space(sizeof(int *) * (dim_x+2));
 749   for(i=0; i< dim_x+2; i++){
 750     access[i] =(int *) space(sizeof(int) * (end_r-beg_r+7));
 751   }
 752   for(i=0;i<end_r -  beg_r +6;i++){
 753     for(j=0;j<dim_x+2;j++){
 754       access[j][i]=INF;
 755     }
 756   }
 757   access[0][0]=dim_x+2;
 758   while(fgets(tmp,sizeof(tmp),in)!=0 && --end_r > -1) /* read a record */
 759   {
 760     float n;
 761     /* int i; */
 762     /* int u; */
 763     beg_r--;
 764     if(beg_r<1){
 765       if(sscanf(tmp,"%d%n",&seq_pos,&temp)==1){ /*  read sequenz pos = 1. spalte */
 766         offset+=temp;
 767         int count;
 768         count=1;
 769         while(sscanf(tmp+offset,"%f%n",&n,&temp)==1){ /* read columns */
 770           offset+=temp;
 771           access[count][-beg_r+5+1]=(int) rint(100*n); /* seq_pos+5 */
 772           /*           printf("%d %d %f\n", count, -beg_r, access[count][-beg_r]); */
 773           count++;
 774         }
 775         offset=0;
 776       }
 777     }
 778   }
 779   fclose(in);
 780   return access;
 781 }
 782
 783
 784
 785
 786
 787
 788
 789 static int ** read_plfold_i(char *fname, const int beg, const int end)
 790 {
 791   FILE *in=fopen(fname,"r");
 792   int i,j;
 793   int **access;
 794   int offset,temp;
 795   temp=0;offset=0;
 796   int seq_pos;
 797   int beg_r, end_r;
 798   beg_r=beg;
 799   end_r=end;
 800
 801   if(in==NULL){
 802     perror(" open error");
 803     exit(EXIT_FAILURE);
 804   }
 805
 806   char tmp[2048]={0x0};
 807   /* char *ptr; */
 808   if(fgets(tmp,sizeof(tmp),in)==0){
 809     perror("Empty File");
 810   }
 811   if(strchr(tmp,'>')){
 812     fprintf(stderr,"file %s is not in RNAplfold format",fname);
 813     exit(EXIT_FAILURE);
 814   }
 815   if(fgets(tmp,sizeof(tmp),in)==0){
 816     perror("No accessibility data");
 817   }
 818   int dim_x;
 819   dim_x=get_max_u(tmp,'\t');
 820   access = (int**) space(sizeof(int *) * (dim_x+2));
 821   for(i=0; i< dim_x+2; i++){
 822     access[i] =(int *) space(sizeof(int) * (end_r-beg_r+1));
 823   }
 824
 825   for(i=0;i<end_r -  beg_r +1;i++){
 826     for(j=0;j<dim_x+2;j++){
 827       access[j][i]=INF;
 828     }
 829   }
 830   access[0][0]=dim_x+2;
 831   while(fgets(tmp,sizeof(tmp),in)!=0 && --end_r > -1) /* read a record */
 832     {
 833       float n;
 834       /* int i; */
 835       /* int u; */
 836       beg_r--;
 837       if(beg_r<1){
 838         if(sscanf(tmp,"%d\t%n",&seq_pos,&temp)==1){
 839           offset+=temp;
 840           int count;
 841           count=1;
 842           while(sscanf(tmp+offset,"%f%n",&n,&temp)==1 ){
 843             offset+=temp;
 844             access[count][seq_pos+5]=(int)  rint( 100 * n); /* round the number */
 845             count++;
 846           }
 847           offset=0;
 848         }
 849       }
 850
 851     }
 852   fclose(in);
 853   return access;
 854 }
 855
 856
 857
 858 static void redraw_output(char *fname, char *output, int plfold_up_flag, char *suffix, char *access, char *sname, char *tname)
 859 {
 860   char *line;
 861   line=NULL;
 862   int count;
 863   short two_seq = 0;
 864   char *results;
 865   char *sequence;
 866   char *query;
 867   char *target;
 868   char *structure;
 869   char *pos;
 870   int posi;
 871   int begin=0, end=0, u=0;
 872   if(output==NULL){
 873     output = (char*) space(sizeof(char)*2);
 874     strcpy(output,".\0");
 875   }
 876   count=0;
 877   if(sname==NULL && tname==NULL){
 878     while((line=get_line(stdin))!=NULL) {
 879       count++;
 880       if(two_seq==0 && *line =='>'){
 881         query=(char*) space(strlen(line)+1);
 882         (void) sscanf(line,"%s",query);
 883         free(line);
 884         line=NULL;
 885         memmove(query, query+1, strlen(query));
 886         two_seq++;
 887       }
 888       else if(two_seq==1 && *line =='>'){
 889         target=(char*) space(strlen(line)+1);
 890         (void) sscanf(line,"%s",target);
 891         free(line);
 892         line=NULL;
 893         memmove(target, target+1, strlen(target));
 894         two_seq++;
 895       }
 896       else if(two_seq == 2) {
 897         int *relative_access;
 898         relative_access=NULL;
 899         printf("%s %s\n", target,query);
 900         if(strchr(line, '(') && strchr(line,'&') && strchr(line, '(')&& strchr(line, ',') && strchr(line,':') && strchr(line, '-')){
 901           results=line;
 902           int length;
 903           pos = strchr(results, ' ');
 904           posi = (int)  (pos - results);
 905           length = posi;
 906           structure = (char *) space((length+1) * sizeof(char));
 907           sscanf(results,"%s",structure); /* parse structure */
 908           char *line2;
 909           if((line2=get_line(stdin))!=NULL){
 910             sequence=(char *) space((length+1)* sizeof(char));
 911             sscanf(line2,"%s",sequence);
 912             if(line2!=NULL){
 913               free(line2);
 914             }
 915           }
 916           else{
 917             printf("your file is fucked up");exit(0);
 918           } /* parse sequence */
 919
 920           sscanf(pos, "%10d,%10d;%10d", &begin,&end,&u); /* parse coordinates */
 921           if(plfold_up_flag){
 922             int **access_s1;
 923             char *file_s1;
 924             /* read_rnaup_file */
 925             file_s1 = (char*) space(sizeof(char) * (strlen(target)+strlen(suffix)+strlen(access)+3));
 926             strcpy(file_s1, access);
 927             strcat(file_s1, "/");
 928             strcat(file_s1, target);
 929             strcat(file_s1, "_");
 930             strcat(file_s1, suffix);
 931             access_s1 = read_rnaup(file_s1, begin, end);
 932             free(file_s1);
 933             relative_access = (int*) space(sizeof(int)*(end-begin+2));
 934             relative_access[0]=access_s1[1][1+5];
 935             int i;
 936             for(i=2;i<(end-begin+2);i++){
 937               relative_access[i-1]=access_s1[i+1][i+5] - access_s1[i][i+4];
 938             }
 939             int l = access_s1[0][0];
 940             while(--l>-1){
 941               free(access_s1[l]);
 942             }
 943             free(access_s1);
 944           }
 945           char *catseq, *catstruct,*output_file;
 946           catseq = (char*) space(strlen(sequence)  *sizeof(char));
 947           catstruct = (char*) space(strlen(structure)  *sizeof(char));
 948           int l1 = strchr(structure, '&')-structure;
 949           strncpy(catseq, sequence, l1);
 950           strcat(catseq, sequence+l1+1);
 951           strncpy(catstruct, structure, l1);
 952           strcat(catstruct, structure+l1+1);
 953           strcat(catseq,"\0");
 954           strcat(catstruct,"\0");
 955
 956           /* printf("%s\n%s\n%s\n%s", catseq,sequence,catstruct,structure); */
 957           cut_point=l1+1;
 958           output_file = (char*) space((strlen(output) + strlen(query)+strlen(target)+50)*sizeof(char));
 959           strcpy(output_file,output);
 960           strcat(output_file, "/");
 961           strcat(output_file,"sno_");
 962           strcat(output_file,query);
 963           strcat(output_file,"_");
 964           strcat(output_file, target);
 965           strcat(output_file,"_u_");
 966           char str[9];
 967           sprintf(str,"%d",u);
 968           strcat(output_file,str);
 969           strcat(output_file,"_");
 970           sprintf(str,"%d",count);
 971           strcat(output_file,str);
 972           strcat(output_file,".ps\0");
 973           PS_rna_plot_snoop_a(catseq, catstruct, output_file,relative_access,NULL);
 974           if(relative_access){
 975             free(relative_access);
 976           }
 977           free(output_file);output_file=NULL;
 978           free(catseq);catseq=NULL;
 979           free(catstruct);catstruct=NULL;
 980           free(structure);
 981           free(sequence);
 982           free(line);
 983           line=NULL;
 984         }
 985         else if(*line=='>'){
 986           free(query);free(target);
 987           two_seq=1;
 988           query = (char*) space(sizeof(char)*(strlen(line)+1));
 989           (void) sscanf(line,"%s",query);
 990           free(line);
 991           memmove(query, query+1, strlen(query));
 992         }
 993       }
 994     }
 995     free(target);
 996     free(query);
 997   }
 998   else if(!(tname==NULL && sname==NULL)) {
 999     FILE* sno, *mrna;
1000     int i;
1001     sno=fopen(sname, "r");
1002     if(sno==NULL){printf("%s: Wrong snoRNA file name\n", sname);}
1003     mrna=fopen(tname, "r");
1004     if(mrna==NULL){printf("%s: Wrong target file name\n", tname);}
1005     char *temp1[MAX_NUM_NAMES], *temp2[MAX_NUM_NAMES], *AS1[MAX_NUM_NAMES], *AS2[MAX_NUM_NAMES], *names1[MAX_NUM_NAMES], *names2[MAX_NUM_NAMES];
1006     int n_seq, n_seq2;
1007     n_seq =  read_clustal(mrna, temp1, names1);
1008     n_seq2 = read_clustal(sno, temp2, names2);
1009     if (n_seq != n_seq2){
1010       for (i=0; temp1[i]; i++) {
1011         free(temp1[i]);
1012         free(temp2[i]);
1013       }
1014       nrerror("unequal number of seqs in alignments");
1015     }
1016     for(i=0;temp1[i];i++){
1017       AS1[i] = (char*) space((strlen(temp1[i])+11)*sizeof(char));
1018       AS2[i] = (char*) space((strlen(temp2[i])+11)*sizeof(char));
1019       strcpy(AS1[i],"NNNNN");
1020       strcpy(AS2[i],"NNNNN");
1021       strcat(AS1[i],temp1[i]);
1022       strcat(AS2[i],temp2[i]);
1023       strcat(AS1[i],"NNNNN\0");
1024       strcat(AS2[i],"NNNNN\0");
1025     }
1026     for (i=0; temp1[i]; i++) {
1027       free(temp1[i]);
1028       free(temp2[i]);
1029     }
1030     AS1[n_seq]=NULL;
1031     AS2[n_seq]=NULL;
1032     int count=0;
1033     while((line=get_line(stdin))!=NULL) {
1034       results=line;
1035       if(strchr(line, '(') && strchr(line,'&') && strchr(line, '(')&& strchr(line, ',') && strchr(line,':') && strchr(line, '-')) {
1036         count++;
1037         int length;
1038         /* int sbegin, send; */
1039         /* int energy; */
1040         snoopT dup;
1041         pos = strchr(results, ' ');
1042         posi = (int)  (pos - results);
1043         length = posi;
1044         dup.structure = (char *) space((length+1) * sizeof(char));
1045         sscanf(results,"%s %10d,%10d;%10d : %10d,%10d (%10f = %10f + %10f + %10f + %10f + 4.1; duplex cov = %10f; stem cov = %10f",
1046                dup.structure,
1047                &begin,
1048                &(dup.i),
1049                &(dup.u),
1050                &(dup.j),
1051                &end,
1052                &(dup.energy),
1053                &(dup.Duplex_El),
1054                &(dup.Duplex_Er),
1055                &(dup.Loop_E),
1056                &(dup.Loop_D),
1057                &(dup.pscd),
1058                &(dup.psct)); /* parse duplex stuff; */
1059         dup.energy*=n_seq;
1060         dup.Duplex_El*=n_seq;
1061         dup.Duplex_Er*=n_seq;
1062         dup.Loop_E*=n_seq;
1063         dup.Loop_D*=n_seq;
1064         aliprint_struc(&dup, (const char**) AS1, (const char**)AS2, names1, names2,count,1);
1065         free(dup.structure);
1066         free(line);
1067       }
1068       else{
1069         free(line);
1070       }
1071     }
1072
1073     for (i=0; AS1[i]; i++) {
1074       free(AS1[i]);
1075       free(AS2[i]);
1076         free(names1[i]);
1077       free(names2[i]);
1078     }
1079     fclose(mrna);
1080     fclose(sno);
1081   }
1082   /*   free(output); */
1083 }