t/SeqIO/Splicedseq.t

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