2 Distances of Secondary Structure Ensembles
3 Peter F Stadler, Ivo L Hofacker, Sebastian Bonhoeffer
13 #include "part_func.h"
14 #include "fold_vars.h"
15 #include "dist_vars.h"
18 #include "read_epars.h"
19 #include "profiledist.h"
20 #include "RNApdist_cmdl.h"
23 #define MAXLENGTH 10000
26 static char rcsid[] = "$Id: RNApdist.c,v 1.8 2002/11/07 12:19:41 ivo Exp $";
28 PRIVATE void command_line(int argc, char *argv[]);
29 PRIVATE void usage(void);
30 PRIVATE void print_aligned_lines(FILE *somewhere);
33 PRIVATE char outfile[FILENAME_MAX_LENGTH];
34 PRIVATE char ruler[] ="....,....1....,....2....,....3....,....4"
35 "....,....5....,....6....,....7....,....8";
36 static int noconv = 0;
38 int main(int argc, char *argv[])
43 int type, length, taxa_list=0;
47 char *line=NULL, fname[FILENAME_MAX_LENGTH], *list_title=NULL;
48 plist *pr_pl, *mfe_pl;
50 pr_pl = mfe_pl = NULL;
52 command_line(argc, argv);
54 if((outfile[0]=='\0')&&(task=='m')&&(edit_backtrack))
55 strcpy(outfile,"backtrack.file");
56 if (outfile[0]!='\0') somewhere = fopen(outfile,"w");
57 if (somewhere==NULL) somewhere = stdout;
58 istty = (isatty(fileno(stdout))&&isatty(fileno(stdin)));
61 if ((istty)&&(n==0)) {
62 printf("\nInput sequence; @ to quit\n");
63 printf("%s\n", ruler);
67 do { /* get sequence to fold */
68 if (line!=NULL) free(line);
70 if ((line=get_line(stdin))==NULL) {type = 999; break;}
71 if (line[0]=='@') type = 999;
74 if (task=='m') taxa_list=1;
78 list_title = strdup(line);
83 if (sscanf(line,">%" XSTR(FILENAME_ID_LENGTH) "s", fname)!=0)
84 strcat(fname, "_dp.ps");
86 printf("%d : %s\n", n+1, line+1);
87 else printf("%s\n",line);
90 if (isalpha(line[0])) {
98 if( (task == 'm')&&(type>800) ) {
100 printf("* END of taxa list\n");
101 printf("> p %d (pdist)\n",n);
102 for (i=1; i<n; i++) {
103 for (j=0; j<i; j++) {
104 printf("%g ",profile_edit_distance(T[i], T[j]));
105 if(edit_backtrack) fprintf(somewhere,"> %d %d\n",i+1,j+1);
106 print_aligned_lines(somewhere);
110 if (type==888) { /* do another distance matrix */
112 printf("%s\n", list_title);
120 if (type == 888) continue;
121 if (outfile[0]!='\0') (void) fclose(somewhere);
122 if (line!= NULL) free(line);
123 return 0; /* finito */
126 length = (int) strlen(line);
127 for (i=0; i<length; i++) {
128 line[i]=toupper(line[i]);
129 if (!noconv && line[i] == 'T') line[i] = 'U';
132 /* init_pf_fold(length); <- obsolete */
133 structure = (char *) space((length+1)*sizeof(char));
134 (void) pf_fold(line,structure);
137 sprintf(fname, "%d_dp.ps", n+1);
139 /* PS_dot_plot(line, fname); <- NOT THREADSAFE and obsolete function! */
141 /* get pairlist of probability matrix */
142 assign_plist_from_pr(&pr_pl, pr, length, 1e-5);
143 /* no previous mfe call thus no mfe structure information known */
144 mfe_pl = (plist *)space(sizeof(plist));
145 mfe_pl[0].i = mfe_pl[0].j = 0;
147 /* call threadsafe dot plot printing function */
148 PS_dot_plot_list(line, fname, pr_pl, mfe_pl, "");
150 T[n] = Make_bp_profile_bppm(pr, length);
151 if((istty)&&(task=='m')) printf("%s\n",structure);
161 dist = profile_edit_distance(T[0],T[1]);
163 print_aligned_lines(somewhere);
171 dist = profile_edit_distance(T[1], T[0]);
173 print_aligned_lines(somewhere);
180 dist = profile_edit_distance(T[1], T[0]);
182 print_aligned_lines(somewhere);
193 nrerror("This can't happen.");
194 } /* END switch task */
195 (void) fflush(stdout);
197 if (line !=NULL) free(line);
201 /* ----------------------------------------------------------------- */
203 PRIVATE void command_line(int argc, char *argv[])
207 char *ns_bases=NULL, *c;
208 char *ParamFile=NULL;
209 struct RNApdist_args_info args_info;
214 #############################################
215 # check the command line parameters
216 #############################################
218 if(RNApdist_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
221 if(args_info.temp_given)
222 temperature = args_info.temp_arg;
224 /* do not take special tetra loop energies into account */
225 if(args_info.noTetra_given)
228 /* set dangle model */
229 if(args_info.dangles_given){
230 dangles = args_info.dangles_arg;
231 if(dangles) dangles = 2;
234 /* set energy model */
235 if(args_info.energyModel_given)
236 energy_set = args_info.energyModel_arg;
238 /* do not allow weak pairs */
239 if(args_info.noLP_given)
242 /* do not allow wobble pairs (GU) */
243 if(args_info.noGU_given)
246 /* do not allow weak closing pairs (AU,GU) */
247 if(args_info.noClosingGU_given)
250 /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
251 if(args_info.noconv_given)
254 /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
255 if(args_info.nsp_given)
256 ns_bases = strdup(args_info.nsp_arg);
258 /* take another energy parameter set */
259 if(args_info.paramFile_given)
260 ParamFile = strdup(args_info.paramFile_arg);
262 if(args_info.compare_given){
263 switch(args_info.compare_arg[0]){
267 case 'c': task=args_info.compare_arg[0];
269 default: RNApdist_cmdline_parser_print_help();
274 if(args_info.backtrack_given){
275 if(strcmp(args_info.backtrack_arg, "none")){
276 strncpy(outfile, args_info.backtrack_arg, FILENAME_MAX_LENGTH-1);
281 /* free allocated memory of command line data structure */
282 RNApdist_cmdline_parser_free (&args_info);
284 /* do some preparations */
286 read_parameter_file(ParamFile);
288 if (ns_bases!=NULL) {
289 nonstandards = space(33);
297 nonstandards[i++]=*c++;
298 nonstandards[i++]=*c;
299 if ((sym)&&(*c!=*(c-1))) {
300 nonstandards[i++]=*c;
301 nonstandards[i++]=*(c-1);
310 /* ------------------------------------------------------------------------- */
312 PRIVATE void usage(void)
314 nrerror("usage: RNApdist [-Xpmfc] [-B [file]] [-T temp] [-4] [-d] [-noGU]\n"
315 " [-noCloseGU] [-noLP] [-e e_set] [-P paramfile] [-nsp pairs]");
318 /*--------------------------------------------------------------------------*/
320 PRIVATE void print_aligned_lines(FILE *somewhere)
323 fprintf(somewhere, "%s\n%s\n", aligned_line[0], aligned_line[1]);
326 /*--------------------------------------------------------------------------*/