Toolbar

Summary

Bio::EnsEMBL::DBSQL::Funcgen::RegulatoryFeatureAdaptor - A database adaptor for fetching and
storing RegulatoryFeature objects.

Package variables

Included modules

Inherit

Bio::EnsEMBL::Funcgen::DBSQL::SetFeatureAdaptor

Synopsis

my $afa = $db->get_RegulatoryFeatureAdaptor();
my $features = $afa->fetch_all_by_Slice($slice);

Description

The RegulatoryFeatureAdaptor is a database adaptor for storing and retrieving
RegulatoryFeature objects.

Methods

Methods description

  Args       : None
  Example    : None
  Description: PROTECTED implementation of superclass abstract method.
               Returns a list of columns to use for queries.
  Returntype : List of strings
  Exceptions : None
  Caller     : Internal
  Status     : At Risk

  Example    : my $regf_featureset = $self->_get_current_FeatureSet;
  Description: Convenience method to get and test the current
  Returntype : Bio::EnsEMBL::Funcgen::FeatureSet
  Exceptions : Throws is FeatureSet is not available
  Caller     : general
  Status     : at risk

  Args       : None
  Example    : None
  Description: PROTECTED implementation of superclass abstract method.
               Returns an additional table joining constraint to use for
			   queries.
  Returntype : List
  Exceptions : None
  Caller     : Internal
  Status     : At Risk

  Arg [1]    : DBI statement handle object
  Example    : None
  Description: PROTECTED implementation of superclass abstract method.
               Creates RegulatoryFeature objects from an executed DBI statement
			   handle.
  Returntype : Listref of Bio::EnsEMBL::RegulatoryFeature objects
  Exceptions : None
  Caller     : Internal
  Status     : At Risk

  Args       : None
  Example    : None
  Description: PROTECTED implementation of superclass abstract method.
               Returns the names and aliases of the tables to use for queries.
  Returntype : List of listrefs of strings
  Exceptions : None
  Caller     : Internal
  Status     : At Risk

  Arg [1]    : Bio::EnsEMBL::Slice
  Example    : my $slice = $sa->fetch_by_region('chromosome', '1');
               my $features = $ofa->fetch_by_Slice($slice);
  Description: Retrieves a list of features on a given slice, specific for the current 
               default RegulatoryFeature set.
  Returntype : Listref of Bio::EnsEMBL::RegulatoryFeature objects
  Exceptions : None
  Caller     : General
  Status     : At Risk

  Arg [1]    : String $stable_id - The stable id of the regulatory feature to retrieve
  Arg [2]    : optional list of FeatureSets
  Example    : my $rf = $rf_adaptor->fetch_by_stable_id('ENSR00000309301');
  Description: Retrieves a regulatory feature via its stable id.
  Returntype : Array ref of Bio::EnsEMBL::Funcgen::RegulatoryFeature objects
  Exceptions : throws if no stable ID provided or FeatureSets aren't valid
               warns if not FeatureSets defined
  Caller     : general
  Status     : at risk

  Arg [1]    : String $stable_id - The stable id of the regulatory feature to retrieve
  Arg [2]    : optional list of FeatureSets
  Example    : my $rf = $rf_adaptor->fetch_by_stable_id('ENSR00000309301');
  Description: Retrieves a regulatory feature via its stable id.
  Returntype : Bio::EnsEMBL::Funcgen::RegulatoryFeature
  Exceptions : none
  Caller     : general
  Status     : Stable

  Args       : List of Bio::EnsEMBL::Funcgen::RegulatoryFeature objects
  Example    : $ofa->store(@features);
  Description: Stores given RegulatoryFeature objects in the database. Should only
               be called once per feature because no checks are made for
			   duplicates. Sets dbID and adaptor on the objects that it stores.
  Returntype : Listref of stored RegulatoryFeatures
  Exceptions : Throws if a list of RegulatoryFeature objects is not provided or if
               the Analysis, CellType and FeatureType objects are not attached or stored.
               Throws if analysis of set and feature do not match
               Warns if RegulatoryFeature already stored in DB and skips store.
  Caller     : General
  Status     : At Risk

Methods code

sub _columns {

  my $self = shift;
  
  return qw(
			rf.regulatory_feature_id rf.seq_region_id
			rf.seq_region_start      rf.seq_region_end
			rf.seq_region_strand     rf.bound_seq_region_start
			rf.bound_seq_region_end  rf.display_label
			rf.feature_type_id       rf.feature_set_id
			rf.stable_id             ra.attribute_feature_id
			ra.attribute_feature_table
	   );

}

