← Back to branch summary

~ubuntu-branches/ubuntu/precise/ncbi-tools6/precise

« back to all changes in this revision

Viewing changes to demo/gi2accn.c

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2005-03-27 12:00:15 UTC
  • mfrom: (2.1.2 hoary)
  • Revision ID: james.westby@ubuntu.com-20050327120015-embhesp32nj73p9r
Tags: 6.1.20041020-3
http://bugs.debian.org/295110
* Fix FTBFS under GCC 4.0 caused by inconsistent use of "static" on
  functions.  (Closes: #295110.)
* Add a watch file, now that we can.  (Upstream's layout needs version=3.)

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/gi2accn.c
 
1
/*   gi2accn.c
 
2
* ===========================================================================
 
3
*
 
4
*                            PUBLIC DOMAIN NOTICE
 
5
*            National Center for Biotechnology Information (NCBI)
 
6
*
 
7
*  This software/database is a "United States Government Work" under the
 
8
*  terms of the United States Copyright Act.  It was written as part of
 
9
*  the author's official duties as a United States Government employee and
 
10
*  thus cannot be copyrighted.  This software/database is freely available
 
11
*  to the public for use. The National Library of Medicine and the U.S.
 
12
*  Government do not place any restriction on its use or reproduction.
 
13
*  We would, however, appreciate having the NCBI and the author cited in
 
14
*  any work or product based on this material
 
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 particular
 
22
*  purpose.
 
23
*
 
24
* ===========================================================================
 
25
*
 
26
* File Name:  gi2accn.c
 
27
*
 
28
* Author:  Jonathan Kans
 
29
*
 
30
* Version Creation Date:   4/15/02
 
31
*
 
32
* $Revision: 6.6 $
 
33
*
 
34
* File Description: 
 
35
*
 
36
* Modifications:  
 
37
* --------------------------------------------------------------------------
 
38
* Date     Name        Description of modification
 
39
* -------  ----------  -----------------------------------------------------
 
40
*
 
41
*
 
42
* ==========================================================================
 
43
*/
 
44
 
 
45
#include <ncbi.h>
 
46
#include <objall.h>
 
47
#include <objsset.h>
 
48
#include <objsub.h>
 
49
#include <objfdef.h>
 
50
#include <sequtil.h>
 
51
#include <gather.h>
 
52
#include <sqnutils.h>
 
53
#include <subutil.h>
 
54
#include <explore.h>
 
55
#include <pmfapi.h>
 
56
#include <ent2api.h>
 
57
 
 
58
static void ConvertGiToAccn (SeqIdPtr sip, Pointer userdata)
 
59
 
 
60
{
 
61
  Int4      gi;
 
62
  SeqIdPtr  newsip;
 
63
 
 
64
  if (sip == NULL) return;
 
65
  if (sip->choice != SEQID_GI) return;
 
66
  gi = sip->data.intvalue;
 
67
  newsip = GetSeqIdForGI (gi);
 
68
  if (newsip == NULL) return;
 
69
  if (newsip->choice == SEQID_GIBBSQ ||
 
70
      newsip->choice == SEQID_GIBBMT ||
 
71
      newsip->choice == SEQID_GI) {
 
72
    SeqIdFree (newsip);
 
73
    return;
 
74
  }
 
75
  SeqIdStripLocus (newsip);
 
76
  sip->choice = newsip->choice;
 
77
  sip->data.ptrvalue = newsip->data.ptrvalue;
 
78
  newsip->choice = SEQID_NOT_SET;
 
79
  newsip->data.ptrvalue = NULL;
 
80
  SeqIdFree (newsip);
 
81
}
 
82
 
 
83
static void UpdateAligns (
 
84
  SeqAlignPtr sap,
 
85
  Pointer userdata
 
86
)
 
87
 
 
88
{
 
89
  VisitSeqIdsInSeqAlign (sap, NULL, ConvertGiToAccn);
 
90
}
 
91
 
 
92
static Boolean IsSipMrna (SeqIdPtr sip, Int4Ptr gilist, Int4 count)
 
93
 
 
94
{
 
95
  BioseqPtr  bsp;
 
96
  Int4       gi = 0;
 
97
  Int4       L, R, mid;
 
98
 
 
99
  if (sip == NULL) return FALSE;
 
100
  bsp = BioseqFind (sip);
 
101
  if (bsp != NULL) return FALSE;
 
102
  if (sip->choice == SEQID_GI) {
 
103
    gi = (Int4) sip->data.intvalue;
 
104
  } else {
 
105
    gi = GetGIForSeqId (sip);
 
106
  }
 
107
  if (gi < 1) return FALSE;
 
108
 
 
109
  L = 0;
 
110
  R = count - 1;
 
111
 
 
112
  while (L < R) {
 
113
    mid = (L + R) / 2;
 
114
    if (gilist [mid] < gi) {
 
115
      L = mid + 1;
 
116
    } else {
 
117
      R = mid;
 
118
    }
 
119
  }
 
120
 
 
121
  if (gilist [R] == gi) return TRUE;
 
122
 
 
123
  return FALSE;
 
124
}
 
125
 
 
126
static Boolean IsMrnaAlignment (SeqAlignPtr align, Int4Ptr gilist, Int4 count)
 
127
 
 
128
{
 
129
  DenseDiagPtr  ddp;
 
130
  DenseSegPtr   dsp;
 
131
  SeqIdPtr      sip;
 
132
  StdSegPtr     ssp;
 
133
  SeqLocPtr     tloc;
 
134
 
 
135
  if (align == NULL) return FALSE;
 
136
  sip = NULL;
 
137
  if (align->segtype == 1) {
 
138
    ddp = (DenseDiagPtr) align->segs;
 
139
    if (ddp != NULL) {
 
140
      for (sip = ddp->id; sip != NULL; sip = sip->next) {
 
141
        if (IsSipMrna (sip, gilist, count)) return TRUE;
 
142
      }
 
143
    }
 
144
  } else if (align->segtype == 2) {
 
145
    dsp = (DenseSegPtr) align->segs;
 
146
    if (dsp != NULL) {
 
147
      for (sip = dsp->ids; sip != NULL; sip = sip->next) {
 
148
        if (IsSipMrna (sip, gilist, count)) return TRUE;
 
149
      }
 
150
    }
 
151
  } else if (align->segtype == 3) {
 
152
    ssp = (StdSegPtr) align->segs;
 
153
    if (ssp != NULL) {
 
154
      for (tloc = ssp->loc; tloc != NULL; tloc = tloc->next) {
 
155
        sip = SeqLocId (tloc);
 
156
        if (IsSipMrna (sip, gilist, count)) return TRUE;
 
157
      }
 
158
    }
 
159
  }
 
160
  return FALSE;
 
161
}
 
162
 
 
163
static SeqAnnotPtr MakeMrnaSeqAnnot (void)
 
164
 
 
165
{
 
166
  AnnotDescrPtr  adp;
 
167
  SeqAnnotPtr    annot;
 
168
  ObjectIdPtr    oip;
 
169
  UserFieldPtr   ufp;
 
170
  UserObjectPtr  uop;
 
171
 
 
172
  annot = SeqAnnotNew ();
 
173
  annot->type = 2;
 
174
 
 
175
  adp = ValNodeNew (NULL);
 
176
  adp->choice = Annot_descr_user;
 
177
  uop = UserObjectNew ();
 
178
  adp->data.ptrvalue = uop;
 
179
  oip = ObjectIdNew ();
 
180
  oip->str = StringSave ("Blast Type");
 
181
  ufp = UserFieldNew ();
 
182
  uop->type = oip;
 
183
  uop->data = ufp;
 
184
  oip = ObjectIdNew ();
 
185
  oip->str = StringSave ("BLASTN - mrna");
 
186
  ufp->label = oip;
 
187
  ufp->choice = 2;
 
188
  ufp->data.intvalue = 1;
 
189
  annot->desc = adp;
 
190
 
 
191
  adp = ValNodeNew (NULL);
 
192
  adp->choice = Annot_descr_user;
 
193
  uop = UserObjectNew ();
 
194
  adp->data.ptrvalue = uop;
 
195
  oip = ObjectIdNew ();
 
196
  oip->str = StringSave ("Hist Seqalign");
 
197
  ufp = UserFieldNew ();
 
198
  uop->type = oip;
 
199
  uop->data = ufp;
 
200
  oip = ObjectIdNew ();
 
201
  oip->str = StringSave ("Hist Seqalign");
 
202
  ufp->label = oip;
 
203
  ufp->choice = 4;
 
204
  ufp->data.boolvalue = TRUE;
 
205
  adp->next = annot->desc;
 
206
  annot->desc = adp;
 
207
 
 
208
  return annot;
 
209
}
 
210
 
 
211
static SeqAnnotPtr ExtractBlastMrna (SeqAlignPtr sap, Pointer PNTR prevlink, Int4Ptr gilist, Int4 count)
 
212
 
 
213
{
 
214
  SeqAnnotPtr  annot = NULL;
 
215
  SeqAlignPtr  next;
 
216
 
 
217
  while (sap != NULL) {
 
218
    next = sap->next;
 
219
 
 
220
    if (IsMrnaAlignment (sap, gilist, count)) {
 
221
      *prevlink = sap->next;
 
222
      sap->next = NULL;
 
223
 
 
224
      if (annot == NULL) {
 
225
        annot = MakeMrnaSeqAnnot ();
 
226
      }
 
227
      if (annot != NULL) {
 
228
        sap->next = annot->data;
 
229
        annot->data = sap;
 
230
      }
 
231
 
 
232
    } else {
 
233
      sap->idx.prevlink = prevlink;
 
234
      prevlink = (Pointer PNTR) &(sap->next);
 
235
    }
 
236
 
 
237
    sap = next;
 
238
  }
 
239
 
 
240
  return annot;
 
241
}
 
242
 
 
243
static void CountGIsInAligns (SeqIdPtr sip, Pointer userdata)
 
244
 
 
245
{
 
246
  Int4Ptr  countp;
 
247
 
 
248
  if (sip == NULL || sip->choice != SEQID_GI || userdata == NULL) return;
 
249
  countp = (Int4Ptr) userdata;
 
250
  (*countp)++;
 
251
}
 
252
 
 
253
typedef struct strgi {
 
254
  Int4     count;
 
255
  Int4Ptr  uidlist;
 
256
} UidData, PNTR UidDataPtr;
 
257
 
 
258
static void StoreGIsFromAligns (SeqIdPtr sip, Pointer userdata)
 
259
 
 
260
{
 
261
  UidDataPtr  udp;
 
262
 
 
263
  if (sip == NULL || sip->choice != SEQID_GI || userdata == NULL) return;
 
264
  udp = (UidDataPtr) userdata;
 
265
  udp->uidlist [udp->count] = (Int4) sip->data.intvalue;
 
266
  (udp->count)++;
 
267
}
 
268
 
 
269
static int LIBCALLBACK SortByInt4 (VoidPtr a, VoidPtr b)
 
270
 
 
271
{
 
272
  if (a == NULL || b == NULL) return 0;
 
273
  return (*(const Int4 *) a) - (*(const Int4 *) b);
 
274
}
 
275
 
 
276
static void FindBlastMrna (SeqAnnotPtr sap, Int4Ptr uidlist, Int4 count)
 
277
 
 
278
{
 
279
  SeqAnnotPtr             annot;
 
280
  Entrez2BooleanReplyPtr  e2br;
 
281
  Entrez2IdListPtr        e2id;
 
282
  Entrez2RequestPtr       e2rq;
 
283
  Entrez2ReplyPtr         e2ry;
 
284
  Int4                    i;
 
285
  Int4Ptr                 gilist;
 
286
  E2ReplyPtr              reply;
 
287
 
 
288
  e2rq = EntrezCreateBooleanRequest (TRUE, FALSE, "nucleotide", NULL, 0, 0, NULL, 0, 0);
 
289
  EntrezAddToBooleanRequest (e2rq, "biomol_mrna [PROP]", 0, NULL, NULL, NULL, 0, 0, NULL, NULL, TRUE, TRUE);
 
290
  EntrezAddToBooleanRequest (e2rq, NULL, ENTREZ_OP_AND, NULL, NULL, NULL, 0, 0, NULL, NULL, TRUE, TRUE);
 
291
  EntrezAddToBooleanRequest (e2rq, NULL, 0, NULL, NULL, NULL, 0, count, uidlist, NULL, TRUE, TRUE);
 
292
  e2ry = EntrezSynchronousQuery (e2rq);
 
293
  e2rq = Entrez2RequestFree (e2rq);
 
294
  if (e2ry == NULL) return;
 
295
  reply = e2ry->reply;
 
296
  if (reply == NULL || reply->choice != E2Reply_eval_boolean) return;
 
297
  e2br = EntrezExtractBooleanReply (e2ry);
 
298
  if (e2br == NULL) return;
 
299
  count = e2br->count;
 
300
  if (count < 1) return;
 
301
  e2id = e2br->uids;
 
302
  if (e2id == NULL || e2id->num < 1 || e2id->uids == NULL) return;
 
303
 
 
304
  gilist = (Int4Ptr) MemNew (sizeof (Int4) * (e2id->num + 1));
 
305
  if (gilist != NULL) {
 
306
 
 
307
    BSSeek (e2id->uids, 0, SEEK_SET);
 
308
    for (i = 0; i < e2id->num; i++) {
 
309
      gilist [i] = Nlm_BSGetUint4 (e2id->uids);
 
310
    }
 
311
    HeapSort (gilist, e2id->num, sizeof (Int4), SortByInt4);
 
312
 
 
313
    annot = ExtractBlastMrna ((SeqAlignPtr) sap->data, (Pointer PNTR) &(sap->data), gilist, e2id->num);
 
314
    if (annot != NULL) {
 
315
      annot->next = sap->next;
 
316
      sap->next = annot;
 
317
    }
 
318
 
 
319
    MemFree (gilist);
 
320
  }
 
321
 
 
322
  Entrez2BooleanReplyFree (e2br);
 
323
}
 
324
 
 
325
static void ProcessBlastNR (SeqAnnotPtr sap)
 
326
 
 
327
{
 
328
  Int4     count = 0;
 
329
  Int4     gi;
 
330
  Int4     i;
 
331
  Int4     j;
 
332
  Int4     last;
 
333
  UidData  ud;
 
334
  Int4Ptr  uidlist;
 
335
 
 
336
  if (sap == NULL || sap->type != 2) return;
 
337
 
 
338
  VisitSeqIdsInSeqAnnot (sap, (Pointer) &count, CountGIsInAligns);
 
339
  if (count < 1) return;
 
340
  uidlist = (Int4Ptr) MemNew (sizeof (Int4) * (count + 1));
 
341
  if (uidlist == NULL) return;
 
342
 
 
343
  ud.count = 0;
 
344
  ud.uidlist = uidlist;
 
345
  VisitSeqIdsInSeqAnnot (sap, (Pointer) &ud, StoreGIsFromAligns);
 
346
 
 
347
  HeapSort (uidlist, count, sizeof (Int4), SortByInt4);
 
348
 
 
349
  /* unique sorted gi list */
 
350
 
 
351
  last = 0;
 
352
  for (i = 0, j = 0; i < count; i++) {
 
353
    gi = uidlist [i];
 
354
    if (gi != last) {
 
355
      uidlist [j] = gi;
 
356
      j++;
 
357
    }
 
358
    last = gi;
 
359
  }
 
360
  count = j;
 
361
 
 
362
  FindBlastMrna (sap, uidlist, count);
 
363
 
 
364
  MemFree (uidlist);
 
365
}
 
366
 
 
367
static void FindBlastNR (SeqAnnotPtr sap, Pointer userdata)
 
368
 
 
369
{
 
370
  AnnotDescrPtr  adp;
 
371
  ObjectIdPtr    oip;
 
372
  UserFieldPtr   ufp;
 
373
  UserObjectPtr  uop;
 
374
 
 
375
  if (sap == NULL || sap->type != 2) return;
 
376
  for (adp = sap->desc; adp != NULL; adp = adp->next) {
 
377
    if (adp->choice != Annot_descr_user) continue;
 
378
    for (uop = adp->data.ptrvalue; uop != NULL; uop = uop->next) {
 
379
      oip = uop->type;
 
380
      if (oip == NULL) continue;
 
381
      if (StringCmp (oip->str, "Blast Type") != 0) continue;
 
382
      ufp = uop->data;
 
383
      if (ufp == NULL) continue;
 
384
      oip = ufp->label;
 
385
      if (oip == NULL) continue;
 
386
      if (StringCmp (oip->str, "BLASTN - nr") != 0) continue;
 
387
      ProcessBlastNR (sap);
 
388
      oip->str = MemFree (oip->str);
 
389
      oip->str = StringSave ("BLASTN - nr minus mrna");
 
390
    }
 
391
  }
 
392
}
 
393
 
 
394
/* Args structure contains command-line arguments */
 
395
 
 
396
#define i_argInputFile     0
 
397
#define o_argOutputFile    1
 
398
#define c_argConvertGIs    2
 
399
#define x_argExtractMrnas  3
 
400
 
 
401
Args myargs [] = {
 
402
  {"Input File", "stdin", NULL, NULL,
 
403
    FALSE, 'i', ARG_FILE_IN, 0.0, 0, NULL},
 
404
  {"Output File", "stdout", NULL, NULL,
 
405
    FALSE, 'o', ARG_FILE_OUT, 0.0, 0, NULL},
 
406
  {"Convert GIs", "T", NULL, NULL,
 
407
    TRUE, 'c', ARG_BOOLEAN, 0.0, 0, NULL},
 
408
  {"Extract mRNAs", "T", NULL, NULL,
 
409
    TRUE, 'x', ARG_BOOLEAN, 0.0, 0, NULL},
 
410
};
 
411
 
 
412
Int2 Main (void)
 
413
 
 
414
{
 
415
  AsnIoPtr     aip;
 
416
  Boolean      convertgis;
 
417
  Boolean      extractmrnas;
 
418
  CharPtr      infile, outfile;
 
419
  SeqEntryPtr  sep;
 
420
 
 
421
  /* standard setup */
 
422
 
 
423
  ErrSetFatalLevel (SEV_MAX);
 
424
  ErrClearOptFlags (EO_SHOW_USERSTR);
 
425
  UseLocalAsnloadDataAndErrMsg ();
 
426
  ErrPathReset ();
 
427
 
 
428
  /* finish resolving internal connections in ASN.1 parse tables */
 
429
 
 
430
  if (! AllObjLoad ()) {
 
431
    Message (MSG_FATAL, "AllObjLoad failed");
 
432
    return 1;
 
433
  }
 
434
  if (! SubmitAsnLoad ()) {
 
435
    Message (MSG_FATAL, "SubmitAsnLoad failed");
 
436
    return 1;
 
437
  }
 
438
  if (! FeatDefSetLoad ()) {
 
439
    Message (MSG_FATAL, "FeatDefSetLoad failed");
 
440
    return 1;
 
441
  }
 
442
  if (! SeqCodeSetLoad ()) {
 
443
    Message (MSG_FATAL, "SeqCodeSetLoad failed");
 
444
    return 1;
 
445
  }
 
446
  if (! GeneticCodeTableLoad ()) {
 
447
    Message (MSG_FATAL, "GeneticCodeTableLoad failed");
 
448
    return 1;
 
449
  }
 
450
 
 
451
  /* process command line arguments */
 
452
 
 
453
  if (! GetArgs ("gi2accn", sizeof (myargs) / sizeof (Args), myargs)) {
 
454
    return 0;
 
455
  }
 
456
 
 
457
  infile = (CharPtr) myargs [i_argInputFile].strvalue;
 
458
  outfile = (CharPtr) myargs [o_argOutputFile].strvalue;
 
459
  convertgis = (Boolean) myargs [c_argConvertGIs].intvalue;
 
460
  extractmrnas = (Boolean) myargs [x_argExtractMrnas].intvalue;
 
461
 
 
462
  aip = AsnIoOpen (infile, "r");
 
463
  if (aip == NULL) {
 
464
    Message (MSG_FATAL, "AsnIoOpen input file failed");
 
465
    return 1;
 
466
  }
 
467
 
 
468
  sep = SeqEntryAsnRead (aip, NULL);
 
469
  AsnIoClose (aip);
 
470
  if (sep == NULL) {
 
471
    Message (MSG_FATAL, "SeqEntryAsnRead failed");
 
472
    return 1;
 
473
  }
 
474
 
 
475
  PubSeqFetchEnable ();
 
476
 
 
477
  LookupFarSeqIDs (sep, FALSE, FALSE, FALSE, TRUE, FALSE);
 
478
 
 
479
  /* Extract mRNA hits from BLAST against nr */
 
480
 
 
481
  if (extractmrnas) {
 
482
    VisitAnnotsInSep (sep, NULL, FindBlastNR);
 
483
    DeleteMarkedObjects (0, OBJ_SEQENTRY, (Pointer) sep);
 
484
  }
 
485
 
 
486
  /* convert gi numbers to accession.version */
 
487
 
 
488
  if (convertgis) {
 
489
    VisitAlignmentsInSep (sep, NULL, UpdateAligns);
 
490
  }
 
491
 
 
492
  PubSeqFetchDisable ();
 
493
 
 
494
  BasicSeqEntryCleanup (sep);
 
495
 
 
496
  aip = AsnIoOpen (outfile, "w");
 
497
  if (aip == NULL) {
 
498
    Message (MSG_FATAL, "AsnIoOpen output file failed");
 
499
    return 1;
 
500
  }
 
501
 
 
502
  SeqEntryAsnWrite (sep, aip, NULL);
 
503
  AsnIoClose (aip);
 
504
 
 
505
  return 0;
 
506
}
 
507