[Bioperl-l] Bio::Coordinate::GeneMapper cds to peptide bug

[Chris Fields]

Aaron,

Do we want to write this up as a set of tests to add to the bioperl test suite?  We can probably add it after the github migration tomorrow.

chris

On May 11, 2010, at 4:26 PM, Aaron Mackey wrote:

> Hi Zerui (and others),
> 
> I've confirmed there seems to be a bug in Bio::Coordinate::GeneMapper,
> specifically in this code:
> 
> lines:
> 1170:               (-start   => int ($loc->start / 3 ) +1,
> 1171:            -end     => int ($loc->end / 3 ) +1,
> 
> both of those lines should look like: int (($loc->start - 1) / 3) + 1
> 
> otherwise for CDS coordinates 1, 2, 3, 4, 5, 6, you get incorrect peptide
> positions 1, 1, 2, 2, 2, 3 (when you instead want 1, 1, 1, 2, 2, 2)
> 
> There is also a problem when mapping exon coordinates that are outside/after
> the CDS region (instead of getting undefined locations, you continue to get
> peptide coordinates, but they are invalid, larger than the protein length).
> 
> Dennis and fringy -- this may affect the SNPtab.pl script I wrote for you,
> as it uses this module to calculate codons for SNPs.
> 
> -Aaron
> 
> P.S. a script the demonstrates the problem:
> 
> use Bio::Coordinate::GeneMapper;
> 
> my $mapper =
>  Bio::Coordinate::GeneMapper
>  ->new( -in => "chr",
>     -out => "propeptide",
>     -exons => [ Bio::Location::Simple
>             ->new( -start => 101,
>                -end   => 109 ),
>             Bio::Location::Simple
>             ->new( -start => 201,
>                -end   => 221 ),
>           ],
>     -cds => Bio::Location::Simple
>             ->new(-start => 101, -end => 209),
>       );
> 
> 
> print join("\t", "chr", "aa"), "\n";
> for my $pos (99..111,199..211) {
>  my $res = $mapper->map(
>    Bio::Location::Simple->new(-start => $pos, -end => $pos, -seq_id => 1));
>  my $start = $res->start; $start = "NA" unless defined $start;
>  my $end = $res->end; $end = "NA" unless defined $end;
>  print join("\t", $pos, $start), "\n";
> }
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l