1 /* emit_main.c
2 * main() for emitting sequences from a stored model
3 * written as a debugging aid
4 *
5 * 1.0: SRE, Tue Jun 15 09:32:43 1993
6 * 2.0: SRE, Fri Sep 10 08:06:37 1993
7 */
8
9 #include <stdlib.h>
10 #include <stdio.h>
11 #include <string.h>
12 #include <time.h>
13
14 #ifdef NEED_GETOPTH
15 #include <getopt.h>
16 #endif
17
18 #include "structs.h"
19 #include "funcs.h"
20 #include "version.h"
21
22 #ifdef MEMDEBUG
23 #include "dbmalloc.h"
24 #endif
25
26 #define OPTIONS "abls:L"
27
28 static char usage[] = "\
29 Usage: covee [-options] <cmfile>\n\
30 where options are:\n\
31 -a : annotate all pairs, not just canonical ones\n\
32 -b : emit single most probable sequence\n\
33 -l : print as mountain landscape\n\
34 -s <seed> : set seed for random()\n\
35 EXPERIMENTAL OPTIONS:\n\
36 -L : calculate expected length distributions for states\n";
37
38
39 static char banner[] = "covee: emit sequences from a covariance model";
40
41 int
main(int argc,char ** argv)42 main(int argc, char **argv)
43 {
44 char *cmfile; /* file to read model from */
45 struct cm_s *cm; /* model */
46 int i; /* counter for sequences */
47 char *seq; /* generated sequence */
48 char *khseq; /* generated structure */
49 struct align_s *ali; /* generated "alignment" */
50
51 int emitnum; /* number of sequences to emit */
52 int seed; /* seed for random number generator */
53 int do_best; /* TRUE if generating only best seq */
54 int do_landscape; /* TRUE if printing as landscape */
55 int watsoncrick; /* TRUE to annotate only canonical pairs */
56 int do_lengths; /* TRUE to do length distributions */
57
58 int optc;
59 extern char *optarg; /* for getopt() */
60 extern int optind; /* for getopt() */
61
62 /***********************************************
63 * Parse command line
64 ***********************************************/
65
66 emitnum = 20;
67 seed = (int) time (0);
68 do_best = FALSE;
69 do_landscape = FALSE;
70 watsoncrick = TRUE;
71 do_lengths = FALSE;
72
73 while ((optc = getopt(argc, argv, OPTIONS)) != -1)
74 switch (optc) {
75
76 case 'a': watsoncrick = FALSE; break;
77 case 'b': do_best = TRUE; break;
78 case 'l': do_landscape = TRUE; break;
79 case 's': seed = atoi(optarg); break;
80 case 'L': do_lengths = TRUE; break;
81
82 default:
83 Die("Error: unrecognized option %c\n", optc);
84 }
85
86 if (argc - optind != 1)
87 Die("Wrong number of arguments.\n%s", usage);
88
89 cmfile = argv[optind];
90 sre_srandom((unsigned) seed);
91
92 /***********************************************
93 * Read in the model
94 ***********************************************/
95
96 if (! ReadCM(cmfile, &cm))
97 Die("Failed to read model from file %s\n", cmfile);
98
99 /***********************************************
100 * Generate sequences from model and print them.
101 ***********************************************/
102
103 puts(banner);
104 printf(" version %s, %s\n\n", RELEASE, RELEASEDATE);
105
106 if (do_lengths)
107 {
108 struct pstate_s *pcm;
109 int statenum;
110 double **lmx;
111 int *min, *max;
112 int y;
113
114 MakePCM(cm, &pcm, &statenum);
115 NormalizePCM(pcm, statenum);
116 LengthDistribution(pcm, statenum, 200, &lmx);
117 LengthBounds(lmx, statenum, 200, 1.0e-6, &min, &max);
118
119 for (y = 0; y < statenum; y++)
120 printf("%4d %4d %4d (%4d) %8.8f %s\n",
121 y, min[y], max[y], max[y]-min[y]+1,
122 (float) (max[y]-min[y]+1) / 200.0,
123 UstatetypeName(pcm[y].statetype));
124
125 Free2DArray(lmx, statenum);
126 free(min);
127 free(max);
128 free(pcm);
129 }
130
131 else if (do_best)
132 {
133 if (! EmitBestSequence(cm, watsoncrick, &ali, &khseq, &seq)) Die("EmitBestSequence() failed");
134
135 if (do_landscape)
136 {
137 if (! PrintAliLandscape(stdout, cm, ali))
138 Warn("PrintAliLandscape failed\n");
139 }
140 else
141 {
142 printf("%s\n", seq);
143 printf("%s\n", khseq);
144 puts("");
145 }
146
147 free(ali);
148 free(khseq);
149 free(seq);
150 }
151
152 else
153 {
154 for (i = 0; i < emitnum; i++)
155 {
156 if (! EmitSequence(cm, watsoncrick, &ali, &khseq, &seq))
157 Die("failed to generate a sequence from the model.");
158
159 if (do_landscape)
160 PrintAliLandscape(stdout, cm, ali);
161 else
162 {
163 printf("seq %2d: %s\n", i, seq);
164 printf(" %s\n", khseq);
165 puts("");
166 }
167
168 free(ali);
169 free(khseq);
170 free(seq);
171 }
172 }
173
174 /***********************************************
175 * Cleanup and exit
176 ***********************************************/
177
178 FreeCM(cm);
179 return 0;
180 }
181
182