Summary

Bio::DB::SeqFeature::Store::GFF3Loader -- GFF3 file loader for Bio::DB::SeqFeature::Store

Package variables

Included modules

Synopsis

  use Bio::DB::SeqFeature::Store;

  # Open the sequence database
  my $db      = Bio::DB::SeqFeature::Store->new( -adaptor => 'DBI::mysql',
                                                 -dsn     => 'dbi:mysql:test',
                                                 -write   => 1 );

  my $loader = Bio::DB::SeqFeature::Store::GFF3Loader->new(-store    => $db,
							   -verbose  => 1,
							   -fast     => 1);

  $loader->load('./my_genome.gff3');

Description

The Bio::DB::SeqFeature::Store::GFF3Loader object parsers GFF3-format
sequence annotation files and loads Bio::DB::SeqFeature::Store
databases. For certain combinations of SeqFeature classes and
SeqFeature::Store databases it features a "fast load" mode which will
greatly accelerate the loading of GFF3 databases by a factor of 5-10.
The GFF3 file format has been extended very slightly to accomodate
Bio::DB::SeqFeature::Store. First, the loader recognizes is a new
directive:

  # #index-subfeatures [0|1]

Note that you can place a space between the two #'s in order to
prevent GFF3 validators from complaining.
If this is true, then subfeatures are indexed (the default) so that
they can be retrieved with a query. See Bio::DB::SeqFeature::Store
for an explanation of this. If false, then subfeatures can only be
accessed through their parent feature. The default is to index all
subfeatures.
Second, the loader recognizes a new attribute tag called index, which
if present, controls indexing of the current feature. Example:

 ctg123	. TF_binding_site 1000 1012 . + . ID=tfbs00001;index=1

You can use this to turn indexing on and off, overriding the default
for a particular feature.

Methods

Methods description

 Title   : new
 Usage   : $loader = Bio::DB::SeqFeature::Store::GFF3Loader->new(@options)
 Function: create a new parser
 Returns : a Bio::DB::SeqFeature::Store::GFF3Loader gff3 parser and loader
 Args    : several - see below
 Status  : public

 Name               Value
 ----               -----

 -store             A writeable Bio::DB::SeqFeature::Store database handle.

 -seqfeature_class  The name of the type of Bio::SeqFeatureI object to create
                      and store in the database (Bio::DB::SeqFeature by default)

 -sf_class          A shorter alias for -seqfeature_class

 -verbose           Send progress information to standard error.

 -fast              If true, activate fast loading (see below)

 -chunk_size        Set the storage chunk size for nucleotide/protein sequences
                       (default 2000 bytes)

 -tmp               Indicate a temporary directory to use when loading non-normalized
                       features.

 Title   : load
 Usage   : $count = $loader->load(@ARGV)
 Function: load the indicated files or filehandles
 Returns : number of feature lines loaded
 Args    : list of files or filehandles
 Status  : public

Methods code

sub new {

  my $self = shift;
  my ($store,$seqfeature_class,$tmpdir,$verbose,$fast,$seq_chunk_size) = rearrange(['STORE',
										    ['SF_CLASS','SEQFEATURE_CLASS'],
										    ['TMP','TMPDIR'],
										    'VERBOSE',
										    'FAST',
										    'CHUNK_SIZE',
										   ],@_);

  $seqfeature_class ||= $self->default_seqfeature_class;
  eval "require $seqfeature_class" unless $seqfeature_class->can('new');
  $self->throw($@) if $@;

  my $normalized = $seqfeature_class->can('subfeatures_are_normalized')
    && $seqfeature_class->subfeatures_are_normalized;

  my $in_table = $seqfeature_class->can('subfeatures_are_stored_in_a_table')
    && $seqfeature_class->subfeatures_are_stored_in_a_table;

  if ($fast) {
    my $canfast = $normalized && $in_table;
    warn <<END unless $canfast;
Only features that support the Bio::DB::SeqFeature::NormalizedTableFeature interface
can be loaded using the -fast method. Reverting to slower feature-by-feature method.
END
    $fast &&= $canfast;
  }

  # try to bring in highres time() function
  eval "require Time::HiRes";

  $tmpdir ||= File::Spec->tmpdir();

  my $tmp_store = Bio::DB::SeqFeature::Store->new(-adaptor  => 'berkeleydb',
						  -temporary=> 1,
						  -dsn      => $tmpdir,
						  -cache    => 1,
						  -write    => 1)
    unless $normalized;

  return bless {
		store            => $store,
		tmp_store        => $tmp_store,
		seqfeature_class => $seqfeature_class,
		fast             => $fast,
		seq_chunk_size   => $seq_chunk_size || DEFAULT_SEQ_CHUNK_SIZE,
		verbose          => $verbose,
		load_data        => {},
		subfeatures_normalized => $normalized,
		subfeatures_in_table   => $in_table,
	       },ref($self) || $self;

}

