← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/tools/gt_sequniq.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/tools/gt_sequniq.c
 
1
/*
 
2
  Copyright (c) 2008-2011 Gordon Gremme <gremme@zbh.uni-hamburg.de>
 
3
  Copyright (c) 2011      Giorgio Gonnella <gonnella@zbh.uni-hamburg.de>
 
4
  Copyright (c) 2008-2011 Center for Bioinformatics, University of Hamburg
 
5
 
 
6
  Permission to use, copy, modify, and distribute this software for any
 
7
  purpose with or without fee is hereby granted, provided that the above
 
8
  copyright notice and this permission notice appear in all copies.
 
9
 
 
10
  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 
11
  WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 
12
  MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 
13
  ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 
14
  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 
15
  ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 
16
  OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 
17
*/
 
18
 
 
19
#include "core/bioseq.h"
 
20
#include "core/fasta.h"
 
21
#include "core/fileutils_api.h"
 
22
#include "core/ma.h"
 
23
#include "core/outputfile.h"
 
24
#include "core/option_api.h"
 
25
#include "core/progressbar.h"
 
26
#include "core/seqiterator_sequence_buffer.h"
 
27
#include "core/string_distri.h"
 
28
#include "core/unused_api.h"
 
29
#include "extended/md5set.h"
 
30
#include "tools/gt_sequniq.h"
 
31
 
 
32
typedef struct {
 
33
  bool seqit, verbose, rev;
 
34
  unsigned long width, nofseqs;
 
35
  GtOutputFileInfo *ofi;
 
36
  GtFile *outfp;
 
37
} GtSequniqArguments;
 
38
 
 
39
static void* gt_sequniq_arguments_new(void)
 
40
{
 
41
  GtSequniqArguments *arguments = gt_calloc((size_t)1, sizeof *arguments);
 
42
  arguments->ofi = gt_outputfileinfo_new();
 
43
  return arguments;
 
44
}
 
45
 
 
46
static void gt_sequniq_arguments_delete(void *tool_arguments)
 
47
{
 
48
  GtSequniqArguments *arguments = tool_arguments;
 
49
  if (!arguments) return;
 
50
  gt_file_delete(arguments->outfp);
 
51
  gt_outputfileinfo_delete(arguments->ofi);
 
52
  gt_free(arguments);
 
53
}
 
54
 
 
55
static GtOptionParser* gt_sequniq_option_parser_new(void *tool_arguments)
 
56
{
 
57
  GtSequniqArguments *arguments = tool_arguments;
 
58
  GtOptionParser *op;
 
59
  GtOption *seqit_option, *verbose_option, *width_option, *rev_option,
 
60
           *nofseqs_option;
 
61
  gt_assert(arguments);
 
62
 
 
63
  op = gt_option_parser_new("[option ...] sequence_file [...] ",
 
64
                            "Filter out repeated sequences in given in given "
 
65
                            "sequence_file(s).");
 
66
 
 
67
  /* -seqit */
 
68
  seqit_option = gt_option_new_bool("seqit", "use sequence iterator",
 
69
                                    &arguments->seqit, false);
 
70
  gt_option_is_development_option(seqit_option);
 
71
  gt_option_parser_add_option(op, seqit_option);
 
72
 
 
73
  /* -nofseqs */
 
74
  nofseqs_option = gt_option_new_ulong("nofseqs", "number of sequences "
 
75
      "(improves efficiency)\ndefault: unspecified",
 
76
      &arguments->nofseqs, 0);
 
77
  gt_option_is_development_option(nofseqs_option);
 
78
  gt_option_hide_default(nofseqs_option);
 
79
  gt_option_parser_add_option(op, nofseqs_option);
 
80
 
 
81
  /* -rev */
 
82
  rev_option = gt_option_new_bool("rev", "filter out also sequences whose "
 
83
      "reverse complement is identical to a sequence already output",
 
84
      &arguments->rev, false);
 
85
  gt_option_parser_add_option(op, rev_option);
 
86
 
 
87
  /* -v */
 
88
  verbose_option = gt_option_new_verbose(&arguments->verbose);
 
89
  gt_option_parser_add_option(op, verbose_option);
 
90
 
 
91
  /* -width */
 
92
  width_option = gt_option_new_width(&arguments->width);
 
93
  gt_option_parser_add_option(op, width_option);
 
94
 
 
95
  gt_outputfile_register_options(op, &arguments->outfp, arguments->ofi);
 
96
 
 
97
  /* option implications */
 
98
  gt_option_imply(verbose_option, seqit_option);
 
99
 
 
100
  gt_option_parser_set_min_args(op, 1U);
 
101
  return op;
 
102
}
 