sub _get_current_FeatureSet {

  my $self = shift;

  my $fset = $self->db->get_FeatureSetAdaptor->fetch_by_name('RegulatoryFeatures');

  if(! defined $fset){
	warn('Could not retrieve current RegulatoryFeatures FeatureSet');
  }
  
  return $fset;

}

sub _left_join {

  my $self = shift;
	
  return (['regulatory_attribute', 'rf.regulatory_feature_id = ra.regulatory_feature_id']);

}

sub _objs_from_sth {

  my ($self, $sth, $mapper, $dest_slice) = @_;

  
	#For EFG this has to use a dest_slice from core/dnaDB whether specified or not.
	#So if it not defined then we need to generate one derived from the species_name and schema_build of the feature we're retrieving.

	# This code is ugly because caching is used to improve speed

	#my $sa = $self->db->get_SliceAdaptor();

	
	my ($sa, $reg_feat);#, $old_cs_id);
	$sa = ($dest_slice) ? $dest_slice->adaptor->db->get_SliceAdaptor() : $self->db->get_SliceAdaptor();
	#don't really need this if we're using DNADBSliceAdaptor?

	#Some of this in now probably overkill as we'll always be using the DNADB as the slice DB
	#Hence it should always be on the same coord system
	#my $aa = $self->db->get_AnalysisAdaptor();
	my $ft_adaptor = $self->db->get_FeatureTypeAdaptor();
	my $fset_adaptor = $self->db->get_FeatureSetAdaptor();
	my (@features, @reg_attrs, $seq_region_id);
	my (%fset_hash, %slice_hash, %sr_name_hash, %sr_cs_hash, %ftype_hash);
	my $skip_feature = 0;

	my %feature_adaptors = (
							'annotated' => $self->db->get_AnnotatedFeatureAdaptor(),
							#external?
						   );
	
	
	my (
	    $dbID,                  $efg_seq_region_id,
	    $seq_region_start,      $seq_region_end,
	    $seq_region_strand,     $bound_seq_region_start,
		$bound_seq_region_end,             $display_label,
		$ftype_id,              $fset_id,
		$stable_id,             $attr_id,
		$attr_type
	);

	$sth->bind_columns(\$
					   dbID,\$                  efg_seq_region_id,\$
					   seq_region_start,\$      seq_region_end,\$
					   seq_region_strand,\$     bound_seq_region_start,\$
					   bound_seq_region_end,\$  display_label,\$
					   ftype_id,\$              fset_id,\$
					   stable_id,\$             attr_id,\$
					   attr_type
					  );


	#This needs doing properly!!
	#my $epsth = $self->prepare("SELECT experiment_id 
    #                                FROM experiment_prediction 
    #                                WHERE regulatory_feature_id = ?");

	my ($asm_cs, $cmp_cs, $asm_cs_name);
	my ($asm_cs_vers, $cmp_cs_name, $cmp_cs_vers);
  
	if ($mapper) {
		$asm_cs      = $mapper->assembled_CoordSystem();
		$cmp_cs      = $mapper->component_CoordSystem();
		$asm_cs_name = $asm_cs->name();
		$asm_cs_vers = $asm_cs->version();
		$cmp_cs_name = $cmp_cs->name();
		$cmp_cs_vers = $cmp_cs->version();
	  }
	
	my ($dest_slice_start, $dest_slice_end);
	my ($dest_slice_strand, $dest_slice_length, $dest_slice_sr_name);

	if ($dest_slice) {
		$dest_slice_start   = $dest_slice->start();
		$dest_slice_end     = $dest_slice->end();
		$dest_slice_strand  = $dest_slice->strand();
		$dest_slice_length  = $dest_slice->length();
		$dest_slice_sr_name = $dest_slice->seq_region_name();
	}

	
  my $slice;
	
  FEATURE: while ( $sth->fetch() ) {

	  if(! $reg_feat || ($reg_feat->dbID != $dbID)){
	
		if($skip_feature){
		  undef $reg_feat;#so we don't duplicate the push for the feature previous to the skip feature
		  $skip_feature = 0;
		}

		if($reg_feat){
		  
		  $reg_feat->regulatory_attributes(\@reg_attrs);# if @reg_attrs;
		  push @features, $reg_feat;
		  undef @reg_attrs;
		}

		#Hack to get not NUL 0 autoinserted values to work
		$bound_seq_region_start = undef if ! $bound_seq_region_start;
		$bound_seq_region_end   = undef if ! $bound_seq_region_end;


	    #Need to build a slice adaptor cache here?
	    #Would only ever want to do this if we enable mapping between assemblies??
	    #Or if we supported the mapping between cs systems for a given schema_build, which would have to be handled by the core api
	    
		#this should only be done once for each regulatory_feature_id
		
		
		#get core seq_region_id
		$seq_region_id = $self->get_core_seq_region_id($efg_seq_region_id);

		if(! $seq_region_id){
		  warn "Cannot get slice for eFG seq_region_id $efg_seq_region_id\n".
			"The region you are using is not present in the current dna DB";
		  next;
		}

	    #if($old_cs_id && ($old_cs_id+ != $cs_id)){
	    #  throw("More than one coord_system for feature query, need to implement SliceAdaptor hash?");
	    #}
	    #$old_cs_id = $cs_id;
	    #Need to make sure we are restricting calls to Experiment and channel(i.e. the same coord_system_id)
	    
		#Get the FeatureSet object
		$fset_hash{$fset_id} = $fset_adaptor->fetch_by_dbID($fset_id) if(! exists $fset_hash{$fset_id});
		
		
	    # Get the slice object
	    $slice = $slice_hash{'ID:'.$seq_region_id};
	    
	    if (!$slice) {
	      $slice                            = $sa->fetch_by_seq_region_id($seq_region_id);
	      $slice_hash{'ID:'.$seq_region_id} = $slice;
	      $sr_name_hash{$seq_region_id}     = $slice->seq_region_name();
	      $sr_cs_hash{$seq_region_id}       = $slice->coord_system();
	    }
	    
	    my $sr_name = $sr_name_hash{$seq_region_id};
	    my $sr_cs   = $sr_cs_hash{$seq_region_id};
	    
	    # Remap the feature coordinates to another coord system if a mapper was provided
	    if ($mapper) {
	      
	      throw("Not yet implmented mapper, check equals are Funcgen calls too!");
	      
	      ($sr_name, $seq_region_start, $seq_region_end, $seq_region_strand)
			= $mapper->fastmap($sr_name, $seq_region_start, $seq_region_end, $seq_region_strand, $sr_cs);
	
	      # Skip features that map to gaps or coord system boundaries
		  if(! defined $sr_name){
			$skip_feature = 1;
			next FEATURE;
		  }
	      
	      # Get a slice in the coord system we just mapped to
	      if ( $asm_cs == $sr_cs || ( $cmp_cs != $sr_cs && $asm_cs->equals($sr_cs) ) ) {
		$slice = $slice_hash{"NAME:$sr_name:$cmp_cs_name:$cmp_cs_vers"}
		  ||= $sa->fetch_by_region($cmp_cs_name, $sr_name, undef, undef, undef, $cmp_cs_vers);
	      } else {
		$slice = $slice_hash{"NAME:$sr_name:$asm_cs_name:$asm_cs_vers"}
		  ||= $sa->fetch_by_region($asm_cs_name, $sr_name, undef, undef, undef, $asm_cs_vers);
	      }
	    }
	    
	    # If a destination slice was provided convert the coords
	    # If the destination slice starts at 1 and is forward strand, nothing needs doing
	    if ($dest_slice) {

	      unless ($dest_slice_start == 1 && $dest_slice_strand == 1) {
			
			#can remove the if $bound_seq_region_start/end once we have updated all reg feature entries and store API

			if ($dest_slice_strand == 1) {
			  $seq_region_start       = $seq_region_start - $dest_slice_start + 1;
			  $seq_region_end         = $seq_region_end   - $dest_slice_start + 1;

			  #if as we never have a seq_region start of 0;
			  $bound_seq_region_start = $bound_seq_region_start - $dest_slice_start + 1 if $bound_seq_region_start;
			  $bound_seq_region_end   = $bound_seq_region_end   - $dest_slice_start + 1 if $bound_seq_region_end;
			  
			} 
			else {
			  my $tmp_seq_region_start       = $seq_region_start;
			  my $tmp_bound_seq_region_start = $bound_seq_region_start;
			  $seq_region_start        = $dest_slice_end - $seq_region_end       + 1;
			  $seq_region_end          = $dest_slice_end - $tmp_seq_region_start + 1;
			  $bound_seq_region_start  = $dest_slice_end - $bound_seq_region_end + 1 if $bound_seq_region_end;
			  $bound_seq_region_end    = $dest_slice_end - $tmp_bound_seq_region_start + 1 if $bound_seq_region_start;
			  $seq_region_strand      *= -1;
			}
	      }
	      
	      # Throw away features off the end of the requested slice
		  #Do not account for bounds here.
	      if ($seq_region_end < 1 || $seq_region_start > $dest_slice_length
			  || ( $dest_slice_sr_name ne $sr_name )){
			$skip_feature = 1;
			next FEATURE;
		  }
	      
	      $slice = $dest_slice;
	    }
	    

		my ($reg_type, $reg_attrs, $ftype);
	   		
		if(defined $ftype_id){
		  $ftype = $ft_adaptor->fetch_by_dbID($ftype_id);
		}
		
		$reg_feat = Bio::EnsEMBL::Funcgen::RegulatoryFeature->new_fast
		  ({
			'start'          => $seq_region_start,
			'end'            => $seq_region_end,
			'bound_start'    => $bound_seq_region_start,
			'bound_end'      => $bound_seq_region_end,
			'strand'         => $seq_region_strand,
			'slice'          => $slice,
			'analysis'       => $fset_hash{$fset_id}->analysis(),
			'adaptor'        => $self,
			'dbID'           => $dbID,
			'display_label'  => $display_label,
			'feature_set'    => $fset_hash{$fset_id},
			'feature_type'   => $ftype,
			#'regulatory_attributes' => $reg_attrs,
			'stable_id'      => $stable_id,
		   });
	  }
	
  
	  #populate attributes array
	  if(defined $attr_id  && ! $skip_feature){

		#These will all be fetched on their native slice, not necessarily the slice we have fetched this
		#reg feature on, hence we need to map the features to the current slice
		#otherwise the bounds may get messed up

		my $attr = $feature_adaptors{$attr_type}->fetch_by_dbID($attr_id);

		#now need to reset start and ends for the current slice
		#This is not redefining the slice, so we may get minus start values
		#grab the seq_region_start/ends here first
		#as resetting directly causes problems
		my $attr_sr_start = $attr->seq_region_start;
		my $attr_sr_end = $attr->seq_region_end;
 		$attr->slice($slice);

		if($slice->strand ==1){
		  $attr->start($attr_sr_start - $slice->start +1);
		  $attr->end($attr_sr_end - $slice->start +1);	
		}else{
		  $attr->start($slice->end - $attr_sr_end +1);
		  $attr->end($slice->end - $attr_sr_start +1);	
		}
		
		push @reg_attrs, $attr;
	  }
	}

  #handle last record
  if($reg_feat){
	$reg_feat->regulatory_attributes(\@reg_attrs);# if(@reg_attrs);
	push @features, $reg_feat;
  }
  

  return\@ features;

}

sub _tables {

  my $self = shift;
	
  return (
		  [ 'regulatory_feature', 'rf' ],
		  [ 'feature_set', 'fs'],
		  [ 'regulatory_attribute', 'ra'],
		 );

}

sub fetch_all_by_Slice {

  my ($self, $slice) = @_;
	
  my $fset = $self->_get_current_FeatureSet;

  return (defined $fset) ? $self->fetch_all_by_Slice_FeatureSets($slice, [$self->_get_current_FeatureSet]) : undef;
}



1;

}

sub fetch_all_by_stable_id_FeatureSets {

  my ($self, $stable_id, @fsets) = @_;

  #Standard implementation exposes logic name as a parameter
  #But it will always be RegulatoryFeature/Build

  throw('Must provide a stable ID') if ! defined $stable_id;
  $stable_id =~ s/ENSR0*//;

  my $constraint = 'rf.stable_id='.$stable_id;


  if(@fsets){#Get current
	
	#need to catch empty array an invalid FeatureSets
	if(scalar(@fsets == 0)){
	  warning("You have not specified any FeatureSets to fetch the RegulatoryFeature from, defaulting to all");
	}
	else{

	  #validate FeatureSets
	  map { 
		if(! (ref($_) && $_->isa('Bio::EnsEMBL::Funcgen::FeatureSet') && $_->dbID)){
		  throw 'You must provide a valid stored Bio::EnsEMBL::Funcgen::FeatureSet';
		}} @fsets;



		 
	  if(scalar(@fsets) == 1){
		$constraint .= ' and rf.feature_set_id='.$fsets[0]->dbID;
	  }else{
		$constraint .= ' AND rf.feature_set_id IN ('.join(', ', (map $_->dbID, @fsets)).')';
	  } 
	}
  }

  return $self->generic_fetch($constraint);

}

