← Back to branch summary

~ubuntu-branches/ubuntu/dapper/ncbi-tools6/dapper

« back to all changes in this revision

Viewing changes to demo/fdselect.c

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2002-04-04 22:13:09 UTC
  • Revision ID: james.westby@ubuntu.com-20020404221309-vfze028rfnlrldct
Tags: upstream-6.1.20011220a
Import upstream version 6.1.20011220a

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/fdselect.c
 
1
/* $Id: fdselect.c,v 6.3 2000/11/22 21:10:12 shavirin Exp $ */
 
2
/*****************************************************************************
 
3
 
 
4
  
 
5
                          PUBLIC DOMAIN NOTICE
 
6
              National Center for Biotechnology Information
 
7
 
 
8
    This software/database is a "United States Government Work" under the
 
9
    terms of the United States Copyright Act.  It was written as part of    
 
10
    the author's official duties as a United States Government employee
 
11
    and thus cannot be copyrighted.  This software/database is freely
 
12
    available to the public for use. The National Library of Medicine and
 
13
    the U.S. Government have not placed any restriction on its use or
 
14
    reproduction.
 
15
 
 
16
    Although all reasonable efforts have been taken to ensure the accuracy
 
17
    and reliability of the software and data, the NLM and the U.S.
 
18
    Government do not and cannot warrant the performance or results that
 
19
    may be obtained by using this software or data. The NLM and the U.S.
 
20
    Government disclaim all warranties, express or implied, including
 
21
    warranties of performance, merchantability or fitness for any
 
22
    particular purpose.
 
23
 
 
24
    Please cite the author in any work or product based on this material.
 
25
 
 
26
   ***************************************************************************
 
27
 
 
28
   File Name:  fdfilter.c
 
29
 
 
30
   Author:  Sergei B. Shavirin
 
31
   
 
32
   Version Creation Date: 09/13/99
 
33
 
 
34
   $Revision: 6.3 $
 
35
 
 
36
   File Description:  Create few subsets of FASTA database.
 
37
 
 
38
   $Log: fdselect.c,v $
 
39
   Revision 6.3  2000/11/22 21:10:12  shavirin
 
40
   Added tax_id parameter into function FDBAddBioseq.
 
41
 
 
42
   Revision 6.2  2000/03/13 18:37:37  madden
 
43
   Added insert_ctrlA Boolean to readdb_get_bioseq_ex
 
44
 
 
45
   Revision 6.1  1999/09/13 16:20:00  shavirin
 
46
   Initial version.
 
47
 
 
48
 
 
49
*****************************************************************************/
 
50
#include <ncbi.h>
 
51
#include <readdb.h>
 
52
#include <ncbiwww.h>
 
53
 
 
54
/* We will use regular WWW encoding of request to make specific
 
55
   filtering using database information file 
 
56
 
 
57
   The most general query string will look like:
 
58
 
 
59
   taxid=555,666,777&owner=5,6,7&div=AAA,BBB,CCC
 
60
 
 
61
   */
 
62
 
 
63
#define TAXID_LABEL  "taxid"
 
64
#define OWNER_LABEL  "owner"
 
65
#define DIV_LABEL    "div"
 
66
#define DBNAME_LABEL "dbname"
 
67
 
 
68
#define MAX_SR_ELEMENTS   36
 
69
#define MAX_FLT_DATABASES 64
 
70
 
 
71
typedef struct _SR_Info
 
72
{
 
73
    Int4 taxid[MAX_SR_ELEMENTS];
 
74
    Int4 owner[MAX_SR_ELEMENTS];
 
75
    Char div[MAX_SR_ELEMENTS][4];
 
76
    Char dbname[64];
 
77
} SR_Info, PNTR SR_InfoPtr;
 
78
 
 
79
typedef struct gilist
 
80
{
 
81
    Int4    count;
 
82
    Int4    allocated;
 
83
    Int4Ptr seq_num;
 
84
} GiList, *GiListPtr;
 
85
 
 
86
typedef struct _NewDBdata
 
87
{
 
88
    CharPtr dbname;
 
89
    GiListPtr glp;
 
90
} NewDBdata, PNTR NewDBdataPtr;
 
