1 /* Last changed Time-stamp: <2007-12-05 13:55:42 ronny> */
3 Ineractive Access to folding Routines
18 #include "energy_par.h"
20 #include "part_func.h"
21 #include "fold_vars.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*); */
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);
37 static char UNUSED rcsid[] = "$Id: RNAfold.c,v 1.23 2007/12/05 13:04:08 ivo Exp $";
39 #define PRIVATE static
40 #define MAX_NUM_NAMES 500
42 /*--------------------------------------------------------------------------*/
44 int main(int argc, char *argv[])
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*/;
52 char *access;char *result_file;
53 char *output_directory;
55 output_directory=NULL;
62 /* long int elapTicks; */
63 /* clock_t Begin, End; */
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*/,
87 redraw /*if used (option I) allow to redraw command line output into ps files */;
97 snoop_subopt_sorted=0;
99 #############################################
100 # check the command line parameters
101 #############################################
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;
167 threshloop=MIN2(threshloop,0);
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"); */
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");
183 redraw_output(result_file, output_directory, plfold_up_flag, suffix, access,NULL,NULL);
186 redraw_output(result_file, output_directory, plfold_up_flag, suffix, access,sname,tname);
191 /* use current ploting function to output the results */
194 /* istty = isatty(fileno(stdout))&&isatty(fileno(stdin)); */
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) {
216 /* skip comment lines and get filenames */
217 while ((*line_s=='*')||(*line_s=='\0')||(*line_s=='>')) {
219 name_s = (char *) space(strlen(line_s)+1);
220 (void) sscanf(line_s,"%s",name_s);
222 if ((line_s = get_line(sno))==NULL) {
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);
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';
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");
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);
247 int dn3=strlen(cstruc)-(length_s-10);
248 strcpy(structure,".....");
249 strcat(structure, cstruc);
251 strcat(structure,".....\0");
255 strcat(structure,".");
257 strcat(structure,"\0");
259 /* Now we fold with constraints the */
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 */
265 if ((line_t = get_line(mrna))==NULL) {
271 /* skip comment lines and get filenames */
272 while ((*line_t=='*')||(*line_t=='\0')||(*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);
278 printf("%s\n", name_t);
280 /* free(string_t); */
282 if ((line_t = get_line(mrna))==NULL) {
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);
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';
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");
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';
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;}
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);
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);
336 if(plfold_up_flag==1){
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);
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);
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);
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]); */
363 free(file_s1);free(string_t);free(line_t);
364 int i = access_s1[0][0];
370 else if(plfold_up_flag==2){
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);
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);
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);
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]); */
400 free(file_s1);free(string_t);free(line_t);
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);
417 length_t= (int) strlen(string_t);
426 snofree_arrays(strlen(string_s)); /* free's base_pair */
427 free(string_s);string_s=NULL;
428 free(name_s);name_s=NULL;
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];
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++) {
446 nrerror("unequal number of seqs in alignments");
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");
458 for (i=0; temp1[i]; i++) {
464 update_fold_params();
465 alisnofold((const char **)AS2, max_asymm, threshloop, min_s2, max_s2, half_stem, max_half_stem);
467 struc = alisnoopfold((const char**) AS1, (const char**) AS2,
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);
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);
489 printf("no target found under the given constraints\n");
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);
499 alisnofree_arrays(strlen(AS2[0]));
501 for (i=0; temp1[i]; i++) {
514 static void print_struc(snoopT *dup, const char *s1, const char *s2, char *name_s, char *name_t, int count, int nice) {
516 l1 = strchr(dup->structure, '&')-dup->structure;
520 psoutput = (char*) space(100*sizeof(char));
521 /* if(dup->i > strlen(s1)-10){ */
525 /* if(dup->i-l1< 0){ */
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");
534 target = (char *) space(l1+1);
535 strncpy(target, (s1+dup->i-l1+5), l1);
539 s4 = (char*) space(sizeof(char) *(n2-9));
540 strncpy(s4, s2+5, n2-10);
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);
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';
559 char str[16];char upos[16];
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);
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);
578 free(temp_seq);free(temp_struc);free(psoutput);free(temp);
586 static void aliprint_struc(snoopT *dup, const char **s1, const char **s2, char **name_t, char **name_s,int count,int nice) {
588 l1 = strchr(dup->structure, '&')-dup->structure;
591 /* if(dup->i > strlen(s1[0])-10){ */
595 /* if(dup->i-l1< 0){ */
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 */
606 for(s=0; s1[s]!=NULL; s++);
609 n1=strlen(s1[0]);n2=strlen(s2[0]);
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';
620 consens = consens_mis((const char**)target);
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 );
630 psoutput = (char*) space(100*sizeof(char));
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';
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';
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]);
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]='-';
667 strcat(psoutput,temp);
668 strcat(psoutput,".ps\0");
669 /* psoutput[strlen(temp)+4+strlen(str)+39]='\0'; */
671 aliPS_color_aln(dup->structure, psoutput, (const char**) target, (const char**) name_t);
674 PS_rna_plot_snoop_a(temp_seq, temp_struct, psoutput, NULL, (const char**) temp_target);
677 for(s=0; s<n_seq; s++){
678 free(temp_target[s]);
685 for(s=0; s<n_seq; s++){
693 static int get_max_u(const char *s, char delim){
697 pch = strchr(s, delim);
700 pch=strchr(pch+1, delim);
708 static int ** read_rnaup(char *fname, const int beg, const int end)
710 FILE *in = fopen(fname,"r"); /* open RNA_up file */
723 perror("RNAup File open error here, Computing next target");
727 char tmp[2048]={0x0};
731 if(fgets(tmp,sizeof(tmp),in)==0){
732 perror("Empty File");
735 fprintf(stderr,"file %s is not in RNAup format\n",fname);
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 */
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));
752 for(i=0;i<end_r - beg_r +6;i++){
753 for(j=0;j<dim_x+2;j++){
757 access[0][0]=dim_x+2;
758 while(fgets(tmp,sizeof(tmp),in)!=0 && --end_r > -1) /* read a record */
765 if(sscanf(tmp,"%d%n",&seq_pos,&temp)==1){ /* read sequenz pos = 1. spalte */
769 while(sscanf(tmp+offset,"%f%n",&n,&temp)==1){ /* read columns */
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]); */
789 static int ** read_plfold_i(char *fname, const int beg, const int end)
791 FILE *in=fopen(fname,"r");
802 perror(" open error");
806 char tmp[2048]={0x0};
808 if(fgets(tmp,sizeof(tmp),in)==0){
809 perror("Empty File");
812 fprintf(stderr,"file %s is not in RNAplfold format",fname);
815 if(fgets(tmp,sizeof(tmp),in)==0){
816 perror("No accessibility data");
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));
825 for(i=0;i<end_r - beg_r +1;i++){
826 for(j=0;j<dim_x+2;j++){
830 access[0][0]=dim_x+2;
831 while(fgets(tmp,sizeof(tmp),in)!=0 && --end_r > -1) /* read a record */
838 if(sscanf(tmp,"%d\t%n",&seq_pos,&temp)==1){
842 while(sscanf(tmp+offset,"%f%n",&n,&temp)==1 ){
844 access[count][seq_pos+5]=(int) rint( 100 * n); /* round the number */
858 static void redraw_output(char *fname, char *output, int plfold_up_flag, char *suffix, char *access, char *sname, char *tname)
871 int begin=0, end=0, u=0;
873 output = (char*) space(sizeof(char)*2);
874 strcpy(output,".\0");
877 if(sname==NULL && tname==NULL){
878 while((line=get_line(stdin))!=NULL) {
880 if(two_seq==0 && *line =='>'){
881 query=(char*) space(strlen(line)+1);
882 (void) sscanf(line,"%s",query);
885 memmove(query, query+1, strlen(query));
888 else if(two_seq==1 && *line =='>'){
889 target=(char*) space(strlen(line)+1);
890 (void) sscanf(line,"%s",target);
893 memmove(target, target+1, strlen(target));
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, '-')){
903 pos = strchr(results, ' ');
904 posi = (int) (pos - results);
906 structure = (char *) space((length+1) * sizeof(char));
907 sscanf(results,"%s",structure); /* parse structure */
909 if((line2=get_line(stdin))!=NULL){
910 sequence=(char *) space((length+1)* sizeof(char));
911 sscanf(line2,"%s",sequence);
917 printf("your file is fucked up");exit(0);
918 } /* parse sequence */
920 sscanf(pos, "%10d,%10d;%10d", &begin,&end,&u); /* parse coordinates */
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);
933 relative_access = (int*) space(sizeof(int)*(end-begin+2));
934 relative_access[0]=access_s1[1][1+5];
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];
939 int l = access_s1[0][0];
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);
954 strcat(catstruct,"\0");
956 /* printf("%s\n%s\n%s\n%s", catseq,sequence,catstruct,structure); */
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_");
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);
975 free(relative_access);
977 free(output_file);output_file=NULL;
978 free(catseq);catseq=NULL;
979 free(catstruct);catstruct=NULL;
986 free(query);free(target);
988 query = (char*) space(sizeof(char)*(strlen(line)+1));
989 (void) sscanf(line,"%s",query);
991 memmove(query, query+1, strlen(query));
998 else if(!(tname==NULL && sname==NULL)) {
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];
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++) {
1014 nrerror("unequal number of seqs in alignments");
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");
1026 for (i=0; temp1[i]; i++) {
1033 while((line=get_line(stdin))!=NULL) {
1035 if(strchr(line, '(') && strchr(line,'&') && strchr(line, '(')&& strchr(line, ',') && strchr(line,':') && strchr(line, '-')) {
1038 /* int sbegin, send; */
1041 pos = strchr(results, ' ');
1042 posi = (int) (pos - results);
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",
1058 &(dup.psct)); /* parse duplex stuff; */
1060 dup.Duplex_El*=n_seq;
1061 dup.Duplex_Er*=n_seq;
1064 aliprint_struc(&dup, (const char**) AS1, (const char**)AS2, names1, names2,count,1);
1065 free(dup.structure);
1073 for (i=0; AS1[i]; i++) {