103
 
 
104
static int gt_sequniq_runner(int argc, const char **argv, int parsed_args,
 
105
                             void *tool_arguments, GtError *err)
 
106
{
 
107
  GtSequniqArguments *arguments = tool_arguments;
 
108
  unsigned long long duplicates = 0, num_of_sequences = 0;
 
109
  int i, had_err = 0;
 
110
  GtMD5Set *md5set;
 
111
 
 
112
  gt_error_check(err);
 
113
  gt_assert(arguments);
 
114
  md5set = gt_md5set_new(arguments->nofseqs);
 
115
  if (!arguments->seqit) {
 
116
    unsigned long j;
 
117
    GtBioseq *bs;
 
118
 
 
119
    for (i = parsed_args; !had_err && i < argc; i++) {
 
120
      if (!(bs = gt_bioseq_new(argv[i], err)))
 
121
        had_err = -1;
 
122
      if (!had_err) {
 
123
        GtMD5SetStatus retval;
 
124
        for (j = 0; j < gt_bioseq_number_of_sequences(bs) && !had_err; j++) {
 
125
          retval = gt_md5set_add_sequence(md5set,
 
126
              gt_bioseq_get_sequence(bs, j),
 
127
              gt_bioseq_get_sequence_length(bs, j),
 
128
              arguments->rev, err);
 
129
          if (retval == GT_MD5SET_NOT_FOUND)
 
130
            gt_fasta_show_entry(gt_bioseq_get_description(bs, j),
 
131
                                gt_bioseq_get_sequence(bs, j),
 
132
                                gt_bioseq_get_sequence_length(bs, j),
 
133
                                arguments->width, arguments->outfp);
 
134
          else if (retval != GT_MD5SET_ERROR)
 
135
            duplicates++;
 
136
          else
 
137
            had_err = -1;
 
138
          num_of_sequences++;
 
139
        }
 
140
        gt_bioseq_delete(bs);
 
141
      }
 
142
    }
 
143
  }
 
144
  else {
 
145
    GtSeqIterator *seqit;
 
146
    GtStrArray *files;
 
147
    off_t totalsize;
 
148
    const GtUchar *sequence;
 
149
    char *desc;
 
150
    unsigned long len;
 
151
 
 
152
    files = gt_str_array_new();
 
153
    for (i = parsed_args; i < argc; i++)
 
154
      gt_str_array_add_cstr(files, argv[i]);
 
155
    totalsize = gt_files_estimate_total_size(files);
 
156
    seqit = gt_seqiterator_sequence_buffer_new(files, err);
 
157
    if (!seqit)
 
158
      had_err = -1;
 
159
    if (!had_err) {
 
160
      if (arguments->verbose) {
 
161
        gt_progressbar_start(gt_seqiterator_getcurrentcounter(seqit,
 
162
                                                            (unsigned long long)
 
163
                                                            totalsize),
 
164
                             (unsigned long long) totalsize);
 
165
      }
 
166
      while (!had_err) {
 
167
        GtMD5SetStatus retval;
 
168
        if ((gt_seqiterator_next(seqit, &sequence, &len, &desc, err)) != 1)
 
169
          break;
 
170
 
 
171
        retval = gt_md5set_add_sequence(md5set, (const char*) sequence, len,
 
172
            arguments->rev, err);
 
173
        if (retval == GT_MD5SET_NOT_FOUND)
 
174
          gt_fasta_show_entry(desc, (const char*) sequence, len,
 
175
                              arguments->width, arguments->outfp);
 
176
        else if (retval != GT_MD5SET_ERROR)
 
177
          duplicates++;
 
178
        else
 
179
          had_err = -1;
 
180
        num_of_sequences++;
 
181
      }
 
182
      if (arguments->verbose)
 
183
        gt_progressbar_stop();
 
184
      gt_seqiterator_delete(seqit);
 
185
    }
 
186
    gt_str_array_delete(files);
 
187
  }
 
188
 
 
189
  /* show statistics */
 
190
  if (!had_err) {
 
191
    fprintf(stderr, "# %lu out of %lu sequences have been removed (%.3f%%)\n",
 
192
            (unsigned long)duplicates, (unsigned long)num_of_sequences,
 
193
            ((double) duplicates / (double)num_of_sequences) * 100.0);
 
194
  }
 
195
 
 
196
  gt_md5set_delete(md5set);
 
197
  return had_err;
 
198
}
 
199
 
 
200
GtTool* gt_sequniq(void)
 
201
{
 
202
  return gt_tool_new(gt_sequniq_arguments_new,
 
203
                     gt_sequniq_arguments_delete,
 
204
                     gt_sequniq_option_parser_new,
 
205
                     NULL,
 
206
                     gt_sequniq_runner);
 
207
}