91
 
 
92
#define NUMARG 8
 
93
 
 
94
Args flt_args[NUMARG] = {
 
95
    { "Title for output database file", 
 
96
      NULL, NULL, NULL, TRUE, 't', ARG_STRING, 0.0, 0, NULL},
 
97
    {"Input file for the filtering (this parameter must be set)",
 
98
     NULL, NULL,NULL,FALSE,'i',ARG_FILE_IN, 0.0,0,NULL},
 
99
    {"Configuration file for database subsets creation",
 
100
     "fdselect.cfg", NULL,NULL,TRUE,'c',ARG_FILE_IN, 0.0,0,NULL},
 
101
    {"Input file with a list of gis",
 
102
     NULL, NULL,NULL,TRUE,'g',ARG_FILE_IN, 0.0,0,NULL},
 
103
    {"Create sparse indexes in the filtered database",
 
104
     "F", NULL,NULL,TRUE,'s',ARG_BOOLEAN, 0.0,0,NULL},
 
105
    {"Query string for creating database subset",
 
106
     NULL, NULL,NULL,TRUE,'q',ARG_STRING, 0.0,0,NULL},
 
107
    {"Input database is protein",
 
108
     "T", NULL,NULL,TRUE,'p',ARG_BOOLEAN, 0.0,0,NULL},
 
109
    {"Logfile name:",
 
110
     "fdselect.log", NULL,NULL,TRUE,'l',ARG_FILE_OUT, 0.0,0,NULL},
 
111
};
 
112
 
 
113
#define FLT_Input  flt_args[1].strvalue
 
114
#define IS_Protein flt_args[6].intvalue
 
115
 
 
116
 
 
117
void FDB_optionsFree(FDB_optionsPtr options)
 
118
{
 
119
    if(options == NULL)
 
120
        return;
 
121
    
 
122
    MemFree(options->db_title);
 
123
    MemFree(options->db_file);
 
124
    MemFree(options->LogFileName);
 
125
    MemFree(options);
 
126
    
 
127
    return;
 
128
}
 
129
 
 
130
Boolean SRReadCharData(CharPtr buffer, CharPtr PNTR div_in)
 
131
{
 
132
    CharPtr tmp, ch, ch2;
 
133
    Int4 j;
 
134
    Char div[MAX_SR_ELEMENTS][4];
 
135
 
 
136
    tmp = StringSave(buffer);
 
137
    MemSet(div, NULLB, sizeof(div));
 
138
 
 
139
    for(ch2 = tmp, j = 0; j < MAX_SR_ELEMENTS; j++) {
 
140
        
 
141
        if((ch = StringChr(ch2, ',')) == NULL) {
 
142
            StringNCpy(div[j], ch2, 3);
 
143
            break;
 
144
        }
 
145
        
 
146
        *ch = NULLB;
 
147
        ch++;
 
148
        StringNCpy(div[j], ch2, 3);
 
149
        ch2 = ch;
 
150
    }
 
151
 
 
152
    MemCpy(div_in, div, sizeof(div));
 
153
 
 
154
    MemFree(tmp);
 
155
    
 
156
    return TRUE;
 
157
}
 
158
 
 
159
Boolean SRReadIntData(CharPtr buffer, Int4Ptr id)
 
160
{
 
161
    CharPtr tmp, ch, ch2;
 
162
    Int4 j;
 
163
    tmp = StringSave(buffer);
 
164
    
 
165
    for(ch2 = tmp, j = 0; j < MAX_SR_ELEMENTS; j++) {
 
166
        
 
167
        if((ch = StringChr(ch2, ',')) == NULL) {
 
168
            id[j] = atol(ch2);
 
169
            id[j+1] = -1; /* Terminating character */
 
170
            break;
 
171
        }
 
172
        
 
173
        *ch = NULLB;
 
174
        ch++;
 
175
        id[j] = atol(ch2);
 
176
        ch2 = ch;
 
177
    }
 
178
 
 
179
    MemFree(tmp);
 
180
    
 
181
    return TRUE;
 
182
}
 
183
 
 
184
SR_InfoPtr SRReadSRInfo(CharPtr buffer)
 
