source: branches/stable/TOOLS/arb_probe_match.cxx

Last change on this file was 18222, checked in by westram, 5 years ago
  • typo in documentation
  • improve CLI help for arb_probe_match
File size: 7.2 KB
Line 
1#include <arbdb.h>
2#include <PT_com.h>
3#include <client.h>
4
5#include <iostream>
6
7using std::cerr;
8using std::endl;
9
10static void help() {
11    cerr <<
12        "Executes a probe match on a running PT server\n"
13        "\n"
14        "Parameters:\n"
15        "\n"
16        " --port <PORT>           port of PT server\n"
17        " --sequence <SEQUENCE>   probe\n"
18        "\n"
19        " --complementFirst       before matching first complement probe (applied before '--alsoRevCompl')\n"
20        " --alsoRevCompl          also report matches of reverse-complement probe\n"
21        "  (Note: specifying both does report complement plus reverse matches)\n"
22        "\n"
23        " --max-hits <NUM>        maximum number of hits (default:1,000,000; 0 means 'unlimited')\n"
24        " --weighted              use weighted matching\n"
25        " --weighted-pos          use weighted matching with pos&strength\n"
26        " --n-matches <NUM1>      accept NUM1 occurrences of 'N' as match (default:2)\n"
27        " --n-match-bound <NUM2>  when NUM2 occurrences of 'N' are encountered => all 'N's count as mismatch (default:5)\n"
28        "                         (inbetween all 'N' above NUM1 count as mismatch)"
29        " --mismatches <NUM>      accept NUM mismatches (default:0)\n"
30
31         << endl;
32}
33
34enum SORT_BY_TYPE {
35    SORT_BY_MISMATCHES = 0,
36    SORT_BY_WEIGHTED_MISMATCHES = 1,
37    SORT_BY_WEIGHTED_MISMATCHES_W_POS_AND_STRENGTH = 2
38};
39
40int ARB_main(int argc, char ** argv) {
41    const char*  port            = NULp;
42    const char*  sequence        = NULp;
43    int          alsoRevCompl    = 0;
44    int          complementFirst = 0;
45    int          max_hits        = 1000000;
46    SORT_BY_TYPE sort_by         = SORT_BY_MISMATCHES;
47    int          n_accept        = 2;
48    int          n_limit         = 5;
49    int          max_mismatches  = 0;
50
51    GB_shell shell;
52
53    if (argc<2) {
54        help();
55        return 0;
56    }
57
58    for (int i = 1; i < argc; i++) {
59        int err = 0;
60        if (!strcmp(argv[i], "--help")) {
61            help();
62            return 0;
63        } else if (!strcmp(argv[i], "--alsoRevCompl")) {
64            alsoRevCompl = 1;
65        } else if (!strcmp(argv[i], "--complementFirst")) {
66            complementFirst = 1;
67        } else if (!strcmp(argv[i], "--weighted")) {
68            sort_by = SORT_BY_WEIGHTED_MISMATCHES;
69        } else if (!strcmp(argv[i], "--weighted-pos")) {
70            sort_by = SORT_BY_WEIGHTED_MISMATCHES_W_POS_AND_STRENGTH;
71        } else if (argc > i+1) {
72            if (!strcmp(argv[i], "--port")) {
73                port = argv[++i];
74            } else if (!strcmp(argv[i], "--sequence")) {
75                sequence = argv[++i];
76            } else if (!strcmp(argv[i], "--max-hits")) {
77                max_hits = atoi(argv[++i]);
78            } else if (!strcmp(argv[i], "--n-matches")) {
79                n_accept = atoi(argv[++i]);
80            } else if (!strcmp(argv[i], "--n-match-bound")) {
81                n_limit = atoi(argv[++i]);
82            } else if (!strcmp(argv[i], "--mismatches")) {
83                max_mismatches = atoi(argv[++i]);
84            } else {
85                err = 1;
86            }
87        } else {
88            err = 1;
89        }
90        if (err) {
91            cerr << "Error: Did not understand argument '" << argv[i]
92                 << "'." << endl << endl;
93            help();
94            return 1;
95        }
96    }
97
98    bool err = false;
99
100    if (!port) {
101        cerr << "need '--port' parameter" << endl;
102        err = true;
103    }
104    if (!sequence) {
105        cerr << "need '--sequence' parameter" << endl;
106        err = true;
107    }
108    if (max_mismatches < 0 || max_mismatches > 10) {
109        cerr << "mismatches must be between 0 and 10" << endl;
110        err = true;
111    }
112    if (err) {
113        help();
114        return 1;
115    }
116
117    T_PT_LOCS  locs;
118    T_PT_MAIN  com;
119    GB_ERROR   error = NULp;
120    aisc_com  *link  = aisc_open(port, com, AISC_MAGIC_NUMBER, &error);
121
122    if (!error) {
123        if (!link) {
124            error = "Cannot contact PT server [1]";
125        }
126        else if (aisc_create(link, PT_MAIN, com, MAIN_LOCS, PT_LOCS, locs, NULp)) {
127            error = "Cannot contact PT server [2]";
128        }
129    }
130
131    if (error) {
132        cerr << "Error in arb_probe_match: " << error << endl;
133        return 1;
134    }
135
136    if (aisc_nput(link, PT_LOCS, locs,
137                  LOCS_MATCH_ALSO_REVCOMP,   (long)alsoRevCompl,
138                  LOCS_COMPLEMENT_FIRST,     (long)complementFirst,
139                  LOCS_MATCH_MAX_HITS,       (long)max_hits,
140                  LOCS_MATCH_SORT_BY,        (long)sort_by,
141                  LOCS_MATCH_N_ACCEPT,       (long)n_accept,
142                  LOCS_MATCH_N_LIMIT,        (long)n_limit,
143                  LOCS_MATCH_MAX_MISMATCHES, (long)max_mismatches,
144                  LOCS_SEARCHMATCH,          sequence,
145                  NULp)) {
146        cerr << "Connection to PT server lost" << endl;
147        return 1;
148    }
149
150    T_PT_MATCHLIST match_list;
151
152    long match_list_cnt    = 0;
153    long matches_truncated = 0;
154
155    bytestring bs;
156    bs.data = NULp;
157
158    char *locs_error = NULp;
159
160    aisc_get(link, PT_LOCS, locs,
161             LOCS_MATCH_LIST,        match_list.as_result_param(),
162             LOCS_MATCH_LIST_CNT,    &match_list_cnt,
163             LOCS_MATCH_STRING,      &bs,
164             LOCS_MATCHES_TRUNCATED, &matches_truncated,
165             LOCS_ERROR,             &locs_error,
166             NULp);
167
168    if (locs_error) {
169        if (*locs_error) cerr << "received error: " << *locs_error << endl;
170        free(locs_error);
171    }
172
173    if (matches_truncated) {
174        std::cout << "match list was truncated!" << endl;
175    }
176
177    std::cout << "acc     \t"
178              << "start\t"
179              << "stop\t"
180              << "pos\t"
181              << "mis\t"
182              << "wmis\t"
183              << "nmis\t"
184              << "dt\t"
185              << "rev\t"
186              << "seq"
187              << std::endl;
188
189    while (match_list.exists()) {
190        char *m_acc, *m_sequence;
191        long m_start, m_stop, m_pos, m_mismatches, m_n_mismatches;
192        long m_reversed;
193        double m_wmismatches, m_dt;
194
195        aisc_get(link, PT_MATCHLIST, match_list,
196                 MATCHLIST_FIELD_ACC,     &m_acc,
197                 MATCHLIST_FIELD_START,   &m_start,
198                 MATCHLIST_FIELD_STOP,    &m_stop,
199                 MATCHLIST_POS,           &m_pos,
200                 MATCHLIST_MISMATCHES,    &m_mismatches,
201                 MATCHLIST_WMISMATCHES,   &m_wmismatches,
202                 MATCHLIST_N_MISMATCHES,  &m_n_mismatches,
203                 MATCHLIST_DT,            &m_dt,
204                 MATCHLIST_OVERLAY,       &m_sequence,
205                 MATCHLIST_REVERSED,      &m_reversed,
206                 MATCHLIST_NEXT,          match_list.as_result_param(),
207                 NULp);
208
209        std::cout << m_acc << "\t"
210                  << m_start << "\t"
211                  << m_stop << "\t"
212                  << m_pos << "\t"
213                  << m_mismatches << "\t"
214                  << m_wmismatches << "\t"
215                  << m_n_mismatches << "\t"
216                  << m_dt << "\t"
217                  << m_reversed << "\t"
218                  << m_sequence << "\t"
219                  << std::endl;
220        fflush(stdout);
221
222        free(m_acc);
223        free(m_sequence);
224    }
225
226    free(bs.data);
227    aisc_close(link, com);
228    return 0;
229}
230
231
Note: See TracBrowser for help on using the repository browser.