Commit 0ae61e75 authored by Ian Longden's avatar Ian Longden
Browse files

mapping fixes to use seq_region_id instead pf names

parent 79a1ae96
......@@ -213,7 +213,6 @@ sub map {
push @tmp, $frm_seq_region_name;
my $seq_region_id = @{$self->adaptor()->seq_regions_to_ids($frm_cs, \@tmp)}[0];
print STDERR "name = $frm_seq_region_name, id = $seq_region_id \n";
#speed critical section:
#try to do simple pointer equality comparisons of the coord system objects
......@@ -241,8 +240,6 @@ sub map {
" of this AssemblyMapper");
}
print STDERR "seq region id is $seq_region_id\n";
print STDERR "mapper is of type ".ref($mapper)."\n";
return $mapper->map_coordinates($seq_region_id, $frm_start, $frm_end,
$frm_strand, $frm);
}
......@@ -482,6 +479,7 @@ sub list_seq_regions {
my($self, $frm_seq_region, $frm_start, $frm_end, $frm_cs) = @_;
#retrieve the seq_region names
my @seq_ids =
$self->list_ids($frm_seq_region,$frm_start,$frm_end,$frm_cs);
......@@ -495,7 +493,7 @@ sub list_seq_regions {
}
#convert them to ids
return @{$self->adaptor()->seq_ids_to_regions($to_cs, \@seq_ids)};
return @{$self->adaptor()->seq_ids_to_regions(\@seq_ids)};
}
#sub list_ids {
......
......@@ -303,9 +303,6 @@ sub map {
#get a list of ranges in the requested region that have not been registered,
#and register them at the same
#print STDERR "frm_start=$frm_start frm_end=$frm_end" .
# "min_start=$min_start min_end=$min_end\n";
my $ranges;
if($is_insert) {
......@@ -332,13 +329,10 @@ sub map {
}
if($fastmap) {
print STDERR "FAST seq region id is $seq_region_id using $frm_seq_region_name\n";
return $mapper->fastmap($frm_seq_region_name, $frm_start, $frm_end,
return $mapper->fastmap($seq_region_id, $frm_start, $frm_end,
$frm_strand, $frm);
}
print STDERR "mapper is of type ".ref($mapper)."\n";
print STDERR "SLOW seq region id is $seq_region_id\n";
return $mapper->map_coordinates($seq_region_id, $frm_start, $frm_end,
$frm_strand, $frm);
}
......
......@@ -195,10 +195,6 @@ sub fetch_all_by_Slice_constraint {
# same seq_region as original slice
my $sr_id = $slice->get_seq_region_id();
print STDERR "sr_id = $sr_id\n";
foreach my $pro (@proj){
print STDERR "projected slice -> ".$pro->to_Slice->get_seq_region_id()."\n";
}
@proj = grep { $_->to_Slice->get_seq_region_id() != $sr_id } @proj;
......@@ -221,8 +217,6 @@ sub fetch_all_by_Slice_constraint {
foreach my $seg (@proj) {
my $offset = $seg->from_start();
my $seg_slice = $seg->to_Slice();
print STDERR "MORE: ".$seg_slice->name."\n";
print STDERR "constraint :- $constraint\n";
my $features = $self->_slice_fetch($seg_slice, $constraint); ## NO RESULTS
# if this was a symlinked slice offset the feature coordinates as needed
......@@ -230,7 +224,6 @@ sub fetch_all_by_Slice_constraint {
FEATURE:
foreach my $f (@$features) {
print STDERR "FEAT: ".ref($f)."\n";
if($offset != 1) {
$f->{'start'} += $offset-1;
......@@ -308,7 +301,7 @@ sub _slice_fetch {
my $asma = $self->db->get_AssemblyMapperAdaptor();
my @features;
warn "jere are $self, @feat_css\n";
# fetch the features from each coordinate system they are stored in
COORD_SYSTEM: foreach my $feat_cs (@feat_css) {
my $mapper;
......@@ -342,16 +335,13 @@ warn "jere are $self, @feat_css\n";
" AND ${tab_syn}.seq_region_start >= $min_start";
}
print STDERR "#constraint $constraint\n";
my $fs = $self->generic_fetch($constraint,undef,$slice);
# features may still have to have coordinates made relative to slice
# start
print STDERR"#number of features is :- ".scalar(@$fs)."\n";
$fs = _remap($fs, $mapper, $slice);
push @features, @$fs;
print STDERR "#number of features is :- ".scalar(@$fs)."\n";
} else {
$mapper = $asma->fetch_by_CoordSystems($slice_cs, $feat_cs);
......@@ -381,14 +371,11 @@ warn "jere are $self, @feat_css\n";
$constraint .= " AND " if($constraint);
$constraint .= "${tab_syn}.seq_region_id IN ($id_str)";
print STDERR ".con: $constraint\n";
my $fs = $self->generic_fetch($constraint, $mapper, $slice);
print STDERR ".number of features is :- ".scalar(@$fs)."\n";
$fs = _remap($fs, $mapper, $slice);
push @features, @$fs;
print STDERR ".number of features is :- ".scalar(@features)."\n";
} else {
# do multiple split queries using start / end constraints
......@@ -411,12 +398,9 @@ warn "jere are $self, @feat_css\n";
" AND ${tab_syn}.seq_region_start >= $min_start";
}
print STDERR "~con: $constraint\n";
my $fs = $self->generic_fetch($constraint,$mapper,$slice);
print STDERR "~number of features is :- ".scalar(@$fs)."\n";
$fs = _remap($fs, $mapper, $slice);
print STDERR "~number after _remap features is :- ".scalar(@$fs)."\n";
push @features, @$fs;
}
......
......@@ -180,11 +180,8 @@ sub fetch_all_by_Transcript {
" AND e.exon_id = et.exon_id";
# fetch all of the exons
print STDERR "\n".$slice."\n\t".$slice->name."\n";
print STDERR "$constraint\n";
my $exons = $self->fetch_all_by_Slice_constraint($slice,$constraint);
print STDERR "ther are ".scalar(@$exons)." exons\n";
#un-override the table definition
$self->{'tables'} = undef;
$self->{'final_clause'} = undef;
......@@ -193,7 +190,6 @@ sub fetch_all_by_Transcript {
if($slice->name() ne $tslice->name()) {
my @out;
foreach my $ex (@$exons) {
print STDERR "$ex being transfered\n";
push @out, $ex->transfer($tslice);
}
$exons = \@out;
......
......@@ -251,7 +251,6 @@ sub _objs_from_sth {
$mapper->fastmap($sr_name, $seq_region_start, $seq_region_end,
$seq_region_strand, $sr_cs);
print STDERR "SR_NAME".$seq_region_id."\n";
#skip features that map to gaps or coord system boundaries
next FEATURE if(!defined($seq_region_id));
......@@ -302,7 +301,6 @@ sub _objs_from_sth {
'display_label' => $display_label,
'score' => $score});
}
print STDERR "FEATURE COUNT $count\n";
return \@features;
}
......
......@@ -84,56 +84,56 @@ sub sql {
## call Bio::EnsEMBL::DBSQL::StatementHandle->dump(0); to end log
## and set $dump to 0
## leave $dump = 1 for continuous log
#my @bind_args=();
#my $dump = 1;
#my $total_time = 0;
my @bind_args=();
my $dump = 0;
my $total_time = 0;
sub dump {
$dump = shift;
if($total_time){
print STDERR "###################Time taken was $total_time\n";
}
$total_time = 0;
}
#sub dump {
# $dump = shift;
# if($total_time){
# print "Time taken was $total_time\n";
# }
# $total_time = 0;
#}
sub bind_param{
my ($self,@args) = @_;
#sub bind_param{
# my ($self,@args) = @_;
$bind_args[$args[0]-1] = $args[1];
$self->SUPER::bind_param(@args);
}
# $bind_args[$args[0]-1] = $args[1];
# $self->SUPER::bind_param(@args);
#}
sub execute {
my ($self,@args) = @_;
#sub execute {
# my ($self,@args) = @_;
if(!$dump){ # skip dumping
return $self->SUPER::execute(@args);
}
my $sql = $self->sql();
# if(!$dump){ # skip dumping
# return $self->SUPER::execute(@args);
# }
# my $sql = $self->sql();
my @chrs = split(//, $sql);
# my @chrs = split(//, $sql);
my $j = 0;
for(my $i =0; $i < @chrs; $i++) {
$chrs[$i] = $bind_args[$j++] if($chrs[$i] eq '?' && defined($bind_args[$j]));
}
# my $j = 0;
# for(my $i =0; $i < @chrs; $i++) {
# $chrs[$i] = $bind_args[$j++] if($chrs[$i] eq '?' && defined($bind_args[$j]));
# }
my $str = join('', @chrs);
# my $str = join('', @chrs);
print STDERR "\nSQL:\n$str\n\n";
# print "\nSQL:\n$str\n\n";
my $time = time;
my $res = $self->SUPER::execute(@args);
$time = time - $time;
# my $time = time;
# my $res = $self->SUPER::execute(@args);
# $time = time - $time;
$total_time += $time;
print "DONE ($time)\n";
return $res;
}
# $total_time += $time;
# print "DONE ($time)\n";
# return $res;
#}
#
# End uncomment
#
sub DESTROY {
my ($obj) = @_;
......
......@@ -1196,7 +1196,6 @@ sub _objs_from_sth {
#
if($dest_mapper) {
print STDERR "NAME: $sr_name\n";
($seq_region_id,$seq_region_start,$seq_region_end,$seq_region_strand) =
$dest_mapper->fastmap($sr_name, $seq_region_start, $seq_region_end,
$seq_region_strand, $sr_cs);
......@@ -1206,13 +1205,9 @@ sub _objs_from_sth {
#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))) {
print STDERR "---$cmp_cs_name, $seq_region_id,undef, undef, undef,
$cmp_cs_vers\n";
$slice = $slice_hash{"ID:".$seq_region_id} ||=
$sa->fetch_by_seq_region_id($seq_region_id);
} else {
print STDERR "----$asm_cs_name, $seq_region_id, undef, undef, undef,
$asm_cs_vers\n";
$slice = $slice_hash{"ID:".$seq_region_id} ||=
$sa->fetch_by_seq_region_id($seq_region_id);
}
......
......@@ -110,7 +110,6 @@ sub fetch_by_Transcript {
# but I guess thats ok ....
foreach my $exon (@{$transcript->get_all_Exons()}) {
print STDERR "exon id $exon";
if($exon->dbID() == $start_exon_id ) {
$start_exon = $exon;
}
......
......@@ -266,6 +266,7 @@ sub fetch_all_by_Slice {
my $slice_end = $slice->end;
my $slice_strand = $slice->strand;
my $slice_seq_region = $slice->seq_region_name;
my $slice_seq_region_id = $slice->get_seq_region_id;
my $coord_system = $slice->coord_system;
if($self->_supported('SLICE')) {
......@@ -309,8 +310,8 @@ sub fetch_all_by_Slice {
foreach my $segment (@{$slice->project($from_coord_system->name,
$from_coord_system->version)}) {
my ($start,$end,$pslice) = @$segment;
$features{$pslice->seq_region_name} ||= [];
push @{$features{$pslice->seq_region_name}},
$features{$pslice->seq_region_name } ||= [];
push @{$features{$pslice->seq_region_name }},
@{&$fetch_method($self, $pslice->seq_region_name,
$pslice->start(),
$pslice->end())};
......@@ -334,17 +335,17 @@ sub fetch_all_by_Slice {
$slice_setter = _guess_slice_setter($feats) if(!$slice_setter);
foreach my $f (@$feats) {
my($sr_name, $start, $end, $strand) =
my($sr_id, $start, $end, $strand) =
$mapper->fastmap($fseq_region,$f->start,$f->end,$f->strand,
$from_coord_system);
#maps to gap
next if(!defined($sr_name));
next if(!defined($sr_id));
#maps to unexpected seq region, probably error in the externally
if($sr_name ne $slice_seq_region) {
warning("Externally created Feature mapped to [$sr_name] " .
"which is not on requested seq_region [$slice_seq_region]");
if($sr_id ne $slice_seq_region_id) {
warning("Externally created Feature mapped to [$sr_id] " .
"which is not on requested seq_region_id [$slice_seq_region_id]");
next;
}
......
......@@ -429,10 +429,7 @@ sub transform {
}
my $projection = $self->project( $cs_name, $cs_version );
print STDERR "projection is".$projection."\n";
foreach my $segment (@$projection) {
print STDERR "feature project ".$segment->[2]."\t".$segment->[0]."\n";
}
if( @$projection != 1 ) {
return undef;
} else {
......
......@@ -172,7 +172,7 @@ sub _objs_from_sth {
my %marker_cache;
my %slice_hash;
my %sr_name_hash;
# my %sr_name_hash;
my %sr_cs_hash;
my %analysis_cache;
my $marker_adp = $self->db->get_MarkerAdaptor;
......@@ -225,12 +225,12 @@ sub _objs_from_sth {
}
#get the slice object
my $slice = $slice_hash{"ID:".$seq_region_id};
my $slice = $slice_hash{$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();
$slice_hash{$seq_region_id} = $slice;
# $sr_name_hash{$seq_region_id} = $slice->seq_region_name();
$sr_cs_hash{$seq_region_id} = $slice->coord_system();
}
......@@ -246,26 +246,18 @@ sub _objs_from_sth {
# if a mapper was provided
#
if($mapper) {
my $sr_name = $sr_name_hash{$seq_region_id};
# my $sr_name = $sr_name_hash{$seq_region_id};
my $sr_cs = $sr_cs_hash{$seq_region_id};
($sr_name,$seq_region_start,$seq_region_end) =
$mapper->fastmap($sr_name, $seq_region_start, $seq_region_end,
0, $sr_cs);
($seq_region_id,$seq_region_start,$seq_region_end) =
$mapper->fastmap($slice->seq_region_name(), $seq_region_start, $seq_region_end, 0, $sr_cs);
#skip features that map to gaps or coord system boundaries
next FEATURE if(!defined($sr_name));
next FEATURE if(!defined($seq_region_id));
#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);
}
$slice = $slice_hash{"$seq_region_id"} ||=
$sa->fetch_by_seq_region_id($seq_region_id);
}
#
......
......@@ -523,16 +523,9 @@ sub add_map_coordinates{
}
my $test = $contig_id;
$test =~ s/\d+//g;
if(length($test) > 0 or $contig_id < 1000){
print STDERR "$contig_id NOT INTEGER\n";
print STDERR stack_trace_dump();
}
$test = $chr_name;
$test =~ s/\d+//g;
if(length($test) > 0 or $chr_name < 1000 ){
print STDERR "$chr_name NOT INTEGER\n";
print STDERR stack_trace_dump();
}
my $from =
Bio::EnsEMBL::Mapper::Unit->new($contig_id, $contig_start, $contig_end);
......
......@@ -706,7 +706,6 @@ sub project {
my $normal_slice_proj =
$slice_adaptor->fetch_normalized_slice_projection($self);
foreach my $segment (@$normal_slice_proj) {
print STDERR "slice project ".$segment->[2]->name."\t".$segment->[2]->coord_system->name."\n";
my $normal_slice = $segment->[2];
$slice_cs = $normal_slice->coord_system();
......@@ -718,11 +717,11 @@ sub project {
my @coords;
if( defined $asm_mapper ) {
print STDERR "HI ".$normal_slice->seq_region_name()."\n";
print STDERR "\t".$normal_slice->start(),
"\t".$normal_slice->end(),
"\t",$normal_slice->strand(),
"\t",$slice_cs->name." -> ".$cs->name."\n";
# print STDERR "HI ".$normal_slice->seq_region_name()."\n";
# print STDERR "\t".$normal_slice->start(),
# "\t".$normal_slice->end(),
# "\t",$normal_slice->strand`(),
# "\t",$slice_cs->name." -> ".$cs->name."\n";
@coords = $asm_mapper->map($normal_slice->seq_region_name(),
$normal_slice->start(),
$normal_slice->end(),
......@@ -734,16 +733,6 @@ sub project {
}
#construct a projection from the mapping results and return it
if($coords[0]->isa('Bio::EnsEMBL::Mapper::Coordinate')) {
print STDERR "asm_mapper is a ".ref($asm_mapper)."\n";
print STDERR "#COORD id is ".$coords[0]->id."\n";
my $test = $coords[0]->id;
$test =~ s/\d+//g;
if(length($test) > 0 ){
print STDERR "NOT INTEGER\n";
print STDERR stack_trace_dump();
}
}
foreach my $coord (@coords) {
my $coord_start = $coord->start();
my $coord_end = $coord->end();
......@@ -764,7 +753,7 @@ sub project {
}
#create slices for the mapped-to coord system
print STDERR "COORD id is ".$coord->id."\n";
# print STDERR "COORD id is ".$coord->id."\n";
# my $slice = $slice_adaptor->fetch_by_region($coord_cs->name(),
# $coord->id(),
# $coord_start,
......@@ -776,8 +765,8 @@ sub project {
$coord->id(),
$coord_start,
$coord_end,
$coord->strand(),);
print STDERR "slice poject ".$slice."\n";
$coord->strand());
my $current_end = $current_start + $length - 1;
push @projection, bless([$current_start, $current_end, $slice],
......
......@@ -328,10 +328,10 @@ sub _list {
my $seq_region_id = @{$self->adaptor()->seq_regions_to_ids($frm_cs, \@tmp)}[0];
if($seq_regions) {
@result = $mapper->list_seq_regions($seq_region_id, $frm_start,
@result = $mapper->list_seq_regions($frm_seq_region_name, $frm_start,
$frm_end, $frm_cs);
} else {
@result = $mapper->list_ids($seq_region_id, $frm_start,
@result = $mapper->list_ids($frm_seq_region_name, $frm_start,
$frm_end, $frm_cs);
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment