SliceAdaptor.pm 18 KB
Newer Older
1 2

#
Simon Potter's avatar
pod  
Simon Potter committed
3
# Ensembl module for Bio::EnsEMBL::DBSQL::SliceAdaptor
4 5 6 7 8 9 10 11 12 13 14
#
# Cared for by Ewan Birney <ensembl-dev@ebi.ac.uk>
#
# Copyright Ewan Birney
#
# You may distribute this module under the same terms as perl itself

# POD documentation - main docs before the code

=head1 NAME

15
Bio::EnsEMBL::DBSQL::SliceAdaptor - Adaptors for slices
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

=head1 SYNOPSIS
  



=head1 DESCRIPTION

Factory for getting out slices of assemblies. WebSlice is the highly
accelerated version for the web site.

=head1 AUTHOR - Ewan Birney

This modules is part of the Ensembl project http://www.ensembl.org

Email ensembl-dev@ebi.ac.uk

Describe contact details here

=head1 APPENDIX

The rest of the documentation details each of the object
methods. Internal methods are usually preceded with a _

=cut


# Let the code begin...


46
package Bio::EnsEMBL::DBSQL::SliceAdaptor;
47 48 49
use vars qw(@ISA);
use strict;

50 51

# Object preamble - inherits from Bio::EnsEMBL::Root
52 53
use Bio::EnsEMBL::DBSQL::BaseAdaptor;
use Bio::EnsEMBL::Slice;
54
use Bio::EnsEMBL::DBSQL::DBAdaptor;
55

56

57
@ISA = ('Bio::EnsEMBL::DBSQL::BaseAdaptor');
58 59


Graham McVicker's avatar
Graham McVicker committed
60
# new is inherited from BaseAdaptor
61

62

63

64
=head2 fetch_by_chr_start_end
65

66 67 68 69 70 71 72 73 74 75 76 77 78
  Arg [1]    : string $chr
               the name of the chromosome to obtain a slice for
  Arg [2]    : int $start
               the start basepair of the slice to obtain in chromosomal 
               coordinates
  Arg [3]    : int $end 
               the end basepair of the slice to obtain in chromosomal 
               coordinates
  Example    : $slice = $slice_adaptor->fetch_by_chr_start_end();
  Description: Creates a slice object on the given chromosome and coordinates.
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : none
  Caller     : general
79 80 81

=cut

Graham McVicker's avatar
Graham McVicker committed
82
sub fetch_by_chr_start_end {
83 84
    my ($self,$chr,$start,$end) = @_;

85 86
    unless($chr) {
      $self->throw("chromosome name argument must be defined and not ''");
87 88
    }

89 90 91 92 93 94 95 96 97
    unless(defined $end) {   # Why defined?  Is '0' a valid end?
      $self->throw("end argument must be defined\n");
    }

    unless(defined $start) {
      $self->throw("start argument must be defined\n");
    }

    if($start > $end) {
98
      $self->throw("start must be less than end: parameters $chr:$start:$end");
99
    }
100
    
101
    my $slice;
102
    my $type = $self->db->assembly_type();
103

Graham McVicker's avatar
Graham McVicker committed
104
    $slice = Bio::EnsEMBL::Slice->new(
105 106 107 108
          -chr_name      => $chr,
          -chr_start     => $start,
          -chr_end       => $end,
          -assembly_type => $type,
109
          -adaptor       => $self
Graham McVicker's avatar
Graham McVicker committed
110
	 );
111 112 113 114 115 116

    return $slice;
}



117
=head2 fetch_by_contig_name
118

Graham McVicker's avatar
Graham McVicker committed
119 120 121 122 123 124 125 126 127 128 129
  Arg [1]    : string $name
               the name of the contig to obtain a slice for
  Arg [2]    : (optional) int $size
               the size of the flanking regions to obtain (aka context size)
  Example    : $slc = $slc_adaptor->fetch_by_contig_name('AB000878.1.1.33983');
  Description: Creates a slice object around the specified contig.  
               If a context size is given, the slice is extended by that 
               number of basepairs on either side of the contig.
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : none
  Caller     : general
130 131 132

=cut

