← Back to branch summary

~ubuntu-branches/ubuntu/quantal/genometools/quantal-backports

« back to all changes in this revision

Viewing changes to src/match/giextract.c

  • Committer: Package Import Robot
  • Author(s): Sascha Steinbiss
  • Date: 2012-07-09 14:10:23 UTC
  • Revision ID: package-import@ubuntu.com-20120709141023-juuu4spm6chqsf9o
Tags: upstream-1.4.1
Import upstream version 1.4.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/match/giextract.c
 
1
/*
 
2
  Copyright (c) 2007-2009 Stefan Kurtz <kurtz@zbh.uni-hamburg.de>
 
3
  Copyright (c) 2007-2009 Center for Bioinformatics, University of Hamburg
 
4
 
 
5
  Permission to use, copy, modify, and distribute this software for any
 
6
  purpose with or without fee is hereby granted, provided that the above
 
7
  copyright notice and this permission notice appear in all copies.
 
8
 
 
9
  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 
10
  WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 
11
  MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 
12
  ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 
13
  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 
14
  ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 
15
  OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 
16
*/
 
17
 
 
18
#include <errno.h>
 
19
#include <stdio.h>
 
20
#include <string.h>
 
21
#include <stdlib.h>
 
22
#ifndef S_SPLINT_S
 
23
#include <ctype.h>
 
24
#include "core/fileutils.h"
 
25
#endif
 
26
#include "core/assert_api.h"
 
27
#include "core/divmodmul.h"
 
28
#include "core/error.h"
 
29
#include "core/fa.h"
 
30
#include "core/fasta.h"
 
31
#include "core/fileutils.h"
 
32
#include "core/format64.h"
 
33
#include "core/ma_api.h"
 
34
#include "core/progressbar.h"
 
35
#include "core/seqiterator_sequence_buffer.h"
 
36
#include "core/xansi_api.h"
 
37
#include "giextract.h"
 
38
#include "echoseq.h"
 
39
 
 
40
#define COMPLETE(VALUE)\
 
41
        ((VALUE)->frompos == 1UL && (VALUE)->topos == 0)
 
42
 
 
43
#define EXTRABUF 128
 
44
 
 
45
#define CHECKPOSITIVE(VAL,FORMAT,WHICH)\
 
46
        if ((VAL) <= 0)\
 
47
        {\
 
48
          gt_error_set(err,"file \"%s\", line " Formatuint64_t \
 
49
                           ": illegal format: %s element = " FORMAT \
 
50
                           " is not a positive integer",\
 
51
                           gt_str_get(fileofkeystoextract),\
 
52
                           PRINTuint64_tcast(linenum+1),\
 
53
                           WHICH,\
 
54
                           VAL);\
 
55
          haserr = true;\
 
56
        }
 
57
 
 
58
typedef struct
 
59
{
 
60
  char *fastakey;
 
61
  unsigned long frompos, topos;
 
62
  bool markhit;
 
63
} Fastakeyquery;
 
64
 
 
65
static int comparefastakeys(const void *a,const void *b)
 
66
{
 
67
  int cmp = strcmp(((Fastakeyquery *) a)->fastakey,
 
68
                   ((Fastakeyquery *) b)->fastakey);
 
69
 
 
70
  if (cmp < 0)
 
71
  {
 
72
    return -1;
 
73
  }
 
74
  if (cmp > 0)
 
75
  {
 
76
    return 1;
 
77
  }
 
78
  if (((Fastakeyquery *) a)->frompos < ((Fastakeyquery *) b)->frompos)
 
79
  {
 
80
    return -1;
 
81
  }
 
82
  if (((Fastakeyquery *) a)->frompos > ((Fastakeyquery *) b)->frompos)
 
83
  {
 
84
    return 1;
 
85
  }
 
86
  if (((Fastakeyquery *) a)->topos < ((Fastakeyquery *) b)->topos)
 
87
  {
 
88
    return -1;
 
89
  }
 
90
  if (((Fastakeyquery *) a)->topos > ((Fastakeyquery *) b)->topos)
 
91
  {
 
92
    return 1;
 
93
  }
 
94
  return 0;
 
95
}
 
96
 
 
97
static unsigned long remdupsfastakeyqueries(Fastakeyquery *fastakeyqueries,
 
98
                                            unsigned long numofqueries,
 
99
                                            bool verbose)
 