185
{
 
186
    WWWInfoPtr info;
 
187
    WWWInfoDataPtr info_data;
 
188
    SR_InfoPtr srip;
 
189
    CharPtr chptr;
 
190
 
 
191
    info_data = (WWWInfoDataPtr) MemNew(sizeof(WWWInfoData));
 
192
    info_data->query = StringSave(buffer);
 
193
    info_data->entries = WWWGetEntries(&info_data->num_entries, 
 
194
                                       info_data->query, FALSE);
 
195
    info = (VoidPtr) info_data;
 
196
 
 
197
    srip = MemNew(sizeof(SR_Info));
 
198
    
 
199
    if((chptr = WWWGetValueByName(info, TAXID_LABEL)) != NULL)
 
200
        SRReadIntData(chptr, srip->taxid);
 
201
    if((chptr = WWWGetValueByName(info, OWNER_LABEL)) != NULL)
 
202
        SRReadIntData(chptr, srip->owner);
 
203
    if((chptr = WWWGetValueByName(info, DIV_LABEL)) != NULL)
 
204
        SRReadCharData(chptr, (CharPtr PNTR) srip->div);
 
205
    if((chptr = WWWGetValueByName(info, DBNAME_LABEL)) != NULL)
 
206
        StringCpy(srip->dbname, chptr);
 
207
 
 
208
    /* If no conditions exists - this is an error */
 
209
    if(srip->taxid[0] <= 0 && srip->owner[0] <= 0 && **srip->div == NULLB) {
 
210
        ErrPostEx(SEV_ERROR, 0, 0, 
 
211
                  "No valid conditions exists in query string");
 
212
        MemFree(srip);
 
213
        srip = NULL;
 
214
    }
 
215
    WWWInfoFree(info);
 
216
    return srip;
 
217
}
 
218
 
 
219
FDB_optionsPtr FDB_CreateCLOptions(Boolean is_prot)
 
220
{
 
221
    FDB_optionsPtr options;
 
222
    Char buffer[128];
 
223
 
 
224
    options = MemNew(sizeof(FDB_options));
 
225
    
 
226
    options->db_title = StringSave(flt_args[0].strvalue);
 
227
 
 
228
    sprintf(buffer, "%s.flt", FLT_Input);
 
229
    options->db_file = StringSave(buffer);
 
230
 
 
231
    options->LogFileName = StringSave(flt_args[7].strvalue);
 
232
    options->is_protein = is_prot;
 
233
    options->parse_mode = TRUE;
 
234
    options->isASN = FALSE;
 
235
    options->asnbin = FALSE;
 
236
    options->is_seqentry = FALSE;
 
237
    options->base_name = NULL;
 
238
    options->dump_info = FALSE;
 
239
    options->sparse_idx = flt_args[4].intvalue;
 
240
 
 
241
    return options;
 
242
}
 
243
 
 
244
#define GI_ALLOC_CHUNK 1024
 
245
 
 
246
SeqIdPtr MySeqIdFree(SeqIdPtr sip)
 
247
{
 
248
    SeqIdPtr sip_tmp;
 
249
    do {
 
250
        sip_tmp = sip->next;
 
251
        SeqIdFree(sip);
 
252
        sip = sip_tmp;
 
253
    } while(sip != NULL);
 
254
    
 
255
    return NULL;
 
256
}
 
257
 
 
258
/* Functions used in filtering by gi number */
 
259
GiListPtr GiListNew(void)
 
260
{
 
261
    GiListPtr glp;
 
262
 
 
263
    glp = MemNew(sizeof(GiList));
 
264
    glp->allocated = GI_ALLOC_CHUNK;
 
265
    glp->seq_num = MemNew(sizeof(Int4) * glp->allocated);
 
266
    glp->count = 0;
 
267
    
 
268
    return glp;
 
269
}
 
270
void GiListFree(GiListPtr glp)
 
271
{
 
272
    if(glp == NULL)
 
273
        return;
 
274
    
 
275
    MemFree(glp->seq_num);
 
276
    MemFree(glp);
 
277
    return;
 
278
}
 
279
 
 
280
Boolean ReadGiList(ReadDBFILEPtr rdfp, GiListPtr glp, CharPtr filename)
 