133
sub fetch_by_contig_name {
Graham McVicker's avatar
Graham McVicker committed
134
   my ($self,$name, $size) = @_;
135 136 137

   if( !defined $size ) {$size=0;}

Graham McVicker's avatar
Graham McVicker committed
138
   my ($chr_name,$start,$end) = $self->_get_chr_start_end_of_contig($name);
139

140 141 142 143 144 145 146 147
   $start -= $size;
   $end += $size;

   if($start < 1) {
     $start  = 1;
   }

   return $self->fetch_by_chr_start_end($chr_name, $start, $end);
Graham McVicker's avatar
Graham McVicker committed
148 149 150 151
 }



Graham McVicker's avatar
Graham McVicker committed
152
=head2 fetch_by_fpc_name
Graham McVicker's avatar
Graham McVicker committed
153

Graham McVicker's avatar
Graham McVicker committed
154 155 156 157 158 159 160
  Arg [1]    : string $fpc_name
  Example    : my $slice = $slice_adaptor->fetch_by_fpc_name('NT_004321');
  Description: Creates a Slice on the region of the assembly where 
               the specified FPC (super) contig lies.
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : none
  Caller     : general
Graham McVicker's avatar
Graham McVicker committed
161 162 163 164 165

=cut

sub fetch_by_fpc_name {
    my ($self,$fpc_name) = @_;
Arne Stabenau's avatar
Arne Stabenau committed
166 167 168
    
    my( $p, $f, $l ) = caller; 
    $self->warn( "$f:$l calls deprecated method fetch_by_fpc_name. Please use fetch_by_supercontig_name instead" );
Graham McVicker's avatar
Graham McVicker committed
169

Arne Stabenau's avatar
Arne Stabenau committed
170 171
    $self->fetch_by_supercontig_name( $fpc_name ); 
}
Graham McVicker's avatar
Graham McVicker committed
172

Arne Stabenau's avatar
Arne Stabenau committed
173 174 175 176 177 178
sub fetch_by_supercontig_name {
  my ($self,$supercontig_name) = @_;
  
  my $assembly_type = $self->db->assembly_type();
  
  my $sth = $self->db->prepare("
179 180
        SELECT chr.name, a.superctg_ori, MIN(a.chr_start), MAX(a.chr_end)
        FROM assembly a, chromosome chr
Arne Stabenau's avatar
Arne Stabenau committed
181 182
        WHERE superctg_name = ?
        AND type = ?
183
        AND chr.chromosome_id = a.chromosome_id
Graham McVicker's avatar
Graham McVicker committed
184 185 186
        GROUP by superctg_name
        ");

Arne Stabenau's avatar
Arne Stabenau committed
187 188 189 190 191 192 193 194 195 196 197 198
  $sth->execute( $supercontig_name, $assembly_type );
  
  my ($chr, $strand, $slice_start, $slice_end) = $sth->fetchrow_array;
  
  my $slice;
  
  $slice = new Bio::EnsEMBL::Slice
    (
     -chr_name => $chr,
     -chr_start =>$slice_start,
     -chr_end => $slice_end,
     -strand => $strand,
199 200
     -assembly_type => $assembly_type,
     -adaptor => $self
Arne Stabenau's avatar
Arne Stabenau committed
201 202 203 204 205
    );
  
  return $slice;
}

Graham McVicker's avatar
Graham McVicker committed
206

Arne Stabenau's avatar
Arne Stabenau committed
207
=head2 list_overlapping_supercontigs
Graham McVicker's avatar
Graham McVicker committed
208

Arne Stabenau's avatar
Arne Stabenau committed
209 210 211 212 213 214 215 216 217
  Arg [1]    : Bio::EnsEMBL::Slice $slice
               overlapping given Sice
  Example    : 
  Description: return the names of the supercontigs that overlap given Slice.  
  Returntype : listref string
  Exceptions : none
  Caller     : general

=cut
Graham McVicker's avatar
Graham McVicker committed
218 219


Arne Stabenau's avatar
Arne Stabenau committed
220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237
sub list_overlapping_supercontigs {
   my ($self,$slice) = @_;
   my $sth = $self->db->prepare( "
      SELECT DISTINCT superctg_name
        FROM assembly a, chromosome c
       WHERE c.chromosome_id = a.chromosome_id 
         AND c.name = ?
         AND a.type = ?
         AND a.chr_end >= ?
         AND a.chr_start <= ?
       " );
   $sth->execute( $slice->chr_name(), $slice->assembly_type(),
		  $slice->chr_start(), $slice->chr_end() );

   my $result = [];
   while( my $aref = $sth->fetchrow_arrayref() ) {
     push( @$result, $aref->[0] );
   }
238

Arne Stabenau's avatar
Arne Stabenau committed
239 240
   return $result;
}
241 242


Web Admin's avatar
Web Admin committed
243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277
=head2 fetch_by_chr_band

 Title   : fetch_by_chr_band
 Usage   :
 Function: create a Slice representing a series of bands
 Example :
 Returns :
 Args    : the band name


=cut

sub fetch_by_chr_band {
    my ($self,$chr,$band) = @_;

    my $type = $self->db->assembly_type();

    my $sth = $self->db->prepare("
        select min(k.chr_start), max(k.chr_end)
          from chromosome as c, karyotype as k
         where c.chromosome_id = k.chromosome_id and c.name=? and k.band like ?
    ");
    $sth->execute( $chr, "$band%" );
    my ( $slice_start, $slice_end) = $sth->fetchrow_array;

    unless( defined($slice_start) ) {
       my $sth = $self->db->prepare("
           select min(k.chr_start), max(k.chr_end)
             from chromosome as c, karyotype as k
            where c.chromosome_id = k.chromosome_id and k.band like ?
       ");
       $sth->execute( "$band%" );
       ( $slice_start, $slice_end) = $sth->fetchrow_array;
    }

278 279 280 281 282 283 284 285 286 287 288
   if(defined $slice_start) {
        return new Bio::EnsEMBL::Slice(
           -chr_name  => $chr,
           -chr_start => $slice_start,
           -chr_end   => $slice_end,
           -strand    => 1,
           -assembly_type => $type
        );
    }

    $self->throw("Band not recognised in database");
Web Admin's avatar
Web Admin committed
289 290 291
}


Graham McVicker's avatar
Graham McVicker committed
292 293
=head2 fetch_by_clone_accession

Graham McVicker's avatar
Graham McVicker committed
294 295 296 297 298 299 300 301 302 303 304
  Arg [1]    : string $clone 
               the embl accession of the clone object to retrieve
  Arg [2]    : (optional) int $size
               the size of the flanking regions to obtain around the clone 
  Example    : $slc = $slc_adaptor->fetch_by_clone_accession('AC000012',1000);
  Description: Creates a Slice around the specified clone.  If a context size 
               is given, the Slice is extended by that number of basepairs on 
               either side of the clone.  Throws if the clone is not golden.
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : thrown if the clone is not in the assembly 
  Caller     : general
305 306 307

=cut

Graham McVicker's avatar
Graham McVicker committed
308
sub fetch_by_clone_accession{
309 310 311
   my ($self,$clone,$size) = @_;

   if( !defined $clone ) {
Graham McVicker's avatar
Graham McVicker committed
312
     $self->throw("Must have clone to fetch Slice of clone");
313 314 315
   }
   if( !defined $size ) {$size=0;}

316
   my $type = $self->db->assembly_type()
317 318 319 320 321
    or $self->throw("No assembly type defined");

   my $sth = $self->db->prepare("SELECT  c.name,
                        a.chr_start,
                        a.chr_end,
322
                        chr.name 
323 324
                    FROM    assembly a, 
                        contig c, 
325 326
                        clone  cl,
                        chromosome chr
327 328 329
                    WHERE c.clone_id = cl.clone_id
                    AND cl.name = '$clone'  
                    AND c.contig_id = a.contig_id 
330 331
                    AND a.type = '$type'
                    AND chr.chromosome_id = a.chromosome_id
332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348
                    ORDER BY a.chr_start"
                    );
   $sth->execute();
 
   my ($contig,$start,$end,$chr_name); 
   my $counter; 
   my $first_start;
   while ( my @row=$sth->fetchrow_array){
       $counter++;
       ($contig,$start,$end,$chr_name)=@row;
       if ($counter==1){$first_start=$start;}      
   }

   if( !defined $contig ) {
       $self->throw("Clone is not on the golden path. Cannot build Slice");
   }
     
349 350 351 352 353 354 355 356
   $first_start -= $size;
   $end += $size;

   if($first_start < 1) {
     $first_start = 1;
   }

   my $slice = $self->fetch_by_chr_start_end($chr_name, $first_start, $end);
357 358 359 360 361
   return $slice;
}



Graham McVicker's avatar
Graham McVicker committed
362
=head2 fetch_by_transcript_stable_id
363

Graham McVicker's avatar
Graham McVicker committed
364 365 366 367 368 369 370 371 372 373 374 375 376 377
  Arg [1]    : string $transcriptid
               The stable id of the transcript around which the slice is 
               desired
  Arg [2]    : (optional) int $size
               The length of the flanking regions the slice should encompass 
               on either side of the transcript (0 by default)
  Example    : $slc = $sa->fetch_by_transcript_stable_id('ENST00000302930',10);
  Description: Creates a slice around the region of the specified transcript. 
               If a context size is given, the slice is extended by that 
               number of basepairs on either side of the 
               transcript.  Throws if the transcript is not golden.
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : none
  Caller     : general
378 379 380

=cut

Graham McVicker's avatar
Graham McVicker committed
381
sub fetch_by_transcript_stable_id{
382 383 384 385 386 387
  my ($self,$transcriptid,$size) = @_;

  # Just get the dbID, then fetch slice by that
  my $ta = $self->db->get_TranscriptAdaptor;
  my $transcript_obj = $ta->fetch_by_stable_id($transcriptid);
  my $dbID = $transcript_obj->dbID;
Graham McVicker's avatar
Graham McVicker committed
388 389
  
  return $self->fetch_by_transcript_id($dbID, $size);
390 391
}

392

Graham McVicker's avatar
Graham McVicker committed
393 394


Graham McVicker's avatar
Graham McVicker committed
395 396
=head2 fetch_by_transcript_id

Graham McVicker's avatar
Graham McVicker committed
397 398 399 400 401 402 403 404 405 406 407 408 409 410
  Arg [1]    : int $transcriptid
               The unique database identifier of the transcript around which 
               the slice is desired
  Arg [2]    : (optional) int $size
               The length of the flanking regions the slice should encompass 
               on either side of the transcript (0 by default)
  Example    : $slc = $sa->fetch_by_transcript_id(24, 1000);
  Description: Creates a slice around the region of the specified transcript. 
               If a context size is given, the slice is extended by that 
               number of basepairs on either side of the 
               transcript. 
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : thrown on incorrect args
  Caller     : general
411 412 413

=cut

Graham McVicker's avatar
Graham McVicker committed
414
sub fetch_by_transcript_id {
415
  my ($self,$transcriptid,$size) = @_;
Graham McVicker's avatar
Graham McVicker committed
416 417

  unless( defined $transcriptid ) {
Graham McVicker's avatar
Graham McVicker committed
418 419
    $self->throw("Must have transcriptid id to fetch Slice of transcript");
  }
Graham McVicker's avatar
Graham McVicker committed
420 421 422

  $size = 0 unless(defined $size);
   
Graham McVicker's avatar
Graham McVicker committed
423 424 425 426 427
  my $ta = $self->db->get_TranscriptAdaptor;
  my $transcript_obj = $ta->fetch_by_dbID($transcriptid);
  
  my %exon_transforms;
  
428
  my $emptyslice;
Graham McVicker's avatar
Graham McVicker committed
429 430 431 432 433 434 435 436 437 438 439
  for my $exon ( @{$transcript_obj->get_all_Exons()} ) {
    $emptyslice = Bio::EnsEMBL::Slice->new( '-empty'   => 1,
					    '-adaptor' => $self,
					    '-ASSEMBLY_TYPE' =>
					    $self->db->assembly_type);     
    my $newExon = $exon->transform( $emptyslice );
    $exon_transforms{ $exon } = $newExon;
  }
  
  $transcript_obj->transform( \%exon_transforms );
  
440 441
  my $start = $transcript_obj->start() - $size;
  my $end = $transcript_obj->end() + $size;
Graham McVicker's avatar
Graham McVicker committed
442
  
443 444 445
  if($start < 1) {
    $start = 1;
  }
446
  
447 448 449
  my $slice = $self->fetch_by_chr_start_end($emptyslice->chr_name,
					    $start, $end);
  return $slice;
450 451
}

452 453


Graham McVicker's avatar
Graham McVicker committed
454
=head2 fetch_by_gene_stable_id
455

Graham McVicker's avatar
Graham McVicker committed
456 457 458 459
  Arg [1]    : string $geneid
               The stable id of the gene around which the slice is 
               desired
  Arg [2]    : (optional) int $size
Graham McVicker's avatar
Graham McVicker committed
460
               The length of the flanking regions the slice should encompass
Graham McVicker's avatar
Graham McVicker committed
461 462
               on either side of the gene (0 by default)
  Example    : $slc = $sa->fetch_by_transcript_stable_id('ENSG00000012123',10);
Graham McVicker's avatar
Graham McVicker committed
463 464 465
  Description: Creates a slice around the region of the specified gene.
               If a context size is given, the slice is extended by that
               number of basepairs on either side of the gene.
Graham McVicker's avatar
Graham McVicker committed
466 467 468
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : none
  Caller     : general
469 470 471

=cut

Graham McVicker's avatar
Graham McVicker committed
472
sub fetch_by_gene_stable_id{
473 474 475 476 477 478 479
   my ($self,$geneid,$size) = @_;

   if( !defined $geneid ) {
       $self->throw("Must have gene id to fetch Slice of gene");
   }
   if( !defined $size ) {$size=0;}

Graham McVicker's avatar
Graham McVicker committed
480
   my ($chr_name,$start,$end) = $self->_get_chr_start_end_of_gene($geneid);
481 482

   if( !defined $start ) {
483
     my $type = $self->db->assembly_type()
Graham McVicker's avatar
Graham McVicker committed
484
       or $self->throw("No assembly type defined");
Graham McVicker's avatar
Graham McVicker committed
485 486
     $self->throw("Gene is not on the golden path '$type'. " .
		  "Cannot build Slice.");
487 488
   }
     
489 490 491 492 493 494 495 496
   $start -= $size;
   $end += $size;
   
   if($start < 1) {
     $start = 1;
   }

   return $self->fetch_by_chr_start_end($chr_name, $start, $end);
497 498 499
}


500

Graham McVicker's avatar
Graham McVicker committed
501
=head2 fetch_by_chr_name
Graham McVicker's avatar
Graham McVicker committed
502

Graham McVicker's avatar
Graham McVicker committed
503 504 505 506 507 508
  Arg [1]    : string $chr_name
  Example    : $slice = $slice_adaptor->fetch_by_chr_name('20'); 
  Description: Retrieves a slice on the region of an entire chromosome
  Returntype : Bio::EnsEMBL::Slice
  Exceptions : thrown if $chr_name arg is not supplied
  Caller     : general
Graham McVicker's avatar
Graham McVicker committed
509 510 511 512 513 514 515 516 517 518 519 520 521 522

=cut

sub fetch_by_chr_name{
   my ($self,$chr_name) = @_;

   unless( $chr_name ) {
       $self->throw("Chromosome name argument required");
   }

   my $chr_start = 1;
   
   #set the end of the slice to the end of the chromosome
   my $ca = $self->db()->get_ChromosomeAdaptor();
523
   my $chromosome = $ca->fetch_by_chr_name($chr_name);
Graham McVicker's avatar
Graham McVicker committed
524 525
   my $chr_end = $chromosome->length();

526 527 528 529 530 531 532 533 534 535 536 537
   my $type = $self->db->assembly_type();

   my $slice = Bio::EnsEMBL::Slice->new
     (
      -chr_name      => $chr_name,
      -chr_start     => 1,
      -chr_end       => $chr_end,
      -assembly_type => $type,
      -adaptor       => $self
     );

   return $slice;
Graham McVicker's avatar
Graham McVicker committed
538 539
}

Graham McVicker's avatar
Graham McVicker committed
540 541


542 543 544 545 546 547 548 549 550 551 552 553
=head2 fetch_by_mapfrag

 Title   : fetch_by_mapfrag
 Usage   : $slice = $slice_adaptor->fetch_by_mapfrag('20');
 Function: Creates a slice of a "mapfrag"
 Returns : Slice object
 Args    : chromosome name


=cut

sub fetch_by_mapfrag{
554
   my ($self,$mymapfrag,$flag,$size) = @_;
555 556 557

   $flag ||= 'fixed-width'; # alt.. 'context'
   $size ||= $flag eq 'fixed-width' ? 200000 : 0;
558
   unless( $mymapfrag ) {
559 560 561 562 563 564 565
       $self->throw("Mapfrag name argument required");
   }

   my( $chr_start,$chr_end);
  
   #set the end of the slice to the end of the chromosome
   my $ca = $self->db()->get_MapFragAdaptor();
566
   my $mapfrag = $ca->fetch_by_synonym($mymapfrag);
567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592
   return undef unless defined $mapfrag;

   if( $flag eq 'fixed-width' ) {
       my $halfsize = int( $size/2 );
       $chr_start = $mapfrag->seq_start - $halfsize;
       $chr_end   = $mapfrag->seq_start + $size - $halfsize;
   } else {
       $chr_start     = $mapfrag->seq_start - $size;
       $chr_end       = $mapfrag->seq_end   + $size;
   }
   my $type = $self->db->assembly_type();

   my $slice = Bio::EnsEMBL::Slice->new
     (
      -chr_name      => $mapfrag->seq,
      -chr_start     => $chr_start,
      -chr_end       => $chr_end,
      -assembly_type => $type,
      -adaptor       => $self
     );

   return $slice;
}



Graham McVicker's avatar
Graham McVicker committed
593 594 595 596 597 598 599 600 601 602 603 604 605 606 607


=head2 _get_chr_start_end_of_contig

 Title   : _get_chr_start_end_of_contig
 Usage   :
 Function: returns the chromosome name, absolute start and absolute end of the 
           specified contig
 Returns : returns chr,start,end
 Args    : contig id

=cut

sub _get_chr_start_end_of_contig {
    my ($self,$contigid) = @_;
608

Graham McVicker's avatar
Graham McVicker committed
609 610 611 612 613 614 615 616 617 618
   if( !defined $contigid ) {
       $self->throw("Must have contig id to fetch Slice of contig");
   }
   
   my $type = $self->db->assembly_type()
    or $self->throw("No assembly type defined");

   my $sth = $self->db->prepare("SELECT  c.name,
                        a.chr_start,
                        a.chr_end,
619 620
                        chr.name 
                    FROM assembly a, contig c, chromosome chr 
Graham McVicker's avatar
Graham McVicker committed
621 622
                    WHERE c.name = '$contigid' 
                    AND c.contig_id = a.contig_id 
623 624
                    AND a.type = '$type'
                    AND chr.chromosome_id = a.chromosome_id"
Graham McVicker's avatar
Graham McVicker committed
625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649
                    );
   $sth->execute();
   my ($contig,$start,$end,$chr_name) = $sth->fetchrow_array;

   if( !defined $contig ) {
     $self->throw("Contig $contigid is not on the golden path of type $type");
   }

   return ($chr_name,$start,$end);
}

=head2 _get_chr_start_end_of_gene

 Title   : get_Gene_chr_bp
 Usage   : 
 Function: 
 Returns :  
 Args    :


=cut


sub _get_chr_start_end_of_gene {
  my ($self,$geneid) =  @_;
650
  
Graham McVicker's avatar
Graham McVicker committed
651 652 653 654 655 656 657 658
  my $type = $self->db->assembly_type()
    or $self->throw("No assembly type defined");
  
  my $sth = $self->db->prepare("SELECT  
   if(a.contig_ori=1,(e.contig_start-a.contig_start+a.chr_start),
                    (a.chr_start+a.contig_end-e.contig_end)),
   if(a.contig_ori=1,(e.contig_end-a.contig_start+a.chr_start),
                    (a.chr_start+a.contig_end-e.contig_start)),
659
     chr.name
Graham McVicker's avatar
Graham McVicker committed
660 661 662 663 664
  
                    FROM    exon e,
                        transcript tr,
                        exon_transcript et,
                        assembly a,
665
                        gene_stable_id gsi,
666
                        chromosome chr
Graham McVicker's avatar
Graham McVicker committed
667 668 669 670 671
                    WHERE e.exon_id=et.exon_id 
                    AND et.transcript_id =tr.transcript_id 
                    AND a.contig_id=e.contig_id 
                    AND a.type = '$type' 
                    AND tr.gene_id = gsi.gene_id
672
                    AND gsi.stable_id = '$geneid'
673
                    AND a.chromosome_id = chr.chromosome_id" 
Graham McVicker's avatar
Graham McVicker committed
674 675
                    );
   $sth->execute();
676

Graham McVicker's avatar
Graham McVicker committed
677 678 679 680 681 682 683 684 685
   my ($start,$end,$chr);
   my @start;
   while ( my @row=$sth->fetchrow_array){
      ($start,$end,$chr)=@row;
       push @start,$start;
       push @start,$end;
   }   
   
   my @start_sorted=sort { $a <=> $b } @start;
686

Graham McVicker's avatar
Graham McVicker committed
687 688
   $start=shift @start_sorted;
   $end=pop @start_sorted;
689

Graham McVicker's avatar
Graham McVicker committed
690
   return ($chr,$start,$end);      
691
}
Graham McVicker's avatar
Graham McVicker committed
692

693
1;
Graham McVicker's avatar
Graham McVicker committed
694 695 696 697 698 699