t/SeqFeature/Generic.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 => 364);
  11
  12     use_ok 'Bio::Seq';
  13     use_ok 'Bio::SeqIO';
  14     use_ok 'Bio::SeqFeature::Generic';
  15 }
  16
  17
  18 my ($feat, $str, $feat2, $pair, @sft);
  19
  20 my $DEBUG = test_debug();
  21
  22 ok $feat = Bio::SeqFeature::Generic->new(
  23     -start => 40,
  24     -end => 80,
  25     -strand => 1,
  26 );
  27 is $feat->primary_tag, '';
  28 is $feat->source_tag, '';
  29 is $feat->display_name, '';
  30
  31 ok $feat = Bio::SeqFeature::Generic->new(
  32     -start => 40,
  33     -end => 80,
  34     -strand => 1,
  35     -primary => 'exon',
  36     -source  => 'internal',
  37     -display_name => 'my exon feature',
  38     -tag => {
  39         silly => 20,
  40         new => 1
  41     }
  42 );
  43
  44 is $feat->start, 40, 'Start of feature location';
  45 is $feat->end, 80, 'End of feature location';
  46 is $feat->primary_tag, 'exon', 'Primary tag';
  47 is $feat->source_tag, 'internal', 'Source tag';
  48 is $feat->display_name, 'my exon feature', 'Display name';
  49 is $feat->phase, undef, 'Undef phase by default';
  50 is $feat->phase(1), 1, 'Phase accessor returns';
  51 is $feat->phase, 1, 'Phase is persistent';
  52
  53 ok $feat->gff_string();
  54
  55 ok $feat2 = Bio::SeqFeature::Generic->new(
  56     -start => 400,
  57     -end => 440,
  58     -strand => 1,
  59     -primary => 'other',
  60     -source  => 'program_a',
  61         -phase => 1,
  62         -tag => {
  63             silly => 20,
  64             new => 1
  65         }
  66 );
  67 is $feat2->phase, 1, 'Set phase from constructor';
  68
  69
  70 # Test attaching a SeqFeature::Generic to a Bio::Seq or SeqFeature::Generic
  71 {
  72     # Make the parent sequence object
  73     my $seq = Bio::Seq->new(
  74         -seq        => 'aaaaggggtttt',
  75         -display_id => 'test',
  76         -alphabet   => 'dna',
  77     );
  78
  79     # Make a SeqFeature
  80     ok my $sf1 = Bio::SeqFeature::Generic->new(
  81         -start  => 4,
  82         -end    => 9,
  83         -strand => 1,
  84     );
  85
  86     # Add the SeqFeature to the parent
  87     ok $seq->add_SeqFeature($sf1);
  88
  89     # Test that it gives the correct sequence
  90     is $sf1->start, 4, 'Start of first seqfeature';
  91     is $sf1->end, 9, 'End of first seqfeature';
  92     is $sf1->strand, 1, 'Strand of first seqfeature';
  93     ok my $sf_seq1 = $sf1->seq;
  94     is $sf_seq1->seq, 'aggggt', 'Sequence of first seqfeature';
  95
  96     # Make a second seqfeature on the opposite strand
  97     ok my $sf2 = Bio::SeqFeature::Generic->new(
  98         -start  => 4,
  99         -end    => 9,
 100         -strand => -1,
 101     );
 102
 103     # Now add the PrimarySeq to the seqfeature before adding it to the parent
 104     ok $sf2->attach_seq($seq->primary_seq);
 105     ok $seq->add_SeqFeature($sf2);
 106
 107     # Test again that we have the correct sequence
 108     is $sf2->start, 4, 'Start of second seqfeature';
 109     is $sf2->end, 9, 'End of second seqfeature';
 110     is $sf2->strand, -1, 'Strand of second seqfeature';
 111     ok my $sf_seq2 = $sf2->seq;
 112     is $sf_seq2->seq, 'acccct', 'Sequence of second seqfeature';
 113 }
 114
 115
 116 # Some tests for bug #947
 117
 118 ok my $sfeat = Bio::SeqFeature::Generic->new(-primary => 'test');
 119
 120 ok $sfeat->add_sub_SeqFeature(
 121     Bio::SeqFeature::Generic->new(
 122         -start => 2,
 123         -end   => 4,
 124         -primary => 'sub1'
 125     ),
 126     'EXPAND'
 127 );
 128
 129 ok $sfeat->add_sub_SeqFeature(
 130     Bio::SeqFeature::Generic->new(
 131         -start => 6,
 132         -end   => 8,
 133         -primary => 'sub2'
 134     ),
 135     'EXPAND'
 136 );
 137
 138 is $sfeat->start, 2, 'sfeat start for EXPAND-ED feature (bug #947)';
 139 is $sfeat->end, 8, 'sfeat end for EXPAND-ED feature (bug #947)';
 140
 141 # Some tests to see if we can set a feature to start at 0
 142 ok $sfeat = Bio::SeqFeature::Generic->new(-start => 0, -end => 0 );
 143
 144 ok defined $sfeat->start;
 145 is $sfeat->start, 0, 'Can create feature starting and ending at 0';
 146 ok defined $sfeat->end;
 147 is $sfeat->end, 0, 'Can create feature starting and ending at 0';
 148
 149
 150 # Test for bug when Locations are not created explicitly
 151
 152 ok my $feat1 = Bio::SeqFeature::Generic->new(
 153     -start => 1,
 154     -end   => 15,
 155     -strand=> 1
 156 );
 157
 158 ok $feat2 = Bio::SeqFeature::Generic->new(
 159     -start => 10,
 160     -end   => 25,
 161     -strand=> 1
 162 );
 163
 164 ok my $overlap = $feat1->location->union($feat2->location);
 165 is $overlap->start, 1;
 166 is $overlap->end,   25;
 167
 168 ok my $intersect = $feat1->location->intersection($feat2->location);
 169 is $intersect->start, 10;
 170 is $intersect->end,   15;
 171
 172
 173 # Now let's test spliced_seq
 174 ok my $seqio = Bio::SeqIO->new(
 175     -file => test_input_file('AY095303S1.gbk'),
 176     -format  => 'genbank'
 177 );
 178 isa_ok $seqio, 'Bio::SeqIO';
 179 ok my $geneseq = $seqio->next_seq;
 180 isa_ok $geneseq, 'Bio::Seq';
 181 ok my ($CDS) = grep { $_->primary_tag eq 'CDS' } $geneseq->get_SeqFeatures;
 182 my $db;
 183
 184 SKIP: {
 185     test_skip(-tests => 5,
 186               -requires_modules => [qw(IO::String
 187                                        LWP::UserAgent
 188                                        HTTP::Request::Common)],
 189               -requires_networking => 1);
 190
 191     use_ok 'Bio::DB::GenBank';
 192
 193     $db = Bio::DB::GenBank->new(-verbose=> $DEBUG);
 194     $CDS->verbose(-1);
 195     my $cdsseq = $CDS->spliced_seq(-db => $db,-nosort => 1);
 196
 197     is $cdsseq->subseq(1,76),
 198        'ATGCAGCCATACGCTTCCGTGAGCGGGCGATGTCTATCTAGACCAGATGCATTGCATGTGATACCGTTTGGGCGAC';
 199     is $cdsseq->translate->subseq(1,100),
 200        'MQPYASVSGRCLSRPDALHVIPFGRPLQAIAGRRFVRCFAKGGQPGDKKKLNVTDKLRLGNTPPTLDVLK'.
 201        'APRPTDAPSAIDDAPSTSGLGLGGGVASPR';
 202     # Test what happens without
 203     $cdsseq = $CDS->spliced_seq(-db => $db,-nosort => 1);
 204     is $cdsseq->subseq(1,76),
 205        'ATGCAGCCATACGCTTCCGTGAGCGGGCGATGTCTATCTAGACCAGATGCATTGCATGTGATACCGTTTGGGCGAC';
 206     is $cdsseq->translate->subseq(1,100),
 207        'MQPYASVSGRCLSRPDALHVIPFGRPLQAIAGRRFVRCFAKGGQPGDKKKLNVTDKLRLGNTPPTLDVLK'.
 208        'APRPTDAPSAIDDAPSTSGLGLGGGVASPR';
 209 }
 210
 211 ok $seqio = Bio::SeqIO->new(
 212     -file => test_input_file('AF032047.gbk'),
 213     -format  => 'genbank'
 214 );
 215 isa_ok $seqio, 'Bio::SeqIO';
 216 ok $geneseq = $seqio->next_seq;
 217 isa_ok $geneseq, 'Bio::Seq';
 218 ok( ($CDS) = grep { $_->primary_tag eq 'CDS' } $geneseq->get_SeqFeatures );
 219 SKIP: {
 220     test_skip(-tests => 2,
 221               -requires_modules => [qw(IO::String
 222                                        LWP::UserAgent
 223                                        HTTP::Request::Common)],
 224               -requires_networking => 1);
 225
 226     my $cdsseq = $CDS->spliced_seq( -db => $db, -nosort => 1);
 227     is $cdsseq->subseq(1,70),
 228        'ATGGCTCGCTTCGTGGTGGTAGCCCTGCTCGCGCTACTCTCTCTGTCTGGCCTGGAGGCTATCCAGCATG';
 229     is $cdsseq->translate->seq,
 230        'MARFVVVALLALLSLSGLEAIQHAPKIQVYSRHPAENGKPNFLNCYVSGFHPSDIEVDLLKNGKKIEKVE'.
 231        'HSDLSFSKDWSFYLLYYTEFTPNEKDEYACRVSHVTFPTPKTVKWDRTM*';
 232 }
 233
 234
 235 # Trans-spliced
 236
 237 ok $seqio = Bio::SeqIO->new(
 238     -format => 'genbank',
 239     -file => test_input_file('NC_001284.gbk')
 240 );
 241 isa_ok $seqio, 'Bio::SeqIO';
 242 ok my $genome = $seqio->next_seq;
 243
 244 for my $cds (grep { $_->primary_tag eq 'CDS' } $genome->get_SeqFeatures) {
 245    ok my $spliced = $cds->spliced_seq(-nosort => 1)->translate->seq;
 246    chop $spliced; # remove stop codon
 247    is $spliced, ($cds->get_tag_values('translation'))[0], 'spliced_seq translation matches expected';
 248 }
 249
 250 # Spliced_seq phase
 251 ok my $seq = Bio::SeqIO->new(
 252     -format => 'fasta',
 253     -file   => test_input_file('sbay_c127.fas')
 254 )->next_seq;
 255
 256 ok my $sf = Bio::SeqFeature::Generic->new(
 257     -verbose => -1,
 258     -start => 263,
 259     -end => 721,
 260     -strand => 1,
 261     -primary => 'splicedgene'
 262 );
 263
 264 ok $sf->attach_seq($seq);
 265
 266 my %phase_check = (
 267     'TTCAATGACT' => 'FNDFYSMGKS',
 268     'TCAATGACTT' => 'SMTSIPWVNQ',
 269     'GTTCAATGAC' => 'VQ*LLFHG*I',
 270 );
 271
 272 for my $phase (-1..3) {
 273     ok my $sfseq = $sf->spliced_seq(-phase => $phase);
 274     ok exists $phase_check{$sfseq->subseq(1,10)};
 275     is $sfseq->translate->subseq(1,10), $phase_check{$sfseq->subseq(1,10)}, 'phase check';
 276 }
 277
 278 # Tags
 279 ok $sf->add_tag_value('note','n1');
 280 ok $sf->add_tag_value('note','n2');
 281 ok $sf->add_tag_value('comment','c1');
 282 is_deeply [sort $sf->get_all_tags()], [sort qw(note comment)] , 'Tags found';
 283 is_deeply [sort $sf->get_tagset_values('note')], [sort qw(n1 n2)] , 'get_tagset_values tag values found';
 284 is_deeply [sort $sf->get_tagset_values(qw(note comment))], [sort qw(c1 n1 n2)] , 'get_tagset_values tag values for multiple tags found';
 285 lives_ok {
 286   is_deeply [sort $sf->get_tag_values('note')], [sort qw(n1 n2)] , 'get_tag_values tag values found';
 287 } 'get_tag_values lives with tag';
 288 lives_ok {
 289   is_deeply [$sf->get_tagset_values('notag') ], [], 'get_tagset_values no tag values found';
 290 } 'get_tagset_values lives with no tag';
 291 throws_ok { $sf->get_tag_values('notag') } qr/tag value that does not exist/, 'get_tag_values throws with no tag';
 292
 293 # Circular sequence SeqFeature tests
 294 $seq = Bio::SeqIO->new(
 295     -format => 'genbank',
 296     -file   => test_input_file('PX1CG.gb')
 297 )->next_seq;
 298
 299 ok $seq->is_circular, 'Phi-X174 genome is circular';
 300
 301 # Retrieving the spliced sequence from any split location requires spliced_seq()
 302
 303 my %sf_data = (
 304     #       start
 305     'A'  => [3981, 136, 1, 1542, 'join(3981..5386,1..136)', 'ATGGTTCGTT'],
 306     'A*' => [4497, 136, 1, 1026, 'join(4497..5386,1..136)', 'ATGAAATCGC'],
 307     'B'  => [5075, 51,  1, 363,  'join(5075..5386,1..51)',  'ATGGAACAAC'],
 308 );
 309
 310 ok my @split_sfs = grep {
 311     $_->location->isa('Bio::Location::SplitLocationI')
 312     } $seq->get_SeqFeatures();
 313
 314 is @split_sfs, 3, 'Only 3 split locations';
 315
 316 for my $sf (@split_sfs) {
 317     isa_ok $sf->location, 'Bio::Location::SplitLocationI';
 318     ok my ($tag) = $sf->get_tag_values('product');
 319     my ($start, $end, $strand, $length, $ftstring, $first_ten) = @{$sf_data{$tag}};
 320
 321     is $sf->location->to_FTstring, $ftstring, 'Feature string';
 322     is $sf->spliced_seq->subseq(1,10), $first_ten, 'First ten nucleotides';
 323     is $sf->strand, $strand, 'Strand';
 324     is $sf->start, $start, 'Start';
 325     is $sf->end, $end, 'End';
 326     is $sf->length, $length, 'Expected length';
 327 }
 328
 329 # spliced_seq() on the reverse strand, bug #88 (github)
 330 $seq = Bio::SeqIO->new( -file => test_input_file('AF222649-rc.gbk') )->next_seq;
 331 # All should start with "ATG"
 332 for my $feat ( $seq->get_SeqFeatures('CDS') ) {
 333     ok $feat->spliced_seq->seq =~ /^ATG/, "Reverse strand is spliced correctly";
 334 }
 335
 336