2 Distances of Secondary Structure Ensembles
3 Peter F Stadler, Ivo L Hofacker, Sebastian Bonhoeffer
13 #include "part_func.h"
15 #include "fold_vars.h"
16 #include "dist_vars.h"
18 #include "profiledist.h"
19 #include "ProfileAln.h"
20 #include "read_epars.h"
22 #include "RNApaln_cmdl.h"
24 #define MAXLENGTH 10000
27 static char rcsid[] = "$Id: RNApaln.c,v 1.3 2005/07/24 08:35:15 ivo Exp $";
29 static double gapo=1.5, gape=0.666, seqw=0.5;
32 PRIVATE void command_line(int argc, char *argv[]);
33 PRIVATE void print_aligned_lines(FILE *somewhere);
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;
41 int main(int argc, char *argv[])
47 int type, length, taxa_list=0;
51 char *line=NULL, fname[FILENAME_MAX_LENGTH], *list_title=NULL;
52 plist *pr_pl, *mfe_pl;
55 command_line(argc, argv);
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)));
64 if ((istty)&&(n==0)) {
65 printf("\nInput sequence; @ to quit\n");
66 printf("%s\n", ruler);
70 do { /* get sequence to fold */
71 if (line!=NULL) free(line);
73 if ((line=get_line(stdin))==NULL) {type = 999; break;}
74 if (line[0]=='@') type = 999;
77 if (task=='m') taxa_list=1;
81 list_title = strdup(line);
86 if (sscanf(line,">%" XSTR(FILENAME_ID_LENGTH) "s", fname)!=0)
87 strcat(fname, "_dp.ps");
89 printf("%d : %s\n", n+1, line+1);
90 else printf("%s\n",line);
93 if (isalpha(line[0])) {
101 if( (task == 'm')&&(type>800) ) {
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);
113 if (type==888) { /* do another distance matrix */
115 printf("%s\n", list_title);
123 if (type == 888) continue;
124 if (outfile[0]!='\0') (void) fclose(somewhere);
125 if (line!= NULL) free(line);
126 return 0; /* finito */
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';
135 pr_pl = mfe_pl = NULL;
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);
151 sprintf(fname, "%d_dp.ps", n+1);
153 /* PS_dot_plot(line, fname); <- NOT THREADSAFE and obsolete function! */
155 /* call threadsafe dot plot printing function */
156 PS_dot_plot_list(line, fname, pr_pl, mfe_pl, "");
158 T[n] = Make_bp_profile_bppm(pr, length);
159 seq[n] = strdup(line);
160 if((istty)&&(task=='m')) printf("%s\n",structure);
171 dist = profile_aln(T[0],seq[0],T[1],seq[1]);
173 print_aligned_lines(somewhere);
176 free(seq[0]); free(seq[1]);
182 dist = profile_aln(T[1], seq[1], T[0], seq[0]);
184 print_aligned_lines(somewhere);
185 free_profile(T[1]); free(seq[1]);
191 dist = profile_aln(T[1], seq[1], T[0],seq[0]);
193 print_aligned_lines(somewhere);
194 free_profile(T[0]); free(seq[0]);
195 T[0] = T[1]; seq[0] = seq[1];
204 nrerror("This can't happen.");
205 } /* END switch task */
206 (void) fflush(stdout);
208 if (line !=NULL) free(line);
212 /* ----------------------------------------------------------------- */
214 PRIVATE void command_line(int argc, char *argv[])
218 char *ns_bases=NULL, *c;
219 char *ParamFile=NULL;
220 struct RNApaln_args_info args_info;
224 if(RNApaln_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
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);
244 if(args_info.mode_given){
245 if(strlen(args_info.mode_arg) != 1){
246 RNApaln_cmdline_parser_print_help(); exit(EXIT_FAILURE);
248 else task = *args_info.mode_arg;
250 case 'p': case 'm': case 'f': case 'c': break;
251 default: RNApaln_cmdline_parser_print_help(); exit(EXIT_FAILURE);
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';
259 else outfile[0] = '\0';
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;
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;
274 /* free allocated memory of command line data structure */
275 RNApaln_cmdline_parser_free (&args_info);
277 /* fprintf(stderr, "%f %f %f %d\n", gapo, gape, seqw, -endgaps); */
278 set_paln_params(gapo, gape, seqw, 1-endgaps);
281 read_parameter_file(ParamFile);
283 if (ns_bases!=NULL) {
284 nonstandards = space(33);
292 nonstandards[i++]=*c++;
293 nonstandards[i++]=*c;
294 if ((sym)&&(*c!=*(c-1))) {
295 nonstandards[i++]=*c;
296 nonstandards[i++]=*(c-1);
304 /*--------------------------------------------------------------------------*/
306 PRIVATE void print_aligned_lines(FILE *somewhere)
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]);
314 /*--------------------------------------------------------------------------*/