t/splicedseq.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use lib 't/lib';
   8     use BioperlTest;
   9
  10     test_begin(-tests => 14);
  11
  12         use_ok('Bio::SeqIO');
  13 }
  14
  15 ok my $str = Bio::SeqIO->new(
  16                         '-file'=> test_input_file('U58726.gb'),
  17                         '-format' => 'GenBank');
  18
  19 my $seq;
  20
  21 ok ( $seq = $str->next_seq() );
  22
  23 # Here is a cute way to verify the sequence by seeing if the
  24 # the translation matches what is annotated in the file -js
  25 foreach my $ft ( grep { $_->primary_tag eq 'CDS'}
  26                  $seq->top_SeqFeatures ) {
  27     if( $ft->has_tag('translation') ) {
  28         my ($translation) = $ft->each_tag_value('translation');
  29         my $t = $ft->spliced_seq(-nosort => 1);
  30         my $pepseq = $t->translate()->seq();
  31         chop($pepseq);# chop is to remove stop codon
  32         is($translation,$pepseq);
  33         }
  34 }
  35
  36 my $stream = Bio::SeqIO->new(-file => test_input_file('M12730.gb'),
  37                               -format => 'genbank');
  38 # Jump down to M12730 which lists CDS join(1959..2355,1..92)
  39 while ($seq->accession ne "M12730") {
  40     $seq = $stream->next_seq;
  41 }
  42 ok(my @features = $seq->get_SeqFeatures(), "get_SeqFeatures()");
  43 my $feat;
  44 foreach my $feat2 ( @features ) {
  45     next unless ($feat2->primary_tag eq "CDS");
  46     my @db_xrefs = $feat2->get_tag_values("db_xref");
  47     if (grep { $_ eq "GI:150830" } @db_xrefs) {
  48        $feat = $feat2;
  49        last;
  50     }
  51 }
  52 my ($protein_seq) = $feat->get_tag_values("translation");
  53 like($protein_seq, qr(^MKERYGTVYKGSQRLIDE.*ANEKQENALYLIIILSRTSIT$),
  54          "protein sequence");
  55 my ($nucleotide_seq) = $feat->spliced_seq(-nosort => 1)->seq;
  56 like($nucleotide_seq, qr(^ATGAAAGAAAGATATGGA.*TCAAGGACTAGTATAACATAA$),
  57          "nucleotide sequence - correct CDS range");
  58 is(length($nucleotide_seq), 489, "nucleotide length");