sub load {

  my $self       = shift;
  my $start      = $self->time();
  my $count = 0;

  for my $file_or_fh (@_) {
    $self->msg("loading $file_or_fh...\n");
    my $fh = $self->open_fh($file_or_fh) or $self->throw("Couldn't open $file_or_fh: $!");
    $count += $self->load_fh($fh);
    $self->msg(sprintf "load time: %5.2fs\n",$self->time()-$start);
  }
  $count;

}

sub store {

 shift->{store}

}

sub tmp_store {

 shift->{tmp_store}

}

sub sfclass {

 shift->{seqfeature_class}

}

sub fast {

 shift->{fast}

}

sub seq_chunk_size {

 shift->{seq_chunk_size}

}

sub verbose {

 shift->{verbose}

}

sub default_seqfeature_class {

  my $self = shift;
  return 'Bio::DB::SeqFeature';

}

sub subfeatures_normalized {

  my $self = shift;
  my $d    = $self->{subfeatures_normalized};
  $self->{subfeatures_normalized} = shift if @_;
  $d;

}

sub subfeatures_in_table {

  my $self = shift;
  my $d    = $self->{subfeatures_in_table};
  $self->{subfeatures_in_table} = shift if @_;
  $d;

}

sub load_fh {

  my $self = shift;
  my $fh   = shift;
  $self->start_load();
  my $count = $self->do_load($fh);
  $self->finish_load();
  $count;

}

sub start_load {

  my $self = shift;
  $self->{load_data}{Parent2Child}     = {};
  $self->{load_data}{Local2GlobalID}   = {};
  $self->{load_data}{TemporaryID}      = "GFFLoad0000000";
  $self->{load_data}{IndexSubfeatures} = 1;
  $self->{load_data}{CurrentFeature}   = undef;
  $self->{load_data}{CurrentID}        = undef;
  $self->store->start_bulk_update() if $self->fast;

}

sub finish_load {

  my $self  = shift;

  $self->msg("Building object tree...");
  my $start = $self->time();
  $self->build_object_tree;
  $self->msg(sprintf "%5.2fs\n",$self->time()-$start);

  if ($self->fast) {
    $self->msg("Loading bulk data into database...");
    $start = $self->time();
    $self->store->finish_bulk_update;
    $self->msg(sprintf "%5.2fs\n",$self->time()-$start);
  }
  eval {$self->store->commit};
  delete $self->{load_data};

}

sub do_load {

  my $self = shift;
  my $fh   = shift;

  my $start = $self->time();
  my $count = 0;
  my $mode  = 'gff';  # or 'fasta'

  while (<$fh>) {
    chomp;

    next unless /^\S/;     # blank line
    $mode = 'gff' if /\t/;  # if it has a tab in it, switch to gff mode

    if (/^\#\s?\#\s*(.+)/) {  ## meta instruction
      $mode = 'gff';
      $self->handle_meta($1);

    } elsif (/^\#/) {
      $mode = 'gff';  # just to be safe
      next;  # comment
    }

    elsif (/^>\s*(\S+)/) { # FASTA lines are coming
      $mode = 'fasta';
      $self->start_or_finish_sequence($1);
    }

    elsif ($mode eq 'fasta') {
      $self->load_sequence($_);
    }

    elsif ($mode eq 'gff') {
      $self->handle_feature($_);
      if (++$count % 1000 == 0) {
	my $now = $self->time();
	my $nl = -t STDOUT && !$ENV{EMACS} ? "\r" : "\n";
	$self->msg(sprintf("%d features loaded in %5.2fs...$nl",$count,$now - $start));
	$start = $now;
      }
    }

    else {
      $self->throw("I don't know what to do with this line:\n$_");
    }
  }
  $self->store_current_feature();      # during fast loading, we will have a feature left at the very end
  $self->start_or_finish_sequence();   # finish any half-loaded sequences
  $self->msg(' 'x80,"\n"); #clear screen
  $count;

}

