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