Merge pull request #328 from at7/add_motif_feature_95

Add motif feature 95
Ensembl · Nov 22, 2018 · 2f676c0 · 2f676c0
2 parents 1bddb64 + 73014ae
commit 2f676c0
Show file tree

Hide file tree

Showing 2 changed files with 56 additions and 50 deletions.
diff --git a/modules/Bio/EnsEMBL/VEP/AnnotationSource/Database/RegFeat.pm b/modules/Bio/EnsEMBL/VEP/AnnotationSource/Database/RegFeat.pm
@@ -181,57 +181,53 @@ sub get_features_by_regions_uncached {
 
     $sub_slice->{coord_system}->{adaptor} ||= $self->get_adaptor('core', 'coordsystem');
     $sub_slice->{adaptor} ||= $self->get_adaptor('core', 'slice');
-# Don't dump MotifFeature for release/94
-#    foreach my $type(qw(RegulatoryFeature MotifFeature)) {
-    foreach my $type(qw(RegulatoryFeature)) {
-      my $features = $self->get_adaptor('funcgen', $type)->fetch_all_by_Slice($sub_slice);
-      next unless defined($features);
-
-      foreach my $rf(@$features) {        
-
-        # weaken a circular ref
-        foreach my $a(@{$rf->{_regulatory_activity} || []}) {
-          weaken($a->{_regulatory_feature}) unless isweak($a->{_regulatory_feature});
-        }
-
-        # cell types
-        if($self->{cell_type} && scalar(@{$self->{cell_type}})) {
-
-          my %cl;
-
-          # get cell type using regulatory_activity objects
-          if($type eq 'RegulatoryFeature') {
-            %cl =
-              map {$_->[0] => $_->[1]}
-              map {$_->[0] =~ s/ /\_/g; $_}
-              map {[$_->get_Epigenome->display_label, $_->activity]}
-              @{$rf->regulatory_activity};
-          }
-
-          # get cell type by fetching regfeats that contain this MotifFeature
-          elsif($type eq 'MotifFeature') {
-            %cl =
-              map {$_->[0] => $_->[1]}
-              map {$_->[0] =~ s/ /\_/g; $_}
-              map {[$_->get_Epigenome->display_label, $_->activity]}
-              map {@{$_->regulatory_activity}}
-              @{$rfa->fetch_all_by_attribute_feature($rf)};
-            $rf->binding_matrix->_elements;
-            $rf->binding_matrix->_min_max_sequence_similarity_score;
-          }
-
-          $rf->{cell_types} = \%cl;
-        }
+    my $type = 'RegulatoryFeature'; 
+    my $features = $self->get_adaptor('funcgen', $type)->fetch_all_by_Slice($sub_slice);
+    next unless defined($features);
+
+    foreach my $rf(@$features) { 
+      # weaken a circular ref
+      foreach my $a(@{$rf->{_regulatory_activity} || []}) {
+        weaken($a->{_regulatory_feature}) unless isweak($a->{_regulatory_feature});
       }
 
-      push @region_features,
-        map { $_->{_vep_feature_type} ||= $type; $_ }
-        map { $_->transfer($sr_slice) }
-        @{$features};
+      # cell types
+      if($self->{cell_type} && scalar(@{$self->{cell_type}})) {
+        my %cl =
+          map {$_->[0] => $_->[1]}
+          map {$_->[0] =~ s/ /\_/g; $_}
+          map {[$_->get_Epigenome->display_label, $_->activity]}
+          @{$rf->regulatory_activity};
+        $rf->{cell_types} = \%cl;
+      }
+    }
+    push @region_features,
+      map { $_->{_vep_feature_type} ||= $type; $_ }
+      map { $_->transfer($sr_slice) }
+      @{$features};
+
+    $type = 'MotifFeature';
+    my @motif_features = ();
+    foreach my $rf (@$features) {
+      push @motif_features, @{$rf->fetch_all_MotifFeatures_with_matching_Peak()};
+    } 
+
+    foreach my $mf(@motif_features) {
+      $mf->binding_matrix->summary_as_hash();
+      if($self->{cell_type} && scalar(@{$self->{cell_type}})) {
+        my %cl =  
+          map {$_->[0] => $_->[1]}
+          map {$_->[0] =~ s/ /\_/g; $_}
+          map { [$_->fetch_PeakCalling->fetch_Epigenome->name, 1] } @{$mf->fetch_all_overlapping_Peaks};
+        $mf->{cell_types} = \%cl;
+      }
     }
+    push @region_features,
+      map { $_->{_vep_feature_type} ||= $type; $_ }
+      map { $_->transfer($sr_slice) }
+      @motif_features;
 
     $cache->{$c}->{$region_start} = \@region_features;
-
     push @return, @region_features;
   }
 

diff --git a/modules/Bio/EnsEMBL/VEP/Pipeline/DumpVEP/Dumper/Regulation.pm b/modules/Bio/EnsEMBL/VEP/Pipeline/DumpVEP/Dumper/Regulation.pm
@@ -36,7 +36,6 @@ use Scalar::Util qw(weaken);
 use Bio::EnsEMBL::VEP::Config;
 use Bio::EnsEMBL::VEP::AnnotationSource::Database::RegFeat;
 use Bio::EnsEMBL::VEP::AnnotationSource::Cache::RegFeat;
-
 use base qw(Bio::EnsEMBL::VEP::Pipeline::DumpVEP::Dumper);
 
 sub run {
@@ -131,15 +130,26 @@ sub clean_regfeat {
     set
     _regulatory_activity
     _regulatory_build
+    overlapping_Peaks
   );
 
   if(defined($rf->{binding_matrix})) {
     $rf->{_variation_effect_feature_cache}->{seq} = $rf->seq;
-
-    foreach my $key(qw(adaptor feature_type analysis dbID associated_transcription_factor_complexes)) {
+    foreach my $key(qw(adaptor feature_type analysis dbID)) {
       delete $rf->{binding_matrix}->{$key};
     }
-    $rf->{binding_matrix}->{associated_transcription_factor_complexes} = [];
+    if (defined($rf->{binding_matrix}->{associated_transcription_factor_complexes})) {
+      foreach my $tfc (@{$rf->{binding_matrix}->{associated_transcription_factor_complexes}}) {
+        foreach my $key(qw(adaptor dbID feature_type)) {
+          delete $tfc->{$key};
+        }
+        foreach my $component (@{$tfc->{components}}) {
+          foreach my $key(qw(adaptor dbID feature_type)) {
+            delete $component->{$key};
+          }
+        }
+      }
+    }
   }
 
   $rf->{feature_type} = $rf->{feature_type}->{so_name} if $rf->{feature_type};