initial import of external_feature parsers

misc-scripts/regulatory_regions/RegulatoryFeatureParser/BaseParser.pm

@@ -4,18 +4,187 @@ use strict;
 
 use Bio::EnsEMBL::DBSQL::DBAdaptor;
 use Bio::EnsEMBL::SimpleFeature;
+use Bio::EnsEMBL::Utils::Argument qw( rearrange );
+use Bio::EnsEMBL::Utils::Exception qw( throw );
+use Bio::EnsEMBL::Funcgen::FeatureSet;
+use Bio::EnsEMBL::Funcgen::FeatureType;
+use Bio::EnsEMBL::Analysis;
+
+# Base functionality for external_feature parsers
+
+sub new {
+  my $caller = shift;
+  my $class = ref($caller) || $caller;
+  my $self = {};
+
+  #my $self = {(
+#			  feature_sets => {(
+#								vista   => {('VISTA enhancer set' => undef)},
+#								cisred  => {(
+#											 'cisRED search regions' => undef,
+#											 'cisRED group motifs'   => undef,
+#											)},
+#								miranda => {('miRanda miRNA' => undef)},
+#							   )},
+#			  )};
+  bless $self, $class;
+
+  #validate and set type, analysis and feature_set here
+  my ($type, $db, $clobber, $archive) = rearrange(['TYPE', 'DB', 'CLOBBER', 'ARCHIVE'], @_);
+  
+  throw('You must define a type of external_feature to import') if(! defined $type);
+
+  if (! ($db && ref($db) &&
+		 $db->isa('Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor'))){
+	throw('You must provide a valid Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor');
+  }
+
+  throw('You can only specify either -clobber or -archive, but not both') if($clobber && $archive);
+
+
+  $self->{'db'} = $db;
+  $self->{type} = $type;
+  $self->{'clobber'} = $clobber if defined $clobber;
+  $self->{'archive'} = $archive if defined $archive;
+  
+  print ":: Parsing and loading $type ExternalFeatures";
+
+  return $self;
+
+}
+
+sub db{
+  my ($self, $db) = @_;
 
+  if($db){
+
+	if(! $db->isa('Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor')){
+	  throw('You must prove a Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor');
+	}
+	
+	$self->{'db'} = $db;
+  }
+
+  return $self->{'db'};
+}
+
+
+sub set_feature_sets{
+  my $self = shift;
+
+  throw('Must provide a set feature_set config hash') if ! defined $self->{'feature_sets'};
+
+
+  my $fset_adaptor = $self->db->get_FeatureSetAdaptor;
+  my $analysis_adaptor = $self->db->get_AnalysisAdaptor;
+  
+  foreach my $fset_name(keys %{$self->{feature_sets}}){
+	
+	#check feature_type is validated?
+
+	my $fset = $fset_adaptor->fetch_by_name($fset_name);
+
+	#what about data sets?
+	#we don't need data sets for external_sets!
+
+	if(defined $fset){
+
+	  if($self->{'clobber'}){
+		#Need to clobber any DBEntries first!!!
+
+		if($self->{'type'} eq 'cisred'){
+		  my @ext_feat_ids = 	@{$self->db->dbc->db_handle->selectall_arrayref('select external_feature_id from external_feature where feature_set_id='.$fset->dbID)};
+		  
+		  if(@ext_feat_ids){
+
+			my ($core_ext_dbid) = @{$self->db->dbc->db_handle->selectall_arrayref('select external_db_id from external_db where db_name="core"')};
+
+			if($core_ext_dbid){
+			  #double table delete?
+			  my $sql = "delete x, ox from object_xref ox, xref x where ox.ensembl_object_type='ExternalFeature' and x.external_db_id=$core_ext_dbid and ox.xref_id=x.xref_id and ensembl_id in(".join(', ', @ext_feat_ids).')';
+			  print ":: Clobbering xrefs for $fset_name\n";
+			  $self->db->dbc->do($sql);
+			}
+		  }
+	
+		  print ":: Clobbering old features for external feature_set:\t$fset_name\n";
+		  my $sql = 'delete from external_feature where feature_set_id='.$fset->dbID;
+		  $self->db->dbc->do($sql);
+		}
+
+	  }
+	  elsif($self->{'archive'}){
+		my $archive_fset =  $fset_adaptor->fetch_by_name($fset_name."_v".$self->{'archive'});
+
+		if(defined $archive_fset){
+		  throw("You are trying to create an archive external feature_set which already exists:\t${fset_name}_v".$self->{archive});
+		}
+
+		my $sql = "UPDATE feature_set set name='$fset_name}_v".$self->{archive}."' where name='$fset_name'";
+		$self->db->dbc->do($sql);
+		undef $fset;
+	  }else{
+		throw("You are trying to create an external feature_set which already exists:\t$fset_name\nMaybe to want to clobber or archive?");
+	  }
+	}
+
+	if(!defined $fset){
+	  #don't need to use RNAFeatureType here as this is the setwide generic feature_type
+	  #or do we have separate tables for external_feature and external_rna_feature?
+
+	  #validate analysis first
+	  my $analysis = $analysis_adaptor->fetch_by_logic_name($self->{'feature_sets'}{$fset_name}{'analysis'}{'-logic_name'});
+	  if(! defined $analysis){
+		
+		print ':: Analysis '.$self->{'feature_sets'}{$fset_name}{'analysis'}{'-logic_name'}.
+		  " not found, storing from config hash\n";
+		
+		$analysis_adaptor->store(Bio::EnsEMBL::Analysis->new
+								 (
+								  %{$self->{'feature_sets'}{$fset_name}{'analysis'}}
+								  #-logic_name    => $self->{'feature_sets'}{$fset_name}{'analysis'}{'logic_name'},
+								  #-description   => $self->{'feature_sets'}{$fset_name}{'analysis'}{'description'},
+								  #-display_label => $self->{'feature_sets'}{$fset_name}{'analysis'}{'display_label'},
+								  #-diplayable    => $self->{'feature_sets'}{$fset_name}{'analysis'}{'displayable'},
+								 ));
+		
+		$analysis = $analysis_adaptor->fetch_by_logic_name($self->{'feature_sets'}{$fset_name}{'analysis'});
+	  }
+
+	  #replace hash config with object
+	  $self->{'feature_sets'}{$fset_name}{'analysis'} = $analysis;
+
+	  $fset = Bio::EnsEMBL::Funcgen::FeatureSet->new(
+													 -name         => $fset_name,
+													 -type         => 'external',
+													 -analysis     => $self->{'feature_sets'}{$fset_name}{'analysis'},
+													 -feature_type => ${$self->{'feature_sets'}{$fset_name}{'feature_type'}},
+													);
+
+	  ($fset) = @{$self->db->get_FeatureSetAdaptor->store($fset)};
+	}
+
+	#Now replace config hash with object
+	$self->{feature_sets}{$fset_name} = $fset;
+  }
+
+  return;
+}
 
-# Base functionality for regualatory feature parsers
 
 
 # --------------------------------------------------------------------------------
 # Delete existing regulatory features etc
 
+#This needn't be done in line any more as we can just make the old set not displayable
+
 sub delete_existing {
 
   my ($db_adaptor, $type) = @_;
 
+
+  throw('delete existing is deprecated');
+
   my $t = lc($type);
 
   print "Deleting existing features & related data for type $type\n";
@@ -52,27 +221,61 @@ sub delete_existing {
 # --------------------------------------------------------------------------------
 # Check that the type specified corresponds to an analysis
 
-sub validate_type {
-
-  my ($db_adaptor, $type) = @_;
+sub validate_and_store_feature_types{
+  my $self = shift;
 
   # LBL enhancers have both positive and negative varieties
-  if (lc($type) eq 'enhancer') {
-    return (validate_type($db_adaptor, 'enhancer_positive') && validate_type($db_adaptor, 'enhancer_negative'));
-  }
 
-  my $sth = $db_adaptor->dbc->prepare("SELECT analysis_id FROM analysis WHERE LOWER(logic_name)=?");
+  #feature type class/naming is logical but not intuitive
+  #+-----------------+-------------------------+----------+----------------------------------------------------+
+  #| feature_type_id | name                    | class    | description                                        |
+  #+-----------------+-------------------------+----------+----------------------------------------------------+
+  #|          398680 | VISTA Enhancer          | Enhancer | Enhancer identified by positive VISTA assay        |
+  #|          398681 | VISTA Target - Negative | Region   | Enhancer negative region identified by VISTA assay |
+  #+-----------------+-------------------------+----------+----------------------------------------------------+
 
-  $sth->execute(lc($type));
-  if ($sth->fetchrow_array()) {
-    print "Type $type is valid\n";
-  } else {
-    print "Type $type is not valid - is there an entry for $type in the analysis table?\n";
-    return 0;
-  }
 
-  return 1;
+  #if (lc($type) eq 'VISTA') {
+  #  return (validate_type($db_adaptor, 'VISTA Enhancer') && validate_type($db_adaptor, 'VISTA Target - Negative'));
+  #}
+
+  #my $sth = $self->db->dbc->prepare("SELECT analysis_id FROM analysis WHERE logic_name=?");
+
+  
+
+  #remove lc as mysql doesn't care about case
+  #$sth->execute($type);
+  #if ($sth->fetchrow_array()) {
+  #  print "Type $type is valid\n";
+  #} else {
+  #  print "Type $type is not valid - is there an entry for $type in the analysis table?\n";
+  #  return 0;
+  #}
+
+  my $ftype_adaptor = $self->db->get_FeatureTypeAdaptor;
+
+  foreach my $ftype_name(keys %{$self->{'feature_types'}}){
+
+	my $ftype = $ftype_adaptor->fetch_by_name($ftype_name);
 
+
+	if(! defined $ftype){
+	  print ":: FeatureType '".$ftype_name."' for external feature_set ".$self->{'type'}." not present\n".
+		":: Storing using type hash definitions\n";
+	
+	  $ftype = Bio::EnsEMBL::Funcgen::FeatureType->new(
+													   -name => $ftype_name,
+													   -class => $self->{'feature_types'}{$ftype_name}{'class'},
+													   -description => $self->{'feature_types'}{$ftype_name}{'description'},
+													  );
+	  ($ftype) = @{$ftype_adaptor->store($ftype)};
+	}
+
+	#Replace hash config with object
+	$self->{'feature_types'}{$ftype_name} = $ftype;
+  }
+
+  return;
 }
 
 
@@ -80,10 +283,11 @@ sub validate_type {
 # Find the maximum existing ID in a table.
 sub find_max_id {
 
-  my ($self, $db_adaptor, $table) = @_;
+  my ($self, $table) = @_;
 
-  my $row = @{$db_adaptor->dbc->db_handle->selectall_arrayref("SELECT MAX(${table}_id) FROM $table")}[0];
+  my $row = @{$self->dbc->db_handle->selectall_arrayref("SELECT MAX(${table}_id) FROM $table")}[0];
   my $max_id = @{$row}[0];
+
   if (!defined $max_id) {
     print "No existing ${table}_ids, will start from 1\n";
     $max_id = 0;
@@ -100,36 +304,46 @@ sub find_max_id {
 # Return hashref keyed on {$type}{$stable_id}
 # Type is always all lower case
 
-sub build_stable_id_cache {
+sub build_display_name_cache {
+
+  my ($self, @types) = @_;
 
-  my ($self, $db_adaptor) = @_;
+  my %display_name_cache;
 
-  my %stable_id_to_internal_id;
+  #validate types here?
 
-  foreach my $type ('gene', 'transcript', 'translation') { # Add exon here if required
+  foreach my $type (@types){#'gene', 'transcript', 'translation') { # Add exon here if required
 
-    print "Caching stable ID -> internal ID links for ${type}s\n";
+	#No display_xref for translations!!
+	throw('Cannot build cache for tranlsations, just use stable_id, check your implementatio') if $type eq 'translation';
+
+    print ":: Caching stable ID ->  display_name links for ${type}s\n";
 
     my $count = 0;
+	my $sql;
 
-    my $sth = $db_adaptor->dbc->prepare("SELECT ${type}_id, stable_id FROM ${type}_stable_id");
-    $sth->execute();
-    my ($internal_id, $stable_id);
-    $sth->bind_columns(\$internal_id, \$stable_id);
+	print "dnadb is ".Data::Dumper::Dumper($self->db->dnadb->dbc);
+	
+	#could do a fetchall_hashref here?
+	my $sth= $self->db->dnadb->dbc->prepare("SELECT t.stable_id, x.display_label FROM s.${type}_stable_id, t.${type} left join x.xref on t.display_xref_id=x.xref_id and s.${type}_id=t.${type}_id");
 
-    while ($sth->fetch) {
+	
 
-      $stable_id_to_internal_id{$type}{$stable_id} = $internal_id;
-      $count++;
+    $sth->execute();
+    my ($stable_id, $display_name);
+    $sth->bind_columns(\$stable_id, \$display_name);
 
+    while ($sth->fetch) {
+	  #removed type
+      $display_name_cache{$stable_id} = $display_name;
+      $count++ if $display_name;#will this be ! NULL?
     }
 
-    print "Got $count $type stable ID -> internal ID mappings\n";
+    print ":: Got $count $type stable ID -> display_name mappings\n";
 
   }
 
-  return \%stable_id_to_internal_id;
-
+  return %display_name_cache;
 }
 
 # --------------------------------------------------------------------------------
@@ -137,9 +351,18 @@ sub build_stable_id_cache {
 
 sub upload_features_and_factors {
 
-  my ($self, $db_adaptor, $objects) = @_;
+  my ($self, $objects) = @_;
+
+
+  throw("upload_features_and_factors is deprecated");
+
+  my $dbc = $self->db->dbc;
+
 
-  my $dbc = $db_adaptor->dbc;
+
+  # we need to convert this to use the external_feature API
+  #this means we need to build external_features in the other modules and test for each RNA_feature_type/motif_feature_type?
+  #keep as feature_types for now, maybe next release?
 
   my $feature_sth = $dbc->prepare("INSERT INTO regulatory_feature (regulatory_feature_id, name, seq_region_id, seq_region_start, seq_region_end, seq_region_strand, analysis_id, regulatory_factor_id) VALUES(?,?,?,?,?,?,?,?)");
   my $factor_sth = $dbc->prepare("INSERT INTO regulatory_factor (regulatory_factor_id, name, type) VALUES(?,?,?)");
@@ -211,65 +434,88 @@ sub upload_features_and_factors {
 # --------------------------------------------------------------------------------
 # Project a feature from one slice to another
 sub project_feature {
-
-  my ($self, $start, $end, $strand, $chr, $slice, $analysis, $new_assembly, $slice_adaptor, $label) = @_;
+  my ($self, $feat, $new_assembly) = @_;
 
   # just use a SimpleFeature for convenience
-  my $feat = Bio::EnsEMBL::SimpleFeature->new
-    (-start    => $start,
-     -end      => $end,
-     -strand   => $strand,
-     -slice    => $slice,
-     -analysis => $analysis,
-     -display_label => $label,
-     -score   => 0);
+  #my $feat = Bio::EnsEMBL::SimpleFeature->new
+  #  (-start    => $start,
+  #   -end      => $end,
+  #   -strand   => $strand,
+  #   -slice    => $slice,
+  #   -analysis => $analysis,
+  #   -display_label => $label,
+  #   -score   => 0);
 
   # project feature to new assembly
   my $feat_slice = $feat->feature_Slice;
   my @segments = @{ $feat_slice->project('chromosome', $new_assembly) };
 
-  next unless (@segments);
 
-  next if (scalar(@segments) > 1);
+  #next what?
+  #next unless (@segments);
+  #next if (scalar(@segments) > 1);
+
+  if(! @segments){
+	print "Failed to project feature:\t".$feat->display_label."\n";
+  }
+  elsif(scalar(@segments) >1){
+	print "Failed to project feature to distinct location:\t".$feat->display_label."\n";
+	return;
+  }
 
   my $proj_slice = $segments[0]->to_Slice;
-  next unless ($feat_slice->length == $proj_slice->length);
+
+  if($feat_slice->length != $proj_slice->length){
+	print "Failed to project feature to comparable length region:\t".$feat->display_label."\n";
+	return;
+  }
 
 
   # everything looks fine, so adjust the coords of the feature
   $feat->start($proj_slice->start);
   $feat->end($proj_slice->end);
   $feat->strand($proj_slice->strand);
-  my $slice_new_asm = $slice_adaptor->fetch_by_region('chromosome', $chr, undef, undef, undef, $new_assembly);
+  my $slice_new_asm = $self->slice_adaptor->fetch_by_region('chromosome', $proj_slice->seq_region_name, undef, undef, undef, $new_assembly);
   $feat->slice($slice_new_asm);
 
   return $feat;
 
 }
 
+sub slice_adaptor{
+  my $self = shift;
+
+  if(! defined $self->{'slice_adaptor'}){
+	$self->{'slice_adaptor'} = $self->db->get_SliceAdaptor;
+  }
+  
+  return $self->{'slice_adaptor'};
+}
+
+
 # --------------------------------------------------------------------------------
 
-sub get_blank_factor_id () {
+#sub get_blank_factor_id () {
 
-  my ($self, $db_adaptor) = @_;
+#  my ($self, $db_adaptor) = @_;
 
-  my $sth = $db_adaptor->dbc->prepare("SELECT regulatory_factor_id FROM regulatory_factor WHERE name=''");
-  $sth->execute();
+#  my $sth = $db_adaptor->dbc->prepare("SELECT regulatory_factor_id FROM regulatory_factor WHERE name=''");
+#  $sth->execute();
 
-  my ($factor_id) = $sth->fetchrow_array();
+#  my ($factor_id) = $sth->fetchrow_array();
 
-  if ($factor_id) {
-    print "Found existing blank factor, id = $factor_id\n";
-  } else {
-     $db_adaptor->dbc->do("INSERT INTO regulatory_factor (name) VALUES ('')");
-     $sth->execute();
-     ($factor_id) = $sth->fetchrow_array();
-     print "Created new blank factor, id = $factor_id\n";
-  }
+#  if ($factor_id) {
+#    print "Found existing blank factor, id = $factor_id\n";
+#  } else {
+#     $db_adaptor->dbc->do("INSERT INTO regulatory_factor (name) VALUES ('')");
+#     $sth->execute();
+#     ($factor_id) = $sth->fetchrow_array();
+#     print "Created new blank factor, id = $factor_id\n";
+#  }
 
-  return $factor_id;
+#  return $factor_id;
 
-}
+#}
 
 
 1;

misc-scripts/regulatory_regions/RegulatoryFeatureParser/cisred.pm

@@ -10,13 +10,21 @@ use File::Basename;
 #
 # http://www.cisred.org/content/databases_methods/human_2/data_files/search_regions.txt
 
+
+#No longer valid urls, now use the following for ensembl formats for all species:
+#http://www.bcgsc.ca/downloads/cisred/temp/cisRED4Ensembl/
+#naminf may not be obvious, may have to cross reference with above previous urls to get build info
+
 # Format of motifs.txt (note group_name often blank)
 
-# name	chromosome	start	end	strand	group_name
-# craHsap1	X	138337029	138337034	1
-# craHsap2	X	138338145	138338150	1
-# craHsap3	X	138338363	138338368	1
-# craHsap4	X	138338388	138338395	1
+#name    chromosome      start   end     strand  group_name      ensembl_gene
+#craHsap1        1       168129978       168129997       -1      1       ENSG00000000457
+#craHsap2        1       168129772       168129781       -1      2       ENSG00000000457
+#craHsap3        1       168129745       168129756       -1      3       ENSG00000000457
+#craHsap4        1       168129746       168129753       -1      4       ENSG00000000457
+#craHsap5        1       168129745       168129752       -1      5       ENSG00000000457
+#craHsap6        1       168129741       168129757       -1      6       ENSG00000000457
+
 
 # Format of search_regions.txt
 # name	chromosome	start	end	strand	ensembl_gene_id
@@ -26,284 +34,301 @@ use File::Basename;
 # 19	1	23602820	23605631	-1	ENSG00000007968
 
 
-use Bio::EnsEMBL::DBSQL::DBAdaptor;
+use RegulatoryFeatureParser::BaseParser;
+use Bio::EnsEMBL::DBEntry;
+use Bio::EnsEMBL::Funcgen::ExternalFeature;
 
 use vars qw(@ISA);
 @ISA = qw(RegulatoryFeatureParser::BaseParser);
 
-# Parse file and return hashref containing:
-#
-# - arrayref of features
-# - arrayref of factors
 
-sub parse {
 
-  my ($self, $db_adaptor, $file, $old_assembly, $new_assembly) = @_;
 
-  my %result;
 
-  my $analysis_adaptor = $db_adaptor->get_AnalysisAdaptor();
-  my $slice_adaptor = $db_adaptor->get_SliceAdaptor();
-
-  my $stable_id_to_internal_id = $self->build_stable_id_cache($db_adaptor);
-
-  print "Parsing $file with cisred parser\n";
-
-  # ----------------------------------------
-  # We need a "blank" factor for those features which aren't assigned factors
-  # Done this way to maintain referential integrity
-
-  my $blank_factor_id = $self->get_blank_factor_id($db_adaptor);
-
-  # ----------------------------------------
-
-  my $feature_internal_id = ($self->find_max_id($db_adaptor, "regulatory_feature")) + 1;
-  my $highest_factor_id = ($self->find_max_id($db_adaptor, "regulatory_factor")) + 1;
-
-  my @features;
-  my @factors;
-  my %factor_ids_by_name; # name -> factor_id
-  my %feature_objects;
+sub new {
+  my $caller = shift;
+  my $class = ref($caller) || $caller;
+
+  my $self = $class->SUPER::new(@_);
+
+  #Set default feature_type and feature_set config
+  $self->{'feature_types'} = {
+							   'cisRED Search Region'   => {
+															class       => 'Region',
+															description => 'cisRED search region',
+												   },
+							   'cisRED Motif' => {
+												  class       => 'Regulatory Motif',
+												  description => 'cisRED group motif',
+												 },
+							 };
+  
+  $self->{feature_sets} = {
+						   'cisRED search regions' => {
+													   feature_type      => \$self->{'feature_types'}{'cisRED Search Region'},
+													   analysis          => 
+													   { 
+														-logic_name    => 'cisRED',
+														-description   => 'cisRED motif search (www.cisred.org)',
+														-display_label => 'cisRED',
+														-displayable   => 1,
+													   },
+													  },
+						   'cisRED group motifs' => {
+													 feature_type      => \$self->{'feature_types'}{'cisRED Motif'},
+													 analysis          => 
+													 { 
+													  -logic_name    => 'cisRED',
+													  -description   => 'cisRED motif search (www.cisred.org)',
+													  -display_label => 'cisRED',
+													  -displayable   => 1,
+													 },
+													},
+						  };
+
+ 
+ 
+  $self->validate_and_store_feature_types;
+  $self->set_feature_sets;
 
-  # ----------------------------------------
-  # Analysis - need one for each type of feature
-  my %analysis;
+  return $self;
+}
 
-  foreach my $anal ("cisRed", "cisred_search") {  # TODO - add other types as necessary
 
-    my $analysis_obj = $analysis_adaptor->fetch_by_logic_name($anal);
 
-    die "Can't get analysis for $anal, skipping" if (!$analysis_obj);
 
-    $analysis{$anal} = $analysis_obj->dbID();
-    print "Analysis ID for $anal is " . $analysis{$anal} . "\n";
 
-  }
+# Parse file and return hashref containing:
+#
+# - arrayref of features
+# - arrayref of factors
 
+sub parse_and_load {
+  my ($self, $file, $old_assembly, $new_assembly) = @_;
+  print ":: Parsing $file with cisRED parser\n";
+
+  my $analysis_adaptor = $self->db->get_AnalysisAdaptor();
+  my %features_by_group; # name -> factor_id
+  my %slice_cache;
+  my $extf_adaptor  = $self->db->get_ExternalFeatureAdaptor;
+  my $dbentry_adaptor = $self->db->get_DBEntryAdaptor;
+  my $ftype_adaptor = $self->db->get_FeatureTypeAdaptor;
+  my $display_name_cache = $self->build_display_name_cache('gene');
   # this object is only used for projection
   my $dummy_analysis = new Bio::EnsEMBL::Analysis(-logic_name => 'CisRedProjection');
 
   # ----------------------------------------
-  # Parse motifs.txt file
-
-  print "Parsing features from $file\n";
+  # We need a "blank" factor for those features which aren't assigned factors
+  # Done this way to maintain referential integrity
+  #my $blank_factor_id = $self->get_blank_factor_id($db_adaptor);
 
+  # Parse motifs.txt file
+  print ":: Parsing cisRED motifs from $file\n";
   my $skipped = 0;
-
-  my $coords_changed = 0;
+  my $skipped_xref = 0;
+  #my $coords_changed = 0;
+  my $cnt = 0;
 
   open (FILE, "<$file") || die "Can't open $file";
   <FILE>; # skip header
-  while (<FILE>) {
 
+  while (<FILE>) {
     next if ($_ =~ /^\s*\#/ || $_ =~ /^\s*$/);
 
-    my %feature;
     # name	chromosome	start	end	strand	group_name   ensembl_gene_id
     my ($motif_name, $chromosome, $start, $end, $strand, $group_name, $gene_id) = split (/\t/);
-    ($gene_id) = $gene_id =~ /(ENS.*G\d{11})/;
-
-    # ----------------------------------------
-    # Feature name & analysis
-
-    $feature{NAME} = $motif_name;
-    $feature{INFLUENCE} = "unknown"; # TODO - what does cisRed store?
-    $feature{ANALYSIS_ID} = $analysis{cisRed};
-
-    # ----------------------------------------
-    # Factor
-
-    # If $group_id is present in %group_sizes we want to create or reuse a factor.
-    # If not, this feature is not associated with any factor.
-    # If the factor is to be created, its name is crtHsapXX where XX is group_id.
-    # TODO - other species prefixes
-
-    $feature{FACTOR_ID} = $blank_factor_id;
-    if ($group_name && $group_name ne '' && $group_name !~ /\s/) {
-      my $factor_id = $factor_ids_by_name{$group_name};
-
-      if (!$factor_id) { # create one
-	my %factor;
-	$factor_id = $highest_factor_id + 1;
-	$factor{INTERNAL_ID} = $factor_id;
-	$factor{NAME} = $group_name;
-	$factor{TYPE} = 'NULL';
-	push @factors, \%factor;
-	$factor_ids_by_name{$factor{NAME}} = $factor_id;
-	$highest_factor_id = $factor_id;
-	#print join("  ", ("Factor: ", $factor{ID}, $factor{NAME}, $factor{TYPE})) . "\n";
-      }
-      $feature{FACTOR_ID} = $factor_id;
-    }
-
-    # ----------------------------------------
-    # Seq_region ID and co-ordinates, projected if necessary
-
-    my $chr_slice = $slice_adaptor->fetch_by_region('chromosome', $chromosome, undef, undef, undef, $old_assembly);
-
-    if (!$chr_slice) {
-      print STDERR "Can't get slice for $chromosome:$start:$end\n";
-      next;
-    }
-
-    my $seq_region_id = $slice_adaptor->get_seq_region_id($chr_slice);
-
-    if (!$seq_region_id) {
-      print STDERR "Can't get seq_region_id for chromosome $chromosome\n";
-      next;
-    }
-
-    $feature{SEQ_REGION_ID} = $seq_region_id;
-
+    #($gene_id) = $gene_id =~ /(ENS.*G\d{11})/;
+
+	if(! exists $slice_cache{$chromosome}){
+	
+	  if($old_assembly){
+		$slice_cache{$chromosome} = $self->slice_adaptor->fetch_by_region('chromosome', 
+																	$chromosome, 
+																	undef, 
+																	undef, 
+																	undef, 
+																	$old_assembly);
+	  }else{
+		$slice_cache{$chromosome} = $self->slice_adaptor->fetch_by_region('chromosome', $chromosome);
+	  }
+
+	  if(! defined 	$slice_cache{$chromosome}){
+		warn "Can't get slice $chromosome for motif $motif_name\n";
+		$skipped++
+		next;
+	  }
+	}
+
+	#get feature_type first
+
+	#we are not maintaining this link in the DB!
+	#Do we need another xref for this or a different table?
+
+	
+	if ($group_name && $group_name ne '' && $group_name !~ /\s/) {
+
+	  if(! exists $features_by_group{$group_name}){
+		$features_by_group{$group_name} = $ftype_adaptor->fetch_by_name('crtHsap'.$group_name);
+
+		if(! defined $features_by_group{$group_name}){
+		  ($features_by_group{$group_name}) = @{$ftype_adaptor->store(Bio::EnsEMBL::Funcgen::FeatureType->new
+																	  (
+																	   -name  => 'crtHsap'.$group_name,
+																	   -class => 'Regulatory Motif',
+																	   -description => 'cisRED group motif',
+																	  ))};
+		}
+	  }
+	}else{
+	  throw('Found cisRED feature $motif_name with no group_name, unable to defined feature_type');
+	}
+
+	my $feature = Bio::EnsEMBL::Funcgen::ExternalFeature->new(
+															  -display_label => $motif_name,
+															  -start         => $start,
+															  -end           => $end,
+															  -strand        => $strand,
+															  -feature_type  => $features_by_group{$group_name}
+															  -feature_set   => $self->{'feature_sets'}{'cisRED group motifs'},
+															  -slice         => $slice_cache{$chromosome},
+															 );
+
+
+	
     # project if necessary
     if ($new_assembly) {
+      $feature = $self->project_feature($feature, $new_assembly);
 
-      #print join("\t", "OLD: ", $start, $end, $strand, $chromosome, $motif_name) . "\n";
-
-      my $projected_feature = $self->project_feature($start, $end, $strand, $chromosome, $chr_slice, $dummy_analysis, $new_assembly, $slice_adaptor, $motif_name);
-
-      $coords_changed++ if ($projected_feature->start() != $start || $projected_feature->end() != $end);
+	  if(! defined $feature){
+		$skipped ++;
+		next;
+	  }
 
-      $start = $projected_feature->start();
-      $end = $projected_feature->end();
-      $strand = $projected_feature->strand();
 
-      $feature{SEQ_REGION_ID} = $slice_adaptor->get_seq_region_id($projected_feature->feature_Slice);
-
-      #print join("\t", "NEW: ", $start, $end, $strand, $chromosome, $motif_name) . "\n";
-
-    }
-
-    $feature{START} = $start;
-    $feature{END} = $end;
-    $feature{STRAND} = $strand;
-
-
-    # ----------------------------------------
-    # Ensembl object - always a gene in cisRed
-
-    my $ensembl_id = $stable_id_to_internal_id->{gene}->{$gene_id};
-
-    if (!$ensembl_id) {
-      print STDERR "Can't get ensembl internal ID for $gene_id, skipping\n";
-      print join("-", $motif_name, $chromosome, $start, $end, $strand, $group_name, $gene_id, "\n");
-      $skipped++;
-      next;
+      #$coords_changed++ if ($projected_feature->start() != $start || $projected_feature->end() != $end);
     }
 
-    $feature{ENSEMBL_TYPE} = "Gene";
-    $feature{ENSEMBL_ID} = $ensembl_id;
-
-    # ----------------------------------------
-    # Feature internal ID
-
-    $feature{INTERNAL_ID} = $feature_internal_id++;
-
-    # ----------------------------------------
-    # Evidence
-
-    $feature{EVIDENCE} = "";
-
-    # ----------------------------------------
-    # Add to object to be returned
-
-    push @features, \%feature;
+	($feature) = @{$extf_adaptor->store($feature)};
+	$cnt++;
+
+	#Build Xref here
+	if($gene_id){
+
+	  #need to take xref core dbname as a parameter
+	  #defaulting to current db at present
+
+	  #should this not have some 'gene' e.g. core_gene?
+
+	  if(! exists $display_name_cache->{$gene_id}){
+		warn "Cannot get ensembl gene id $gene_id for motif $motif_name\n";
+		$skipped_xref++;
+		next;
+	  }
+
+	  my $dbentry = Bio::EnsEMBL::DBEntry->new(
+											   -dbname                 => 'core',
+											   -release                => $self->db->dnadb->dbname,
+											   -status                 => 'KNOWNXREF',
+											   -display_label_linkable => 1,
+											   -db_display_name        => $self->db->dnadb->dbname,
+											   -type                   => 'MISC',
+											   -dbprimary_acc          => $gene_id,
+											   -display_name           => $display_name_cache->{$gene_id},
+											   -info_type              => 'MISC'
+											   #could have version here if we use the correct dnadb to build the cache
+											  );
+	  $dbentry_adaptor->store($dbentry);  
+	}
+  }
+	
+  close FILE;
 
-    #print "Feature: ";
-    #foreach my $field (keys %feature) {
-    #	print $field . ": " . $feature{$field} . " ";
-    #}
-    #print "\n";
+  print ":: Stored $cnt cisRED ExternalFeature motif\n";
+  print ":: Skipped $skipped cisRED ExternalFeature motif imports\n";
+  print ":: Skipped an additional $skipped_xref DBEntry imports\n";
 
-  }
 
-  close FILE;
 
   # ----------------------------------------
   # Search regions 
   # read search_regions.txt from same location as $file
-  my $search_regions_file = dirname($file) . "/search_regions.txt";
 
-  my $skipped_sr = 0;
-
-  my @search_regions;
-  print "Parsing search regions from $search_regions_file\n";
+  my $search_regions_file = dirname($file) . "/search_regions.txt";
+	$skipped = 0;
+	$cnt = 0;
+	
+	print ":: Parsing cisRED search regions from $search_regions_file\n";
   open (SEARCH_REGIONS, "<$search_regions_file") || die "Can't open $search_regions_file";
   <SEARCH_REGIONS>; # skip header
   while (<SEARCH_REGIONS>) {
     chomp;
     my ($id, $chromosome, $start, $end, $strand, $ensembl_gene_id) = split;
-    my $gene_id = $stable_id_to_internal_id->{gene}->{$ensembl_gene_id};
+    my $gene_id = $display_name_cache->{gene}->{$ensembl_gene_id};
+	my $name = "CisRed_Search_$id";
+
     if (!$gene_id) {
       warn("Can't get internal ID for $ensembl_gene_id\n");
-      $skipped_sr++;
-      next;
-    }
-    my $sr_chr_slice = $slice_adaptor->fetch_by_region('chromosome', $chromosome, undef, undef, undef, $old_assembly);
-    if (!$sr_chr_slice) {
-      print STDERR "Can't get slice for $chromosome:$start:$end\n";
-      next;
-    }
-    my $sr_seq_region_id = $slice_adaptor->get_seq_region_id($sr_chr_slice);
-    if (!$sr_seq_region_id) {
-      print STDERR "Can't get seq_region_id for chromosome $chromosome\n";
+      $skipped++;
       next;
     }
 
-    my $name = "CisRed_Search_$id";
+	if(! exists $slice_cache{$chromosome}){
+	  
+	  if($old_assembly){
+		$slice_cache{$chromosome} = $self->slice_adaptor->fetch_by_region('chromosome', 
+																		  $chromosome, 
+																		  undef, 
+																		  undef, 
+																		  undef, 
+																		  $old_assembly);
+	  }else{
+		$slice_cache{$chromosome} = $self->slice_adaptor->fetch_by_region('chromosome', $chromosome);
+	  }
+
+	  if(! defined 	$slice_cache{$chromosome}){
+		warn "Can't get slice $chromosome for for search region $name\n";
+		next;
+	  }
+	}
+
+	my $search_feature = Bio::EnsEMBL::Funcgen::ExternalFeature->new
+	  (
+	   -display_label => $name,
+	   -start         => $start,
+	   -end           => $end,
+	   -strand        => $strand,
+	   -feature_type  => $self->{'feature_sets'}{'cisRED group motifs'}->feature_type,
+	   -feature_set   => $self->{'feature_sets'}{'cisRED group motifs'},
+	   -slice         => $slice_cache{$chromosome},
+	  );
+																
+
 
     # project if necessary
     if ($new_assembly) {
-
-      #print join("\t", "OLD: ", $start, $end, $strand, $chromosome, $name) . "\n";
-
-      my $projected_region = $self->project_feature($start, $end, $strand, $chromosome, $sr_chr_slice, $dummy_analysis, $new_assembly, $slice_adaptor, "CisRed_Search_$id");
-
-      $start = $projected_region->start();
-      $end = $projected_region->end();
-      $strand = $projected_region->strand();
-
-      #print join("\t", "NEW: ", $start, $end, $strand, $chromosome, $name) . "\n";
-
-    }
-
-    my %search_region;
-    $search_region{NAME} = $name;
-    $search_region{SEQ_REGION_ID} = $sr_seq_region_id;
-    $search_region{START} = $start;
-    $search_region{END} = $end;
-    $search_region{STRAND} = ($strand =~ /\+/ ? 1 : -1);
-    $search_region{ENSEMBL_OBJECT_TYPE} = 'Gene';
-    $search_region{ENSEMBL_OBJECT_ID} = $gene_id;
-    $search_region{ANALYSIS_ID} = $analysis{cisred_search};
-    push @search_regions, \%search_region;
-
+      $search_feature = $self->project_feature($search_feature);
+	  
+	  if(! defined $search_feature){
+		$skipped ++;
+		next;
+	  }
+	}
+
+	$extf_adaptor->store($search_feature);
   }
-  close(SEARCH_REGIONS);
-
-
-  # ----------------------------------------
 
-  $result{FEATURES} = \@features;
-  $result{FACTORS} = \@factors;
-  $result{SEARCH_REGIONS} = \@search_regions;
-
-  print "Parsed " . scalar(@{$result{FEATURES}}) . " features, " . scalar(@{$result{FACTORS}}) . " factors and " . scalar(@{$result{SEARCH_REGIONS}}) . " search regions\n";
+  close(SEARCH_REGIONS);
 
-  print "Skipped $skipped features and $skipped_sr search regions due to missing stable ID - internal ID mappings\n";
+  
+  print ":: Stored $cnt cisRED search region ExternalFeatures\n";
+  print ":: Skipped $skipped cisRED search region ExternalFeatures\n";
 
-  print "$coords_changed features had their co-ordinates changed as a result of assembly mapping.\n" if ($new_assembly);
+  #print "$coords_changed features had their co-ordinates changed as a result of assembly mapping.\n" if ($new_assembly);
 
-  return \%result;
+  return;
 
 }
 
-sub new {
-
-  my $self = {};
-  bless $self, "RegulatoryFeatureParser::cisred";
-  return $self;
 
-}
 
 1;

misc-scripts/regulatory_regions/RegulatoryFeatureParser/miranda.pm

@@ -8,7 +8,9 @@ use strict;
 #  Similarity	hsa-miR-23a	miRanda	miRNA_target	1	919787	919807	+	.	71	transcript id "ENST00000310998"
 
 use RegulatoryFeatureParser::BaseParser;
-use Bio::EnsEMBL::DBSQL::DBAdaptor;
+use Bio::EnsEMBL::DBEntry;
+use Bio::EnsEMBL::Funcgen::ExternalFeature;
+
 
 use vars qw(@ISA);
 @ISA = qw(RegulatoryFeatureParser::BaseParser);
@@ -18,44 +20,94 @@ use vars qw(@ISA);
 # - arrayref of features
 # - arrayref of factors
 
-sub parse {
 
-  my ($self, $db_adaptor, $file, $old_assembly, $new_assembly) = @_;
 
-  my %result;
 
-  print "Parsing $file with miranda parser\n";
+sub new {
+  my $caller = shift;
+  my $class = ref($caller) || $caller;
+
+  my $self = $class->SUPER::new(@_);
+
+  #Set default feature_type and feature_set config
+  $self->{'feature_types'} = {
+							  'miRanda'   => {
+														   class       => 'RNA',
+														   description => 'miRanda microRNA',
+														  },
+							 };
+  $self->{feature_sets} = {
+						   'miRanda miRNA' => {
+											   feature_type      => \$self->{'feature_types'}{'cisRED Search Region'},
+											   analysis          => 
+											   { 
+												-logic_name    => 'miRanda',
+												-description   => 'miRanda microRNA target prediction (http://cbio.mskcc.org/mirnaviewer/)',
+												-display_label => 'miRanda',
+												-displayable   => 1,
+											   },
+											  },						   
+						  };
+
+ 
+ 
+  $self->validate_and_store_feature_types;
+  $self->set_feature_sets;
 
-  my $feature_internal_id = ($self->find_max_id($db_adaptor, "regulatory_feature")) + 1;
-  my $highest_factor_id = ($self->find_max_id($db_adaptor, "regulatory_factor")) + 1;
+  return $self;
+}
 
-  my $analysis_adaptor = $db_adaptor->get_AnalysisAdaptor();
-  my $slice_adaptor = $db_adaptor->get_SliceAdaptor();
 
-  my @features;
-  my @factors;
-  my %factor_ids_by_name; # name -> factor_id
-  my %feature_objects;
-  my %anal;
 
-  # TODO - regulatory_factor_coding
 
-  my $stable_id_to_internal_id = $self->build_stable_id_cache($db_adaptor);
+sub parse_and_load{
 
+  my ($self, $file, $old_assembly, $new_assembly) = @_;
+
+  print ":: Parsing miRanda data from:\t$file\n";
+
+  my $analysis_adaptor = $db_adaptor->get_AnalysisAdaptor();
+  my %features_by_name; # name -> feature_type
+  my %slice_cache;
+  my $display_name_cache = $self->build_display_name_cache('gene');
   # this object is only used for projection
-  my $dummy_analysis = new Bio::EnsEMBL::Analysis(-logic_name => 'CisRedProjection');
+  my $dummy_analysis = new Bio::EnsEMBL::Analysis(-logic_name => 'miRandaProjection');
+  my $skipped = 0;
+
 
   open (FILE, "<$file") || die "Can't open $file";
 
   while (<FILE>) {
-
-    next if ($_ =~ /^\s*\#/ || $_ =~ /^\s*$/);
-
-    my %feature;
+	next if ($_ =~ /^\s*\#/ || $_ =~ /^\s*$/);
 
     my ($group, $seq, $method, $feature, $chr, $start, $end, $str, $phase, $score, $pvalue, $type, $id_ignore, $id) = split;
     my $strand = ($str =~ /\+/ ? 1 : -1);
     $id =~ s/[\"\']//g;  # strip quotes
+	$id .= ':'.$seq;
+
+	if(! exists $slice_cache{$chromosome}){
+	
+	  if($old_assembly){
+		$slice_cache{$chr} = $self->slice_adaptor->fetch_by_region('chromosome', 
+																	$chr, 
+																	undef, 
+																	undef, 
+																	undef, 
+																	$old_assembly);
+	  }else{
+		$slice_cache{$chr} = $self->slice_adaptor->fetch_by_region('chromosome', $chr);
+	  }
+
+	  if(! defined 	$slice_cache{$chr}){
+		warn "Can't get slice $chr for sequence $id\n";
+		$skipped++;
+		next;
+	  }
+	}
+
+
+
+
 
     # ----------------------------------------
     # Feature name
@@ -63,12 +115,25 @@ sub parse {
     # For miRNA_target, individual features don't have a unique name, so create
     # a composite one. Also set influence.
 
-    $feature{NAME} = $id .":" . $seq;
+  
     $feature{INFLUENCE} = "negative";
 
     # ----------------------------------------
     # Factor
 
+	 if(! exists $features_by_group{$group_name}){
+		$features_by_group{$group_name} = $ftype_adaptor->fetch_by_name('crtHsap'.$group_name);
+
+		if(! defined $features_by_group{$group_name}){
+		  ($features_by_group{$group_name}) = @{$ftype_adaptor->store(Bio::EnsEMBL::Funcgen::FeatureType->new
+																	  (
+																	   -name  => 'crtHsap'.$group_name,
+																	   -class => 'Regulatory Motif',
+																	   -description => 'cisRED group motif',
+																	  ))};
+		}
+	  }
+
     # $seq is the name of a factor - if it's already there, find its ID, otherwise add it
     my $factor_id = $factor_ids_by_name{$seq};
     if (!$factor_id) {