281
{
 
282
    FILE *fd;
 
283
    Int4 gi, retvalue, seqnum;
 
284
 
 
285
    if((fd = FileOpen(filename, "r")) == NULL)
 
286
        return FALSE;
 
287
 
 
288
    while((retvalue = fscanf(fd, "%d", &gi)) != EOF) {
 
289
        if(retvalue == 0) continue;
 
290
        
 
291
        if(glp->count >= glp->allocated) {
 
292
            glp->allocated += GI_ALLOC_CHUNK;
 
293
            glp->seq_num = Realloc(glp->seq_num, 
 
294
                                   sizeof(Int4) * glp->allocated);
 
295
        }
 
296
 
 
297
        seqnum = readdb_gi2seq(rdfp, gi);
 
298
            
 
299
        if(seqnum < 0) {
 
300
            ErrPostEx(SEV_WARNING, 0,0, "Gi %d is not found", gi);
 
301
            continue;
 
302
        }
 
303
 
 
304
        glp->seq_num[glp->count] = seqnum;
 
305
        glp->count++;
 
306
    }
 
307
 
 
308
    FileClose(fd);
 
309
 
 
310
    return TRUE;
 
311
}
 
312
/* Here we will check, that data[2] == tax_id, 
 
313
   data[3] == owner, div=div */
 
314
Boolean CheckSRCondition(SR_InfoPtr srip, Int4Ptr data, CharPtr div)
 
315
{
 
316
    Int4    i;
 
317
    CharPtr chptr;
 
318
    Boolean  cond_ok = FALSE;
 
319
 
 
320
    /* checking tax_id */
 
321
 
 
322
    if(srip->taxid[0] > 0) { /* At least one element exists */
 
323
        cond_ok = FALSE;
 
324
        for(i = 0; srip->taxid[i] > 0 && i < MAX_SR_ELEMENTS; i++) {
 
325
            if(data[2] == srip->taxid[i]) {
 
326
                cond_ok = TRUE;
 
327
                break;
 
328
            }
 
329
        }
 
330
        if(cond_ok == FALSE)
 
331
            return FALSE;
 
332
    }
 
333
 
 
334
    /* checking owner */
 
335
 
 
336
    if(srip->owner[0] > 0) { /* At least one element exists */
 
337
        cond_ok = FALSE;
 
338
        for(i = 0; srip->owner[i] > 0 && i < MAX_SR_ELEMENTS; i++) {
 
339
            if(data[3] == srip->owner[i]) {
 
340
                cond_ok = TRUE;
 
341
                break;
 
342
            }
 
343
        }
 
344
        if(cond_ok == FALSE)
 
345
            return FALSE;
 
346
    }
 
347
 
 
348
    /* checking division */
 
349
 
 
350
    if(*srip->div[0] != NULLB) { /* At least one element exists */
 
351
        cond_ok = FALSE;
 
352
        for(i = 0; *srip->div[i] != NULLB && i < MAX_SR_ELEMENTS; i++) {
 
353
            if(!StringCmp(div, srip->div[i])) {
 
354
                cond_ok = TRUE;
 
355
                break;
 
356
            }
 
357
        }
 
358
        if(cond_ok == FALSE)
 
359
            return FALSE;
 
360
    }
 
361
    
 
362
    return TRUE;
 
363
}
 
364
 
 
365
Boolean FDGetGiListByQuery(ReadDBFILEPtr rdfp, SR_InfoPtr srip, 
 
366
                           GiListPtr glp, CharPtr filename)
 