sub fetch_by_stable_id {

  my ($self, $stable_id) = @_;

  my $fset = $self->_get_current_FeatureSet;

  return (defined $fset) ? $self->fetch_all_by_stable_id_FeatureSets($stable_id, $self->_get_current_FeatureSet)->[0] : undef;

}

sub store {

  my ($self, @rfs) = @_;
	
  if (scalar(@rfs) == 0) {
	throw('Must call store with a list of RegulatoryFeature objects');
  }

  my $sth = $self->prepare("
		INSERT INTO regulatory_feature (
			seq_region_id,         seq_region_start,
			seq_region_end,        bound_seq_region_start,
			bound_seq_region_end,  seq_region_strand,
            display_label,         feature_type_id,
            feature_set_id,        stable_id
		) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?)");
  
  my $sth2 = $self->prepare("
		INSERT INTO regulatory_attribute (
              regulatory_feature_id, attribute_feature_id, attribute_feature_table
		) VALUES (?, ?, ?)");
  
  my $db = $self->db();
  
  foreach my $rf (@rfs) {
	
	if( ! ref $rf || ! $rf->isa('Bio::EnsEMBL::Funcgen::RegulatoryFeature') ) {
	  throw('Feature must be an RegulatoryFeature object');
	}
	
	if ( $rf->is_stored($db) ) {
	  #does not accomodate adding Feature to >1 feature_set
	  warning('RegulatoryFeature [' . $rf->dbID() . '] is already stored in the database');
	  next;
	}
	
	#Have to do this for Analysis separately due to inheritance
	if ( ! $rf->analysis->is_stored($db)) {
	  throw('A stored Bio::EnsEMBL::Analysis must be attached to the RegulatoryFeature objects to be stored.');
	}
	
	if (! $rf->feature_set->is_stored($db)) {
	  throw('A stored Bio::EnsEMBL::Funcgen::FeatureSet must be attached to the RegulatoryFeature objects to be stored.');
	}
	
	if (! $rf->feature_type->is_stored($db)) {
	  throw('A stored Bio::EnsEMBL::Funcgen::FeatureType must be attached to the RegulatoryFeature objects to be stored.');
	}
	  


	#sanity check analysis matches feature_set analysis
	if($rf->analysis->dbID() != $rf->feature_set->analysis->dbID()){
	  throw("RegulatoryFeature analysis(".$rf->analysis->logic_name().") does not match FeatureSet analysis(".$rf->feature_set->analysis->logic_name().")\n".
			"Cannot store mixed analysis sets");
	}

	#Complex analysis to be stored as one in analysis table, or have feature_set_prediciton link table?
	#Or only have single analysis feature which can contribute to multi analysis "regulons"
	#Or can we have multiple entries in feature_set with the same id but different analyses?
	#This would still not be specific for each feature, nor would the regulatory_feature analysis_id
	#reflect all the combined analyses.  Maybe just the one which contributed most?
	
	
	my $seq_region_id;
	($rf, $seq_region_id) = $self->_pre_store($rf);
	$rf->adaptor($self);#Set adaptor first to allow attr feature retreival for bounds
	#This is only required when storing, 
	
	$sth->bind_param(1, $seq_region_id,             SQL_INTEGER);
	$sth->bind_param(2, $rf->start(),               SQL_INTEGER);
	$sth->bind_param(3, $rf->end(),                 SQL_INTEGER);
	$sth->bind_param(4, $rf->bound_start(),         SQL_INTEGER);
	$sth->bind_param(5, $rf->bound_end(),           SQL_INTEGER);
	$sth->bind_param(6, $rf->strand(),              SQL_TINYINT);
	$sth->bind_param(7, $rf->{'display_label'},     SQL_VARCHAR);#Direct access so we always store the binary string
	$sth->bind_param(8, $rf->feature_type->dbID(),  SQL_INTEGER);
	$sth->bind_param(9, $rf->feature_set->dbID(),   SQL_INTEGER);
	$sth->bind_param(10, $rf->{'stable_id'},        SQL_INTEGER);
	
	$sth->execute();
	$rf->dbID( $sth->{'mysql_insertid'} );

	my $table_type;
	my %attrs = %{$rf->_attribute_cache()};

	foreach my $table(keys %attrs){
	  ($table_type = $table)  =~ s/_feature//;

	  foreach my $id(keys %{$attrs{$table}}){
		$sth2->bind_param(1, $rf->dbID,   SQL_INTEGER);
		$sth2->bind_param(2, $id,         SQL_INTEGER);
		$sth2->bind_param(3, $table_type, SQL_VARCHAR);
		$sth2->execute();
	  }
	}
  }
  
  return\@ rfs;

}

General documentation

This module was created by Nathan Johnson.
This module is part of the Ensembl project: /

Post comments or questions to the Ensembl development list: ensembl-dev@ebi.ac.uk