FastaSequence.pm now also overrides Slice::subseq

Ensembl · Aug 7, 2017 · d896527 · d896527
1 parent 44453ce
commit d896527
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 4 deletions.
diff --git a/modules/Bio/EnsEMBL/Variation/Utils/FastaSequence.pm b/modules/Bio/EnsEMBL/Variation/Utils/FastaSequence.pm
@@ -214,17 +214,23 @@ sub setup_fasta {
   # dont re-redefine
   unless($Bio::EnsEMBL::Slice::_fasta_redefined) {
 
-    # "back up" the original seq method
+    # "back up" the original seq and subseq method
     # that way we can fall back to it later
     no warnings 'redefine';
     *Bio::EnsEMBL::Slice::_fasta_old_db_seq = \&Bio::EnsEMBL::Slice::seq unless $offline;
 
-    # redefine Slice's seq() method
+    no warnings 'redefine';
+    *Bio::EnsEMBL::Slice::_fasta_old_db_subseq = \&Bio::EnsEMBL::Slice::subseq unless $offline;
+
+    # redefine Slice's seq() and subseq() methods
     # to a new method that uses the FASTA sequence
     # and some nice caching
     no warnings 'redefine';
     *Bio::EnsEMBL::Slice::seq = _new_slice_seq();
 
+    no warnings 'redefine';
+    *Bio::EnsEMBL::Slice::subseq = _new_slice_seq();
+
     # this method does the actual fetching
     # separate so it can easily use Bio::DB::Fasta or Faidx
     *Bio::EnsEMBL::Slice::_raw_seq = _raw_seq();
@@ -262,6 +268,9 @@ sub revert_fasta {
   no warnings 'redefine';
   *Bio::EnsEMBL::Slice::seq = \&Bio::EnsEMBL::Slice::_fasta_old_db_seq;
 
+  no warnings 'redefine';
+  *Bio::EnsEMBL::Slice::subseq = \&Bio::EnsEMBL::Slice::_fasta_old_db_subseq;
+
   undef $Bio::EnsEMBL::Slice::fasta_db;
   undef $Bio::EnsEMBL::Slice::_fasta_PARs;
   undef $Bio::EnsEMBL::Slice::_fasta_synonyms;
@@ -349,10 +358,13 @@ sub _get_fasta_db {
 # sequence is cached in 1MB range around the current slice
 sub _new_slice_seq {
   return sub {
-    my $self = shift;
+    my ($self, $start, $end, $strand) = @_;
     my ($seq, $length) = ('', 0);
 
-    my ($sr_name, $start, $end, $strand) = ($self->seq_region_name, $self->start, $self->end, $self->strand);
+    $start = $start ? ($self->start + $start) - 1 : $self->start;
+    $end   = $end   ? ($self->start + $end) - 1   : $self->end;
+    $strand = defined($strand) ? $strand * $self->strand : $self->strand;
+    my $sr_name = $self->seq_region_name;
 
     # indels
     return "" if $start > $end;

diff --git a/modules/t/fastaSequence.t b/modules/t/fastaSequence.t
@@ -88,6 +88,17 @@ clear_fasta_cache();
 is($slice1->expand(5, 5)->invert->seq, $seq1, "rev - expand single bp slice 5' and 3' and invert after cache clear");
 
 
+# test subseq
+clear_fasta_cache();
+$slice1 = $sa->fetch_by_region('chromosome', 21, 25606450, 25606460);
+is($slice1->subseq, 'CATGGCACAAC', 'subseq - no params');
+is($slice1->subseq(1, 5), 'CATGG', 'subseq 1');
+is($slice1->subseq(3, 5), 'TGG', 'subseq 2');
+is($slice1->subseq(3, 5, -1), 'CCA', 'subseq rev');
+is($slice1->subseq(5, 4), '', 'subseq e > s');
+is($slice1->subseq(-1, 5), 'CGCATGG', "subseq overlap 5'");
+
+
 # test going off ends
 $slice1 = $sa->fetch_by_region('chromosome', 21, 1, 10);
 is($slice1->seq, 'N' x 10, "start of chrom");