root/trunk/AWTC/AWTC_next_neighbours.cxx

Revision 8744, 17.2 KB (checked in by westram, 46 hours ago)
  • AISC clients
    • client-handles for server-side objects are distinguishable types now (esp. can no longer be confounded w/o compiler messages)
    • main server object has to be passed to aisc_close
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1// ================================================================ //
2//                                                                  //
3//   File      : AWTC_next_neighbours.cxx                           //
4//   Purpose   :                                                    //
5//                                                                  //
6//   Institute of Microbiology (Technical University Munich)        //
7//   http://www.arb-home.de/                                        //
8//                                                                  //
9// ================================================================ //
10
11#include "awtc_next_neighbours.hxx"
12
13#include <servercntrl.h>
14#include <PT_com.h>
15#include <client.h>
16#include <aw_window.hxx>
17#include <aw_root.hxx>
18#include <arbdbt.h>
19#include <arb_strbuf.h>
20
21#include <climits>
22
23struct PT_FF_comImpl {
24    aisc_com  *link;
25    T_PT_MAIN  com;
26    T_PT_LOCS  locs;
27
28    PT_FF_comImpl()
29        : link(NULL)
30    {
31        ff_assert(!com.exists());
32        ff_assert(!locs.exists());
33    }
34};
35
36// -------------------
37//      FamilyList
38
39FamilyList::FamilyList()
40    : next(NULL),
41      name(NULL),
42      matches(0),
43      rel_matches(0.0)
44{}
45
46FamilyList::~FamilyList() {
47    free(name);
48    delete next;
49}
50
51FamilyList *FamilyList::insertSortedBy_matches(FamilyList *other) {
52    ff_assert(next == NULL); // only insert unlinked instace!
53
54    if (!other) {
55        return this;
56    }
57
58    if (matches >= other->matches) {
59        next = other;
60        return this;
61    }
62
63    // insert into other
64    FamilyList *rest = other->next;
65    other->next      = rest ? insertSortedBy_matches(rest) : this;
66
67    return other;
68}
69
70FamilyList *FamilyList::insertSortedBy_rel_matches(FamilyList *other) {
71    ff_assert(next == NULL); // only insert unlinked instace!
72
73    if (rel_matches >= other->rel_matches) {
74        next = other;
75        return this;
76    }
77
78    // insert into other
79    FamilyList *rest = other->next;
80    other->next      = rest ? insertSortedBy_rel_matches(rest) : this;
81
82    return other;
83
84}
85
86// ---------------------
87//      FamilyFinder
88
89FamilyFinder::FamilyFinder(bool rel_matches_, RelativeScoreScaling scaling_)
90    : rel_matches(rel_matches_),
91      scaling(scaling_), 
92      family_list(NULL),
93      hits_truncated(false),
94      real_hits(-1),
95      range(-1, -1)
96{
97}
98
99FamilyFinder::~FamilyFinder() {
100    delete_family_list();
101}
102
103void FamilyFinder::delete_family_list() {
104    delete family_list;
105    family_list = NULL;
106}
107
108// ------------------------
109//      PT_FamilyFinder
110
111PT_FamilyFinder::PT_FamilyFinder(GBDATA *gb_main_, int server_id_, int oligo_len_, int mismatches_, bool fast_flag_, bool rel_matches_, RelativeScoreScaling scaling_)
112    : FamilyFinder(rel_matches_, scaling_), 
113      gb_main(gb_main_),
114      server_id(server_id_),
115      oligo_len(oligo_len_),
116      mismatches(mismatches_),
117      fast_flag(fast_flag_)
118{
119    // 'mismatches'  = the number of allowed mismatches
120    // 'fast_flag'   = 0 -> do complete search, 1 -> search only oligos starting with 'A'
121    // 'rel_matches' = 0 -> score is number of oligo-hits, 1 -> score is relative to longer sequence (target or source) * 10
122
123    ci = new PT_FF_comImpl;
124}
125
126PT_FamilyFinder::~PT_FamilyFinder() {
127    delete_family_list();
128    close();
129    delete ci;
130}
131
132GB_ERROR PT_FamilyFinder::init_communication() {
133    const char *user = "PT_FamilyFinder";
134
135    ff_assert(!ci->locs.exists());
136   
137    // connect PT server
138    if (aisc_create(ci->link, PT_MAIN, ci->com,
139                    MAIN_LOCS, PT_LOCS, ci->locs,
140                    LOCS_USER, user,
141                    NULL)) {
142        return GB_export_error("Cannot initialize communication");
143    }
144    return 0;
145}
146
147
148GB_ERROR PT_FamilyFinder::open(const char *servername) {
149    GB_ERROR error = 0;
150   
151    ff_assert(!ci->com.exists() && !ci->link);
152
153    if (arb_look_and_start_server(AISC_MAGIC_NUMBER, servername)) {
154        error = "Cannot contact PT  server";
155    }
156    else {
157        const char *socketid = GBS_read_arb_tcp(servername);
158        if (!socketid) error = GB_await_error();
159        else {
160            ci->link = aisc_open(socketid, ci->com, AISC_MAGIC_NUMBER);
161            if (!ci->link) error = "Cannot contact PT server [1]";
162            else if (init_communication()) error = "Cannot contact PT server [2]";
163        }
164    }
165
166    ff_assert(error || (ci->com.exists() && ci->link));
167   
168    return error;
169}
170
171void PT_FamilyFinder::close() {
172    if (ci->link) aisc_close(ci->link, ci->com);
173    ci->link = 0;
174    ff_assert(!ci->com.exists());
175    ci->locs.clear();
176}
177
178GB_ERROR PT_FamilyFinder::retrieve_family(const char *sequence, FF_complement compl_mode, int max_results, double min_score) {
179    delete_family_list();
180
181    char     *compressed_sequence = GB_command_interpreter(gb_main, sequence, "|keep(acgtunACGTUN)", 0, 0);
182    GB_ERROR  error               = NULL;
183
184    if      (!compressed_sequence)    error = GB_await_error();
185    else if (!compressed_sequence[0]) error = "No data in sequence(-region)";
186    else {
187        bytestring bs;
188        bs.data = compressed_sequence;
189        bs.size = strlen(bs.data)+1;
190
191        /* Start find_family() at the PT_SERVER
192         *
193         * Here we have to make a loop, until the match count of the
194         * first member is big enough
195         */
196
197        ff_assert(!range.is_empty());
198         
199        // create and init family finder object
200        T_PT_FAMILYFINDER ffinder;
201        if (aisc_create(ci->link, PT_LOCS, ci->locs,
202                        LOCS_FFINDER, PT_FAMILYFINDER, ffinder,
203                        FAMILYFINDER_PROBE_LEN,       long(oligo_len),          // oligo length (12 hardcoded till July 2008)
204                        FAMILYFINDER_MISMATCH_NUMBER, long(mismatches),         // number of mismatches (0 hardcoded till July 2008)
205                        FAMILYFINDER_FIND_TYPE,       long(fast_flag),          // 0: complete search, 1: quick search (only search oligos starting with 'A')
206                        FAMILYFINDER_SORT_TYPE,       long(uses_rel_matches()), // 0: matches, 1: relative matches (0 hardcoded till July 2008)
207                        FAMILYFINDER_REL_SCORING,     long(get_scaling()),      // scaling of relative scores
208                        FAMILYFINDER_SORT_MAX,        long(max_results),        // speed up family sorting (only sort retrieved results)
209                        FAMILYFINDER_MIN_SCORE,       min_score,                // limit hits by score
210                        FAMILYFINDER_COMPLEMENT,      long(compl_mode),         // any combination of: 1 = forward, 2 = reverse, 4 = reverse-complement, 8 = complement (1 hardcoded in PT-Server till July 2008)
211                        FAMILYFINDER_RANGE_STARTPOS,  long(range.start()),
212                        FAMILYFINDER_RANGE_ENDPOS,    long(range.is_limited() ? range.end() : -1),
213                        FAMILYFINDER_FIND_FAMILY,     &bs,                      // RPC (has to be last parameter!)
214                        NULL))
215        {
216            error = "Communication error with PT server ('retrieve_family')";
217        }
218        else {
219            char *ff_error;
220
221            // Read family list
222            T_PT_FAMILYLIST f_list;
223            aisc_get(ci->link, PT_FAMILYFINDER, ffinder,
224                     FAMILYFINDER_FAMILY_LIST,      f_list.as_result_param(),
225                     FAMILYFINDER_FAMILY_LIST_SIZE, &real_hits,
226                     FAMILYFINDER_ERROR,            &ff_error,
227                     NULL);
228
229            if (ff_error[0]) {
230                error = GBS_global_string("PTSERVER: %s", ff_error);
231            }
232            else {
233                hits_truncated = false;
234                if (max_results<1) max_results = INT_MAX;
235
236                FamilyList *tail = NULL;
237                while (f_list.exists()) {
238                    if (max_results == 0) {
239                        hits_truncated = true;
240                        break;
241                    }
242                    max_results--;
243
244                    FamilyList *fl = new FamilyList();
245
246                    (tail ? tail->next : family_list) = fl;
247                    tail                              = fl;
248                    fl->next                          = NULL;
249
250                    aisc_get(ci->link, PT_FAMILYLIST, f_list,
251                             FAMILYLIST_NAME,        &fl->name,
252                             FAMILYLIST_MATCHES,     &fl->matches,
253                             FAMILYLIST_REL_MATCHES, &fl->rel_matches,
254                             FAMILYLIST_NEXT,        f_list.as_result_param(),
255                             NULL);
256                }
257            }
258            free(ff_error);
259        }
260
261        free(compressed_sequence);
262    }
263    return error;
264}
265
266GB_ERROR PT_FamilyFinder::searchFamily(const char *sequence, FF_complement compl_mode, int max_results, double min_score) {
267    // searches the PT-server for species related to 'sequence'.
268    //
269    // relation-score is calculated by fragmenting the sequence into oligos of length 'oligo_len' and
270    // then summarizing the number of hits.
271    //
272    // 'max_results' limits the length of the generated result list (low scores deleted first)
273    //               if < 1 -> don't limit
274    //
275    // 'min_score' limits the results by score (use 0 for unlimited results)
276    //
277    // When using restrict_2_region(), only pass the corresponding part via 'sequence' (not the full alignment)
278
279    GB_ERROR error = range.is_empty() ? "Specified range is empty" : NULL;
280    if (!error) {
281        error             = open(GBS_ptserver_tag(server_id));
282        if (!error) error = retrieve_family(sequence, compl_mode, max_results, min_score);
283        close();
284    }
285    return error;
286}
287
288const char *PT_FamilyFinder::results2string() {
289    GBS_strstruct *out = GBS_stropen(1000);
290    for (FamilyList *fl = family_list; fl; fl = fl->next) {
291        GBS_strnprintf(out, 100, "%s/%li/%3.5f,", fl->name, fl->matches, fl->rel_matches*100);
292    }
293    GBS_str_cut_tail(out, 1);
294    RETURN_LOCAL_ALLOC(GBS_strclose(out));
295}
296
297void AWTC_create_common_next_neighbour_vars(AW_root *aw_root) {
298    static bool created = false;
299    if (!created) {
300        aw_root->awar_int(AWAR_NN_OLIGO_LEN,   12);
301        aw_root->awar_int(AWAR_NN_MISMATCHES,  0);
302        aw_root->awar_int(AWAR_NN_FAST_MODE,   0);
303        aw_root->awar_int(AWAR_NN_REL_MATCHES, 1);
304        aw_root->awar_int(AWAR_NN_REL_SCALING, RSS_BOTH_MIN);
305
306        created = true;
307    }
308}
309
310void AWTC_create_common_next_neighbour_fields(AW_window *aws) {
311    // used in several figs:
312    // - ad_spec_nn.fig
313    // - ad_spec_nnm.fig
314    // - faligner/family_settings.fig
315
316    aws->at("oligo_len");
317    aws->create_input_field(AWAR_NN_OLIGO_LEN, 3);
318
319    aws->at("mismatches");
320    aws->create_input_field(AWAR_NN_MISMATCHES, 3);
321
322    aws->at("mode");
323    aws->create_option_menu(AWAR_NN_FAST_MODE, 0, 0);
324    aws->insert_default_option("Complete", "", 0);
325    aws->insert_option("Quick", "", 1);
326    aws->update_option_menu();
327
328    aws->at("score");
329    aws->create_option_menu(AWAR_NN_REL_MATCHES, 0, 0);
330    aws->insert_option("absolute", "", 0);
331    aws->insert_default_option("relative", "", 1);
332    aws->update_option_menu();
333
334    aws->at("scaling");
335    aws->create_option_menu(AWAR_NN_REL_SCALING, 0, 0);
336    aws->insert_option        ("to source POC",  "", RSS_SOURCE);
337    aws->insert_option        ("to target POC",  "", RSS_TARGET);
338    aws->insert_default_option("to maximum POC", "", RSS_BOTH_MAX);
339    aws->insert_option        ("to minimum POC", "", RSS_BOTH_MIN);
340    aws->update_option_menu();
341}
342
343// --------------------------------------------------------------------------------
344
345#ifdef UNIT_TESTS
346
347#include <test_unit.h>
348
349class ff_tester {
350    GBDATA *gb_main;
351public:
352    bool    relativeMatches;
353    bool    fastMode;
354    bool    partial;
355    bool    shortOligo;
356    double  min_score;
357
358    RelativeScoreScaling scaling;
359
360    ff_tester(GBDATA *gb_main_)
361        : gb_main(gb_main_),
362          relativeMatches(false),
363          fastMode(false),
364          partial(false),
365          shortOligo(false),
366          min_score(0.0),
367          scaling(RSS_BOTH_MAX)
368    {}
369
370    const char *get_result(GB_ERROR& error) {
371        int             oligoLen = shortOligo ? (partial ? 3 : 6) : (partial ? 10 : 18);
372        PT_FamilyFinder ff(gb_main, TEST_SERVER_ID, oligoLen, 1, fastMode, relativeMatches, scaling);
373       
374        const char *sequence;
375        if (partial) {
376            ff.restrict_2_region(PosRange(39, 91)); // alignment range of bases 30-69 of sequence of 'LgtLytic'
377            sequence = "UCUAGCUUGCUAGACGGGUGGCGAG" "GGUAACCGUAGGGGA"; // bases 30-54 of sequence of 'LgtLytic' + 15 bases from 'DcdNodos' (outside region)
378        }
379        else {
380            // sequence of 'LgtLytic' in TEST_pt.arb:
381            sequence = "AGAGUUUGAUCAAGUCGAACGGCAGCACAGUCUAGCUUGCUAGACGGGUGGCGAGUGGCGAACGGACUUGGGGAAACUCAAGCUAAUACCGCAUAAUCAUGACUGGGGUGAAGUCGUAACAAGGUAGCCGUAGGGGAACCUGCGGCUGGAUCACCUCCUN";
382        }
383
384        error = ff.searchFamily(sequence, FF_FORWARD, 4, min_score);
385        return ff.results2string();
386    }
387};
388
389
390#define TEST_RELATIVES_COMMON(tester,expctd) \
391        GB_ERROR    error;                                      \
392        const char *result   = tester.get_result(error);        \
393        const char *expected = expctd;                          \
394        TEST_ASSERT_NO_ERROR(error);                            \
395   
396#define TEST_ASSERT_RELATIVES(tester,expctd) do {               \
397        TEST_RELATIVES_COMMON(tester,expctd);                   \
398        TEST_ASSERT_EQUAL(result, expected);                    \
399    } while(0)
400
401#define TEST_ASSERT_REL__BROK(tester,expctd) do {               \
402        TEST_RELATIVES_COMMON(tester,expctd);                   \
403        TEST_ASSERT_EQUAL__BROKEN(result, expected);            \
404    } while(0)
405   
406void TEST_SLOW_PT_FamilyFinder() {
407    GB_shell shell;
408    TEST_SETUP_GLOBAL_ENVIRONMENT("ptserver");
409
410    GBDATA *gb_main = GB_open("no.arb", "c");
411
412    // check some error cases
413    {
414        PT_FamilyFinder ffe(gb_main, TEST_SERVER_ID, 0, 1, 0, 0, RSS_BOTH_MAX);
415        TEST_ASSERT_CONTAINS(ffe.searchFamily("whatever", FF_FORWARD, 4, 0.0), "minimum oligo length is 1");
416    }
417   
418    ff_tester test(gb_main);
419
420    ff_tester ______RESET = test; TEST_ASSERT_RELATIVES(test, "LgtLytic/142/97.93103,HllHalod/62/43.05556,AclPleur/59/38.06452,PtVVVulg/52/34.21053");
421    test.partial          = true; TEST_ASSERT_RELATIVES(test, "LgtLytic/18/11.76471,VblVulni/5/3.24675,VbhChole/4/2.59740,DcdNodos/4/2.59740");
422    test.shortOligo       = true; TEST_ASSERT_RELATIVES(test, "PtVVVulg/38/22.75449,AclPleur/38/22.35294,VbhChole/38/23.60248,VblVulni/38/23.60248");
423    test.relativeMatches  = true; TEST_ASSERT_RELATIVES(test, "DsssDesu/38/38.77551,CltBotul/38/34.23423,PsAAAA00/38/32.75862,Bl0LLL00/38/25.67568");
424    test.min_score        = 32.6; TEST_ASSERT_RELATIVES(test, "DsssDesu/38/38.77551,CltBotul/38/34.23423,PsAAAA00/38/32.75862");
425    test                  = ______RESET;
426    test.shortOligo       = true; TEST_ASSERT_RELATIVES(test, "LgtLytic/153/97.45223,AclPleur/138/82.63473,VbhChole/133/84.17722,VblVulni/133/84.17722");
427    test.relativeMatches  = true; TEST_ASSERT_RELATIVES(test, "LgtLytic/153/97.45223,HllHalod/133/85.25641,VbhChole/133/84.17722,VblVulni/133/84.17722");
428    test.fastMode         = true; TEST_ASSERT_RELATIVES(test, "LgtLytic/42/26.75159,VblVulni/37/23.41772,HllHalod/36/23.07692,Stsssola/36/23.07692");
429    test.min_score        = 26.7; TEST_ASSERT_RELATIVES(test, "LgtLytic/42/26.75159");
430    test.min_score        = 26.8; TEST_ASSERT_RELATIVES(test, "");
431    test                  = ______RESET;
432    test.fastMode         = true; TEST_ASSERT_RELATIVES(test, "LgtLytic/40/27.58621,HllHalod/18/12.50000,AclPleur/17/10.96774,PtVVVulg/15/9.86842");
433    test.min_score        = 17.0; TEST_ASSERT_RELATIVES(test, "LgtLytic/40/27.58621,HllHalod/18/12.50000,AclPleur/17/10.96774");
434    test.min_score        = 17.5; TEST_ASSERT_RELATIVES(test, "LgtLytic/40/27.58621,HllHalod/18/12.50000");
435    test                  = ______RESET;
436
437    test.shortOligo      = true;
438    test.relativeMatches = true;
439    test.scaling         = RSS_BOTH_MAX; TEST_ASSERT_RELATIVES(test, "LgtLytic/153/97.45223,HllHalod/133/85.25641,VbhChole/133/84.17722,VblVulni/133/84.17722");
440    test.scaling         = RSS_BOTH_MIN; TEST_ASSERT_RELATIVES(test, "LgtLytic/153/98.07692,AclPleur/138/88.46154,DsssDesu/84/88.42105,CltBotul/95/87.96296");
441    test.scaling         = RSS_TARGET;   TEST_ASSERT_RELATIVES(test, "LgtLytic/153/97.45223,DsssDesu/84/88.42105,CltBotul/95/87.96296,PsAAAA00/97/85.84071");
442    test.scaling         = RSS_SOURCE;   TEST_ASSERT_RELATIVES(test, "LgtLytic/153/98.07692,AclPleur/138/88.46154,VbhChole/133/85.25641,VblVulni/133/85.25641");
443    test.partial         = true;
444    test.shortOligo      = false;
445    test.scaling         = RSS_BOTH_MAX; TEST_ASSERT_RELATIVES(test, "LgtLytic/18/11.76471,VblVulni/5/3.24675,VbhChole/4/2.59740,DcdNodos/4/2.59740");
446    test.scaling         = RSS_BOTH_MIN; TEST_ASSERT_RELATIVES(test, "LgtLytic/18/56.25000,VblVulni/5/15.62500,VbhChole/4/12.50000,DcdNodos/4/12.50000");
447    test.scaling         = RSS_TARGET;   TEST_ASSERT_RELATIVES(test, "LgtLytic/18/11.76471,VblVulni/5/3.24675,VbhChole/4/2.59740,DcdNodos/4/2.59740");
448    test.scaling         = RSS_SOURCE;   TEST_ASSERT_RELATIVES(test, "LgtLytic/18/56.25000,VblVulni/5/15.62500,VbhChole/4/12.50000,DcdNodos/4/12.50000");
449    test                 = ______RESET;
450
451    GB_close(gb_main);
452}
453
454#endif
Note: See TracBrowser for help on using the browser.