367
{ 
 
368
    FILE *fd;
 
369
    Int4 length, data[7];
 
370
    Char div[32];
 
371
    Char buffer[1024];
 
372
    
 
373
    if((fd = FileOpen(filename, "r")) == NULL) {
 
374
        ErrPostEx(SEV_ERROR, 0,0, "Unable to open input info file");
 
375
        return FALSE;
 
376
    }
 
377
    
 
378
    length = sizeof(buffer);
 
379
    while(fgets(buffer, length, fd) != NULL) {
 
380
        sscanf(buffer, "%d %d %d %d %s %d %d %d",
 
381
               &data[0], &data[1], &data[2], &data[3],
 
382
               div, &data[4], &data[5], &data[6]);
 
383
        
 
384
        /* If line agree with condition seq_num added to the list */
 
385
        if(CheckSRCondition(srip, data, div)) {
 
386
            if(glp->count >= glp->allocated) {
 
387
                glp->allocated += GI_ALLOC_CHUNK;
 
388
                glp->seq_num = Realloc(glp->seq_num, 
 
389
                                       sizeof(Int4) * glp->allocated);
 
390
            }
 
391
            glp->seq_num[glp->count] = data[0];
 
392
            glp->count++;
 
393
        }
 
394
    }
 
395
    
 
396
    FileClose(fd);
 
397
    
 
398
    return TRUE;
 
399
}
 
400
 
 
401
VoidPtr NewDBThread(VoidPtr data)
 
402
{
 
403
    NewDBdataPtr ndbp;
 
404
    ReadDBFILEPtr rdfp;
 
405
    FDB_optionsPtr options;
 
406
    FormatDBPtr fdbp;
 
407
    BioseqPtr bsp;
 
408
    Int4 i, seqnum;
 
409
 
 
410
    if((ndbp = data) == NULL)
 
411
        return NULL;
 
412
 
 
413
    if((rdfp = readdb_new (FLT_Input, IS_Protein)) == NULL) {
 
414
        ErrPostEx(SEV_ERROR, 0, 0, 
 
415
                  "Failure to intialise database %s", FLT_Input);
 
416
        return NULL;
 
417
    }
 
418
    
 
419
    /* ---------------------------------------------- */
 
420
    /* ------- Initializing thread options  --------- */
 
421
    /* ---------------------------------------------- */
 
422
 
 
423
    if((options = FDB_CreateCLOptions(IS_Protein)) == NULL)
 
424
        return NULL;
 
425
    
 
426
    MemFree(options->db_file);
 
427
    options->db_file = StringSave(ndbp->dbname);
 
428
    
 
429
    /* ---------------------------------------------- */
 
430
    /* ----- Initializing formatdb structure  ------- */
 
431
    /* ---------------------------------------------- */
 
432
    
 
433
    if ((fdbp = FormatDBInit(options)) == NULL)
 
434
        return NULL;
 
435
    
 
436
    /* ---------------------------------------------- */
 
437
    /* ---------------- Main loop ------------------- */
 
438
    /* ---------------------------------------------- */
 
439
    
 
440
    for(i = 0; i < ndbp->glp->count; i++) {
 
441
        
 
442
        if((seqnum = ndbp->glp->seq_num[i]) == -1)
 
443
            continue;
 
444
        
 
445
        bsp = readdb_get_bioseq_ex(rdfp, seqnum, FALSE, FALSE);        
 
446
        FDBAddBioseq(fdbp, bsp, 0);
 
447
 
 
448
        SeqIdSetFree(bsp->id);
 
449
        BioseqFreeComponents(bsp);
 
450
        MemFree(bsp);
 
451
 
 
452
        /* BioseqFree(bsp); */
 
453
    }
 
454
 
 
455
    GiListFree(ndbp->glp);
 
456
    
 
457
    if(FormatDBClose(fdbp))
 
458
        return NULL;
 
459
    
 
460
    readdb_destruct(rdfp);
 
461
    FDB_optionsFree(options);
 
462
    
 
463
    return NULL;
 
464
}
 
465
 
 
466
NewDBdataPtr GetNdprByString(CharPtr string, ReadDBFILEPtr rdfp, 
 
467
                             CharPtr indname)
 
468
{
 
469
    NewDBdataPtr ndbp;
 
470
    SR_InfoPtr srip;
 
471
    
 
472
    if((srip = SRReadSRInfo(string)) == NULL)
 
473
        return NULL;
 
474
    
 
475
    ndbp = MemNew(sizeof(NewDBdata));
 
476
    ndbp->glp = GiListNew();
 
477
    
 
478
    if(!FDGetGiListByQuery(rdfp, srip, ndbp->glp, indname))
 
479
        return NULL;
 
480
    
 
481
    /* If database name set - replacing default one */
 
482
    if(*srip->dbname != NULLB) {
 
483
        ndbp->dbname = StringSave(srip->dbname);
 
484
    }
 
485
    MemFree(srip);
 
486
    return ndbp;
 
487
}
 