sub handle_meta {

  my $self = shift;
  my $instruction = shift;

  if ($instruction =~ /sequence-region\s+(.+)\s+(-?\d+)\s+(-?\d+)/i) {
    my $feature = $self->sfclass->new(-name        => $1,
				      -seq_id      => $1,
				      -start       => $2,
				      -end         => $3,
				      -primary_tag => 'region');
    $self->store->store($feature);
    return;
  }

  if ($instruction =~/index-subfeatures\s+(\S+)/i) {
    $self->{load_data}{IndexSubfeatures} = $1;
    $self->store->index_subfeatures($1);
    return;
  }

}

sub handle_feature {

  my $self     = shift;
  my $gff_line = shift;
  my $ld       = $self->{load_data};

  my @columns = map {$_ eq '.' ? undef : $_ } split /\t/,$gff_line;
  return unless @columns >= 8;
  my ($refname,$source,$method,$start,$end, $score,$strand,$phase,$attributes)      = @columns;
  $strand = $Strandedness{$strand||0};

  my ($reserved,$unreserved) = $self->parse_attributes($attributes);

  my $name        = ($reserved->{Name}   && $reserved->{Name}[0]);

  my $has_loadid  = defined $reserved->{ID}[0];

  my $feature_id  = $reserved->{ID}[0] || $ld->{TemporaryID}++;
  my @parent_ids  = @{$reserved->{Parent}} if $reserved->{Parent};

  my $index_it = $ld->{IndexSubfeatures};
  if (exists $reserved->{Index} || exists $reserved->{index}) {
    $index_it = $reserved->{Index}[0] || $reserved->{index}[0];
  }

  # Everything in the unreserved hash becomes an attribute, so we copy
  # some attributes over
  $unreserved->{Note}   = $reserved->{Note}   if exists $reserved->{Note};
  $unreserved->{Alias}  = $reserved->{Alias}  if exists $reserved->{Alias};
  $unreserved->{Target} = $reserved->{Target} if exists $reserved->{Target};
  $unreserved->{Gap}    = $reserved->{Gap}    if exists $reserved->{Gap};
  $unreserved->{load_id}= $reserved->{ID}     if exists $reserved->{ID};

  # TEMPORARY HACKS TO SIMPLIFY DEBUGGING
  push @{$unreserved->{Alias}},$feature_id  if $has_loadid;
  $unreserved->{parent_id} =\@ parent_ids   if @parent_ids;

  # POSSIBLY A PERMANENT HACK -- TARGETS BECOME ALIASES
  # THIS IS TO ALLOW FOR TARGET-BASED LOOKUPS
  if (exists $reserved->{Target}) {
    my %aliases = map {$_=>1} @{$unreserved->{Alias}};
    for my $t (@{$reserved->{Target}}) {
      (my $tc = $t) =~ s/\s+.*$//;  # get rid of coordinates
      $name ||= $tc;
      push @{$unreserved->{Alias}},$tc unless $name eq $tc || $aliases{$tc};
    }
  }

  my @args = (-display_name => $name,
	      -seq_id       => $refname,
	      -start        => $start,
	      -end          => $end,
	      -strand       => $strand || 0,
	      -score        => $score,
	      -phase        => $phase,
	      -primary_tag  => $method || 'feature',
	      -source       => $source,
	      -tag          => $unreserved,
	      -attributes   => $unreserved,
	     );

  # Here's where we handle feature lines that have the same ID (multiple locations, not
  # parent/child relationships)

  my $old_feat;

  # Current feature is the same as the previous feature, which hasn't yet been loaded
  if (defined $ld->{CurrentID} && $ld->{CurrentID} eq $feature_id) {
    $old_feat = $ld->{CurrentFeature};
  }

  # Current feature is the same as a feature that was loaded earlier
  elsif (my $id = $self->{load_data}{Local2GlobalID}{$feature_id}) {
    $old_feat = $self->fetch($feature_id)
      or $self->warn(<<END);
ID=$feature_id has been used more than once, but it cannot be found in the database.
This can happen if you have specified fast loading, but features sharing the same ID
are not contiguous in the GFF file. This will be loaded as a separate feature.
Line $.: "$_"
END
  }

  # contiguous feature, so add a segment
  if (defined $old_feat) {
    $self->add_segment($old_feat,$self->sfclass->new(@args));
    return;
  }

  # we get here if this is a new feature
  # first of all, store the current feature if it is there
  $self->store_current_feature() if defined $ld->{CurrentID};

  # now create the new feature
  # (index top-level features only if policy asks us to)
  my $feature = $self->sfclass->new(@args);
  $feature->object_store($self->store) if $feature->can('object_store');  # for lazy table features
  $ld->{CurrentFeature} = $feature;
  $ld->{CurrentID}      = $feature_id;

  my $top_level = !@parent_ids;
  my $has_id    = defined $reserved->{ID}[0];
  $index_it   ||= $top_level;

  $ld->{IndexIt}{$feature_id}++    if $index_it;
  $ld->{TopLevel}{$feature_id}++   if !$self->{fast} && $top_level;  # need to track top level features

  # remember parentage
  for my $parent (@parent_ids) {
    push @{$ld->{Parent2Child}{$parent}},$feature_id;
  }

}

