1 /* Last changed Time-stamp: <2005-07-23 16:50:24 ivo> */
3 Compute duplex structure of two RNA strands
15 #include "fold_vars.h"
18 #include "read_epars.h"
21 #include "RNAaliduplex_cmdl.h"
24 static char rcsid[] = "$Id: RNAaliduplex.c,v 1.1 2007/08/26 10:08:44 ivo Exp $";
26 PRIVATE void print_struc(duplexT const *dup);
28 #define MAX_NUM_NAMES 500
30 int main(int argc, char *argv[]){
31 struct RNAaliduplex_args_info args_info;
32 char *AS1[MAX_NUM_NAMES], *AS2[MAX_NUM_NAMES], *names1[MAX_NUM_NAMES], *names2[MAX_NUM_NAMES];
33 char *ParamFile=NULL, *c, *ns_bases=NULL;
34 int i, r, n_seq, n_seq2;
35 int istty, delta=-1, sym;
38 FILE *file1=NULL, *file2=NULL; /* input alignments */
43 #############################################
44 # check the command line parameters
45 #############################################
47 if(RNAaliduplex_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
49 if(args_info.temp_given) temperature = args_info.temp_arg;
50 /* do not take special tetra loop energies into account */
51 if(args_info.noTetra_given) tetra_loop=0;
52 /* set dangle model */
53 if(args_info.dangles_given){
54 if((args_info.dangles_arg < 0) || (args_info.dangles_arg > 3))
55 warn_user("required dangle model not implemented, falling back to default dangles=2");
57 dangles = args_info.dangles_arg;
59 /* do not allow weak pairs */
60 if(args_info.noLP_given) noLonelyPairs = 1;
61 /* do not allow wobble pairs (GU) */
62 if(args_info.noGU_given) noGU = 1;
63 /* do not allow weak closing pairs (AU,GU) */
64 if(args_info.noClosingGU_given) no_closingGU = 1;
65 /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
66 if(args_info.noconv_given) noconv = 1;
67 /* take another energy parameter set */
68 if(args_info.paramFile_given) ParamFile = strdup(args_info.paramFile_arg);
69 /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
70 if(args_info.nsp_given) ns_bases = strdup(args_info.nsp_arg);
72 if(args_info.deltaEnergy_given) delta = (int) (0.1+args_info.deltaEnergy_arg*100);
74 if(args_info.sorted_given) subopt_sorted = 1;
76 /* check unnamed options a.k.a. filenames of input alignments */
77 if(args_info.inputs_num == 2){
78 file1 = fopen(args_info.inputs[0], "r");
80 fprintf(stderr, "can't open %s\n", args_info.inputs[0]);
82 file2 = fopen(args_info.inputs[1], "r");
84 fprintf(stderr, "can't open %s\n", args_info.inputs[1]);
88 RNAaliduplex_cmdline_parser_print_help();
92 if (!(file1 && file2)){
93 fprintf(stderr, "No input files");
94 RNAaliduplex_cmdline_parser_print_help();
98 n_seq = read_clustal(file1, AS1, names1);
99 n_seq2 = read_clustal(file2, AS2, names2);
102 if(n_seq != n_seq2) nrerror("unequal number of seqs in alignments");
105 /* free allocated memory of command line data structure */
106 RNAaliduplex_cmdline_parser_free (&args_info);
109 #############################################
111 #############################################
115 if (ParamFile != NULL)
116 read_parameter_file(ParamFile);
118 if (ns_bases != NULL) {
119 nonstandards = space(33);
127 nonstandards[i++]=*c++;
128 nonstandards[i++]=*c;
129 if ((sym)&&(*c!=*(c-1))) {
130 nonstandards[i++]=*c;
131 nonstandards[i++]=*(c-1);
138 istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
141 #############################################
143 #############################################
145 update_fold_params();
149 subopt = aliduplex_subopt((const char **)AS1, (const char **)AS2, delta, 5);
150 for (sub=subopt; sub->i >0; sub++) {
152 free(sub->structure);
156 mfe = aliduplexfold((const char **)AS1, (const char **)AS2);
160 for (i=0; AS1[i]; i++) {
161 free(AS1[i]); free(names1[i]);
162 free(AS2[i]); free(names2[i]);
167 PRIVATE void print_struc(duplexT const *dup) {
169 l1 = strchr(dup->structure, '&')-dup->structure;
170 printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n", dup->structure, dup->i+1-l1,
171 dup->i, dup->j, dup->j+strlen(dup->structure)-l1-2, dup->energy);