Commit 5bd6ce87 authored by Ian Longden's avatar Ian Longden
Browse files

fixes for LRGSlices and new methods to find nearest gene to a feature (this...

fixes for LRGSlices and new methods to find nearest gene to a feature (this will change greatly as this is a very simple case)
parent 9c9a1ce1
...@@ -1872,6 +1872,76 @@ sub fetch_all_by_transcript_supporting_evidence { ...@@ -1872,6 +1872,76 @@ sub fetch_all_by_transcript_supporting_evidence {
return \@genes; return \@genes;
} }
=head2 fetch_nearest_Gene_by_Feature
Arg [1] : Feature object
Example : $genes = $gene_adaptor->fetch_nearest_Gene_by_Feature($feat);
Description: Gets the nearest gene to the feature
Returntype : Listref of Bio::EnsEMBL::Gene
Caller : general
Status : UnStable
=cut
sub fetch_nearest_Gene_by_Feature{
my $self = shift;
my $feat = shift;
my $stranded = shift;
my $min_dist = 0;
my $gene_id = 0;
my $overlapping = $feat->get_overlapping_Genes();
if(defined(@{$overlapping}[0])){
# print "found overlapping\n";
return @{$overlapping}[0];
# }
# else{
# print "No overlapping trying nearest\n";
# }
my $seq_region_id = $feat->slice->adaptor->get_seq_region_id($feat->slice);
my $start = ($feat->start + $feat->slice->start) -1;
my $end = ($feat->end + $feat->slice->start) -1;
my $sql1 = "select g.gene_id, (? - g.seq_region_end) as 'dist' from gene g where ";
if($stranded){
$sql1 .= "g.seq_region_strand = ".$feat->strand." and ";
}
$sql1 .= "seq_region_id = ? and g.seq_region_end < ? order by dist limit 1";
my $sql1_sth = $self->prepare($sql1) || die "Could not prepare $sql1";
$sql1_sth->execute($start, $seq_region_id, $start) || die "Could not execute sql";
$sql1_sth->bind_columns(\$gene_id, \$min_dist) || die "Could mot bin columns";
$sql1_sth->fetch() || die "Could not fetch result";
$sql1_sth->finish();
my $sql2 = "select g.gene_id, (g.seq_region_start - ?) as 'dist' from gene g where ";
if($stranded){
$sql2 .= "g.seq_region_strand = ".$feat->strand." and ";
}
$sql2 .= "seq_region_id = ? and g.seq_region_start > ? order by dist limit 1";
my $sql2_sth = $self->prepare($sql2) || die "could not prepare $sql2";
my ($tmp_min_dist, $tmp_gene_id);
$sql2_sth->execute($end, $seq_region_id, $end) || die "Could not execute sql";
$sql2_sth->bind_columns(\$tmp_gene_id, \$tmp_min_dist) || die "Could mot bin columns";
$sql2_sth->fetch() || die "Could not fetch result";
$sql2_sth->finish();
if($tmp_min_dist < $min_dist || !$gene_id){
$gene_id = $tmp_gene_id
}
if($gene_id){
return $self->fetch_by_dbID($gene_id);
}
return undef;
}
########################## ##########################
# # # #
......
...@@ -116,7 +116,7 @@ sub new { ...@@ -116,7 +116,7 @@ sub new {
rearrange(['START','END','STRAND','SLICE','ANALYSIS', 'SEQNAME', rearrange(['START','END','STRAND','SLICE','ANALYSIS', 'SEQNAME',
'DBID', 'ADAPTOR'], @_); 'DBID', 'ADAPTOR'], @_);
if($slice) { if($slice) {
if(!ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) { if(!ref($slice) || !($slice->isa('Bio::EnsEMBL::Slice') or $slice->isa('Bio::EnsEMBL::LRGSlice')) ) {
throw('-SLICE argument must be a Bio::EnsEMBL::Slice not '.$slice); throw('-SLICE argument must be a Bio::EnsEMBL::Slice not '.$slice);
} }
} }
...@@ -365,7 +365,7 @@ sub slice { ...@@ -365,7 +365,7 @@ sub slice {
if(@_) { if(@_) {
my $sl = shift; my $sl = shift;
if(defined($sl) && (!ref($sl) || !$sl->isa('Bio::EnsEMBL::Slice'))) { if(defined($sl) && (!ref($sl) || !($sl->isa('Bio::EnsEMBL::Slice') or $sl->isa('Bio::EnsEMBL::LRGSlice')) )) {
throw('slice argument must be a Bio::EnsEMBL::Slice'); throw('slice argument must be a Bio::EnsEMBL::Slice');
} }
...@@ -462,8 +462,8 @@ sub transform { ...@@ -462,8 +462,8 @@ sub transform {
my $projection = $self->project( $cs_name, $cs_version ); my $projection = $self->project( $cs_name, $cs_version );
if( @$projection != 1 and !defined($to_slice)) { if( @$projection != 1 and !defined($to_slice)) {
warn "MORE than one projection and NO slice specified "; # warn "MORE than one projection and NO slice specified ";
warn "from ".$self->slice->name." to $cs_name, $cs_version\n"; # warn "from ".$self->slice->name." to $cs_name, $cs_version\n";
return undef; return undef;
} }
my $index = 0; my $index = 0;
...@@ -1203,6 +1203,47 @@ sub overlaps { ...@@ -1203,6 +1203,47 @@ sub overlaps {
} }
=head2 get_overlapping_Genes
Description: Get all the genes that overlap this feature.
Returntype : list ref of Bio::EnsEMBL::Gene
Caller : general
Status : UnStable
=cut
sub get_overlapping_Genes{
my $self = shift;
my $slice = $self->feature_Slice;
return $slice->get_all_Genes();
}
# query for absolute nearest.
# select x.display_label, g.gene_id, g.seq_region_start, ABS(cast((32921638 - g.seq_region_end) as signed)) as 'dist' from gene g, xref x where g.display_xref_id = x.xref_id and seq_region_id = 27513 order by ABS(cast((32921638 - g.seq_region_end) as signed)) limit 10;
=head2 get_nearest_Gene
Description: Get all the nearest gene to the feature
Returntype : Bio::EnsEMBL::Gene
Caller : general
Status : UnStable
=cut
sub get_nearest_Gene {
my $self = shift;
my $stranded = shift;
my $ga = Bio::EnsEMBL::Registry->get_adaptor($self->adaptor->db->species,"core","Gene");
return $ga->fetch_nearest_Gene_by_Feature($self, $stranded);
}
############################################## ##############################################
......
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