sub store_current_feature {

  my $self    = shift;

  my $ld   = $self->{load_data};
  defined $ld->{CurrentFeature} or return;
  my $f    = $ld->{CurrentFeature};

  my $normalized = $self->subfeatures_normalized;
  my $indexed    = $ld->{IndexIt}{$ld->{CurrentID}};

  # logic is as follows:
  # 1. If the feature is an indexed feature, then we store it into the main database
  #    so that it can be searched. It doesn't matter whether it is a top-level feature
  #    or a subfeature.
  # 2. If the feature class is normalized, but not indexed, then we store it into the
  #    main database using the "no_index" method. This will make it accessible to
  #    queries on the top level parent, but it won't come up by itself in range or
  #    attribute searches.
  # 3. Otherwise, this is an unindexed subfeature; we store it in the temporary database
  #    until the object build step, at which point it gets integrated into its object tree
  #    and copied into the main database.

  if ($indexed) {
    $self->store->store($f);
  }

  elsif ($normalized) {
    $self->store->store_noindex($f)
  }

  else {
    $self->tmp_store->store_noindex($f)
  }
	
  my $id        = $f->primary_id;    # assigned by store()
  $ld->{Local2GlobalID}{$ld->{CurrentID}} = $id;

  undef $ld->{IndexIt}{$ld->{CurrentID}} if $normalized;  # no need to remember this
  undef $ld->{CurrentID};
  undef $ld->{CurrentFeature};

}

sub build_object_tree {

  my $self = shift;
  $self->subfeatures_in_table ? $self->build_object_tree_in_tables : $self->build_object_tree_in_features;

}

sub build_object_tree_in_tables {

  my $self = shift;
  my $store = $self->store;
  my $ld    = $self->{load_data};

  while (my ($load_id,$children) = each %{$ld->{Parent2Child}}) {
    my $parent_id = $ld->{Local2GlobalID}{$load_id} or $self->throw("$load_id doesn't have a primary id");
    my @children  = map {$ld->{Local2GlobalID}{$_}} @$children;

    # this updates the table that keeps track of parent/child relationships,
    # but does not update the parent object -- so (start,end) had better be right!!!
    $store->add_SeqFeature($parent_id,@children);
  }

}

sub build_object_tree_in_features {

  my $self  = shift;
  my $store      = $self->store;
  my $tmp        = $self->tmp_store;
  my $ld         = $self->{load_data};
  my $normalized = $self->subfeatures_normalized;

  while (my ($load_id) = each %{$ld->{TopLevel}}) {
    my $feature  = $self->fetch($load_id)
      or $self->throw("$load_id (id=$ld->{Local2GlobalID}{$load_id}) should have a database entry, but doesn't");
    $self->attach_children($store,$ld,$load_id,$feature);
    $feature->primary_id(undef) unless $ld->{IndexIt}{$load_id};  # Indexed objects are updated, not created anew
    $store->store($feature);
  }

}

sub attach_children {

  my $self = shift;
  my ($store,$ld,$load_id,$feature)  = @_;

  my $children   = $ld->{Parent2Child}{$load_id} or return;
  for my $child_id (@$children) {
    my $child = $self->fetch($child_id)
      or $self->throw("$child_id should have a database entry, but doesn't");
    $self->attach_children($store,$ld,$child_id,$child);   # recursive call
    $feature->add_SeqFeature($child);
  }

}