100
{
 
101
  if (numofqueries == 0)
 
102
  {
 
103
    return 0;
 
104
  } else
 
105
  {
 
106
    Fastakeyquery *storeptr, *readptr;
 
107
    unsigned long newnumofqueries;
 
108
 
 
109
    for (storeptr = fastakeyqueries, readptr = fastakeyqueries+1;
 
110
         readptr < fastakeyqueries + numofqueries;
 
111
         readptr++)
 
112
    {
 
113
      if (strcmp(storeptr->fastakey,readptr->fastakey) != 0 ||
 
114
          storeptr->frompos != readptr->frompos ||
 
115
          storeptr->topos != readptr->topos)
 
116
      {
 
117
        storeptr++;
 
118
        if (storeptr != readptr)
 
119
        {
 
120
          size_t len;
 
121
 
 
122
          storeptr->frompos = readptr->frompos;
 
123
          storeptr->topos = readptr->topos;
 
124
          len = strlen(readptr->fastakey);
 
125
          storeptr->fastakey = gt_realloc(storeptr->fastakey,
 
126
                                          sizeof (char) * (len+1));
 
127
          strcpy(storeptr->fastakey,readptr->fastakey);
 
128
        }
 
129
      }
 
130
    }
 
131
    newnumofqueries = (unsigned long) (storeptr - fastakeyqueries + 1);
 
132
    if (newnumofqueries < numofqueries)
 
133
    {
 
134
      if (verbose)
 
135
      {
 
136
        printf("# removed %lu duplicate queries\n",
 
137
                numofqueries - newnumofqueries);
 
138
      }
 
139
      for (storeptr = fastakeyqueries + newnumofqueries;
 
140
           storeptr < fastakeyqueries + numofqueries;
 
141
           storeptr++)
 
142
      {
 
143
        gt_free(storeptr->fastakey);
 
144
        storeptr->fastakey = NULL;
 
145
      }
 
146
    }
 
147
    return newnumofqueries;
 
148
  }
 
149
}
 
150
 
 
151
static void fastakeyqueries_delete(Fastakeyquery *fastakeyqueries,
 
152
                                   unsigned long numofqueries)
 
153
{
 
154
  if (fastakeyqueries != NULL)
 
155
  {
 
156
    unsigned long idx;
 
157
 
 
158
    for (idx=0; idx<numofqueries; idx++)
 
159
    {
 
160
      gt_free(fastakeyqueries[idx].fastakey);
 
161
    }
 
162
    gt_free(fastakeyqueries);
 
163
  }
 
164
}
 
165
 
 
166
static int extractkeyfromcurrentline(Fastakeyquery *fastakeyptr,
 
167
                                     unsigned long keysize,
 
168
                                     const GtStr *currentline,
 
169
                                     uint64_t linenum,
 
170
                                     const GtStr *fileofkeystoextract,
 
171
                                     GtError *err)
 
172
{
 
173
  char *lineptr = gt_str_get(currentline);
 
174
  long readlongfrompos, readlongtopos;
 
175
  size_t idx;
 
176
  bool haserr = false;
 
177
 
 
178
  for (idx = 0; lineptr[idx] != '\0' && !isspace(lineptr[idx]); idx++)
 
179
    /* nothing */ ;
 
180
  if (keysize == 0)
 
181
  {
 
182
    fastakeyptr->fastakey = gt_malloc(sizeof (char) * (idx+1));
 
183
  } else
 
184
  {
 
185
    if (idx != (size_t) keysize)
 
186
    {
 
187
      gt_error_set(err,"key \"%*.*s\" is not of size %lu",(int) idx,
 
188
                   (int) idx,lineptr,keysize);
 
189
      haserr = true;
 
190
    }
 
191
  }
 
192
  if (!haserr)
 
193
  {
 
194
    strncpy(fastakeyptr->fastakey,lineptr,idx);
 
195
    fastakeyptr->fastakey[idx] = '\0';
 
196
    fastakeyptr->frompos = 1UL;
 
197
    fastakeyptr->topos = 0;
 
198
    if (sscanf(lineptr+idx,"%ld %ld\n",&readlongfrompos,&readlongtopos) == 2)
 
199
    {
 
200
      CHECKPOSITIVE(readlongfrompos,"%ld","second");
 
201
      if (!haserr)
 
202
      {
 
203
        fastakeyptr->frompos = (unsigned long) readlongfrompos;
 
204
        CHECKPOSITIVE(readlongtopos,"%ld","third");
 
205
      }
 
206
      if (!haserr)
 
207
      {
 
208
        fastakeyptr->topos = (unsigned long) readlongtopos;
 
209
      }
 
210
    }
 
211
  }
 
212
  if (!haserr)
 
213
  {
 
214
    fastakeyptr->markhit = false;
 
215
    if (!COMPLETE(fastakeyptr) && fastakeyptr->frompos > fastakeyptr->topos)
 
216
    {
 
217
      gt_error_set(err, "file \"%s\", line " Formatuint64_t
 
218
                        "illegal format: second value "
 
219
                        "%lu is larger than third value %lu",
 
220
                        gt_str_get(fileofkeystoextract),
 
221
                        PRINTuint64_tcast(linenum+1),
 
222
                        fastakeyptr->frompos,
 
223
                        fastakeyptr->topos);
 
224
      haserr = true;
 
225
    }
 
226
  }
 
227
  return haserr ? -1 : 0;
 
228
}
 
229
 
 
230
static Fastakeyquery *readfileofkeystoextract(bool verbose,
 
231
                                              unsigned long *numofqueries,
 
232
                                              const GtStr *fileofkeystoextract,
 
233
                                              GtError *err)
 
234
{
 
235
  FILE *fp;
 
236
  GtStr *currentline;
 
237
  bool haserr = false;
 
238
  uint64_t linenum;
 
239
  Fastakeyquery *fastakeyqueries;
 
240
#undef SKDEBUG
 
241
#ifdef SKDEBUG
 
242
  unsigned long i;
 
243
#endif
 
244
 
 
245
  gt_error_check(err);
 
246
  *numofqueries = gt_file_number_of_lines(gt_str_get(fileofkeystoextract));
 
247
  if (*numofqueries == 0)
 
248
  {
 
249
    gt_error_set(err,"empty file \"%s\" not allowed",
 
250
                 gt_str_get(fileofkeystoextract));
 
251
    return NULL;
 
252
  }
 
253
  fp = gt_fa_fopen(gt_str_get(fileofkeystoextract),"r",err);
 
254
  if (fp == NULL)
 
255
  {
 
256
    return NULL;
 
257
  }
 
258
  if (verbose)
 
259
  {
 
260
    printf("# opened keyfile \"%s\"\n",gt_str_get(fileofkeystoextract));
 
261
  }
 
262
  fastakeyqueries = gt_malloc(sizeof (*fastakeyqueries) * (*numofqueries));
 
263
  currentline = gt_str_new();
 
264
  for (linenum = 0; gt_str_read_next_line(currentline, fp) != EOF; linenum++)
 
265
  {
 
266
    if (extractkeyfromcurrentline(fastakeyqueries + linenum,
 
267
                                  0,
 
268
                                  currentline,
 
269
                                  linenum,
 
270
                                  fileofkeystoextract,
 
271
                                  err) != 0)
 
272
    {
 
273
      haserr = true;
 
274
      break;
 
275
    }
 
276
    gt_str_reset(currentline);
 
277
  }
 
278
  gt_str_delete(currentline);
 
279
  gt_fa_fclose(fp);
 
280
  if (haserr)
 
281
  {
 
282
    fastakeyqueries_delete(fastakeyqueries,*numofqueries);
 
283
    return NULL;
 
284
  }
 
285
  qsort(fastakeyqueries,(size_t) *numofqueries,sizeof (*fastakeyqueries),
 
286
        comparefastakeys);
 
287
  if (verbose)
 
288
  {
 
289
    printf("# %lu fastakey-queries successfully parsed and sorted\n",
 
290
            *numofqueries);
 
291
  }
 
292
  *numofqueries = remdupsfastakeyqueries(fastakeyqueries,*numofqueries,verbose);
 
293
#ifdef SKDEBUG
 
294
  for (i=0; i<*numofqueries; i++)
 
295
  {
 
296
    printf("%lu %s\n",i,fastakeyqueries[i].fastakey);
 
297
  }
 
298
#endif
 
299
  return fastakeyqueries;
 
300
}
 
301
 
 
302
static unsigned long searchdesinfastakeyqueries(const char *extractkey,
 
303
                                                const Fastakeyquery
 
304
                                                  *fastakeyqueries,
 
305
                                                unsigned long numofqueries)
 
306
{
 
307
  const Fastakeyquery *leftptr, *rightptr, *midptr;
 
308
  int cmp;
 
309
 
 
310
  leftptr = fastakeyqueries;
 
311
  rightptr = fastakeyqueries + numofqueries - 1;
 
312
  while (leftptr <= rightptr)
 
313
  {
 
314
    midptr = leftptr + GT_DIV2((unsigned long) (rightptr-leftptr));
 
315
    cmp = strcmp(extractkey,midptr->fastakey);
 
316
    if (cmp == 0)
 
317
    {
 
318
      if (midptr > fastakeyqueries &&
 
319
          strcmp(extractkey,(midptr-1)->fastakey) == 0)
 
320
      {
 
321
        rightptr = midptr - 1;
 
322
      } else
 
323
      {
 
324
        return (unsigned long) (midptr - fastakeyqueries);
 
325
      }
 
326
    } else
 
327
    {
 
328
      if (cmp < 0)
 
329
      {
 
330
        rightptr = midptr-1;
 
331
      } else
 
332
      {
 
333
        leftptr = midptr + 1;
 
334
      }
 
335
    }
 
336
  }
 
337
  return numofqueries;
 
338
}
 
339
 
 
340
static void outputnonmarked(const Fastakeyquery *fastakeyqueries,
 
341
                            unsigned long numofqueries)
 
342
{
 
343
  unsigned long idx, countmissing = 0;
 
344
 
 
345
  for (idx=0; idx<numofqueries; idx++)
 
346
  {
 
347
    if (!fastakeyqueries[idx].markhit)
 
348
    {
 
349
      printf("unsatisfied %s",fastakeyqueries[idx].fastakey);
 
350
      if (COMPLETE(fastakeyqueries + idx))
 
351
      {
 
352
        printf(" complete\n");
 
353
      } else
 
354
      {
 
355
        printf(" %lu %lu\n",fastakeyqueries[idx].frompos,
 
356
                            fastakeyqueries[idx].topos);
 
357
      }
 
358
      countmissing++;
 
359
    }
 
360
  }
 
361
  printf("# number of unsatified fastakey-queries: %lu\n",countmissing);
 
362
}
 
363
 
 
364
static const char *desc2key(unsigned long *keylen,const char *desc,
 
365
                            GtError *err)
 
366
{
 
367
  unsigned long i, firstpipe = 0, secondpipe = 0;
 
368
 
 
369
  gt_error_check(err);
 
370
  for (i=0; desc[i] != '\0'; i++)
 
371
  {
 
372
    if (desc[i] == '|')
 
373
    {
 
374
      if (firstpipe > 0)
 
375
      {
 
376
        gt_assert(i>0);
 
377
        secondpipe = i;
 
378
        break;
 
379
      }
 
380
      gt_assert(i>0);
 
381
      firstpipe = i;
 
382
    }
 
383
  }
 
384
  if (firstpipe == 0 || secondpipe == 0)
 
385
  {
 
386
    gt_error_set(err,"Cannot find key in description \"%s\"",desc);
 
387
    return NULL;
 
388
  }
 
389
  gt_assert(firstpipe < secondpipe);
 
390
  *keylen = secondpipe - firstpipe - 1;
 
391
  return desc + firstpipe + 1;
 
392
}
 
393
 
 
394
static int giextract_encodedseq2fasta(FILE *fpout,
 
395
                                      const GtEncseq *encseq,
 
396
                                      unsigned long seqnum,
 
397
                                      const Fastakeyquery *fastakeyquery,
 
398
                                      unsigned long linewidth,
 
399
                                      GT_UNUSED GtError *err)
 
400
{
 
401
  const char *desc;
 
402
  unsigned long desclen;
 
403
  bool haserr = false;
 
404
 
 
405
  desc = gt_encseq_description(encseq, &desclen, seqnum);
 
406
  gt_xfputc('>',fpout);
 
407
  if (fastakeyquery != NULL && !COMPLETE(fastakeyquery))
 
408
  {
 
409
    printf("%s %lu %lu ",fastakeyquery->fastakey,
 
410
                         fastakeyquery->frompos,
 
411
                         fastakeyquery->topos);
 
412
  }
 
413
  gt_xfwrite(desc,sizeof *desc,(size_t) desclen,fpout);
 
414
  if (!haserr)
 
415
  {
 
416
    unsigned long frompos, topos, seqstartpos, seqlength ;
 
417
 
 
418
    gt_xfputc('\n',fpout);
 
419
    seqstartpos = gt_encseq_seqstartpos(encseq, seqnum);
 
420
    seqlength = gt_encseq_seqlength(encseq, seqnum);
 
421
    if (fastakeyquery != NULL && !COMPLETE(fastakeyquery))
 
422
    {
 
423
      frompos = fastakeyquery->frompos-1;
 
424
      topos = fastakeyquery->topos - fastakeyquery->frompos + 1;
 
425
    } else
 
426
    {
 
427
      frompos = 0;
 
428
      topos = seqlength;
 
429
    }
 
430
    gt_encseq2symbolstring(fpout,
 
431
                           encseq,
 
432
                           GT_READMODE_FORWARD,
 
433
                           seqstartpos + frompos,
 
434
                           topos,
 
435
                           linewidth);
 
436
  }
 
437
  return haserr ? -1 : 0;
 
438
}
 
439
 
 
440
#define MAXFIXEDKEYSIZE 11
 
441
 
 
442
typedef struct
 
443
{
 
444
  char key[MAXFIXEDKEYSIZE+1];
 
445
  unsigned long seqnum;
 
446
} Fixedsizekey;
 
447
 
 
448
static int compareFixedkeys(const void *a,const void *b)
 
449
{
 
450
  const Fixedsizekey *fa = a;
 
451
  const Fixedsizekey *fb = b;
 
452
  return strcmp(((const Fixedsizekey *) fa)->key,
 
453
                ((const Fixedsizekey *) fb)->key);
 
454
}
 
455
 
 
456
#define GT_KEYSTABFILESUFFIX ".kys"
 
457
 
 
458
int gt_extractkeysfromdesfile(const char *indexname,
 
459
                              bool sortkeys,
 
460
                              GtLogger *logger,
 
461
                              GtError *err)
 
462
{
 
463
  FILE *fpin, *fpout = NULL;
 
464
  GtStr *line = NULL;
 
465
  const char *keyptr;
 
466
  unsigned long keylen, constantkeylen = 0, linenum;/* incorrectorder = 0;*/
 
467
  bool haserr = false, firstdesc = true;
 
468
  char *previouskey = NULL;
 
469
  Fixedsizekey *keytab = NULL, *keytabptr = NULL;
 
470
  GtEncseq *encseq = NULL;
 
471
  unsigned long numofentries = 0;
 
472
  const unsigned long linewidth = 60UL;
 
473
 
 
474
  fpin = gt_fa_fopen_with_suffix(indexname,GT_DESTABFILESUFFIX,"rb",err);
 
475
  if (fpin == NULL)
 
476
  {
 
477
    return -1;
 
478
  }
 
479
  if (!sortkeys)
 
480
  {
 
481
    fpout = gt_fa_fopen_with_suffix(indexname,GT_KEYSTABFILESUFFIX,"wb",err);
 
482
    if (fpout == NULL)
 
483
    {
 
484
      haserr = true;
 
485
    }
 
486
  }
 
487
  if (!haserr)
 
488
  {
 
489
    line = gt_str_new();
 
490
  }
 
491
  for (linenum = 0; !haserr && gt_str_read_next_line(line, fpin) != EOF;
 
492
       linenum++)
 
493
  {
 
494
    keyptr = desc2key(&keylen,gt_str_get(line),err);
 
495
    if (keyptr == NULL)
 
496
    {
 
497
      haserr = true;
 
498
      break;
 
499
    }
 
500
    if (keylen == 0)
 
501
    {
 
502
      gt_error_set(err,"key of length 0 in \"%s\" not expected",
 
503
                   gt_str_get(line));
 
504
      haserr = true;
 
505
      break;
 
506
    }
 
507
    if (firstdesc)
 
508
    {
 
509
      if (keylen > (unsigned long) CHAR_MAX)
 
510
      {
 
511
        gt_error_set(err,"key \"%*.*s\" of length %lu not allowed; "
 
512
                         "no key must be larger than %d",
 
513
                          (int) keylen,(int) keylen,keyptr,keylen,CHAR_MAX);
 
514
        haserr = true;
 
515
        break;
 
516
      }
 
517
      constantkeylen = keylen;
 
518
      previouskey = gt_malloc(sizeof (char) * (constantkeylen+1));
 
519
      firstdesc = false;
 
520
      if (!sortkeys)
 
521
      {
 
522
        gt_xfputc((char) constantkeylen,fpout);
 
523
      } else
 
524
      {
 
525
        GtEncseqLoader *el;
 
526
        if (constantkeylen > (unsigned long) MAXFIXEDKEYSIZE)
 
527
        {
 
528
          gt_error_set(err,"key \"%*.*s\" of length %lu not allowed; "
 
529
                           "no key must be larger than %d",
 
530
                            (int) keylen,(int) keylen,keyptr,keylen,
 
531
                            MAXFIXEDKEYSIZE);
 
532
          haserr = true;
 
533
          break;
 
534
        }
 
535
        el = gt_encseq_loader_new();
 
536
        gt_encseq_loader_set_logger(el, logger);
 
537
        encseq = gt_encseq_loader_load(el, indexname, err);
 
538
        gt_encseq_loader_delete(el);
 
539
        if (encseq == NULL)
 
540
        {
 
541
          haserr = true;
 
542
          break;
 
543
        }
 
544
        numofentries = gt_encseq_num_of_sequences(encseq);
 
545
        gt_assert(numofentries > 0);
 
546
        keytab = gt_malloc(sizeof (*keytab) * numofentries);
 
547
        keytabptr = keytab;
 
548
      }
 
549
    } else
 
550
    {
 
551
      if (constantkeylen != keylen)
 
552
      {
 
553
        gt_error_set(err,"key \"%*.*s\" of length %lu: all keys must be of "
 
554
                         "the same length which for all previously seen "
 
555
                         "headers is %lu",
 
556
                         (int) keylen,(int) keylen,keyptr,keylen,
 
557
                         constantkeylen);
 
558
        haserr = true;
 
559
        break;
 
560
      }
 
561
      gt_assert(previouskey != NULL);
 
562
      if (!sortkeys && strncmp(previouskey,keyptr,(size_t) constantkeylen) >= 0)
 
563
      {
 
564
        gt_error_set(err,"previous key \"%s\" is not lexicographically smaller "
 
565
                         "than current key \"%*.*s\"",
 
566
                         previouskey,(int) keylen,(int) keylen,keyptr);
 
567
        haserr = true;
 
568
        break;
 
569
        /*
 
570
        printf("previous key \"%s\" (no %lu) is lexicographically larger "
 
571
               "than current key \"%*.*s\"\n",
 
572
               previouskey,linenum,(int) keylen,(int) keylen,keyptr);
 
573
        incorrectorder++;
 
574
        */
 
575
      }
 
576
    }
 
577
    if (!sortkeys)
 
578
    {
 
579
      gt_xfwrite(keyptr,sizeof *keyptr,(size_t) keylen,fpout);
 
580
      gt_xfputc('\0',fpout);
 
581
    } else
 
582
    {
 
583
      gt_assert(keytabptr != NULL);
 
584
      strncpy(keytabptr->key,keyptr,(size_t) constantkeylen);
 
585
      keytabptr->key[constantkeylen] = '\0';
 
586
      keytabptr->seqnum = linenum;
 
587
      keytabptr++;
 
588
    }
 
589
    strncpy(previouskey,keyptr,(size_t) constantkeylen);
 
590
    previouskey[constantkeylen] = '\0';
 
591
    gt_str_reset(line);
 
592
  }
 
593
  if (!haserr)
 
594
  {
 
595
    gt_logger_log(logger,"number of keys of length %lu = %lu",
 
596
                constantkeylen,linenum);
 
597
    /*
 
598
    gt_logger_log(logger,"number of incorrectly ordered keys = %lu",
 
599
                incorrectorder);
 
600
    */
 
601
  }
 
602
  gt_str_delete(line);
 
603
  gt_fa_fclose(fpin);
 
604
  gt_fa_fclose(fpout);
 
605
  gt_free(previouskey);
 
606
  if (!haserr && sortkeys)
 
607
  {
 
608
    gt_assert(keytabptr != NULL);
 
609
    gt_assert(numofentries > 0);
 
610
    gt_assert(keytabptr == keytab + numofentries);
 
611
    qsort(keytab,(size_t) numofentries,sizeof (*keytab),compareFixedkeys);
 
612
    gt_assert(keytabptr != NULL);
 
613
    for (keytabptr = keytab; !haserr && keytabptr < keytab + numofentries;
 
614
         keytabptr++)
 
615
    {
 
616
      if (giextract_encodedseq2fasta(stdout,
 
617
                                     encseq,
 
618
                                     keytabptr->seqnum,
 
619
                                     NULL,
 
620
                                     linewidth,
 
621
                                     err) != 0)
 
622
      {
 
623
        haserr = true;
 
624
        break;
 
625
      }
 
626
    }
 
627
  }
 
628
  if (encseq != NULL)
 
629
  {
 
630
    gt_encseq_delete(encseq);
 
631
    encseq = NULL;
 
632
  }
 
633
  gt_free(keytab);
 
634
  return haserr ? -1 : 0;
 
635
}
 
636
 
 
637
bool gt_deskeysfileexists(const char *indexname)
 
638
{
 
639
  return gt_file_with_suffix_exists(indexname,GT_KEYSTABFILESUFFIX);
 
640
}
 
641
 
 
642
static unsigned long searchfastaqueryindes(const char *extractkey,
 
643
                                           const char *keytab,
 
644
                                           unsigned long numofkeys,
 
645
                                           unsigned long keysize)
 
646
{
 
647
  unsigned long left = 0, right = numofkeys - 1, mid;
 
648
  int cmp;
 
649
 
 
650
  while (left <= right)
 
651
  {
 
652
    mid = left + GT_DIV2((unsigned long) (right-left));
 
653
    cmp = strcmp(extractkey,keytab + 1UL + mid * (keysize+1));
 
654
    if (cmp < 0)
 
655
    {
 
656
      gt_assert(mid > 0);
 
657
      right = mid-1;
 
658
    } else
 
659
    {
 
660
      if (cmp > 0)
 
661
      {
 
662
        left = mid+1;
 
663
      } else
 
664
      {
 
665
        gt_assert(mid < numofkeys);
 
666
        return mid;
 
667
      }
 
668
    }
 
669
  }
 
670
  return numofkeys;
 
671
}
 
672
 
 
673
static int itersearchoverallkeys(const GtEncseq *encseq,
 
674
                                 const char *keytab,
 
675
                                 unsigned long numofkeys,
 
676
                                 unsigned long keysize,
 
677
                                 const GtStr *fileofkeystoextract,
 
678
                                 unsigned long linewidth,
 
679
                                 GtError *err)
 
680
{
 
681
  FILE *fp;
 
682
  GtStr *currentline;
 
683
  uint64_t linenum;
 
684
  unsigned long seqnum, countmissing = 0;
 
685
  bool haserr = false;
 
686
  Fastakeyquery fastakeyquery;
 
687
 
 
688
  if (linewidth == 0)
 
689
  {
 
690
    gt_error_set(err,"use option width to specify line width for formatting");
 
691
    return -1;
 
692
  }
 
693
  fp = gt_fa_fopen(gt_str_get(fileofkeystoextract),"r",err);
 
694
  if (fp == NULL)
 
695
  {
 
696
    return -1;
 
697
  }
 
698
  currentline = gt_str_new();
 
699
  fastakeyquery.fastakey = gt_malloc(sizeof (char) * (keysize+1));
 
700
  for (linenum = 0; gt_str_read_next_line(currentline, fp) != EOF; linenum++)
 
701
  {
 
702
    if (extractkeyfromcurrentline(&fastakeyquery,
 
703
                                  keysize,
 
704
                                  currentline,
 
705
                                  linenum,
 
706
                                  fileofkeystoextract,
 
707
                                  err) != 0)
 
708
    {
 
709
      haserr = true;
 
710
      break;
 
711
    }
 
712
    seqnum = searchfastaqueryindes(fastakeyquery.fastakey,keytab,numofkeys,
 
713
                                   keysize);
 
714
    if (seqnum < numofkeys)
 
715
    {
 
716
      if (giextract_encodedseq2fasta(stdout,
 
717
                                     encseq,
 
718
                                     seqnum,
 
719
                                     &fastakeyquery,
 
720
                                     linewidth,
 
721
                                     err) != 0)
 
722
      {
 
723
        haserr = true;
 
724
        break;
 
725
      }
 
726
    } else
 
727
    {
 
728
      countmissing++;
 
729
    }
 
730
    gt_str_reset(currentline);
 
731
  }
 
732
  if (!haserr && countmissing > 0)
 
733
  {
 
734
    printf("# number of unsatified fastakey-queries: %lu\n",countmissing);
 
735
  }
 
736
  gt_str_delete(currentline);
 
737
  gt_fa_fclose(fp);
 
738
  gt_free(fastakeyquery.fastakey);
 
739
  return haserr ? - 1 : 0;
 
740
}
 
