~ubuntu-branches/ubuntu/oneiric/bioperl/oneiric

« back to all changes in this revision

Viewing changes to Bio/DB/SeqFeature/Store/GFF3Loader.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2007-09-21 22:52:22 UTC
  • mfrom: (1.2.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20070921225222-tt20m2yy6ycuy2d8
Tags: 1.5.2.102-1
* Developer release.
* Upgraded source package to debhelper 5 and standards-version 3.7.2.
* Added libmodule-build-perl and libtest-harness-perl to
  build-depends-indep.
* Disabled automatic CRAN download.
* Using quilt instead of .diff.gz to manage modifications.
* Updated Recommends list for the binary package.
* Moved the "production-quality" scripts to /usr/bin/.
* New maintainer: Debian-Med packaging team mailing list.
* New uploaders: Charles Plessy and Steffen Moeller.
* Updated Depends, Recommends and Suggests.
* Imported in Debian-Med's SVN repository on Alioth.
* Executing the regression tests during package building.
* Moved the Homepage: field out from the package's description.
* Updated watch file.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
package Bio::DB::SeqFeature::Store::GFF3Loader;
 
2
 
 
3
# $Id: GFF3Loader.pm,v 1.15.4.3 2006/10/02 23:10:17 sendu Exp $
 
4
 
 
5
=head1 NAME
 
6
 
 
7
Bio::DB::SeqFeature::Store::GFF3Loader -- GFF3 file loader for Bio::DB::SeqFeature::Store
 
8
 
 
9
=head1 SYNOPSIS
 
10
 
 
11
  use Bio::DB::SeqFeature::Store;
 
12
 
 
13
  # Open the sequence database
 
14
  my $db      = Bio::DB::SeqFeature::Store->new( -adaptor => 'DBI::mysql',
 
15
                                                 -dsn     => 'dbi:mysql:test',
 
16
                                                 -write   => 1 );
 
17
 
 
18
  my $loader = Bio::DB::SeqFeature::Store::GFF3Loader->new(-store    => $db,
 
19
                                                           -verbose  => 1,
 
20
                                                           -fast     => 1);
 
21
 
 
22
  $loader->load('./my_genome.gff3');
 
23
 
 
24
 
 
25
=head1 DESCRIPTION
 
26
 
 
27
The Bio::DB::SeqFeature::Store::GFF3Loader object parsers GFF3-format
 
28
sequence annotation files and loads Bio::DB::SeqFeature::Store
 
29
databases. For certain combinations of SeqFeature classes and
 
30
SeqFeature::Store databases it features a "fast load" mode which will
 
31
greatly accelerate the loading of GFF3 databases by a factor of 5-10.
 
32
 
 
33
The GFF3 file format has been extended very slightly to accomodate
 
34
Bio::DB::SeqFeature::Store. First, the loader recognizes is a new
 
35
directive:
 
36
 
 
37
  # #index-subfeatures [0|1]
 
38
 
 
39
Note that you can place a space between the two #'s in order to
 
40
prevent GFF3 validators from complaining.
 
41
 
 
42
If this is true, then subfeatures are indexed (the default) so that
 
43
they can be retrieved with a query. See L<Bio::DB::SeqFeature::Store>
 
44
for an explanation of this. If false, then subfeatures can only be
 
45
accessed through their parent feature. The default is to index all
 
46
subfeatures.
 
47
 
 
48
Second, the loader recognizes a new attribute tag called index, which
 
49
if present, controls indexing of the current feature. Example:
 
50
 
 
51
 ctg123 . TF_binding_site 1000 1012 . + . ID=tfbs00001;index=1
 
52
 
 
53
You can use this to turn indexing on and off, overriding the default
 
54
for a particular feature.
 
55
 
 
56
=cut
 
57
 
 
58
 
 
59
# load utility - incrementally load the store based on GFF3 file
 
60
#
 
61
# two modes:
 
62
#   slow mode -- features can occur in any order in the GFF3 file
 
63
#   fast mode -- all features with same ID must be contiguous in GFF3 file
 
64
 
 
65
use strict;
 
66
use Carp 'croak';
 
67
use IO::File;
 
68
use Bio::DB::GFF::Util::Rearrange;
 
69
use Bio::DB::SeqFeature::Store;
 
70
use File::Spec;
 
71
use base 'Bio::Root::Root';
 
72
 
 
73
use constant DEFAULT_SEQ_CHUNK_SIZE => 2000;
 
74
 
 
75
my %Special_attributes =(
 
76
                         Gap    => 1, Target => 1,
 
77
                         Parent => 1, Name   => 1,
 
78
                         Alias  => 1, ID     => 1,
 
79
                         index  => 1, Index  => 1,
 
80
                        );
 
81
my %Strandedness = ( '+'  => 1,
 
82
                     '-'  => -1,
 
83
                     '.'  => 0,
 
84
                     ''   => 0,
 
85
                     0    => 0,
 
86
                     1    => 1,
 
87
                     -1   => -1,
 
88
                     +1   => 1,
 
89
                     undef => 0,
 
90
                   );
 
91
 
 
92
=head2 new
 
93
 
 
94
 Title   : new
 
95
 Usage   : $loader = Bio::DB::SeqFeature::Store::GFF3Loader->new(@options)
 
96
 Function: create a new parser
 
97
 Returns : a Bio::DB::SeqFeature::Store::GFF3Loader gff3 parser and loader
 
98
 Args    : several - see below
 
99
 Status  : public
 
100
 
 
101
This method creates a new GFF3 loader and establishes its connection
 
102
with a Bio::DB::SeqFeature::Store database. Arguments are -name=E<gt>$value
 
103
pairs as described in this table:
 
104
 
 
105
 Name               Value
 
106
 ----               -----
 
107
 
 
108
 -store             A writeable Bio::DB::SeqFeature::Store database handle.
 
109
 
 
110
 -seqfeature_class  The name of the type of Bio::SeqFeatureI object to create
 
111
                      and store in the database (Bio::DB::SeqFeature by default)
 
112
 
 
113
 -sf_class          A shorter alias for -seqfeature_class
 
114
 
 
115
 -verbose           Send progress information to standard error.
 
116
 
 
117
 -fast              If true, activate fast loading (see below)
 
118
 
 
119
 -chunk_size        Set the storage chunk size for nucleotide/protein sequences
 
120
                       (default 2000 bytes)
 
121
 
 
122
 -tmp               Indicate a temporary directory to use when loading non-normalized
 
123
                       features.
 
124
 
 
125
When you call new(), a connection to a Bio::DB::SeqFeature::Store
 
126
database should already have been established and the database
 
127
initialized (if appropriate).
 
128
 
 
129
Some combinations of Bio::SeqFeatures and Bio::DB::SeqFeature::Store
 
130
databases support a fast loading mode. Currently the only reliable
 
131
implementation of fast loading is the combination of DBI::mysql with
 
132
Bio::DB::SeqFeature. The other important restriction on fast loading
 
133
is the requirement that a feature that contains subfeatures must occur
 
134
in the GFF3 file before any of its subfeatures. Otherwise the
 
135
subfeatures that occurred before the parent feature will not be
 
136
attached to the parent correctly. This restriction does not apply to
 
137
normal (slow) loading.
 
138
 
 
139
If you use an unnormalized feature class, such as
 
140
Bio::SeqFeature::Generic, then the loader needs to create a temporary
 
141
database in which to cache features until all their parts and subparts
 
142
have been seen. This temporary databases uses the "bdb" adaptor. The
 
143
-tmp option specifies the directory in which that database will be
 
144
created. If not present, it defaults to the system default tmp
 
145
directory specified by File::Spec-E<gt>tmpdir().
 
146
 
 
147
The -chunk_size option allows you to tune the representation of
 
148
DNA/Protein sequence in the Store database. By default, sequences are
 
149
split into 2000 base/residue chunks and then reassembled as
 
150
needed. This avoids the problem of pulling a whole chromosome into
 
151
memory in order to fetch a short subsequence from somewhere in the
 
152
middle. Depending on your usage patterns, you may wish to tune this
 
153
parameter using a chunk size that is larger or smaller than the
 
154
default.
 
155
 
 
156
=cut
 
157
 
 
158
sub new {
 
159
  my $self = shift;
 
160
  my ($store,$seqfeature_class,$tmpdir,$verbose,$fast,$seq_chunk_size) = rearrange(['STORE',
 
161
                                                                                    ['SF_CLASS','SEQFEATURE_CLASS'],
 
162
                                                                                    ['TMP','TMPDIR'],
 
163
                                                                                    'VERBOSE',
 
164
                                                                                    'FAST',
 
165
                                                                                    'CHUNK_SIZE',
 
166
                                                                                   ],@_);
 
167
 
 
168
  $seqfeature_class ||= $self->default_seqfeature_class;
 
169
  eval "require $seqfeature_class" unless $seqfeature_class->can('new');
 
170
  $self->throw($@) if $@;
 
171
 
 
172
  my $normalized = $seqfeature_class->can('subfeatures_are_normalized')
 
173
    && $seqfeature_class->subfeatures_are_normalized;
 
174
 
 
175
  my $in_table = $seqfeature_class->can('subfeatures_are_stored_in_a_table')
 
176
    && $seqfeature_class->subfeatures_are_stored_in_a_table;
 
177
 
 
178
  if ($fast) {
 
179
    my $canfast = $normalized && $in_table;
 
180
    warn <<END unless $canfast;
 
181
Only features that support the Bio::DB::SeqFeature::NormalizedTableFeature interface
 
182
can be loaded using the -fast method. Reverting to slower feature-by-feature method.
 
183
END
 
184
    $fast &&= $canfast;
 
185
  }
 
186
 
 
187
  # try to bring in highres time() function
 
188
  eval "require Time::HiRes";
 
189
 
 
190
  $tmpdir ||= File::Spec->tmpdir();
 
191
 
 
192
  my $tmp_store = Bio::DB::SeqFeature::Store->new(-adaptor  => 'berkeleydb',
 
193
                                                  -temporary=> 1,
 
194
                                                  -dsn      => $tmpdir,
 
195
                                                  -cache    => 1,
 
196
                                                  -write    => 1)
 
197
    unless $normalized;
 
198
 
 
199
  return bless {
 
200
                store            => $store,
 
201
                tmp_store        => $tmp_store,
 
202
                seqfeature_class => $seqfeature_class,
 
203
                fast             => $fast,
 
204
                seq_chunk_size   => $seq_chunk_size || DEFAULT_SEQ_CHUNK_SIZE,
 
205
                verbose          => $verbose,
 
206
                load_data        => {},
 
207
                subfeatures_normalized => $normalized,
 
208
                subfeatures_in_table   => $in_table,
 
209
               },ref($self) || $self;
 
210
}
 
211
 
 
212
=head2 load
 
213
 
 
214
 Title   : load
 
215
 Usage   : $count = $loader->load(@ARGV)
 
216
 Function: load the indicated files or filehandles
 
217
 Returns : number of feature lines loaded
 
218
 Args    : list of files or filehandles
 
219
 Status  : public
 
220
 
 
221
Once the loader is created, invoke its load() method with a list of
 
222
GFF3 or FASTA file paths or previously-opened filehandles in order to
 
223
load them into the database. Compressed files ending with .gz, .Z and
 
224
.bz2 are automatically recognized and uncompressed on the fly. Paths
 
225
beginning with http: or ftp: are treated as URLs and opened using the
 
226
LWP GET program (which must be on your path).
 
227
 
 
228
FASTA files are recognized by their initial "E<gt>" character. Do not feed
 
229
the loader a file that is neither GFF3 nor FASTA; I don't know what
 
230
will happen, but it will probably not be what you expect.
 
231
 
 
232
=cut
 
233
 
 
234
sub load {
 
235
  my $self       = shift;
 
236
  my $start      = $self->time();
 
237
  my $count = 0;
 
238
 
 
239
  for my $file_or_fh (@_) {
 
240
    $self->msg("loading $file_or_fh...\n");
 
241
    my $fh = $self->open_fh($file_or_fh) or $self->throw("Couldn't open $file_or_fh: $!");
 
242
    $count += $self->load_fh($fh);
 
243
    $self->msg(sprintf "load time: %5.2fs\n",$self->time()-$start);
 
244
  }
 
245
  $count;
 
246
}
 
247
 
 
248
=head2 accessors
 
249
 
 
250
The following read-only accessors return values passed or created during new():
 
251
 
 
252
 store()          the long-term Bio::DB::SeqFeature::Store object
 
253
 
 
254
 tmp_store()      the temporary Bio::DB::SeqFeature::Store object used
 
255
                    during loading
 
256
 
 
257
 sfclass()        the Bio::SeqFeatureI class
 
258
 
 
259
 fast()           whether fast loading is active
 
260
 
 
261
 seq_chunk_size() the sequence chunk size
 
262
 
 
263
 verbose()        verbose progress messages
 
264
 
 
265
=cut
 
266
 
 
267
sub store          { shift->{store}            }
 
268
sub tmp_store      { shift->{tmp_store}        }
 
269
sub sfclass        { shift->{seqfeature_class} }
 
270
sub fast           { shift->{fast}             }
 
271
sub seq_chunk_size { shift->{seq_chunk_size}             }
 
272
sub verbose        { shift->{verbose}          }
 
273
 
 
274
=head2 Internal Methods
 
275
 
 
276
The following methods are used internally and may be overidden by
 
277
subclasses.
 
278
 
 
279
=over 4
 
280
 
 
281
=item default_seqfeature_class
 
282
 
 
283
  $class = $loader->default_seqfeature_class
 
284
 
 
285
Return the default SeqFeatureI class (Bio::DB::SeqFeature).
 
286
 
 
287
=cut
 
288
 
 
289
sub default_seqfeature_class {
 
290
  my $self = shift;
 
291
  return 'Bio::DB::SeqFeature';
 
292
}
 
293
 
 
294
=item subfeatures_normalized
 
295
 
 
296
  $flag = $loader->subfeatures_normalized([$new_flag])
 
297
 
 
298
Get or set a flag that indicates that the subfeatures are
 
299
normalized. This is deduced from the SeqFeature class information.
 
300
 
 
301
=cut
 
302
 
 
303
sub subfeatures_normalized {
 
304
  my $self = shift;
 
305
  my $d    = $self->{subfeatures_normalized};
 
306
  $self->{subfeatures_normalized} = shift if @_;
 
307
  $d;
 
308
}
 
309
 
 
310
=item subfeatures_in_table
 
311
 
 
312
  $flag = $loader->subfeatures_in_table([$new_flag])
 
313
 
 
314
Get or set a flag that indicates that feature/subfeature relationships
 
315
are stored in a table. This is deduced from the SeqFeature class and
 
316
Store information.
 
317
 
 
318
=cut
 
319
 
 
320
sub subfeatures_in_table {
 
321
  my $self = shift;
 
322
  my $d    = $self->{subfeatures_in_table};
 
323
  $self->{subfeatures_in_table} = shift if @_;
 
324
  $d;
 
325
}
 
326
 
 
327
=item load_fh
 
328
 
 
329
  $count = $loader->load_fh($filehandle)
 
330
 
 
331
Load the GFF3 data at the other end of the filehandle and return true
 
332
if successful. Internally, load_fh() invokes:
 
333
 
 
334
  start_load();
 
335
  do_load($filehandle);
 
336
  finish_load();
 
337
 
 
338
=cut
 
339
 
 
340
sub load_fh {
 
341
  my $self = shift;
 
342
  my $fh   = shift;
 
343
  $self->start_load();
 
344
  my $count = $self->do_load($fh);
 
345
  $self->finish_load();
 
346
  $count;
 
347
}
 
348
 
 
349
 
 
350
=item start_load, finish_load
 
351
 
 
352
These methods are called at the start and end of a filehandle load.
 
353
 
 
354
=cut
 
355
 
 
356
sub start_load {
 
357
  my $self = shift;
 
358
  $self->{load_data}{Parent2Child}     = {};
 
359
  $self->{load_data}{Local2GlobalID}   = {};
 
360
  $self->{load_data}{TemporaryID}      = "GFFLoad0000000";
 
361
  $self->{load_data}{IndexSubfeatures} = 1;
 
362
  $self->{load_data}{CurrentFeature}   = undef;
 
363
  $self->{load_data}{CurrentID}        = undef;
 
364
  $self->store->start_bulk_update() if $self->fast;
 
365
}
 
366
 
 
367
sub finish_load {
 
368
  my $self  = shift;
 
369
 
 
370
  $self->msg("Building object tree...");
 
371
  my $start = $self->time();
 
372
  $self->build_object_tree;
 
373
  $self->msg(sprintf "%5.2fs\n",$self->time()-$start);
 
374
 
 
375
  if ($self->fast) {
 
376
    $self->msg("Loading bulk data into database...");
 
377
    $start = $self->time();
 
378
    $self->store->finish_bulk_update;
 
379
    $self->msg(sprintf "%5.2fs\n",$self->time()-$start);
 
380
  }
 
381
  eval {$self->store->commit};
 
382
  delete $self->{load_data};
 
383
}
 
384
 
 
385
=item do_load
 
386
 
 
387
  $count = $loader->do_load($fh)
 
388
 
 
389
This is called by load_fh() to load the GFF3 file's filehandle and
 
390
return the number of lines loaded.
 
391
 
 
392
=cut
 
393
 
 
394
sub do_load {
 
395
  my $self = shift;
 
396
  my $fh   = shift;
 
397
 
 
398
  my $start = $self->time();
 
399
  my $count = 0;
 
400
  my $mode  = 'gff';  # or 'fasta'
 
401
 
 
402
  while (<$fh>) {
 
403
    chomp;
 
404
 
 
405
    next unless /^\S/;     # blank line
 
406
    $mode = 'gff' if /\t/;  # if it has a tab in it, switch to gff mode
 
407
 
 
408
    if (/^\#\s?\#\s*(.+)/) {  ## meta instruction
 
409
      $mode = 'gff';
 
410
      $self->handle_meta($1);
 
411
 
 
412
    } elsif (/^\#/) {
 
413
      $mode = 'gff';  # just to be safe
 
414
      next;  # comment
 
415
    }
 
416
 
 
417
    elsif (/^>\s*(\S+)/) { # FASTA lines are coming
 
418
      $mode = 'fasta';
 
419
      $self->start_or_finish_sequence($1);
 
420
    }
 
421
 
 
422
    elsif ($mode eq 'fasta') {
 
423
      $self->load_sequence($_);
 
424
    }
 
425
 
 
426
    elsif ($mode eq 'gff') {
 
427
      $self->handle_feature($_);
 
428
      if (++$count % 1000 == 0) {
 
429
        my $now = $self->time();
 
430
        my $nl = -t STDOUT && !$ENV{EMACS} ? "\r" : "\n";
 
431
        $self->msg(sprintf("%d features loaded in %5.2fs...$nl",$count,$now - $start));
 
432
        $start = $now;
 
433
      }
 
434
    }
 
435
 
 
436
    else {
 
437
      $self->throw("I don't know what to do with this line:\n$_");
 
438
    }
 
439
  }
 
440
  $self->store_current_feature();      # during fast loading, we will have a feature left at the very end
 
441
  $self->start_or_finish_sequence();   # finish any half-loaded sequences
 
442
  $self->msg(' 'x80,"\n"); #clear screen
 
443
  $count;
 
444
}
 
445
 
 
446
=item handle_meta
 
447
 
 
448
  $loader->handle_meta($meta_directive)
 
449
 
 
450
This method is called to handle meta-directives such as
 
451
##sequence-region. The method will receive the directive with the
 
452
initial ## stripped off.
 
453
 
 
454
=cut
 
455
 
 
456
sub handle_meta {
 
457
  my $self = shift;
 
458
  my $instruction = shift;
 
459
 
 
460
  if ($instruction =~ /sequence-region\s+(.+)\s+(-?\d+)\s+(-?\d+)/i) {
 
461
    my $feature = $self->sfclass->new(-name        => $1,
 
462
                                      -seq_id      => $1,
 
463
                                      -start       => $2,
 
464
                                      -end         => $3,
 
465
                                      -primary_tag => 'region');
 
466
    $self->store->store($feature);
 
467
    return;
 
468
  }
 
469
 
 
470
  if ($instruction =~/index-subfeatures\s+(\S+)/i) {
 
471
    $self->{load_data}{IndexSubfeatures} = $1;
 
472
    $self->store->index_subfeatures($1);
 
473
    return;
 
474
  }
 
475
}
 
476
 
 
477
=item handle_feature
 
478
 
 
479
  $loader->handle_feature($gff3_line)
 
480
 
 
481
This method is called to process a single GFF3 line. It manipulates
 
482
information stored a data structure called $self-E<gt>{load_data}.
 
483
 
 
484
=cut
 
485
 
 
486
sub handle_feature {
 
487
  my $self     = shift;
 
488
  my $gff_line = shift;
 
489
  my $ld       = $self->{load_data};
 
490
 
 
491
  my @columns = map {$_ eq '.' ? undef : $_ } split /\t/,$gff_line;
 
492
  return unless @columns >= 8;
 
493
  my ($refname,$source,$method,$start,$end, $score,$strand,$phase,$attributes)      = @columns;
 
494
  $strand = $Strandedness{$strand||0};
 
495
 
 
496
  my ($reserved,$unreserved) = $self->parse_attributes($attributes);
 
497
 
 
498
  my $name        = ($reserved->{Name}   && $reserved->{Name}[0]);
 
499
 
 
500
  my $has_loadid  = defined $reserved->{ID}[0];
 
501
 
 
502
  my $feature_id  = $reserved->{ID}[0] || $ld->{TemporaryID}++;
 
503
  my @parent_ids  = @{$reserved->{Parent}} if $reserved->{Parent};
 
504
 
 
505
  my $index_it = $ld->{IndexSubfeatures};
 
506
  if (exists $reserved->{Index} || exists $reserved->{index}) {
 
507
    $index_it = $reserved->{Index}[0] || $reserved->{index}[0];
 
508
  }
 
509
 
 
510
  # Everything in the unreserved hash becomes an attribute, so we copy
 
511
  # some attributes over
 
512
  $unreserved->{Note}   = $reserved->{Note}   if exists $reserved->{Note};
 
513
  $unreserved->{Alias}  = $reserved->{Alias}  if exists $reserved->{Alias};
 
514
  $unreserved->{Target} = $reserved->{Target} if exists $reserved->{Target};
 
515
  $unreserved->{Gap}    = $reserved->{Gap}    if exists $reserved->{Gap};
 
516
  $unreserved->{load_id}= $reserved->{ID}     if exists $reserved->{ID};
 
517
 
 
518
  # TEMPORARY HACKS TO SIMPLIFY DEBUGGING
 
519
  push @{$unreserved->{Alias}},$feature_id  if $has_loadid;
 
520
  $unreserved->{parent_id} = \@parent_ids   if @parent_ids;
 
521
 
 
522
  # POSSIBLY A PERMANENT HACK -- TARGETS BECOME ALIASES
 
523
  # THIS IS TO ALLOW FOR TARGET-BASED LOOKUPS
 
524
  if (exists $reserved->{Target}) {
 
525
    my %aliases = map {$_=>1} @{$unreserved->{Alias}};
 
526
    for my $t (@{$reserved->{Target}}) {
 
527
      (my $tc = $t) =~ s/\s+.*$//;  # get rid of coordinates
 
528
      $name ||= $tc;
 
529
      push @{$unreserved->{Alias}},$tc unless $name eq $tc || $aliases{$tc};
 
530
    }
 
531
  }
 
532
 
 
533
  my @args = (-display_name => $name,
 
534
              -seq_id       => $refname,
 
535
              -start        => $start,
 
536
              -end          => $end,
 
537
              -strand       => $strand || 0,
 
538
              -score        => $score,
 
539
              -phase        => $phase,
 
540
              -primary_tag  => $method || 'feature',
 
541
              -source       => $source,
 
542
              -tag          => $unreserved,
 
543
              -attributes   => $unreserved,
 
544
             );
 
545
 
 
546
  # Here's where we handle feature lines that have the same ID (multiple locations, not
 
547
  # parent/child relationships)
 
548
 
 
549
  my $old_feat;
 
550
 
 
551
  # Current feature is the same as the previous feature, which hasn't yet been loaded
 
552
  if (defined $ld->{CurrentID} && $ld->{CurrentID} eq $feature_id) {
 
553
    $old_feat = $ld->{CurrentFeature};
 
554
  }
 
555
 
 
556
  # Current feature is the same as a feature that was loaded earlier
 
557
  elsif (my $id = $self->{load_data}{Local2GlobalID}{$feature_id}) {
 
558
    $old_feat = $self->fetch($feature_id)
 
559
      or $self->warn(<<END);
 
560
ID=$feature_id has been used more than once, but it cannot be found in the database.
 
561
This can happen if you have specified fast loading, but features sharing the same ID
 
562
are not contiguous in the GFF file. This will be loaded as a separate feature.
 
563
Line $.: "$_"
 
564
END
 
565
  }
 
566
 
 
567
  # contiguous feature, so add a segment
 
568
  if (defined $old_feat) {
 
569
    $self->add_segment($old_feat,$self->sfclass->new(@args));
 
570
    return;
 
571
  }
 
572
 
 
573
  # we get here if this is a new feature
 
574
  # first of all, store the current feature if it is there
 
575
  $self->store_current_feature() if defined $ld->{CurrentID};
 
576
 
 
577
  # now create the new feature
 
578
  # (index top-level features only if policy asks us to)
 
579
  my $feature = $self->sfclass->new(@args);
 
580
  $feature->object_store($self->store) if $feature->can('object_store');  # for lazy table features
 
581
  $ld->{CurrentFeature} = $feature;
 
582
  $ld->{CurrentID}      = $feature_id;
 
583
 
 
584
  my $top_level = !@parent_ids;
 
585
  my $has_id    = defined $reserved->{ID}[0];
 
586
  $index_it   ||= $top_level;
 
587
 
 
588
  $ld->{IndexIt}{$feature_id}++    if $index_it;
 
589
  $ld->{TopLevel}{$feature_id}++   if !$self->{fast} && $top_level;  # need to track top level features
 
590
 
 
591
  # remember parentage
 
592
  for my $parent (@parent_ids) {
 
593
    push @{$ld->{Parent2Child}{$parent}},$feature_id;
 
594
  }
 
595
 
 
596
}
 
597
 
 
598
=item store_current_feature
 
599
 
 
600
  $loader->store_current_feature()
 
601
 
 
602
This method is called to store the currently active feature in the
 
603
database. It uses a data structure stored in $self-E<gt>{load_data}.
 
604
 
 
605
=cut
 
606
 
 
607
 
 
608
sub store_current_feature {
 
609
  my $self    = shift;
 
610
 
 
611
  my $ld   = $self->{load_data};
 
612
  defined $ld->{CurrentFeature} or return;
 
613
  my $f    = $ld->{CurrentFeature};
 
614
 
 
615
  my $normalized = $self->subfeatures_normalized;
 
616
  my $indexed    = $ld->{IndexIt}{$ld->{CurrentID}};
 
617
 
 
618
  # logic is as follows:
 
619
  # 1. If the feature is an indexed feature, then we store it into the main database
 
620
  #    so that it can be searched. It doesn't matter whether it is a top-level feature
 
621
  #    or a subfeature.
 
622
  # 2. If the feature class is normalized, but not indexed, then we store it into the
 
623
  #    main database using the "no_index" method. This will make it accessible to
 
624
  #    queries on the top level parent, but it won't come up by itself in range or
 
625
  #    attribute searches.
 
626
  # 3. Otherwise, this is an unindexed subfeature; we store it in the temporary database
 
627
  #    until the object build step, at which point it gets integrated into its object tree
 
628
  #    and copied into the main database.
 
629
 
 
630
  if ($indexed) {
 
631
    $self->store->store($f);
 
632
  }
 
633
 
 
634
  elsif ($normalized) {
 
635
    $self->store->store_noindex($f)
 
636
  }
 
637
 
 
638
  else {
 
639
    $self->tmp_store->store_noindex($f)
 
640
  }
 
641
        
 
642
  my $id        = $f->primary_id;    # assigned by store()
 
643
  $ld->{Local2GlobalID}{$ld->{CurrentID}} = $id;
 
644
 
 
645
  undef $ld->{IndexIt}{$ld->{CurrentID}} if $normalized;  # no need to remember this
 
646
  undef $ld->{CurrentID};
 
647
  undef $ld->{CurrentFeature};
 
648
}
 
649
 
 
650
=item build_object_tree
 
651
 
 
652
 $loader->build_object_tree()
 
653
 
 
654
This method gathers together features and subfeatures and builds the graph that connects them.
 
655
 
 
656
=cut
 
657
 
 
658
###
 
659
# put objects together
 
660
#
 
661
sub build_object_tree {
 
662
  my $self = shift;
 
663
  $self->subfeatures_in_table ? $self->build_object_tree_in_tables : $self->build_object_tree_in_features;
 
664
}
 
665
 
 
666
=item build_object_tree_in_tables
 
667
 
 
668
 $loader->build_object_tree_in_tables()
 
669
 
 
670
This method gathers together features and subfeatures and builds the
 
671
graph that connects them, assuming that parent/child relationships
 
672
will be stored in a database table.
 
673
 
 
674
=cut
 
675
 
 
676
sub build_object_tree_in_tables {
 
677
  my $self = shift;
 
678
  my $store = $self->store;
 
679
  my $ld    = $self->{load_data};
 
680
 
 
681
  while (my ($load_id,$children) = each %{$ld->{Parent2Child}}) {
 
682
    my $parent_id = $ld->{Local2GlobalID}{$load_id} or $self->throw("$load_id doesn't have a primary id");
 
683
    my @children  = map {$ld->{Local2GlobalID}{$_}} @$children;
 
684
 
 
685
    # this updates the table that keeps track of parent/child relationships,
 
686
    # but does not update the parent object -- so (start,end) had better be right!!!
 
687
    $store->add_SeqFeature($parent_id,@children);
 
688
  }
 
689
 
 
690
}
 
691
 
 
692
=item build_object_tree_in_features
 
693
 
 
694
 $loader->build_object_tree_in_features()
 
695
 
 
696
This method gathers together features and subfeatures and builds the
 
697
graph that connects them, assuming that parent/child relationships are
 
698
stored in the seqfeature objects themselves.
 
699
 
 
700
=cut
 
701
 
 
702
sub build_object_tree_in_features {
 
703
  my $self  = shift;
 
704
  my $store      = $self->store;
 
705
  my $tmp        = $self->tmp_store;
 
706
  my $ld         = $self->{load_data};
 
707
  my $normalized = $self->subfeatures_normalized;
 
708
 
 
709
  while (my ($load_id) = each %{$ld->{TopLevel}}) {
 
710
    my $feature  = $self->fetch($load_id)
 
711
      or $self->throw("$load_id (id=$ld->{Local2GlobalID}{$load_id}) should have a database entry, but doesn't");
 
712
    $self->attach_children($store,$ld,$load_id,$feature);
 
713
    $feature->primary_id(undef) unless $ld->{IndexIt}{$load_id};  # Indexed objects are updated, not created anew
 
714
    $store->store($feature);
 
715
  }
 
716
 
 
717
}
 
718
 
 
719
=item attach_children
 
720
 
 
721
 $loader->attach_children($store,$load_data,$load_id,$feature)
 
722
 
 
723
This recursively adds children to features and their subfeatures. It
 
724
is called when subfeatures are directly contained within other
 
725
features, rather than stored in a relational table.
 
726
 
 
727
=cut
 
728
 
 
729
sub attach_children {
 
730
  my $self = shift;
 
731
  my ($store,$ld,$load_id,$feature)  = @_;
 
732
 
 
733
  my $children   = $ld->{Parent2Child}{$load_id} or return;
 
734
  for my $child_id (@$children) {
 
735
    my $child = $self->fetch($child_id)
 
736
      or $self->throw("$child_id should have a database entry, but doesn't");
 
737
    $self->attach_children($store,$ld,$child_id,$child);   # recursive call
 
738
    $feature->add_SeqFeature($child);
 
739
  }
 
740
}
 
741
 
 
742
=item fetch
 
743
 
 
744
 my $feature = $loader->fetch($load_id)
 
745
 
 
746
Given a load ID (from the ID= attribute) this method returns the
 
747
feature from the temporary database or the permanent one, depending on
 
748
where it is stored.
 
749
 
 
750
=cut
 
751
 
 
752
sub fetch {
 
753
  my $self    = shift;
 
754
  my $load_id = shift;
 
755
  my $ld      = $self->{load_data};
 
756
  my $id      = $ld->{Local2GlobalID}{$load_id};
 
757
 
 
758
  return
 
759
    $self->subfeatures_normalized || $ld->{IndexIt}{$load_id}
 
760
      ? $self->store->fetch($id)
 
761
      : $self->tmp_store->fetch($id);
 
762
}
 
763
 
 
764
=item add_segment
 
765
 
 
766
 $loader->add_segment($parent,$child)
 
767
 
 
768
This method is used to add a split location to the parent.
 
769
 
 
770
=cut
 
771
 
 
772
sub add_segment {
 
773
  my $self = shift;
 
774
  my ($parent,$child) = @_;
 
775
 
 
776
  if ($parent->can('add_segment')) { # probably a lazy table feature
 
777
    my $segment_count =  $parent->can('denormalized_segment_count') ? $parent->denormalized_segment_count
 
778
                       : $parent->can('denormalized_segments ')     ? $parent->denormalized_segments
 
779
                       : $parent->can('segments')                   ? $parent->segments
 
780
                       : 0;
 
781
    unless ($segment_count) {  # convert into a segmented object
 
782
      my $segment;
 
783
      if ($parent->can('clone')) {
 
784
        $segment = $parent->clone;
 
785
      } else {
 
786
        my %clone   = %$parent;
 
787
        $segment = bless \%clone,ref $parent;
 
788
      }
 
789
      delete $segment->{segments};
 
790
      eval {$segment->object_store(undef) };
 
791
      $segment->primary_id(undef);
 
792
 
 
793
      # this updates the object and expands its start and end positions without writing
 
794
      # the segments into the database as individual objects
 
795
      $parent->add_segment($segment);
 
796
    }
 
797
    $parent->add_segment($child);
 
798
    1; # for debugging
 
799
  }
 
800
 
 
801
  # a conventional Bio::SeqFeature::Generic object - create a split location
 
802
  else {
 
803
    my $current_location = $parent->location;
 
804
    if ($current_location->can('add_sub_Location')) {
 
805
      $current_location->add_sub_Location($child->location);
 
806
    } else {
 
807
      eval "require Bio::Location::Split" unless Bio::Location::Split->can('add_sub_Location');
 
808
      my $new_location = Bio::Location::Split->new();
 
809
      $new_location->add_sub_Location($current_location);
 
810
      $new_location->add_sub_Location($child->location);
 
811
      $parent->location($new_location);
 
812
    }
 
813
  }
 
814
}
 
815
 
 
816
=item parse_attributes
 
817
 
 
818
 ($reserved,$unreserved) = $loader->parse_attributes($attribute_line)
 
819
 
 
820
This method parses the information contained in the $attribute_line
 
821
into two hashrefs, one containing the values of reserved attribute
 
822
tags (e.g. ID) and the other containing the values of unreserved ones.
 
823
 
 
824
=cut
 
825
 
 
826
sub parse_attributes {
 
827
  my $self = shift;
 
828
  my $att  = shift;
 
829
  my @pairs =  map {my ($name,$value) = split /=/; [unescape($name) => unescape($value)] } split /;/,$att;
 
830
  my (%reserved,%unreserved);
 
831
  foreach (@pairs) {
 
832
    my $tag    = $_->[0];
 
833
    my @values = split /,/,$_->[1];
 
834
 
 
835
    if ($Special_attributes{$tag}) {  # reserved attribute
 
836
      push @{$reserved{$tag}},@values;
 
837
    } else {
 
838
      push @{$unreserved{$tag}},@values
 
839
    }
 
840
  }
 
841
  return (\%reserved,\%unreserved);
 
842
}
 
843
 
 
844
=item start_or_finish_sequence
 
845
 
 
846
  $loader->start_or_finish_sequence('Chr9')
 
847
 
 
848
This method is called at the beginning and end of a fasta section.
 
849
 
 
850
=cut
 
851
 
 
852
 
 
853
# this gets called at the beginning and end of a fasta section
 
854
sub start_or_finish_sequence {
 
855
  my $self  = shift;
 
856
  my $seqid = shift;
 
857
  if (my $sl    = $self->{fasta_load}) {
 
858
    if (defined $sl->{seqid}) {
 
859
      $self->store->insert_sequence($sl->{seqid},$sl->{sequence},$sl->{offset});
 
860
      delete $self->{fasta_load};
 
861
    }
 
862
  }
 
863
  if (defined $seqid) {
 
864
    $self->{fasta_load} = {seqid  => $seqid,
 
865
                           offset => 0,
 
866
                           sequence => ''};
 
867
  }
 
868
}
 
869
 
 
870
=item load_sequence
 
871
 
 
872
  $loader->load_sequence('gatttcccaaa')
 
873
 
 
874
This method is called to load some amount of sequence after
 
875
start_or_finish_sequence() is first called.
 
876
 
 
877
=cut
 
878
 
 
879
sub load_sequence {
 
880
  my $self = shift;
 
881
  my $seq  = shift;
 
882
  my $sl   = $self->{fasta_load} or return;
 
883
  my $cs   = $self->seq_chunk_size;
 
884
  $sl->{sequence} .= $seq;
 
885
  while (length $sl->{sequence} >= $cs) {
 
886
    my $chunk = substr($sl->{sequence},0,$cs);
 
887
    $self->store->insert_sequence($sl->{seqid},$chunk,$sl->{offset});
 
888
    $sl->{offset} += length $chunk;
 
889
    substr($sl->{sequence},0,$cs) = '';
 
890
  }
 
891
}
 
892
 
 
893
=item open_fh
 
894
 
 
895
 my $io_file = $loader->open_fh($filehandle_or_path)
 
896
 
 
897
This method opens up the indicated file or pipe, using some intelligence to recognized compressed files and URLs and doing 
 
898
the right thing.
 
899
 
 
900
=cut
 
901
 
 
902
 
 
903
sub open_fh {
 
904
  my $self  = shift;
 
905
  my $thing = shift;
 
906
 
 
907
  no strict 'refs';
 
908
 
 
909
  return $thing                                  if defined fileno($thing);
 
910
  return IO::File->new("gunzip -c $thing |")     if $thing =~ /\.gz$/;
 
911
  return IO::File->new("uncompress -c $thing |") if $thing =~ /\.Z$/;
 
912
  return IO::File->new("bunzip2 -c $thing |")    if $thing =~ /\.bz2$/;
 
913
  return IO::File->new("GET $thing |")           if $thing =~ /^(http|ftp):/;
 
914
  return IO::File->new($thing);
 
915
}
 
916
 
 
917
sub msg {
 
918
  my $self = shift;
 
919
  my @msg  = @_;
 
920
  return unless $self->verbose;
 
921
  print STDERR @msg;
 
922
}
 
923
 
 
924
=item time
 
925
 
 
926
 my $time = $loader->time
 
927
 
 
928
This method returns the current time in seconds, using Time::HiRes if available.
 
929
 
 
930
=cut
 
931
 
 
932
sub time {
 
933
  return Time::HiRes::time() if Time::HiRes->can('time');
 
934
  return time();
 
935
}
 
936
 
 
937
=item escape
 
938
 
 
939
 my $unescaped = GFF3Loader::unescape($escaped)
 
940
 
 
941
This is an internal utility.  It is the same as CGI::Util::unescape,
 
942
but don't change pluses into spaces and ignores unicode escapes.
 
943
 
 
944
=cut
 
945
 
 
946
sub unescape {
 
947
  my $todecode = shift;
 
948
  $todecode =~ s/%([0-9a-fA-F]{2})/chr hex($1)/ge;
 
949
  return $todecode;
 
950
}
 
951
 
 
952
1;
 
953
__END__
 
954
 
 
955
=back
 
956
 
 
957
=head1 BUGS
 
958
 
 
959
This is an early version, so there are certainly some bugs. Please
 
960
use the BioPerl bug tracking system to report bugs.
 
961
 
 
962
=head1 SEE ALSO
 
963
 
 
964
L<bioperl>,
 
965
L<Bio::DB::SeqFeature::Store>,
 
966
L<Bio::DB::SeqFeature::Segment>,
 
967
L<Bio::DB::SeqFeature::NormalizedFeature>,
 
968
L<Bio::DB::SeqFeature::GFF3Loader>,
 
969
L<Bio::DB::SeqFeature::Store::DBI::mysql>,
 
970
L<Bio::DB::SeqFeature::Store::bdb>
 
971
 
 
972
=head1 AUTHOR
 
973
 
 
974
Lincoln Stein E<lt>lstein@cshl.orgE<gt>.
 
975
 
 
976
Copyright (c) 2006 Cold Spring Harbor Laboratory.
 
977
 
 
978
This library is free software; you can redistribute it and/or modify
 
979
it under the same terms as Perl itself.
 
980
 
 
981
=cut
 
982
 
 
983