source: tags/arb-6.0/TOOLS/arb_primer.cxx

Last change on this file was 11021, checked in by westram, 9 years ago
  • remove class init via memset (arb_prm_struct + AW_xfig)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.7 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : arb_primer.cxx                                    //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include <arbdbt.h>
12#include <arb_strarray.h>
13
14#define ADD_LEN        10
15#define PRM_BUFFERSIZE 256
16
17
18struct arb_prm_struct : virtual Noncopyable {
19    ConstStrArray   alignment_names;
20    int             al_len;
21    int             max_name;
22    GBDATA         *gb_main;
23    char            buffer[PRM_BUFFERSIZE];
24    const char     *source;
25    int             prmanz;
26    int             prmlen;
27    int             prmsmin;
28    char          **data;
29    int             sp_count;
30    int             key_cnt;
31    int             one_key_cnt;
32    int             reduce;
33    FILE           *out;
34    char           *outname;
35
36    arb_prm_struct()
37        : al_len(0),
38          max_name(0),
39          gb_main(NULL),
40          source(NULL),
41          prmanz(0),
42          prmlen(0),
43          prmsmin(0),
44          data(NULL),
45          sp_count(0),
46          key_cnt(0),
47          one_key_cnt(0),
48          reduce(0),
49          out(NULL),
50          outname(NULL)
51    {
52        memset(buffer, 0, sizeof(buffer));
53    }
54    ~arb_prm_struct() {
55        if (data) {
56            for (int i = 0; i<sp_count; ++i) free(data[i]);
57            free(data);
58        }
59        free(outname);
60    }
61
62};
63static arb_prm_struct aprm;
64
65inline int getNumFromStdin() {
66    // returns entered number (or 0)
67    int i = 0;
68    if (fgets(aprm.buffer, PRM_BUFFERSIZE, stdin)) {
69        i = atoi(aprm.buffer);
70    }
71    return i;
72}
73
74static GB_ERROR arb_prm_menu() {
75    printf(" Please select an Alignment:\n");
76    int i;
77    for (i=1; aprm.alignment_names[i-1]; ++i) {
78        printf("%i: %s\n", i, aprm.alignment_names[i-1]);
79    }
80    aprm.max_name = i;
81
82    GB_ERROR error = NULL;
83    i = getNumFromStdin();
84    if ((i<1) || (i>=aprm.max_name)) {
85        error = GBS_global_string("selection %i out of range", i);
86    }
87    else {
88        aprm.source = aprm.alignment_names[i-1];
89
90        printf("This module will search for primers for all positions.\n"
91               "   The best result is one primer for all (marked) taxa , the worst case\n"
92               "   are n primers for n taxa.\n"
93               "   Please specify the maximum number of primers:\n"
94            );
95        aprm.prmanz = getNumFromStdin();
96
97        printf("Select minimum length of a primer, the maximum will be (minimum + %i)\n", ADD_LEN);
98        i = getNumFromStdin();
99        if ((i<4) || (i>30)) {
100            error = GBS_global_string("selection %i out of range", i);
101        }
102        else {
103            aprm.prmlen = i;
104
105            printf("There may be short sequences or/and deletes in full sequences\n"
106                   "   So a primer normally does not match all sequences\n"
107                   "   Specify minimum percentage of species (0-100 %%):\n");
108            i = getNumFromStdin();
109            if ((i<1) || (i>100)) {
110                error = GBS_global_string("selection %i out of range", i);
111            }
112            else {
113                aprm.prmsmin = i;
114
115                printf("Write output to file (enter \"\" to write to screen)\n");
116                if (fgets(aprm.buffer, PRM_BUFFERSIZE, stdin)) {
117                    char *lf      = strchr(aprm.buffer, '\n');
118                    if (lf) lf[0] = 0; // remove linefeed from filename
119                    aprm.outname  = strdup(aprm.buffer);
120                }
121                else {
122                    aprm.outname  = strdup("");
123                }
124            }
125        }
126    }
127    return error;
128}
129
130static GB_ERROR arb_prm_read(int /* prmanz */) {
131    GBDATA *gb_presets = GBT_get_presets(aprm.gb_main);
132    {
133        GBDATA *gb_source = GB_find_string(gb_presets, "alignment_name", aprm.source, GB_IGNORE_CASE, SEARCH_GRANDCHILD);
134        GBDATA *gb_len    = GB_brother(gb_source, "alignment_len");
135        aprm.al_len       = GB_read_int(gb_len);
136    }
137
138    int sp_count = GBT_count_marked_species(aprm.gb_main);
139    aprm.data    = (char **)calloc(sp_count, sizeof(char *));
140
141    sp_count = 0;
142    for (GBDATA *gb_species = GBT_first_marked_species(aprm.gb_main);
143         gb_species;
144         gb_species = GBT_next_marked_species(gb_species))
145    {
146        GBDATA *gb_source = GB_entry(gb_species, aprm.source);
147        if (gb_source) {
148            GBDATA *gb_source_data = GB_entry(gb_source, "data");
149            if (gb_source_data) {
150                const char *hdata = GB_read_char_pntr(gb_source_data);
151
152                if (!hdata) {
153                    GB_print_error();
154                }
155                else {
156                    char *data             = (char *)calloc(sizeof(char), aprm.al_len+1);
157                    aprm.data[sp_count ++] = data;
158                   
159                    if (sp_count % 50 == 0) printf("Reading taxa %i\n", sp_count);
160
161                    int size = GB_read_string_count(gb_source_data);
162                    int i;
163                    for (i=0; i<size; i++) {
164                        char c = hdata[i];
165                        if ((c>='a') && (c<='z')) {
166                            data[i] = c-'a'+'A';
167                        }
168                        else {
169                            data[i] = c;
170                        }
171                    }
172                    for (; i<aprm.al_len; i++) {
173                        data[i] = '.';
174                    }
175                    data[i] = 0;
176                }
177            }
178        }
179    }
180    printf("%i taxa read\n", sp_count);
181    aprm.sp_count = sp_count;
182    if (sp_count == 0) {
183        return "No marked taxa found";
184    }
185    return NULL;
186}
187
188static long arb_count_keys(const char * /* key */, long val, void *)
189{
190    if (val >1) {
191        aprm.key_cnt++;
192    }
193    else {
194        aprm.one_key_cnt++;
195    }
196    return val;
197}
198
199static long arb_print_primer(const char *key, long val, void *)
200{
201    if (val <= 1) return val;
202    int gc = 0;
203    const char *p;
204    for (p = key; *p; p++) {
205        if (*p == 'G' || *p == 'C') gc++;
206    }
207    fprintf(aprm.out, "  %s matching %4li taxa     GC = %3i%%\n",
208            key, val, 100*gc/(int)strlen(key));
209    return val;
210}
211
212#define is_base(c) (((c>='a') && (c<='z')) || ((c>='A')&&(c<='Z')))
213
214static int primer_print(char *dest, char * source, int size)
215{
216    char c;
217    c = *(source++);
218    if (!is_base(c)) return 1;
219    while (size) {
220        while (!is_base(c)) {
221            c = *(source++);
222            if (!c) return 1;
223        }
224        if (c == 'N' || c == 'n') return 1;
225        *(dest++) = c;
226        size--;
227        if (!c) return 1;
228        c = 0;
229    }
230    *dest = 0;
231    return 0;
232}
233
234
235static long arb_reduce_primer_len(const char *key, long val, void *cl_hash) {
236    GB_HASH* hash = (GB_HASH*)cl_hash;
237    char     buffer[256];
238    int      size = strlen(key)-aprm.reduce;
239
240    strncpy(buffer, key, size);
241    buffer[size] = 0;
242    val += GBS_read_hash(hash, buffer);
243    GBS_write_hash(hash, buffer, val);
244    return val;
245}
246
247static void arb_prm_primer(int /* prmanz */)
248{
249    GB_HASH *mhash;
250    int      sp;
251    char    *buffer;
252    int      pos;
253    int      prmlen;
254    int      pspecies;
255    int      cutoff_cnt;
256    int     *best_primer_cnt;
257    int     *best_primer_new;
258    int     *best_primer_swap;
259
260    prmlen = aprm.prmlen + ADD_LEN + 1;
261
262    buffer = (char *) calloc(sizeof(char), prmlen + 1);
263    best_primer_cnt = (int *)calloc(prmlen+1, sizeof(int));
264    best_primer_new = (int *)calloc(prmlen+1, sizeof(int));
265
266    for (pos = 0; pos < aprm.al_len; pos++) {
267        prmlen = aprm.prmlen + ADD_LEN;
268        mhash = GBS_create_hash(1024, GB_MIND_CASE);
269        pspecies = 0;
270        if (pos % 50 == 0) printf("Pos. %i (%i)\n", pos, aprm.al_len);
271        cutoff_cnt = aprm.prmanz+1;
272        for (sp = 0; sp < aprm.sp_count; sp++) {    // build initial hash table
273            if (!primer_print(buffer, aprm.data[sp] + pos, prmlen)) {
274                GBS_incr_hash(mhash, buffer);
275                pspecies++;
276            }
277        }
278        if (pspecies*100 >= aprm.prmsmin * aprm.sp_count) {     // reduce primer length
279            for (; prmlen >= aprm.prmlen; prmlen-=aprm.reduce) {
280                GB_HASH *hash = GBS_create_hash(aprm.prmanz, GB_MIND_CASE);
281
282                aprm.key_cnt = 0;
283                aprm.one_key_cnt = 0;
284                GBS_hash_do_loop(mhash, arb_count_keys, NULL);
285                if ((aprm.key_cnt + aprm.one_key_cnt < cutoff_cnt) &&
286                    //  (aprm.key_cnt > aprm.one_key_cnt) &&
287                    (aprm.key_cnt<best_primer_cnt[prmlen+1])) {
288                    fprintf(aprm.out, "%3i primer found len %3i(of %4i taxa) for position %i\n", aprm.key_cnt, prmlen, pspecies, pos);
289                    GBS_hash_do_loop(mhash, arb_print_primer, NULL);
290                    fprintf(aprm.out, "\n\n");
291                    cutoff_cnt = aprm.key_cnt;
292                }
293                best_primer_new[prmlen] = aprm.key_cnt;
294                aprm.reduce = 1;
295                while (aprm.key_cnt > aprm.prmanz*4) {
296                    aprm.key_cnt/=4;
297                    aprm.reduce++;
298                }
299                GBS_hash_do_loop(mhash, arb_reduce_primer_len, hash);
300                GBS_free_hash(mhash);
301                mhash = hash;
302            }
303        }
304        else {
305            for (; prmlen>0; prmlen--) best_primer_new[prmlen] = aprm.prmanz+1;
306        }
307        GBS_free_hash(mhash);
308        best_primer_swap = best_primer_new;
309        best_primer_new = best_primer_cnt;
310        best_primer_cnt = best_primer_swap;
311        mhash = 0;
312    }
313
314    free(best_primer_new);
315    free(best_primer_cnt);
316    free(buffer);
317}
318
319int ARB_main(int argc, char *argv[]) {
320    const char *path  = NULL;
321
322    while (argc >= 2) {
323        if (strcmp(argv[1], "--help") == 0) {
324            fprintf(stderr,
325                    "Usage: arb_primer [dbname]\n"
326                    "Searches sequencing primers\n");
327            return EXIT_FAILURE;
328        }
329        path = argv[1];
330        argv++; argc--;
331    }
332
333    if (!path) path = ":";
334
335    GB_ERROR error = NULL;
336    GB_shell shell;
337    aprm.gb_main   = GB_open(path, "r");
338    if (!aprm.gb_main) {
339        error = GBS_global_string("Can't open db '%s' (Reason: %s)", path, GB_await_error());
340    }
341    else {
342        GB_begin_transaction(aprm.gb_main);
343        GBT_get_alignment_names(aprm.alignment_names, aprm.gb_main);
344        GB_commit_transaction(aprm.gb_main);
345
346        error = arb_prm_menu();
347
348        if (!error) {
349            GB_begin_transaction(aprm.gb_main);
350            error = arb_prm_read(aprm.prmanz);
351            if (!error) {
352                GB_commit_transaction(aprm.gb_main);
353                if (strlen(aprm.outname)) {
354                    aprm.out = fopen(aprm.outname, "w");
355                    if (!aprm.out) {
356                        error = GB_IO_error("writing", aprm.outname);
357                    }
358                    else {
359                        arb_prm_primer(aprm.prmanz);
360                        fclose(aprm.out);
361                    }
362                }
363                else {
364                    aprm.out = stdout;
365                    arb_prm_primer(aprm.prmanz);
366                }
367            }
368            else {
369                GB_abort_transaction(aprm.gb_main);
370            }
371        }
372        GB_close(aprm.gb_main);
373    }
374
375    if (error) {
376        fprintf(stderr, "Error in arb_primer: %s\n", error);
377        return EXIT_FAILURE;
378    }
379    return EXIT_SUCCESS;
380}
Note: See TracBrowser for help on using the repository browser.