sub add_segment {

  my $self = shift;
  my ($parent,$child) = @_;

  if ($parent->can('add_segment')) { # probably a lazy table feature
    my $segment_count =  $parent->can('denormalized_segment_count') ? $parent->denormalized_segment_count
                       : $parent->can('denormalized_segments ')     ? $parent->denormalized_segments
		       : $parent->can('segments')                   ? $parent->segments
		       : 0;
    unless ($segment_count) {  # convert into a segmented object
      my $segment;
      if ($parent->can('clone')) {
	$segment = $parent->clone;
      } else {
	my %clone   = %$parent;
	$segment = bless\% clone,ref $parent;
      }
      delete $segment->{segments};
      eval {$segment->object_store(undef) };
      $segment->primary_id(undef);

      # this updates the object and expands its start and end positions without writing
      # the segments into the database as individual objects
      $parent->add_segment($segment);
    }
    $parent->add_segment($child);
    1; # for debugging
  }

  # a conventional Bio::SeqFeature::Generic object - create a split location
  else {
    my $current_location = $parent->location;
    if ($current_location->can('add_sub_Location')) {
      $current_location->add_sub_Location($child->location);
    } else {
      eval "require Bio::Location::Split" unless Bio::Location::Split->can('add_sub_Location');
      my $new_location = Bio::Location::Split->new();
      $new_location->add_sub_Location($current_location);
      $new_location->add_sub_Location($child->location);
      $parent->location($new_location);
    }
  }

}

sub parse_attributes {

  my $self = shift;
  my $att  = shift;
  my @pairs =  map {my ($name,$value) = split /=/; [unescape($name) => unescape($value)] } split /;/,$att;
  my (%reserved,%unreserved);
  foreach (@pairs) {
    my $tag    = $_->[0];
    my @values = split /,/,$_->[1];

    if ($Special_attributes{$tag}) {  # reserved attribute
      push @{$reserved{$tag}},@values;
    } else {
      push @{$unreserved{$tag}},@values
    }
  }
  return (\%reserved,\%unreserved);

}

sub start_or_finish_sequence {

  my $self  = shift;
  my $seqid = shift;
  if (my $sl    = $self->{fasta_load}) {
    if (defined $sl->{seqid}) {
      $self->store->insert_sequence($sl->{seqid},$sl->{sequence},$sl->{offset});
      delete $self->{fasta_load};
    }
  }
  if (defined $seqid) {
    $self->{fasta_load} = {seqid  => $seqid,
			   offset => 0,
			   sequence => ''};
  }

}

sub load_sequence {

  my $self = shift;
  my $seq  = shift;
  my $sl   = $self->{fasta_load} or return;
  my $cs   = $self->seq_chunk_size;
  $sl->{sequence} .= $seq;
  while (length $sl->{sequence} >= $cs) {
    my $chunk = substr($sl->{sequence},0,$cs);
    $self->store->insert_sequence($sl->{seqid},$chunk,$sl->{offset});
    $sl->{offset} += length $chunk;
    substr($sl->{sequence},0,$cs) = '';
  }

}

sub msg {

  my $self = shift;
  my @msg  = @_;
  return unless $self->verbose;
  print STDERR @msg;

}

General documentation

The following read-only accessors return values passed or created during new():

 store()          the long-term Bio::DB::SeqFeature::Store object

 tmp_store()      the temporary Bio::DB::SeqFeature::Store object used
                    during loading

 sfclass()        the Bio::SeqFeatureI class

 fast()           whether fast loading is active

 seq_chunk_size() the sequence chunk size

 verbose()        verbose progress messages

The following methods are used internally and may be overidden by
subclasses.
    
    Return the default SeqFeatureI class (Bio::DB::SeqFeature).

This is an early version, so there are certainly some bugs. Please
use the BioPerl bug tracking system to report bugs.

bioperl,
Bio::DB::SeqFeature::Store,
Bio::DB::SeqFeature::Segment,
Bio::DB::SeqFeature::NormalizedFeature,
Bio::DB::SeqFeature::GFF3Loader,
Bio::DB::SeqFeature::Store::DBI::mysql,
Bio::DB::SeqFeature::Store::bdb

Lincoln Stein <lstein@cshl.org>.
Copyright (c) 2006 Cold Spring Harbor Laboratory.
This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself.