[Bioperl-l] Bio::RangeI intersection and Bio::DB::GFF

[Brian Osborne]

Marco,

Odd, because the intersection() code is quite simple and it's clear how it
should behave. What version of Bioperl are you using? I'm looking at the
latest, in bioperl-live...

Brian O.


On 5/2/06 4:32 PM, "Marco Blanchette" <mblanche at berkeley.edu> wrote:

> Brian--
> 
> Even when both elements of intersection() are from the negative strand, the
> return object is from the positive strand and $overlap is actually the
> revervese complement of the intersection between the 2 exons. Here is part
> of the output from the script below:
> 
> ===
> ex1     Strand: -1
> CTTTTTTCCACATACGTCGTCAACGTGATTCGACCTTTTCCGGTTTATTAGTTGAACATGGCAGTCGGCAAAAATA
> AAGGTCTTTCCAAGGGTGGTAAGAAGGGCGG
> TAAGAAGAAGGTGGTGGACCCGTTTTCTCGCAAGGACTG
> ex2     Strand: -1
> CTTTTTTCCACATACGTCGTCAACGTGATTCGACCTTTTCCGGTTTATTAGTTGAACATGGCAGTCGGCAAAAATA
> AAGGTCTTTCCAAGGGTGGTAAGAAGGGCGG
> TAAGAAGAAGGTGGTGGACCCGTTTTCTCGCAAGGACTGGTACGATGTCAAAGCTCCGAATATGTTTCAAACCCGT
> CAAATCG
> overlap Strand: 1
> CAGTCCTTGCGAGAAAACGGGTCCACCACCTTCTTCTTACCGCCCTTCTTACCACCCTTGGAAAGACCTTTATTTT
> TGCCGACTGCCATGTTCAACTAATAAACCGG
> AAAAGGTCGAATCACGTTGACGACGTATGTGGAAAAAAG
> ...
> 
> If both are from the positive strand, the return object is positive as in:
> 
> ===
> ex1     Strand: 1
> CAACGCAGACGTGGTACGGCGTTTTAAATCTGATAACATTTTGAACCGGGAATTATTTTAGAGTACCATTCTTTGT
> TTTGTGCCTGTTTCAGTATAAATTAATTATG
> CGCCTGATTTAAAGTACAAAATGTGTAAATATATCACCTTACCGTCGCGGGTGCACCCAATTGTGCTTTGATGAAT
> AAATATACATATATGCAACATATATAACTTC
> CTGTGTTAGTATAAGTGTATGTCAGCCAAAAACAAATATATATATGAGTGTTTATCGGCATTCGTGTGCTGGCAGA
> GCAGCGATCAAAGCTGCGTTCGGTACTCGTT
> GACTGGCCCAAGAATGAATTCTCGTGCAAGTGTGTTGATAAAAAGTATACGTATGTAT
> ex2     Strand: 1
> ATCGACAGTTGCCATCGTCGTTATTCCAGCACTAATTTAAAAAAAATTCGATCAACGCAGACGTG
> overlap Strand: 1
> CAACGCAGACGTG
> 
> Is there something I am missing? Here is the script generating the output
> 
> Many thanks all...
> 
> Marco
> 
> 
> use strict;
> use warnings;
> use Bio::DB::GFF;
> 
> MAIN:{
> 
>     my $db = Bio::DB::GFF->new( -adaptor => 'dbi::mysql',
>                                 -dsn =>
> 'dbi:mysql:database=dmel_43_LS;host=riolab.net',
>                                 -user => 'guest');
>     my $test_db = $db->segment('4');
>     
>     # Load up the exons into $exons_p
>     for my $gene ($test_db->features(-types => 'gene')){
> 
>         my $exons_p = extractExons($gene);
> 
>         cluster($exons_p) unless ($#{$exons_p} == -1);
> 
>     }
> }
> 
> sub extractExons {
>     my $gene = shift;
>     my %ex_list;
>     my @tcs = $gene->features(    -type =>'processed_transcript',
>                                     -attributes =>{Gene => $gene->group});
>                  
>     for my $tc (@tcs){
>         my @exons = $tc->features (-type => 'exon',
>                                      -attributes => {Parent => $tc->group}
> );        
>     
>         for (@exons){
>             my $ex_id    = $_->id;
>             $ex_list{$ex_id} = $_ unless (exists $ex_list{$ex_id});
> 
>         }
>     
>     }
>     my @values = values %ex_list;
>     return(\@values);
> }
> 
> sub cluster {
>     my $exons_p = shift;
>     
>     for (my $s = 0; $s <= $#{$exons_p}; $s++){
>         for (my $t = $s+1; $t <= $#{$exons_p}; $t++){
>             my $exon1 = $exons_p->[$s];
>             my $exon2 = $exons_p->[$t];
>             
>             if (!($exon1->equals($exon2)) && $exon1->overlaps($exon2)){
>             
>                 my $overlap = $exon1->intersection($exon2);
>                 
>                 print "===\n";;
>                 print     "ex1\tStrand: ", $exon1->strand, "\n",
>                         $exon1->seq, "\n";
>                 print     "ex2\tStrand: ", $exon2->strand, "\n",
>                         $exon2->seq, "\n";
>                 print "overlap\tStrand: ", $overlap->strand, "\n",
>                         $overlap->seq, "\n";
>             }
>         }
>     }
> }
> 
> On 5/2/06 13:17, "Brian Osborne" <osborne1 at optonline.net> wrote:
> 
>> Marco,
>> 
>> Yes, this is how intersection() is supposed to work. If both of the Range
>> objects have the same strand then the strand information is returned as part
>> of the result but if they aren't on the same strand then no strand
>> information is returned.
>> 
>> Brian O.
>> 
>> 
>> On 5/2/06 3:30 PM, "Marco Blanchette" <mblanche at berkeley.edu> wrote:
>> 
>>> Dear all--
>>> 
>>> I have been trying to use the intersection function to extract overlapping
>>> region from alternatively spliced exons as in the following script. The
>>> returned object from the 'my $overlap = $exon1->intersection($exon2);' is
>>> actually loosing the strand of $exon1 if $exon1 is from the negative strand.
>>> Is this behavior expected? Should I check the strand of $exon1 before
>>> working on the object return by any Bio::RangeI function?
>>> 
>>> Many thanks 
>>> 
>>> #!/usr/bin/perl
>>> use strict;
>>> use warnings;
>>> use Bio::DB::GFF;
>>> 
>>> MAIN:{
>>> 
>>>     my $db = Bio::DB::GFF->new( -adaptor => 'dbi::mysql',
>>>                                 -dsn =>
>>> 'dbi:mysql:database=dmel_43_LS;host=riolab.net',
>>>                                 -user => 'guest');
>>>     my $test_db = $db->segment('4');
>>>     
>>>     # Load up the exons into $exons_p
>>>     for my $gene ($test_db->features(-types => 'gene')){
>>> 
>>>         my $exons_p = extractExons($gene);
>>> 
>>>         cluster($exons_p) unless ($#{$exons_p} == -1);
>>> 
>>>     }
>>> }
>>> 
>>> sub extractExons {
>>>     my $gene = shift;
>>>     my %ex_list;
>>>     my @tcs = $gene->features(    -type =>'processed_transcript',
>>>                                     -attributes =>{Gene => $gene->group});
>>>                
>>>     for my $tc (@tcs){
>>>         my @exons = $tc->features (-type => 'exon',
>>>                                      -attributes => {Parent => $tc->group}
>>> );        
>>>     
>>>         for (@exons){
>>>             my $ex_id    = $_->id;
>>>             $ex_list{$ex_id} = $_ unless (exists $ex_list{$ex_id});
>>> 
>>>         }
>>>     
>>>     }
>>>     my @values = values %ex_list;
>>>     return(\@values);
>>> }
>>> 
>>> sub cluster {
>>>     my $exons_p = shift;
>>>     
>>>     for (my $s = 0; $s <= $#{$exons_p}; $s++){
>>>         for (my $t = $s+1; $t <= $#{$exons_p}; $t++){
>>>             my $exon1 = $exons_p->[$s];
>>>             my $exon2 = $exons_p->[$t];
>>>             
>>>             if (!($exon1->equals($exon2)) && $exon1->overlaps($exon2)){
>>>             
>>>                 my $overlap = $exon1->intersection($exon2);
>>>                
>>>                 print "===\n";;
>>>                 print "ex1\n", $exon1->seq, "\n";
>>>                 print "ex2\n", $exon2->seq, "\n";
>>>                 print "overlap\n", $overlap->seq, "\n";
>>>             }
>>>         }
>>>     }
>>> }
>>> ______________________________
>>> Marco Blanchette, Ph.D.
>>> 
>>> mblanche at uclink.berkeley.edu
>>> 
>>> Donald C. Rio's lab
>>> Department of Molecular and Cell Biology
>>> 16 Barker Hall
>>> University of California
>>> Berkeley, CA 94720-3204
>>> 
>>> Tel: (510) 642-1084
>>> Cell: (510) 847-0996
>>> Fax: (510) 642-6062
>> 
>> 
> 
> ______________________________
> Marco Blanchette, Ph.D.
> 
> mblanche at uclink.berkeley.edu
> 
> Donald C. Rio's lab
> Department of Molecular and Cell Biology
> 16 Barker Hall
> University of California
> Berkeley, CA 94720-3204
> 
> Tel: (510) 642-1084
> Cell: (510) 847-0996
> Fax: (510) 642-6062