t/Seq/PrimarySeq.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5 use Data::Dumper;
   6
   7 BEGIN {
   8     use lib '.';
   9     use Bio::Root::Test;
  10
  11     test_begin( -tests => 155 );
  12
  13     use_ok('Bio::PrimarySeq');
  14     use_ok('Bio::Location::Simple');
  15     use_ok('Bio::Location::Fuzzy');
  16     use_ok('Bio::Location::Split');
  17 }
  18
  19
  20 # Bare object
  21 ok my $seq = Bio::PrimarySeq->new(), 'Bare object';
  22 isa_ok $seq, 'Bio::PrimarySeqI';
  23 is $seq->id, undef;
  24 is $seq->seq, undef;
  25 is $seq->length, 0;
  26 is $seq->alphabet, undef;
  27
  28
  29 # Empty sequence
  30 ok $seq = Bio::PrimarySeq->new( -seq => '');
  31 is $seq->seq, '';
  32 is $seq->length, 0;
  33 is $seq->alphabet, undef;
  34
  35
  36 # Basic tests
  37 ok $seq = Bio::PrimarySeq->new(
  38     '-seq'              => 'TTGGTGGCGTCAACT',
  39     '-display_id'       => 'new-id',
  40     '-alphabet'         => 'dna',
  41     '-accession_number' => 'X677667',
  42     '-desc'             => 'Sample Bio::Seq object'
  43 );
  44 ok defined $seq;
  45 is $seq->accession_number(), 'X677667';
  46 is $seq->seq(),              'TTGGTGGCGTCAACT';
  47 is $seq->display_id(),       'new-id';
  48 is $seq->alphabet(),         'dna';
  49 is $seq->is_circular(),      undef;
  50 ok $seq->is_circular(1);
  51 is $seq->is_circular(0), 0;
  52
  53 # check IdentifiableI and DescribableI interfaces
  54 isa_ok $seq, 'Bio::IdentifiableI';
  55 isa_ok $seq, 'Bio::DescribableI';
  56
  57 # make sure all methods are implemented
  58 is $seq->authority("bioperl.org"), "bioperl.org";
  59 is $seq->namespace("t"),           "t";
  60 is $seq->namespace, "t";
  61 is $seq->version(0), 0;
  62 is $seq->lsid_string(),      "bioperl.org:t:X677667";
  63 is $seq->namespace_string(), "t:X677667.0";
  64 $seq->version(47);
  65 is $seq->version, 47;
  66 is $seq->namespace_string(), "t:X677667.47";
  67 is $seq->description(),      'Sample Bio::Seq object';
  68 is $seq->display_name(),     "new-id";
  69
  70 my $location = Bio::Location::Simple->new(
  71     '-start'  => 2,
  72     '-end'    => 5,
  73     '-strand' => -1
  74 );
  75 is( $seq->subseq($location), 'ACCA' );
  76
  77 my $splitlocation = Bio::Location::Split->new();
  78 $splitlocation->add_sub_Location(
  79     Bio::Location::Simple->new(
  80         '-start'  => 1,
  81         '-end'    => 4,
  82         '-strand' => 1
  83     )
  84 );
  85
  86 $splitlocation->add_sub_Location(
  87     Bio::Location::Simple->new(
  88         '-start'  => 7,
  89         '-end'    => 12,
  90         '-strand' => -1
  91     )
  92 );
  93
  94 is( $seq->subseq($splitlocation), 'TTGGTGACGC' );
  95
  96 my $fuzzy = Bio::Location::Fuzzy->new(
  97     -start  => '<3',
  98     -end    => '8',
  99     -strand => 1
 100 );
 101
 102 is( $seq->subseq($fuzzy), 'GGTGGC' );
 103
 104 my $trunc = $seq->trunc( 1, 4 );
 105 isa_ok $trunc, 'Bio::PrimarySeqI';
 106 is $trunc->seq(), 'TTGG' or diag( "Expecting TTGG. Got " . $trunc->seq() );
 107
 108 $trunc = $seq->trunc($splitlocation);
 109 isa_ok( $trunc, 'Bio::PrimarySeqI' );
 110 is( $trunc->seq(), 'TTGGTGACGC' );
 111
 112 $trunc = $seq->trunc($fuzzy);
 113 isa_ok( $trunc, 'Bio::PrimarySeqI' );
 114 is( $trunc->seq(), 'GGTGGC' );
 115
 116 my $rev = $seq->revcom();
 117 isa_ok( $rev, 'Bio::PrimarySeqI' );
 118
 119 is $rev->seq(), 'AGTTGACGCCACCAA'
 120   or diag( 'revcom() failed, was ' . $rev->seq() );
 121
 122 is $rev->display_id, 'new-id';
 123 is( $rev->alphabet(),    'dna', 'alphabet copied through revcom' );
 124 TODO: {
 125     local $TODO =
 126       'all attributes of primaryseqs are not currently copied through revcoms';
 127     is( $rev->namespace, 't', 'namespace copied through revcom' );
 128     is( $rev->namespace_string(),
 129         "t:X677667.47", 'namespace_string copied through revcom' );
 130     is( $rev->is_circular(), 0,     'is_circular copied through revcom' );
 131 }
 132
 133 #
 134 # Translate
 135 #
 136
 137 my $aa = $seq->translate();    # TTG GTG GCG TCA ACT
 138 is $aa->seq, 'LVAST', "Translation: " . $aa->seq;
 139
 140 # tests for non-standard initiator codon coding for
 141 # M by making translate() look for an initiator codon and
 142 # terminator codon ("complete", the 5th argument below)
 143 $seq->seq('TTGGTGGCGTCAACTTAA');    # TTG GTG GCG TCA ACT TAA
 144 $aa = $seq->translate( undef, undef, undef, undef, 1 );
 145 is $aa->seq, 'MVAST', "Translation: " . $aa->seq;
 146
 147 # same test as previous, but using named parameter
 148 $aa = $seq->translate( -complete => 1 );
 149 is $aa->seq, 'MVAST', "Translation: " . $aa->seq;
 150
 151 # find ORF, ignore codons outside the ORF or CDS
 152 $seq->seq('TTTTATGGTGGCGTCAACTTAATTT');    # ATG GTG GCG TCA ACT
 153 $aa = $seq->translate( -orf => 1 );
 154 is $aa->seq, 'MVAST*', "Translation: " . $aa->seq;
 155
 156 # smallest possible ORF
 157 $seq->seq("ggggggatgtagcccc");             # atg tga
 158 $aa = $seq->translate( -orf => 1 );
 159 is $aa->seq, 'M*', "Translation: " . $aa->seq;
 160
 161 # same as previous but complete, so * is removed
 162 $aa = $seq->translate(
 163     -orf      => 1,
 164     -complete => 1
 165 );
 166 is $aa->seq, 'M', "Translation: " . $aa->seq;
 167
 168 # ORF without termination codon
 169 # should warn, let's change it into throw for testing
 170 $seq->verbose(2);
 171 $seq->seq("ggggggatgtggcccc");    # atg tgg ccc
 172 eval { $seq->translate( -orf => 1 ); };
 173 like( $@, qr/\batgtggccc\b/i );
 174 $seq->verbose(-1);
 175 $aa = $seq->translate( -orf => 1 );
 176 is $aa->seq, 'MWP', "Translation: MWP";
 177 $seq->verbose(0);
 178
 179 # use non-standard codon table where terminator is read as Q
 180 $seq->seq('ATGGTGGCGTCAACTTAG');    # ATG GTG GCG TCA ACT TAG
 181 $aa = $seq->translate( -codontable_id => 6 );
 182 is $aa->seq, 'MVASTQ' or diag( "Translation: " . $aa->seq );
 183
 184 # insert an odd character instead of terminating with *
 185 $aa = $seq->translate( -terminator => 'X' );
 186 is $aa->seq, 'MVASTX' or diag( "Translation: " . $aa->seq );
 187
 188 # change frame from default
 189 $aa = $seq->translate( -frame => 1 );    # TGG TGG CGT CAA CTT AG
 190 is $aa->seq, 'WWRQL' or diag( "Translation: " . $aa->seq );
 191
 192 $aa = $seq->translate( -frame => 2 );    # GGT GGC GTC AAC TTA G
 193 is $aa->seq, 'GGVNL' or diag( "Translation: " . $aa->seq );
 194
 195 # TTG is initiator in Standard codon table? Afraid so.
 196 $seq->seq("ggggggttgtagcccc");           # ttg tag
 197 $aa = $seq->translate( -orf => 1 );
 198 is $aa->seq, 'L*' or diag( "Translation: " . $aa->seq );
 199
 200 # Replace L at 1st position with M by setting complete to 1
 201 $seq->seq("ggggggttgtagcccc");           # ttg tag
 202 $aa = $seq->translate(
 203     -orf      => 1,
 204     -complete => 1
 205 );
 206 is $aa->seq, 'M' or diag( "Translation: " . $aa->seq );
 207
 208 # Ignore non-ATG initiators (e.g. TTG) in codon table
 209 $seq->seq("ggggggttgatgtagcccc");        # atg tag
 210 $aa = $seq->translate(
 211     -orf      => 1,
 212     -start    => "atg",
 213     -complete => 1
 214 );
 215 is $aa->seq, 'M' or diag( "Translation: " . $aa->seq );
 216
 217 # test for character '?' in the sequence string
 218 is $seq->seq('TTGGTGGCG?CAACT'), 'TTGGTGGCG?CAACT';
 219
 220 # test for some aliases
 221 $seq = Bio::PrimarySeq->new(
 222     -id          => 'aliasid',
 223     -description => 'Alias desc'
 224 );
 225 is( $seq->description, 'Alias desc' );
 226 is( $seq->display_id,  'aliasid' );
 227
 228 # Test alphabet
 229
 230 ok $seq->seq('actgx');
 231 is $seq->alphabet, 'protein', 'Alphabet';
 232 ok $seq->seq('actge');
 233 is $seq->alphabet, 'protein';
 234 ok $seq->seq('actgf');
 235 is $seq->alphabet, 'protein';
 236 ok $seq->seq('actgi');
 237 is $seq->alphabet, 'protein';
 238 ok $seq->seq('actgj');
 239 is $seq->alphabet, 'protein';
 240 ok $seq->seq('actgl');
 241 is $seq->alphabet, 'protein';
 242 ok $seq->seq('actgo');
 243 is $seq->alphabet, 'protein';
 244 ok $seq->seq('actgp');
 245 is $seq->alphabet, 'protein';
 246 ok $seq->seq('actgq');
 247 is $seq->alphabet, 'protein';
 248 ok $seq->seq('actgz');
 249 is $seq->alphabet, 'protein';
 250 ok $seq->seq('actgn');
 251 is $seq->alphabet, 'dna';
 252 ok $seq->seq('acugn');
 253 is $seq->alphabet, 'rna';
 254 ok $seq->seq('bdhkm');
 255 is $seq->alphabet, 'protein';
 256 ok $seq->seq('rsvwx');
 257 is $seq->alphabet, 'protein';
 258 ok $seq->seq('AAACTYAAAAGAATTGRCGG'); # valid degenerate DNA PCR primer sequence (90% ACGTN)
 259 is $seq->alphabet, 'dna';
 260 ok $seq->seq('AAACTYAAAKGAATTGRCGG'); # another primer previously detected as protein (85% ACGTN)
 261 is $seq->alphabet, 'dna';
 262 ok $seq->seq('YWACTYAAAKGARTTGRCGG'); # 70% ACGTN. Everything <= 70% is considered a protein
 263 is $seq->alphabet, 'protein';
 264 ok $seq->seq('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'); # Bug 2438
 265 is $seq->alphabet, 'protein', 'Bug 2438';
 266 ok $seq->seq('CAGTCXXXXXXXXXXXXXXXXXXXXXXXXXXXCAGCG');
 267 is $seq->alphabet, 'protein';
 268
 269 ok $seq->seq('actgn', 'protein'); # accept specified alphabet, no matter what
 270 is $seq->alphabet, 'protein';
 271 ok $seq->seq('bdhkm', 'dna');
 272 is $seq->alphabet, 'dna';
 273
 274
 275 # Bug #2864:
 276
 277 $seq = Bio::PrimarySeq->new( -display_id => 0, -seq => 'GATC' );
 278
 279 is $seq->display_id, 0, "Bug #2864";
 280
 281 # Test that the check for terminators inside the translated protein
 282 # works when the terminator isn't '*':
 283
 284 $seq = Bio::PrimarySeq->new(-seq=>'ATGCTCTAAGCAGGGTAA'); # ML*AG*
 285 eval { $aa = $seq->translate(-complete=>1, -throw=>1, -terminator=>'#') };
 286 my $error = $@;
 287 ok $error =~ /\QTerminator codon inside CDS!\E/, 'Terminator + inside sequence';
 288
 289 $seq = Bio::PrimarySeq->new(-seq=>'ATGCTCGCAGGGTAA'); # MLAG*
 290 $aa = $seq->translate(-complete=>1, -throw=>1, -terminator=>'#');
 291 is $aa->seq, 'MLAG';
 292
 293
 294 # Test length method
 295 ok $seq = Bio::PrimarySeq->new(), 'Length method';
 296 is $seq->length, 0;
 297 ok $seq->length(123);
 298 is $seq->length, 123;
 299
 300 ok $seq = Bio::PrimarySeq->new( -seq => 'ATGCTCTAAGCAGGGTAA' );
 301 is $seq->length, 18;
 302 ok $seq->seq('ATGCTCTAAG');
 303 is $seq->length, 10;
 304 is $seq->seq(undef), undef;
 305 is $seq->length, 0;
 306
 307 ok $seq = Bio::PrimarySeq->new( -length => 123 );
 308 is $seq->length, 123;
 309
 310 ok $seq = Bio::PrimarySeq->new( -seq => 'ATGCTCTAAGCAGGGTAA' );
 311 is $seq->length, 18;
 312 ok $seq->length( $seq->length ); # save memory by removing seq
 313 is $seq->seq( undef ), undef;    # ... but keeping a record of length
 314 is $seq->length, 18;
 315 is $seq->seq, undef;
 316 ok $seq->seq('ACGT');
 317 is $seq->length, 4; # manually-specified length changed when sequence is changed
 318
 319 throws_ok { $seq->length(666); } qr/.+/; # Cannot lie about length
 320
 321
 322 # Sequence validation method
 323 is $seq->validate_seq( undef    ), 1;
 324 is $seq->validate_seq( ''       ), 1;
 325 is $seq->validate_seq( 'acgt'   ), 1;
 326 is $seq->validate_seq( 'ACGT'   ), 1;
 327 is $seq->validate_seq( 'XFRH'   ), 1;
 328 is $seq->validate_seq( '-~'     ), 1; # gap symbols
 329 is $seq->validate_seq( '-.*?=~' ), 1; # other valid symbols
 330 is $seq->validate_seq( 'AAAA$'  ), 0;
 331 is $seq->validate_seq( 'tt&tt'  ), 0;
 332
 333
 334 # Test direct option (no sequence validation)
 335 throws_ok { $seq = Bio::PrimarySeq->new(-seq => 'A\T$AGQ+T'); } qr/.+/, 'Validation';
 336 ok $seq = Bio::PrimarySeq->new( -seq => 'A\T$AGQ+T', -direct => 1 );
 337 is $seq->seq, 'A\T$AGQ+T';
 338
 339
 340 # Set a sequence by reference
 341 my $string = 'AAAACCCCGGGGTTTT';
 342 ok $seq = Bio::PrimarySeq->new( -ref_to_seq => \$string );
 343 is $seq->seq, 'AAAACCCCGGGGTTTT';
 344
 345
 346 # Test internal PrimarySeqI _find_orfs function and translate( -orf => 'longest' )
 347 {
 348     my @tests = (
 349         #tiny test
 350         ['TTTTATGGTGGCGTCAACTTAATTT',
 351          [[4,22,18,1]],
 352         ],
 353
 354         #bigger test (this is a tomato unigene)
 355         ['GAAGGCTGGTTCTGAGTTGGATCTATGTTTGATGAAGGGAAGTAGACCGGAGGTCTTGCATCAGCAATATTAGTACCAAATCCAGGTGGAGGCGCATCCTGTCTCCGTTGCATTTCAACTTTCATTTCAGCAATCTGTTGCATCAGTTGCATGATCAATTCATTCTGTTCCACTACAGTGGGCTGAGCGACCACAACGTCAGTAAGACGCCCTTCGTCATTGTTGTCTCCCATAACTGTTTTTCCTTTATCTGAATTTGATCGAGGGAAGGAATCTGTAGGACCTTTCGATCTGGTGAAGTAAGGATGATCTGCCAGCTTTATTGACACAGATCAGTAAAAAGGTACCTGAAAGGTAAAAACAACTCAAAGGCAAATTTGTTAGTGCATATCCAGAGTACAAAATGCTTAATATCGCACATAAAACCGATAAACACACAAGTCGTTTTGTTTGAGGATATCTTAACCCACGAATAAGGACGGATATATATTTTGAACAAACAGGAATTTGTTTGTTTGGCGTTATCTTGGGAAATCTG',
 356          [[98,254,156,2],[347,476,129,2],[219,303,84,0],[16,73,57,1],[403,454,51,1],[310,358,48,1],[235,280,45,1],[491,536,45,2],[150,186,36,0],[507,537,30,0],[5,32,27,2],[511,538,27,1],[24,45,21,0],[305,326,21,2],[450,465,15,0]],
 357         ],
 358
 359
 360        );
 361     foreach my $test (@tests) {
 362         my ($test_seq, $orfs) = @$test;
 363         my @orfs = Bio::PrimarySeqI::_find_orfs_nucleotide(
 364             undef,
 365             $test_seq,
 366             Bio::Tools::CodonTable->new,
 367             undef,
 368            ); # ATG GTG GCG TCA ACT
 369         is_deeply( \@orfs, $orfs, '_find_orfs 1')
 370             or diag "for $test_seq, _find_orfs returned:\n"
 371                     .Dumper([map [@$_], @orfs]);
 372
 373         is_deeply( $orfs->[0],
 374                    (sort {$b->[2] <=> $a->[2]} @$orfs)[0],
 375                    'orfs are sorted by descending length'
 376                   );
 377
 378         # make sure we get the same sequence by taking the longest orf
 379         # nucleotide from the test data and translating it, as by
 380         # calling translate with -orf => 'longest'
 381         is(
 382             Bio::PrimarySeq
 383               ->new( -seq => $test_seq, -id => 'fake_id' )
 384               ->translate( -orf => 'longest' )
 385               ->seq,
 386
 387             Bio::PrimarySeq
 388               ->new( -seq => substr( $test_seq, $orfs->[0][0], $orfs->[0][2] ),
 389                      -id => 'foo'
 390                     )
 391               ->translate
 392               ->seq,
 393             'got correct -orf => "longest" seq',
 394            );
 395     }
 396 }