488
void TrimSpaces(CharPtr buffer)
 
489
{
 
490
    CharPtr chptr;
 
491
    Int4 i;
 
492
    
 
493
    if(buffer == NULL)
 
494
        return;
 
495
    
 
496
    for(chptr = buffer, i = 0; *chptr != NULLB; chptr++) {
 
497
        if(IS_WHITESP(*chptr))
 
498
            continue;
 
499
        buffer[i] = *chptr;
 
500
        i++;
 
501
    }
 
502
    buffer[i] = NULLB;
 
503
    return;
 
504
}
 
505
 
 
506
Int2 Main(void)
 
507
{
 
508
    Int4 i, glp_count, buflen;
 
509
    NewDBdataPtr ndbp[MAX_FLT_DATABASES];
 
510
    Char indname[128], buffer[1024];
 
511
    FILE *fd;
 
512
    ReadDBFILEPtr rdfp;
 
513
    
 
514
    if ( !GetArgs ("fdselect", NUMARG, flt_args) )
 
515
        return NULL;
 
516
    
 
517
    if ( !ErrSetLog (flt_args[7].strvalue) ) { /* Logfile */
 
518
        ErrShow();
 
519
    } else {
 
520
        ErrSetOpts (ERR_CONTINUE, ERR_LOG_ON);
 
521
    }
 
522
 
 
523
    /* ---------------------------------------------- */
 
524
    /* ------ Initializing readdb structures -------- */
 
525
    /* ---------------------------------------------- */
 
526
 
 
527
    if((rdfp = readdb_new (FLT_Input, IS_Protein)) == NULL) {
 
528
        ErrPostEx(SEV_ERROR, 0, 0, 
 
529
                  "Failure to intialise database %s", FLT_Input);
 
530
        return NULL;
 
531
    }
 
532
        
 
533
    sprintf(indname, "%s.%cdi", FLT_Input, IS_Protein? 'p' : 'n');
 
534
    
 
535
    if(flt_args[3].strvalue != NULL) {/* This is list of gis */
 
536
        
 
537
        ndbp[0] = MemNew(sizeof(NewDBdata));
 
538
        ndbp[0]->glp = GiListNew();
 
539
        glp_count = 1;
 
540
        
 
541
        if(!ReadGiList(rdfp, ndbp[0]->glp, flt_args[2].strvalue))
 
542
            return -1;
 
543
        
 
544
    } else if(flt_args[5].strvalue != NULL) {/* This is string for gis */
 
545
        
 
546
        TrimSpaces(flt_args[5].strvalue);
 
547
        
 
548
        ndbp[0] = GetNdprByString(flt_args[5].strvalue, rdfp, indname);
 
549
        glp_count = 1;
 
550
 
 
551
    } else {
 
552
 
 
553
        /* Default behaviour is to read and interprete every line of the
 
554
           configuration file */
 
555
        
 
556
        if((fd = FileOpen(flt_args[2].strvalue, "r")) == NULL) {
 
557
            ErrPostEx(SEV_ERROR, 0, 0, "Failure to open config file %s",
 
558
                      flt_args[2].strvalue);
 
559
        }
 
560
 
 
561
        buflen = sizeof(buffer);
 
562
        for(glp_count = 0; fgets(buffer, buflen, fd) != NULL && 
 
563
                glp_count < MAX_FLT_DATABASES;) {
 
564
            
 
565
            if(*buffer == '#')
 
566
                continue;
 
567
 
 
568
            TrimSpaces(buffer);
 
569
            
 
570
            ndbp[glp_count] = GetNdprByString(buffer, rdfp, indname);
 
571
            glp_count++;
 
572
        }
 
573
    }
 
574
    
 
575
    readdb_destruct(rdfp);
 
576
 
 
577
    for(i = 0; i < glp_count; i++) {
 
578
        NlmThreadCreate(NewDBThread, (VoidPtr) ndbp[i]);
 
579
    }
 
580
 
 
581
    NlmThreadJoinAll();
 
582
 
 
583
    return 0;
 
584
}
 
585