741
 
 
742
static int readkeysize(const char *indexname,GtError *err)
 
743
{
 
744
  FILE *fp;
 
745
  bool haserr = false;
 
746
  char cc;
 
747
 
 
748
  gt_error_check(err);
 
749
  fp = gt_fa_fopen_with_suffix(indexname,GT_KEYSTABFILESUFFIX,"rb",err);
 
750
  if (fp == NULL)
 
751
  {
 
752
    haserr = true;
 
753
  }
 
754
  if (!haserr)
 
755
  {
 
756
    GT_UNUSED size_t ret;
 
757
 
 
758
    ret = fread(&cc,sizeof cc, (size_t) 1, fp);
 
759
    if (ferror(fp))
 
760
    {
 
761
      gt_error_set(err,"error when trying to read first byte of file %s%s: %s",
 
762
                   indexname,GT_KEYSTABFILESUFFIX,strerror(errno));
 
763
      haserr = true;
 
764
    }
 
765
  }
 
766
  gt_assert(cc >= 0);
 
767
  gt_fa_xfclose(fp);
 
768
  return haserr ? -1 : (int) cc;
 
769
}
 
770
 
 
771
int gt_extractkeysfromfastaindex(const char *indexname,
 
772
                                 const GtStr *fileofkeystoextract,
 
773
                                 unsigned long linewidth,GtError *err)
 
774
{
 
775
  GtEncseq *encseq = NULL;
 
776
  GtEncseqLoader *el = NULL;
 
777
  bool haserr = false;
 
778
  unsigned long numofdbsequences = 0, keysize = 0;
 
779
 
 
780
  el = gt_encseq_loader_new();
 
781
  encseq = gt_encseq_loader_load(el, indexname, err);
 
782
  gt_encseq_loader_delete(el);
 
783
  if (encseq == NULL)
 
784
  {
 
785
    haserr = true;
 
786
  }
 
787
  if (!haserr)
 
788
  {
 
789
    int retval;
 
790
 
 
791
    numofdbsequences = gt_encseq_num_of_sequences(encseq);
 
792
    retval = readkeysize(indexname,err);
 
793
    if (retval < 0)
 
794
    {
 
795
      haserr = true;
 
796
    }
 
797
    keysize = (unsigned long) retval;
 
798
  }
 
799
  if (!haserr)
 
800
  {
 
801
    char *keytab;
 
802
    unsigned long keytablength;
 
803
 
 
804
    keytablength = 1UL + numofdbsequences * (keysize+1);
 
805
    keytab = gt_fa_mmap_check_size_with_suffix(indexname,
 
806
                                               GT_KEYSTABFILESUFFIX,
 
807
                                               keytablength,
 
808
                                               sizeof (GtUchar),
 
809
                                               err);
 
810
    if (keytab == NULL)
 
811
    {
 
812
      haserr = true;
 
813
    } else
 
814
    {
 
815
      if (itersearchoverallkeys(encseq,keytab,numofdbsequences,
 
816
                                keysize,fileofkeystoextract,
 
817
                                linewidth,err) != 0)
 
818
      {
 
819
        haserr = true;
 
820
      }
 
821
    }
 
822
    gt_fa_xmunmap(keytab);
 
823
  }
 
824
  if (encseq != NULL)
 
825
  {
 
826
    gt_encseq_delete(encseq);
 
827
    encseq = NULL;
 
828
  }
 
829
  return haserr ? -1 : 0;
 
830
}
 
831
 
 
832
int gt_extractkeysfromfastafile(bool verbose,
 
833
                                GtFile *outfp,
 
834
                                unsigned long width,
 
835
                                const GtStr *fileofkeystoextract,
 
836
                                GtStrArray *referencefiletab,
 
837
                                GtError *err)
 
838
{
 
839
  GtSeqIterator *seqit;
 
840
  const GtUchar *sequence;
 
841
  char *desc, *headerbufferspace = NULL, *keyspace = NULL;
 
842
  const char *keyptr;
 
843
  unsigned long allockeyspace = 0, len, keylen, numofqueries, keyposition,
 
844
                countmarkhit = 0;
 
845
  int had_err = 0;
 
846
  off_t totalsize;
 
847
  Fastakeyquery *fastakeyqueries;
 
848
  size_t headerbuffersize = 0, headerlength;
 
849
 
 
850
  gt_error_check(err);
 
851
  fastakeyqueries = readfileofkeystoextract(verbose,&numofqueries,
 
852
                                            fileofkeystoextract,err);
 
853
  if (fastakeyqueries == NULL)
 
854
  {
 
855
    return -1;
 
856
  }
 
857
  totalsize = gt_files_estimate_total_size(referencefiletab);
 
858
  if (verbose)
 
859
  {
 
860
    printf("# estimated total size is " Formatuint64_t "\n",
 
861
            PRINTuint64_tcast(totalsize));
 
862
  }
 
863
  seqit = gt_seqiterator_sequence_buffer_new(referencefiletab, err);
 
864
  if (!seqit)
 
865
  {
 
866
    had_err = -1;
 
867
  }
 
868
  if (!had_err && verbose)
 
869
  {
 
870
    gt_progressbar_start(gt_seqiterator_getcurrentcounter(seqit,
 
871
                                                          (unsigned long long)
 
872
                                                          totalsize),
 
873
                                                          (unsigned long long)
 
874
                                                          totalsize);
 
875
  }
 
876
  while (had_err != -1 && countmarkhit < numofqueries)
 
877
  {
 
878
    had_err = gt_seqiterator_next(seqit, &sequence, &len, &desc, err);
 
879
    if (had_err != 1)
 
880
    {
 
881
      break;
 
882
    }
 
883
    keyptr = desc2key(&keylen,desc,err);
 
884
    if (keyptr == NULL)
 
885
    {
 
886
      had_err = -1;
 
887
    } else
 
888
    {
 
889
      if (allockeyspace < keylen)
 
890
      {
 
891
        keyspace = gt_realloc(keyspace,sizeof (*keyspace) * (keylen+1));
 
892
        allockeyspace = keylen;
 
893
      }
 
894
      gt_assert(keyspace != NULL);
 
895
      strncpy(keyspace,keyptr,(size_t) keylen);
 
896
      keyspace[keylen] = '\0';
 
897
      keyposition = searchdesinfastakeyqueries(keyspace,fastakeyqueries,
 
898
                                               numofqueries);
 
899
      if (keyposition < numofqueries)
 
900
      {
 
901
        while (keyposition < numofqueries &&
 
902
               strcmp(fastakeyqueries[keyposition].fastakey,keyspace) == 0)
 
903
        {
 
904
#ifndef NDEBUG
 
905
          if (fastakeyqueries[keyposition].markhit)
 
906
          {
 
907
            fprintf(stderr,"key %s was already found before\n",
 
908
                     fastakeyqueries[keyposition].fastakey);
 
909
            exit(GT_EXIT_PROGRAMMING_ERROR);
 
910
          }
 
911
#endif
 
912
          headerlength = strlen(desc);
 
913
          if (headerbuffersize < headerlength + EXTRABUF + 1)
 
914
          {
 
915
            headerbuffersize = headerlength + EXTRABUF + 1;
 
916
            headerbufferspace = gt_realloc(headerbufferspace,
 
917
                                           sizeof (*headerbufferspace)
 
918
                                           * headerbuffersize);
 
919
          }
 
920
          if (COMPLETE(fastakeyqueries + keyposition))
 
921
          {
 
922
            /*
 
923
            (void) snprintf(headerbufferspace,headerbuffersize,
 
924
                            "%*.*s complete %s",
 
925
                            (int) keylen,(int) keylen,keyspace,
 
926
                            desc);
 
927
            */
 
928
            gt_fasta_show_entry(desc, (const char *) sequence, len, width,
 
929
                                outfp);
 
930
          } else
 
931
          {
 
932
            (void) snprintf(headerbufferspace,headerbuffersize,
 
933
                            "%*.*s %lu %lu %s",
 
934
                            (int) keylen,(int) keylen,keyspace,
 
935
                            fastakeyqueries[keyposition].frompos,
 
936
                            fastakeyqueries[keyposition].topos,
 
937
                            desc);
 
938
            gt_fasta_show_entry(headerbufferspace,
 
939
                                (const char *)
 
940
                                (sequence+fastakeyqueries[keyposition].
 
941
                                                          frompos - 1),
 
942
                                fastakeyqueries[keyposition].topos -
 
943
                                fastakeyqueries[keyposition].frompos+1,
 
944
                                width, outfp);
 
945
          }
 
946
          fastakeyqueries[keyposition].markhit = true;
 
947
          countmarkhit++;
 
948
          keyposition++;
 
949
        }
 
950
      }
 
951
#ifdef SKDEBUG
 
952
      printf("%s 1 %lu\n",keyspace, len);
 
953
#endif
 
954
    }
 
955
  }
 
956
  gt_free(headerbufferspace);
 
957
  gt_free(keyspace);
 
958
  if (verbose)
 
959
  {
 
960
    gt_progressbar_stop();
 
961
  }
 
962
  if (verbose)
 
963
  {
 
964
    outputnonmarked(fastakeyqueries,numofqueries);
 
965
  }
 
966
  fastakeyqueries_delete(fastakeyqueries,numofqueries);
 
967
  gt_seqiterator_delete(seqit);
 
968
  return had_err;
 
969
}