t/SeqIO/Splicedseq.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use Bio::Root::Test;
   8
   9     test_begin(-tests => 27);
  10
  11     use_ok('Bio::Seq');
  12     use_ok('Bio::SeqIO');
  13 }
  14
  15 ok my $str = Bio::SeqIO->new(-file   => test_input_file('U58726.gb'),
  16                              -format => 'GenBank');
  17 my $seq;
  18 ok ( $seq = $str->next_seq() );
  19
  20 # Here is a cute way to verify the sequence by seeing if the
  21 # the translation matches what is annotated in the file -js
  22 foreach my $ft ( grep { $_->primary_tag eq 'CDS'}
  23                  $seq->top_SeqFeatures ) {
  24     if( $ft->has_tag('translation') ) {
  25         my ($translation) = $ft->get_tag_values('translation');
  26         my $t = $ft->spliced_seq(-nosort => 1);
  27         my $pepseq = $t->translate()->seq();
  28         chop($pepseq); # chop is to remove stop codon
  29         is($translation, $pepseq);
  30     }
  31 }
  32
  33 my $stream = Bio::SeqIO->new(-file   => test_input_file('M12730.gb'),
  34                              -format => 'genbank');
  35 # Jump down to M12730 which lists CDS join(1959..2355,1..92)
  36 while ($seq->accession ne "M12730") {
  37     $seq = $stream->next_seq;
  38 }
  39 ok(my @features = $seq->get_SeqFeatures(), "get_SeqFeatures()");
  40 my $feat;
  41 foreach my $feat2 ( @features ) {
  42     next unless ($feat2->primary_tag eq "CDS");
  43     my @db_xrefs = $feat2->get_tag_values("db_xref");
  44     if (grep { $_ eq "GI:150830" } @db_xrefs) {
  45        $feat = $feat2;
  46        last;
  47     }
  48 }
  49 my ($protein_seq) = $feat->get_tag_values("translation");
  50 like($protein_seq, qr(^MKERYGTVYKGSQRLIDE.*ANEKQENALYLIIILSRTSIT$),
  51                    "protein sequence");
  52 my ($nucleotide_seq) = $feat->spliced_seq(-nosort => 1)->seq;
  53 like($nucleotide_seq, qr(^ATGAAAGAAAGATATGGA.*TCAAGGACTAGTATAACATAA$),
  54                      "nucleotide sequence - correct CDS range");
  55 is(length($nucleotide_seq), 489, "nucleotide length");
  56
  57 #  Test for Fix spliced seq #72
  58 my $str2 = Bio::SeqIO->new(-file   => test_input_file('AF032047.gbk'),
  59                            -format => 'GenBank');
  60 my @feats = $str2-> next_seq -> get_SeqFeatures;
  61 # feat[1] has 2 exons from remote sequence AF032048.1
  62 my $len_nodb;
  63 warnings_like { $len_nodb = length($feats[1]->spliced_seq()->seq); }
  64               [ {carped => qr/cannot get remote location for/},
  65                 {carped => qr/cannot get remote location for/}
  66                ],
  67               "appropriate warning if db not provided for remote sequence";
  68 ok($len_nodb == 374, "correct number of Ns added if remote sequence not provided");
  69
  70 # Test for cut by origin features
  71 my $seq_obj = Bio::Seq->new(-display_id => 'NC_008309',
  72                             -seq        => 'AAAAACCCCCGGGGGTTTTT');
  73 $seq_obj->is_circular(1);
  74 my $loc_obj = Bio::Factory::FTLocationFactory->from_string('join(16..20,1..2)');
  75 my $cut_feat = Bio::SeqFeature::Generic->new(-primary_tag => 'CDS',
  76                                              -location    => $loc_obj,
  77                                              -tag => { locus_tag  => 'HS_1792',
  78                                                        product    => 'hypothetical protein',
  79                                                        protein_id => 'YP_718205.1',
  80                                                       } );
  81 $seq_obj->add_SeqFeature($cut_feat);
  82 is $cut_feat->seq->seq,         'TTTTTAA', 'cut by origin sequence using $feat->seq';
  83 is $cut_feat->spliced_seq->seq, 'TTTTTAA', 'cut by origin sequence using $feat->spliced_seq';
  84 is $cut_feat->start,             16,       'cut by origin start using $feat->start';
  85 is $cut_feat->end,               2,        'cut by origin end using $feat->end';
  86 is $cut_feat->location->start,   16,       'cut by origin start using $feat->location->start';
  87 is $cut_feat->location->end,     2,        'cut by origin end using $feat->location->end';
  88
  89 SKIP: {
  90     test_skip(-tests => 3,
  91               -requires_modules    => [qw(Bio::DB::GenBank
  92                                           LWP::UserAgent )],
  93               -requires_networking => 1);
  94     my $db_in;
  95     eval {
  96         ok $db_in = Bio::DB::GenBank->new();
  97         my $seq_obj = $db_in->get_Seq_by_id('AF032048.1');
  98     };
  99     if ($@) {
 100         print "$@\n";
 101         skip  "Warning: Problem accessing GenBank entry AF032048.1 "
 102             . "to test spliced_seq on remote DBs", 2;
 103     }
 104
 105     my $len_w_db;
 106     warning_is { $len_w_db = length($feats[1]->spliced_seq(-db => $db_in)->seq) }
 107                [],
 108                "no warnings if GenBank db provided for remote sequence";
 109     ok($len_w_db == 374, "correct length if remote sequence is provided")
 110 }