From cbfc45d33433db4bae7a9c8042f726d965dcad43 Mon Sep 17 00:00:00 2001
From: Ewan Birney <birney@sanger.ac.uk>
Date: Tue, 16 Apr 2002 15:19:15 +0000
Subject: [PATCH] commits for slices for laura - laura - clone.t is now
failing. Unclear if this is me or you
---
modules/Bio/EnsEMBL/BaseAlignFeature.pm | 315 +-----------------
.../EnsEMBL/DBSQL/AssemblyMapperAdaptor.pm | 37 +-
modules/Bio/EnsEMBL/DBSQL/DBAdaptor.pm | 88 +----
.../EnsEMBL/DBSQL/DnaAlignFeatureAdaptor.pm | 79 +++--
modules/Bio/EnsEMBL/DBSQL/SliceAdaptor.pm | 18 +-
modules/Bio/EnsEMBL/FeaturePair.pm | 5 +
modules/Bio/EnsEMBL/Mapper.pm | 8 +-
modules/Bio/EnsEMBL/Slice.pm | 50 ++-
modules/t/dnaalignfeature.t | 8 +-
modules/t/slice.t | 110 ++++++
10 files changed, 274 insertions(+), 444 deletions(-)
create mode 100644 modules/t/slice.t
diff --git a/modules/Bio/EnsEMBL/BaseAlignFeature.pm b/modules/Bio/EnsEMBL/BaseAlignFeature.pm
index 994e8b9e7e..c4bcc19e2d 100644
--- a/modules/Bio/EnsEMBL/BaseAlignFeature.pm
+++ b/modules/Bio/EnsEMBL/BaseAlignFeature.pm
@@ -174,319 +174,6 @@ sub cigar_string {
}
-
-=head2 _generic_parse_cigar
-
- Arg 1 : int query_unit,
- either 1 (peptide) or 3 (dna)
-
- Arg 2 : int hit_unit,
- either 1 (peptide) or 3 (dna)
-
- a DNA,DNA alignment should use 1,1 for the query_unit and hit_unit
- Function : Converts the cigar_string contained by the module into
- an array of ungapped Bio::EnsEMBL::FeaturePair.
-
- See sub cigar_string for an explanation of what that is.
-
- Returntype: list of Bio::EnsEMBL::FeaturePairs
- Exceptions: Mainly internal, eg, no cigar line, wrong arguments or
- miswritten query_unit or hit_unit
-
- Caller : _parse_cigar in derived class
-
-=cut
-
-
-sub _generic_parse_cigar2 {
- my ( $self, $query_unit, $hit_unit ) = @_;
-
- my $string = $self->cigar_string;
-
- if (!defined($string)) {
- $self->throw("No cigar string defined in object. This should be caught by the cigar_string method and never happen");
- }
-
- my @pieces = split(/:/,$string);
-
- my @features;
-
- foreach my $piece (@pieces) {
-
- my $start1;
- my $start2;
- my $length;
-
- if ($piece =~ /(\S+)\,(\S+)\,(\S+)/) {
- $start1 = $1;
- $start2 = $2;
- $length = $3;
- } else {
- $self->throw("Can't parse cigar string element [$piece]. Invalid format. Should be start,end,length*strand");
- }
-
- my $strand = 1;
-
- if ($length == 0) {
- $self->throw("Length of piece is 0 - impossible");
- }
- if ($length < 0) {
- $strand = -1;
- }
-
- $length = abs($length);
-
- my $feature1 = new Bio::EnsEMBL::SeqFeature();
-
- $feature1->start($start1);
- $feature1->end ($start1 + $length - 1);
- $feature1->strand($strand);
- $feature1->score($self->score);
- $feature1->seqname($self->seqname);
- $feature1->phase($self->phase);
- $feature1->p_value($self->p_value);
- $feature1->percent_id($self->percent_id);
-
- my $feature2 = new Bio::EnsEMBL::SeqFeature();
-
- my $mapped_length;
-
- # explicit if statements to avoid rounding problems
- # and make sure we have sane coordinate systems
- if( $query_unit == 1 && $hit_unit == 3 ) {
- $mapped_length = $length*3;
- } elsif( $query_unit == 3 && $hit_unit == 1 ) {
- $mapped_length = $length / 3;
- } elsif ( $query_unit == 1 && $hit_unit == 1 ) {
- $mapped_length = $length;
- } else {
- $self->throw("Internal error $query_unit $hit_unit, currently only allowing 1 or 3 ");
- }
-
- if( int($mapped_length) != $mapped_length ) {
- $self->throw("Internal error with mismapped length of hit, query $query_unit, hit $hit_unit, length $length");
- }
-
-
- if ($strand == 1) {
- $feature2->start($start2);
- $feature2->end ($start2 + $strand*($mapped_length - 1));
- } else {
- $feature2->end ($start2);
- $feature2->start($start2 + $strand*($mapped_length - 1));
- }
- $feature2->strand($strand);
- $feature2->score($self->score);
- $feature2->seqname($self->hseqname);
- $feature2->phase($self->phase);
- $feature2->p_value($self->p_value);
- $feature2->percent_id($self->percent_id);
-
- my $fp = new Bio::EnsEMBL::FeaturePair(-feature1 => $feature1,
- -feature2 => $feature2);
-
-
- $fp->analysis($self->analysis);
-
- push(@features,$fp);
- }
- return @features;
-}
-
-
-=head2 _generic_parse_features
-
- Arg 1 : listref Bio::EnsEMBL::FeaturePair $features
-
- Arg 2 : int query_unit,
- either 1 (peptide) or 3 (dna)
-
- Arg 3 : int hit_unit,
- either 1 (peptide) or 3 (dna)
-
- a DNA,DNA alignment should use 1,1 for the query_unit and hit_unit
-
- Usage : Internal method - not used.
-
- Function : Converts an array of FeaturePairs into a gapped feature with
- a cigar string describing the
-
- See sub cigar_string for an explanation of what that is.
-
- Exception: If the argument passed is not an array reference
-
- All the features must have arisen from the same source
- i.e. a blast HSP or some other alignment. Thus
- exceptions are thrown when the scores,percent ids,p_values
- seqnames , hseqnames and strands differ amongst the input
- features.
-
- All the features must not overlap in order to provide a
- sane gapped alignment. An exception is thrown if they do.
-
- If any element of the array is not a Bio::EnsEMBL::FeaturePair
-
- If there are no elements in the array
-
- If the hit length is not exactly 3 times the query length
-
- Caller : Called internally to the module by the constructor
-
-=cut
-
-sub _generic_parse_features2 {
- my ($self,$features, $query_unit, $hit_unit ) = @_;
-
- if (ref($features) ne "ARRAY") {
- $self->throw("features must be an array reference not a [" . ref($features) . "]");
- }
-
- #print ::LOG "Enter new ",ref( $self ), " with ",scalar( @$features ), " features.\n";
-
- my @f = sort {$a->start <=> $b->start} @$features;
-
- if (scalar(@f) == 0) {
- $self->throw("No features in the array to parse");
- }
- my $hstrand = $f[0]->hstrand;
- my $strand = $f[0]->strand;
- my $name = $f[0]->seqname;
- my $hname = $f[0]->hseqname;
- my $score = $f[0]->score;
- my $percent = $f[0]->percent_id;
- my $pvalue = $f[0]->p_value;
- my $analysis = $f[0]->analysis;
- my $phase = $f[0]->phase;
-
- # implicit starnd 1 for peptide sequences
- ( defined $strand ) || ( $strand = 1 );
- ( defined $hstrand ) || ( $hstrand = 1 );
- my $ori = $strand * $hstrand;
-
-
- my $prev1;
- my $prev2;
-
- my $string;
-
-
- my $f1start = $f[0]->start;
- my $f1end = $f[$#f]->end;
-
- my $f2start;
- my $f2end;
-
- if ( $hstrand == 1 ) {
- $f2start = $f[0]->hstart;
- $f2end = $f[$#f]->hend;
- } else {
- $f2end = $f[0]->hend;
- $f2start = $f[$#f]->hstart;
- }
-
- foreach my $f (@f) {
- if (!$f->isa("Bio::EnsEMBL::FeaturePair")) {
- $self->throw("Array element [$f] is not a Bio::EnsEMBL::FeaturePair");
- }
-
- ##print STDERR "Processing " . $f->gffstring . "\n";
- if( defined $f->hstrand() ) {
- if ($f->hstrand != $hstrand) {
- $self->throw("Inconsistent hstrands in feature array");
- }
- }
- if( defined $f->strand() ) {
- if ($f->strand != $strand) {
- $self->throw("Inconsistent strands in feature array");
- }
- }
-
- if ($name ne $f->seqname) {
- $self->throw("Inconsistent names in feature array [$name - " . $f->seqname . "]");
- }
- if ($hname ne $f->hseqname) {
- $self->throw("Inconsistent names in feature array [$hname - " . $f->hseqname . "]");
- }
- if ($score ne $f->score) {
- $self->throw("Inconsisent scores in feature array [$score - " . $f->score . "]");
- }
- if ($percent ne $f->percent_id) {
- $self->throw("Inconsistent pids in feature array [$percent - " . $f->percent_id . "]");
- }
-
- my $start1 = $f->start;
- my $start2 = $f->hstart();
-
- if (defined($prev2)) {
- if ( $ori == 1 ) {
- if ($f->hstart < $prev2) {
- $self->throw("Inconsistent coordinates in feature array " . $f->hstart . " < $prev2");
- }
- } else {
- if ($f->hend > $prev2) {
- $self->throw("Inconsistent coordinates in feature array " . $f->hend . " > $prev2");
- }
- }
-
- }
- my $length = ($f->end - $f->start + 1);
- my $hlength = ($f->hend - $f->hstart + 1);
-
- # using multiplication to avoid rounding errors, hence the
- # switch from query to hit for the ratios
-
- if( $length * $hit_unit != $hlength * $query_unit ) {
- $self->throw( "Feature lengths not comparable Lengths:".$length." ".$hlength." Ratios:".$query_unit." ". $hit_unit );
- }
-
- # now make length have ori sign
- $length *= $ori;
-
- $string = $string . $start1 . "," . $start2 . "," . $length . ":";
-
- if( $ori == 1 ) {
- $prev2 = $f->hend;
- } else {
- $prev2 = $f->hstart;
- }
-
- }
-
- $string =~ s/:$//;
-
- my $feature1 = new Bio::EnsEMBL::SeqFeature();
-
- $feature1->start($f1start);
- $feature1->end ($f1end);
- $feature1->strand($strand);
- $feature1->score($score);
- $feature1->percent_id($percent);
- $feature1->p_value($pvalue);
- $feature1->seqname($name);
- $feature1->percent_id($percent);
- $feature1->p_value($pvalue);
- $feature1->phase($phase);
- $feature1->analysis($analysis);
-
- my $feature2 = new Bio::EnsEMBL::SeqFeature();
-
- $feature2->start($f2start);
- $feature2->end ($f2end);
- $feature2->strand($hstrand);
- $feature2->score($score);
- $feature2->seqname($hname);
- $feature2->percent_id($percent);
- $feature2->p_value($pvalue);
- $feature2->phase($phase);
- $feature2->analysis($analysis);
-
- $self->feature1($feature1);
- $self->feature2($feature2);
- $self->cigar_string($string);
-
- #print ::LOG "Exit with cigar $string.\n";
-}
-
=head2 ungapped_features
Arg : None.
@@ -516,7 +203,7 @@ sub ungapped_features {
}
-=head2 _generic_parse_cigar2
+=head2 _generic_parse_cigar
Arg 1 : int query_unit,
either 1 (peptide) or 3 (dna)
diff --git a/modules/Bio/EnsEMBL/DBSQL/AssemblyMapperAdaptor.pm b/modules/Bio/EnsEMBL/DBSQL/AssemblyMapperAdaptor.pm
index 03e01374c2..fb0f8fcb3c 100644
--- a/modules/Bio/EnsEMBL/DBSQL/AssemblyMapperAdaptor.pm
+++ b/modules/Bio/EnsEMBL/DBSQL/AssemblyMapperAdaptor.pm
@@ -118,6 +118,7 @@ sub register_region{
$self->throw("$assmapper is not a Bio::EnsEMBL::AssemblyMapper")
unless $assmapper->isa("Bio::EnsEMBL::AssemblyMapper");
+
# This is to reassure people that I've got this piece of
# SQL correct ;) EB
#
@@ -132,8 +133,12 @@ sub register_region{
# SCP: Those not clauses are a bit confusing, so I've
# ditched them :-) "NOT a > b" equiv to "a <= b"
+ # No. SCP - not right! The NOT has a bracket'd AND and so
+ # one can't just expand it like this EB. Replaced it and
+ # if *anyone* wants to change this, talk to me or James G first!
- my $sth = $self->prepare(qq{
+
+ my $select = qq{
select
ass.contig_start,
ass.contig_end,
@@ -148,24 +153,28 @@ sub register_region{
where
chr.name = '$chr_name' and
ass.chromosome_id = chr.chromosome_id and
- $end >= ass.chr_start and
- $start <= ass.chr_end and
+ NOT (
+ $start >= ass.chr_end and
+ $end <= ass.chr_start) and
ass.type = '$type'
- });
+ };
+ #print STDERR "Going to select $select\n";
+ my $sth = $self->prepare($select);
+
$sth->execute();
-
+
while( my $arrayref = $sth->fetchrow_arrayref ) {
- my ($contig_start,$contig_end,$contig_id,$contig_ori,$chr_name,$chr_start,$chr_end) = @$arrayref;
- if( $assmapper->_have_registered_contig($contig_id) == 0 ) {
- $assmapper->_register_contig($contig_id);
- $assmapper->_mapper->add_map_coordinates(
- $contig_id, $contig_start, $contig_end, $contig_ori,
- $chr_name, $chr_start, $chr_end
- );
- }
+ my ($contig_start,$contig_end,$contig_id,$contig_ori,$chr_name,$chr_start,$chr_end) = @$arrayref;
+ if( $assmapper->_have_registered_contig($contig_id) == 0 ) {
+ $assmapper->_register_contig($contig_id);
+ $assmapper->_mapper->add_map_coordinates(
+ $contig_id, $contig_start, $contig_end, $contig_ori,
+ $chr_name, $chr_start, $chr_end
+ );
+ }
}
-
+
}
diff --git a/modules/Bio/EnsEMBL/DBSQL/DBAdaptor.pm b/modules/Bio/EnsEMBL/DBSQL/DBAdaptor.pm
index 3494c6e7b2..4b35a1d141 100755
--- a/modules/Bio/EnsEMBL/DBSQL/DBAdaptor.pm
+++ b/modules/Bio/EnsEMBL/DBSQL/DBAdaptor.pm
@@ -2227,82 +2227,6 @@ sub perl_only_contigs{
}
-=head2 delete_Clone
-
- Title : delete_Clone
- Usage : $obj->delete_Clone($clone_id)
- Function: Deletes clone, including contigs, but not its genes
- Example :
- Returns :
- Args :
-
-
-=cut
-
-sub delete_Clone{
- my ($self,$clone_id) = @_;
-
- #$self->warn("Obj->delete_Clone is a deprecated method!
-#Calling Clone->delete instead!");
-
- (ref($clone_id)) && $self->throw ("Passing an object reference instead of a variable\n");
-
- my $clone = new Bio::EnsEMBL::DBSQL::Clone( -id => $clone_id,
- -dbobj => $self );
-
- return $clone->delete();
-}
-
-=head2 cloneid_to_geneid
-
- Title : cloneid_to_geneid
- Usage :
- Function:
- Example :
- Returns :
- Args :
-
-
-=cut
-
-sub cloneid_to_geneid{
- my ($self,$cloneid) = @_;
-
- $self->warn("Obj->cloneid_to_geneid is a deprecated method!
-Calling Clone->get_all_geneid instead!");
-
- (ref($cloneid)) && $self->throw ("Passing an object reference instead of a variable!\n");
-
- my $clone = new Bio::EnsEMBL::DBSQL::Clone( -id => $cloneid,
- -dbobj => $self );
-
- return $clone->get_all_my_geneid();
-}
-
-
-=head2 gene_Obj
-
- Title : gene_Obj
- Usage : my $geneobj = $db->gene_Obj
- Function: Returns the gene object database handle
- Example :
- Returns : Bio::EnsEMBL::DB::Gene_ObjI
- Args :
-
-=cut
-
-sub gene_Obj {
- my( $self ) = @_;
-
- my( $go );
- unless ($go = $self->{'_gene_obj'}) {
- require Bio::EnsEMBL::DBSQL::Gene_Obj;
- $go = Bio::EnsEMBL::DBSQL::Gene_Obj->new($self);
- $self->{'_gene_obj'} = $go;
- }
- return $go;
-}
-
sub get_GeneAdaptor {
my( $self ) = @_;
@@ -2401,6 +2325,18 @@ sub get_RawContigAdaptor {
}
+sub get_SliceAdaptor {
+ my( $self ) = @_;
+
+ my( $rca );
+ unless ($rca = $self->{'_slice_adaptor'}) {
+ require Bio::EnsEMBL::DBSQL::SliceAdaptor;
+ $a = Bio::EnsEMBL::DBSQL::SliceAdaptor->new($self);
+ $self->{'_slice_adaptor'} = $a;
+ }
+ return $a;
+}
+
=head2 get_AnalysisAdaptor
Title : get_AnalysisAdaptor
diff --git a/modules/Bio/EnsEMBL/DBSQL/DnaAlignFeatureAdaptor.pm b/modules/Bio/EnsEMBL/DBSQL/DnaAlignFeatureAdaptor.pm
index b740f997be..f5cd27a4b4 100644
--- a/modules/Bio/EnsEMBL/DBSQL/DnaAlignFeatureAdaptor.pm
+++ b/modules/Bio/EnsEMBL/DBSQL/DnaAlignFeatureAdaptor.pm
@@ -91,22 +91,8 @@ sub fetch_by_dbID{
}
my $contig = $self->db->get_RawContigAdaptor->fetch_by_dbID($contig_id);
-
- my $out = Bio::EnsEMBL::FeatureFactory->new_feature_pair();;
- $out->start($start);
- $out->end($end);
- $out->strand($strand);
- $out->score($score);
- $out->hstart($hstart);
- $out->hend($hend);
- $out->hstrand($hstrand);
- $out->hseqname($hname);
- $out->cigar($cigar);
-
- $out->analysis($self->db->get_AnalysisAdaptor->fetch_by_dbID($analysis_id));
-
- $out->seqname($contig->id);
- $out->attach_seq($contig->seq);
+ my $analysis = $self->db->get_AnalysisAdaptor->fetch_by_dbID($analysis_id);
+ my $out= $self->_new_feature($start,$end,$strand,$score,$hstart,$hend,$hstrand,$hname,$cigar,$analysis,$contig->name,$contig->seq);
return $out;
@@ -148,25 +134,15 @@ sub fetch_by_contig_id{
}
- my $out = Bio::EnsEMBL::FeatureFactory->new_feature_pair();;
- $out->start($start);
- $out->end($end);
- $out->strand($strand);
- $out->score($score);
- $out->hstart($hstart);
- $out->hend($hend);
- $out->hstrand($hstrand);
- $out->hseqname($hname);
- $out->cigar($cigar);
+ my $out= $self->_new_feature($start,$end,$strand,$score,$hstart,$hend,$hstrand,$hname,$cigar,$ana{$analysis_id},$contig->name,$contig->seq);
- $out->analysis($ana{$analysis_id});
- $out->seqname($contig->name);
- $out->attach_seq($contig->seq);
push(@f,$out);
}
-
return @f;
}
+
+
+
sub fetch_by_contig_id_and_logic_name{
my($self, $cid, $logic_name) = @_;
@@ -280,6 +256,13 @@ sub fetch_by_contig_id_and_dbname{
}
+sub fetch_by_Slice_and_score {
+ my ($self,$slice,$score) = @_;
+
+ $self->warn("*** CODE GOING TO BE FIXED - RETURNING 0***");
+ return ();
+}
+
=head2 fetch_by_assembly_location
Title : fetch_by_assembly_location
@@ -299,6 +282,10 @@ sub fetch_by_assembly_location{
$self->throw("Assembly location must be start,end,chr,type");
}
+
+ $self->warn("*** CODE GOING TO BE FIXED - RETURNING 0***");
+ return ();
+
if( $start !~ /^\d/ || $end !~ /^\d/ ) {
$self->throw("start/end must be numbers not $start,$end (have you typed the location in the right way around - start,end,chromosome,type");
}
@@ -805,6 +792,38 @@ sub fetch_featurepair_list_by_assembly_location_nad_logic_name{
}
+
+sub _new_feature {
+ my ($self,$start,$end,$strand,$score,$hstart,$hend,$hstrand,$hseqname,$cigar,$analysis,$seqname,$seq) = @_;
+
+ if( !defined $seq ) {
+ $self->throw("Internal error - wrong number of arguments to new_feature");
+ }
+
+ my $f1 = Bio::EnsEMBL::SeqFeature->new();
+ my $f2 = Bio::EnsEMBL::SeqFeature->new();
+
+ $f1->start($start);
+ $f1->end($end);
+ $f1->strand($strand);
+ $f1->score($score);
+ $f1->seqname($seqname);
+ $f1->attach_seq($seq);
+
+ $f2->start($hstart);
+ $f2->end($hend);
+ $f2->strand($hstrand);
+ $f2->seqname($hseqname);
+
+ $f1->analysis($analysis);
+ $f2->analysis($analysis);
+
+
+ my $out = Bio::EnsEMBL::DnaDnaAlignFeature->new( -cigar_string => $cigar, -feature1 => $f1, -feature2 => $f2);
+
+ return $out;
+}
+
1;
diff --git a/modules/Bio/EnsEMBL/DBSQL/SliceAdaptor.pm b/modules/Bio/EnsEMBL/DBSQL/SliceAdaptor.pm
index a897b0200b..aa919b447c 100644
--- a/modules/Bio/EnsEMBL/DBSQL/SliceAdaptor.pm
+++ b/modules/Bio/EnsEMBL/DBSQL/SliceAdaptor.pm
@@ -47,10 +47,14 @@ package Bio::EnsEMBL::DBSQL::SliceAdaptor;
use vars qw(@ISA);
use strict;
-# Object preamble - inherits from Bio::Root::RootI
+use Bio::EnsEMBL::DBSQL::BaseAdaptor;
+use Bio::EnsEMBL::Slice;
+@ISA = ('Bio::EnsEMBL::DBSQL::BaseAdaptor');
+# new is inherieted from BaseAdaptor
+
=head2 new_slice
Title : new_slice
@@ -64,10 +68,16 @@ use strict;
=cut
sub new_slice{
- my ($chr,$start,$end,$type) = @_;
+ my ($self,$chr,$start,$end,$type) = @_;
- die "Not implemented new slice yet";
+ my $slice = Bio::EnsEMBL::Slice->new( -chr_name => $chr,
+ -chr_start => $start,
+ -chr_end => $end,
+ -assembly_type => $type,
+ -adaptor => $self);
+
+ return $slice;
}
@@ -84,7 +94,7 @@ sub new_slice{
=cut
sub new_web_slice{
- my ($chr,$start,$end,$type) = @_;
+ my ($self,$chr,$start,$end,$type) = @_;
die "Not implemented new slice yet";
diff --git a/modules/Bio/EnsEMBL/FeaturePair.pm b/modules/Bio/EnsEMBL/FeaturePair.pm
index 38d964cc7e..8a4d7a0cd4 100755
--- a/modules/Bio/EnsEMBL/FeaturePair.pm
+++ b/modules/Bio/EnsEMBL/FeaturePair.pm
@@ -108,6 +108,11 @@ sub new {
FEATURE2
)],@args);
+
+ #if( !defined $feature1 || !defined $feature2 ) {
+ # $self->throw("Have not defined either feature1 or feature2. You can define empty feature1 and feature2 if so desired");
+ #}
+
# Store the features in the object
$feature1 && $self->feature1($feature1);
diff --git a/modules/Bio/EnsEMBL/Mapper.pm b/modules/Bio/EnsEMBL/Mapper.pm
index c9c881eaea..d896224a3d 100644
--- a/modules/Bio/EnsEMBL/Mapper.pm
+++ b/modules/Bio/EnsEMBL/Mapper.pm
@@ -268,6 +268,7 @@ sub map_coordinates{
sub add_map_coordinates{
my ($self, $contig_id, $contig_start, $contig_end, $contig_ori, $chr_name, $chr_start, $chr_end) = @_;
+
if( !defined $contig_ori ) {
$self->throw("Need 7 arguments!");
}
@@ -334,10 +335,15 @@ sub add_map_coordinates{
sub list_pairs{
my ($self, $id, $start, $end, $type) = @_;
+
if( !defined $type ) {
$self->throw("Must start,end,id,type as coordinates");
}
+ if( $start > $end ) {
+ $self->throw("Start is greater than end for id $id, start $start, end $end\n");
+ }
+
# perhaps a little paranoid/excessive
if( $self->_is_sorted == 0 ) {
$self->_sort();
@@ -349,7 +355,6 @@ sub list_pairs{
if( $type eq $self->to ) {
$self_func = \&Bio::EnsEMBL::Mapper::Pair::to;
-
if( !defined $self->{'_pair_hash_to'}->{$id} ) {
return ();
}
@@ -365,6 +370,7 @@ sub list_pairs{
my @output;
foreach my $p ( @list ) {
+
if( &$self_func($p)->end < $start ) {
next;
}
diff --git a/modules/Bio/EnsEMBL/Slice.pm b/modules/Bio/EnsEMBL/Slice.pm
index cf06c04be6..5741aa36e8 100644
--- a/modules/Bio/EnsEMBL/Slice.pm
+++ b/modules/Bio/EnsEMBL/Slice.pm
@@ -61,7 +61,19 @@ sub new {
my $self = {};
bless $self,$class;
-
+
+ my ($chr,$start,$end,$type,$adaptor) = $self->_rearrange([qw(CHR_NAME CHR_START CHR_END ASSEMBLY_TYPE ADAPTOR)],@args);
+
+ if( !defined $chr || !defined $start || !defined $end || !defined $type ) {
+ $self->throw("Do not have all the parameters for slice");
+ }
+
+ $self->chr_name($chr);
+ $self->chr_start($start);
+ $self->chr_end($end);
+ $self->assembly_type($type);
+ $self->adaptor($adaptor);
+
# set stuff in self from @args
return $self;
}
@@ -88,9 +100,17 @@ are the methods to implement
=cut
sub get_all_SimilarityFeatures_above_score{
- my ($self,@args) = @_;
+ my ($self,$score) = @_;
- $self->throw("Ewan has not implemented this function! Complain!!!!");
+ my @out;
+
+ if( !defined $score ) {
+ $self->throw("No defined score.");
+ }
+
+ push(@out,$self->adaptor->db->get_DnaAlignFeatureAdaptor->fetch_by_Slice_and_score($self,$score));
+
+ return @out;
}
@@ -392,5 +412,29 @@ sub assembly_type{
}
+
+=head2 adaptor
+
+ Title : adaptor
+ Usage : $obj->adaptor($newval)
+ Function:
+ Example :
+ Returns : value of adaptor
+ Args : newvalue (optional)
+
+
+=cut
+
+sub adaptor{
+ my ($self,$value) = @_;
+ if( defined $value) {
+ $self->{'adaptor'} = $value;
+ }
+ return $self->{'adaptor'};
+
+}
+
+
+
1;
diff --git a/modules/t/dnaalignfeature.t b/modules/t/dnaalignfeature.t
index e99d356d1d..1be6fe0c0b 100644
--- a/modules/t/dnaalignfeature.t
+++ b/modules/t/dnaalignfeature.t
@@ -4,7 +4,7 @@
BEGIN { $| = 1;
use Test;
- plan tests => 7;
+ plan tests => 8;
}
my $loaded = 0;
@@ -95,5 +95,9 @@ ok($outf);
ok($outf->start == 5);
ok($outf->end == 14);
-#ok( scalar($outf->ungapped_features) == 2);
+ok( scalar($outf->ungapped_features) == 2);
+
+($outf) = $dna_f_ad->fetch_by_assembly_location(1,400,1,'NCBI_28');
+#ok($outf);
+#ok( scalar($outf->ungapped_features) == 2);
diff --git a/modules/t/slice.t b/modules/t/slice.t
new file mode 100644
index 0000000000..21066492c3
--- /dev/null
+++ b/modules/t/slice.t
@@ -0,0 +1,110 @@
+
+ use lib 't';
+
+BEGIN { $| = 1;
+ use Test;
+ plan tests => 4;
+}
+
+my $loaded = 0;
+END {print "not ok 1\n" unless $loaded;}
+
+use EnsTestDB;
+use Bio::EnsEMBL::DBLoader;
+
+
+$loaded = 1;
+
+ok(1);
+
+# Database will be dropped when this
+# object goes out of scope
+my $ens_test = EnsTestDB->new;
+
+$ens_test->do_sql_file("t/minidatabase.dump");
+
+ok($ens_test);
+
+
+my $db = $ens_test->get_DBSQL_Obj;
+
+$sla= $db->get_SliceAdaptor();
+
+$slice = $sla->new_slice('1',4,400,'NCBI_28');
+
+ok($slice);
+
+&write_feature();
+
+ok(1);
+
+
+($outf) = $slice->get_all_SimilarityFeatures_above_score(5);
+
+#ok($outf);
+
+
+
+sub write_feature {
+
+$dna_f_ad = $db->get_DnaAlignFeatureAdaptor();
+
+
+
+
+$feature1 = new Bio::EnsEMBL::SeqFeature();
+$feature1->start(5);
+$feature1->end (7);
+$feature1->strand(1);
+$feature1->score(10);
+$feature1->seqname(1);
+#$feature1->analysis($self->analysis);
+
+$feature2 = new Bio::EnsEMBL::SeqFeature();
+$feature2->start (105);
+$feature2->end (107);
+$feature2->strand (1);
+$feature2->score (10);
+$feature2->seqname("dummy-hid");
+
+$fp = new Bio::EnsEMBL::FeaturePair(-feature1 => $feature1,
+ -feature2 => $feature2);
+
+push(@feats,$fp);
+
+
+$feature1 = new Bio::EnsEMBL::SeqFeature();
+$feature1->start(10);
+$feature1->end (14);
+$feature1->strand(1);
+$feature1->score(10);
+$feature1->seqname(1);
+
+#$feature1->analysis($self->analysis);
+
+$feature2 = new Bio::EnsEMBL::SeqFeature();
+$feature2->start (106);
+$feature2->end (110);
+$feature2->strand (1);
+$feature2->score (10);
+$feature2->seqname('dummy-hid');
+
+$fp2 = new Bio::EnsEMBL::FeaturePair(-feature1 => $feature1,
+ -feature2 => $feature2);
+
+
+push(@feats,$fp2);
+
+
+$dnaf = Bio::EnsEMBL::DnaDnaAlignFeature->new( -features => \@feats );
+
+
+$dnaf->seqname(1);
+$dnaf->hseqname('dummy-hid');
+
+$dnaf->analysis($db->get_AnalysisAdaptor->fetch_by_logic_name("dummy-blast"));
+
+
+$dna_f_ad->store(1,$dnaf);
+
+}
--
GitLab