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"
17 #include "read_epars.h"
20 #include "RNAduplex_cmdl.h"
22 PRIVATE void print_struc(duplexT const *dup);
24 static char rcsid[] = "$Id: RNAduplex.c,v 1.5 2007/08/26 09:41:12 ivo Exp $";
26 /*--------------------------------------------------------------------------*/
28 int main(int argc, char *argv[]){
29 struct RNAduplex_args_info args_info;
30 char fname[FILENAME_MAX_LENGTH], *input_string, *s1, *s2, *orig_s1, *orig_s2, *c, *ParamFile=NULL, *ns_bases=NULL;
31 unsigned int input_type;
36 s1 = s2 = orig_s1 = orig_s2 = NULL;
41 #############################################
42 # check the command line parameters
43 #############################################
45 if(RNAduplex_cmdline_parser (argc, argv, &args_info) != 0) exit(1);
47 if(args_info.temp_given) temperature = args_info.temp_arg;
48 /* do not take special tetra loop energies into account */
49 if(args_info.noTetra_given) tetra_loop=0;
50 /* set dangle model */
51 if(args_info.dangles_given){
52 if((args_info.dangles_arg < 0) || (args_info.dangles_arg > 3))
53 warn_user("required dangle model not implemented, falling back to default dangles=2");
55 dangles = args_info.dangles_arg;
57 /* do not allow weak pairs */
58 if(args_info.noLP_given) noLonelyPairs = 1;
59 /* do not allow wobble pairs (GU) */
60 if(args_info.noGU_given) noGU = 1;
61 /* do not allow weak closing pairs (AU,GU) */
62 if(args_info.noClosingGU_given) no_closingGU = 1;
63 /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
64 if(args_info.noconv_given) noconv = 1;
65 /* take another energy parameter set */
66 if(args_info.paramFile_given) ParamFile = strdup(args_info.paramFile_arg);
67 /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
68 if(args_info.nsp_given) ns_bases = strdup(args_info.nsp_arg);
70 if(args_info.deltaEnergy_given) delta = (int) (0.1+args_info.deltaEnergy_arg*100);
72 if(args_info.sorted_given) subopt_sorted = 1;
74 /* free allocated memory of command line data structure */
75 RNAduplex_cmdline_parser_free (&args_info);
78 #############################################
80 #############################################
82 if (ParamFile != NULL)
83 read_parameter_file(ParamFile);
85 if (ns_bases != NULL) {
86 nonstandards = space(33);
94 nonstandards[i++]=*c++;
96 if ((sym)&&(*c!=*(c-1))) {
98 nonstandards[i++]=*(c-1);
104 istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
107 #############################################
108 # main loop: continue until end of file
109 #############################################
112 duplexT mfe, *subopt;
114 ########################################################
115 # handle user input from 'stdin'
116 ########################################################
118 if(istty) print_tty_input_seq_str("Input two sequences (one line each)");
120 /* extract filename from fasta header if available */
122 while((input_type = get_input_line(&input_string, 0)) == VRNA_INPUT_FASTA_HEADER){
123 printf(">%s\n", input_string);
124 (void) sscanf(input_string, "%" XSTR(FILENAME_ID_LENGTH) "s", fname);
128 /* break on any error, EOF or quit request */
129 if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)){ break;}
130 /* else assume a proper sequence of letters of a certain alphabet (RNA, DNA, etc.) */
132 s1 = strdup(input_string);
136 /* get second sequence */
137 while((input_type = get_input_line(&input_string, 0)) == VRNA_INPUT_FASTA_HEADER){
138 printf(">%s\n", input_string);
141 /* break on any error, EOF or quit request */
142 if(input_type & (VRNA_INPUT_QUIT | VRNA_INPUT_ERROR)){ break;}
143 /* else assume a proper sequence of letters of a certain alphabet (RNA, DNA, etc.) */
145 s2 = strdup(input_string);
149 /* convert DNA alphabet to RNA if not explicitely switched off */
154 /* store case-unmodified sequence */
155 orig_s1 = strdup(s1);
156 orig_s2 = strdup(s2);
157 /* convert sequence to uppercase letters only */
161 if (istty) printf("lengths = %d,%d\n", strlen(s1), strlen(s2));
164 ########################################################
165 # done with 'stdin' handling, now init everything properly
166 ########################################################
168 update_fold_params();
171 ########################################################
172 # begin actual computations
173 ########################################################
177 subopt = duplex_subopt(s1, s2, delta, 5);
178 for (sub=subopt; sub->i >0; sub++) {
180 free(sub->structure);
185 mfe = duplexfold(s1, s2);
189 (void) fflush(stdout);
194 s1 = s2 = orig_s1 = orig_s2 = NULL;
199 PRIVATE void print_struc(duplexT const *dup) {
201 l1 = strchr(dup->structure, '&')-dup->structure;
202 printf("%s %3d,%-3d : %3d,%-3d (%5.2f)\n", dup->structure, dup->i+1-l1,
203 dup->i, dup->j, dup->j+strlen(dup->structure)-l1-2, dup->energy);