Bio::Tools::CodonTable and Bio::Tools::IUPAC: prepare with dzil.
[bioperl-live.git] / t / Seq / Seq.t
blobcaa3fa196b581f0805474f02713ebca893c003e7
1 # -*-Perl-*- Test Harness script for Bioperl
2 # $Id$
4 use strict;
6 BEGIN {
7     use Bio::Root::Test;
9     test_begin(-tests => 76);
11     use_ok('Bio::Seq');
12     use_ok('Bio::Seq::RichSeq');
13     use_ok('Bio::SeqFeature::Generic');
14     use_ok('Bio::Species');
15     use_ok('Bio::Annotation::SimpleValue');
18 ok my $seq = Bio::Seq->new(-seq=>'ACTGTGGCGTCAACT',
19                         -desc=>'Sample Bio::Seq object',
20                         -alphabet => 'dna',
21                         -is_circular => 1
22                        );
23 isa_ok($seq,"Bio::AnnotatableI");
25 ok $seq->is_circular;
26 ok not $seq->is_circular(0);
27 ok not $seq->is_circular;
29 my $trunc = $seq->trunc(1,4);
30 is $trunc->length, 4, 'truncated sequence length';
32 is $trunc->seq, 'ACTG', 'truncated sequence string';
34 # test ability to get str function
35 is $seq->seq(),  'ACTGTGGCGTCAACT' ;
37 ok $seq = Bio::Seq->new(-seq=>'actgtggcgtcaact',
38                         -desc=>'Sample Bio::Seq object',
39                         -display_id => 'something',
40                         -accession_number => 'accnum',
41                         -alphabet => 'dna' );
43 is uc $seq->alphabet, 'DNA' , 'alphabet';
45 # basic methods
47 is $seq->id(), 'something',  "id";
48 is $seq->accession_number, 'accnum', "accession number";
49 is $seq->subseq(5, 9),  'tggcg', "subseq";
51 # check IdentifiableI and DescribableI interfaces
52 isa_ok $seq, 'Bio::IdentifiableI';
53 isa_ok $seq, 'Bio::DescribableI';
54 # make sure all methods are implemented
55 is $seq->authority("bioperl.org"), "bioperl.org";
56 is $seq->namespace("t"), "t";
57 is $seq->version(0), 0;
58 is $seq->lsid_string(), "bioperl.org:t:accnum";
59 is $seq->namespace_string(), "t:accnum.0";
60 is $seq->description(), 'Sample Bio::Seq object';
61 is $seq->display_name(), "something";
63 # check that feature accession works regardless of lazy things going on
64 is scalar($seq->top_SeqFeatures()), 0;
65 is scalar($seq->flush_SeqFeatures()), 0;
67 my $newfeat = Bio::SeqFeature::Generic->new( -start => 10,
68                                              -end => 12,
69                                              -primary => 'silly',
70                                              -source => 'stuff');
73 $seq->add_SeqFeature($newfeat);
74 is $seq->feature_count, 1;
76 my $species = Bio::Species->new
77     (-verbose => 1,
78      -classification => [ qw( sapiens Homo Hominidae
79                               Catarrhini Primates Eutheria
80                               Mammalia Vertebrata Chordata
81                               Metazoa Eukaryota )]);
82 $seq->species($species);
83 is $seq->species->binomial, 'Homo sapiens';
84 $seq->annotation->add_Annotation('description',
85                  Bio::Annotation::SimpleValue->new(-value => 'desc-here'));
86 my ($descr) = $seq->annotation->get_Annotations('description');
87 is $descr->value(), 'desc-here';
88 is $descr->tagname(), 'description';
91 #  translation tests
94 my $trans = $seq->translate();
95 is  $trans->seq(), 'TVAST' , 'translated sequence';
97 # unambiguous two character codons like 'ACN' and 'GTN' should give out an amino
98 # acid ...with the addendum that there should be no assumption by the method
99 # to complete the codon unless specified, using the -complete_codons flag.
101 $seq->seq('ACTGTGGCGTCAACN');
102 $trans = $seq->translate();
103 is $trans->seq(), 'TVAST', 'translated sequence with explicit unambiguous codons';
105 $seq->seq('ACTGTGGCGTCAAC');
106 $trans = $seq->translate();
107 is $trans->seq(), 'TVAS', 'translated sequence with unknown unambiguous codons';
109 $seq->seq('ACTGTGGCGTCAAC');
110 $trans = $seq->translate(-complete_codons => 1);
111 is $trans->seq(), 'TVAST', 'translated sequence with unknown unambiguous codons, completed';
113 $seq->seq('ACTGTGGCGTCAACA');
114 $trans = $seq->translate();
115 is $trans->seq(), 'TVAST', 'translated sequence with unambiguous codons';
117 $seq->seq('ACTGTGGCGTCAACAG');
118 $trans = $seq->translate();
119 is $trans->seq(), 'TVAST', 'translated sequence with unambiguous codons';
121 $seq->seq('ACTGTGGCGTCAACAGT');
122 $trans = $seq->translate(-complete_codons => 1);
123 is $trans->seq(), 'TVASTV', 'translated sequence with unknown unambiguous codons, completed';
125 $seq->seq('ACTGTGGCGTCAACAGTA');
126 $trans = $seq->translate();
127 is $trans->seq(), 'TVASTV', 'translated sequence with unambiguous codons';
129 $seq->seq('AC');
130 is $seq->translate(-complete_codons => 1)->seq , 'T', 'translated sequence with unknown unambiguous codons, completed';
132 #difference between the default and full CDS translation
134 $seq->seq('atgtggtaa');
135 $trans = $seq->translate();
136 is $trans->seq(), 'MW*' , 'translated sequence with stop';
138 $seq->seq('atgtggtaa');
139 $trans = $seq->translate(undef,undef,undef,undef,1);
140 is $trans->seq(), 'MW', 'translated sequence';
142 #frame
143 my $string;
144 my @frames = (0, 1, 2);
145 foreach my $frame (@frames) {
146     $string .= $seq->translate(undef, undef, $frame)->seq;
147     $string .= $seq->revcom->translate(undef, undef, $frame)->seq;
149 is $string, 'MW*LPHCGYHVVTT';
151 #Translating with all codon tables using method defaults
152 $string = '';
153 my @codontables = qw(0 1 2 3 4 5 6 9 10 11 12 13 14 16 21
154     22 23 24 25 26 27 28 29 30 31);
155 foreach my $ct (@codontables) {
156     $string .= $seq->translate(undef, undef, undef, $ct)->seq;
158 is $string, 'MW*MW*MW*MW*MW*MW*MWQMW*MW*MW*MW*MW*MWYMW*MW*MW*MW*MW*MW*MW*MWQMWQMWYMWEMWE';
160 # CDS translation set to throw an exception for internal stop codons
161 $seq->seq('atgtggtaataa');
162 eval {
163     $seq->translate(undef, undef, undef, undef, 'CDS' , 'throw');
165 like ($@, qr/EX/);
167 $seq->seq('atgtggtaataa');
168 is( $seq->translate('J', '-',)->seq, 'MWJJ');
170 # tests for RichSeq
171 ok my $richseq = Bio::Seq::RichSeq->new( -seq => 'atgtggtaataa',
172                                       -accession_number => 'AC123',
173                                       -alphabet => 'rna',
174                                       -molecule => 'mRNA',
175                                       -id => 'id1',
176                                       -dates => [ '2001/1/1' ],
177                                       -pid => '887821',
178                                       -keywords => 'JUNK1;JUNK2',
179                                       -division => 'Fungi',
180                                       -secondary_accessions => 'AC1152' );
182 is ($richseq->seq, 'atgtggtaataa');
183 is ($richseq->display_id, 'id1');
184 is (($richseq->get_dates)[0], '2001/1/1');
185 is (($richseq->get_secondary_accessions)[0], 'AC1152');
186 is ($richseq->accession_number, 'AC123');
187 is ($richseq->alphabet, 'rna');
188 is ($richseq->molecule, 'mRNA');
189 is ($richseq->pid, 887821);
190 is ($richseq->division, 'Fungi');
191 is ($richseq->keywords, 'JUNK1; JUNK2');
192 $richseq->seq_version('2');
193 is ($richseq->seq_version, 2);
195 # Test adding a feature to a RichSeq type, then
196 # trunc() and see if the feature vanishes (we shouldn't
197 # be using clone() for RichSeq types)
198 $richseq->add_SeqFeature($newfeat);
199 is $richseq->feature_count, 1;
200 my $newrichseq = $richseq->trunc(1,5);
201 is $newrichseq->feature_count, 0, "Don't use clone for trunc of Bio::Seq::RichSeq";
202 is $newrichseq->length, 5;
204 # tests for subtle misbehaviors
205 $seq = Bio::Seq->new(-primary_id => 'blah', -accession_number => 'foo');
206 is ($seq->accession_number, $seq->primary_seq->accession_number);
207 is ($seq->primary_id, $seq->primary_seq->primary_id);
208 $seq->accession_number('blurb');
209 $seq->primary_id('bar');
210 is ($seq->accession_number, $seq->primary_seq->accession_number);
211 is ($seq->primary_id, $seq->primary_seq->primary_id);
214 # Bug #2864:
216 $seq = Bio::Seq->new(-display_id => 0, -seq => 'GATC');
218 is $seq->display_id, 0, "Bug #2864";
220 # transcribe/rev_transcribe
222 $seq = Bio::Seq->new( -id => 'seq1', -alphabet=>'dna',
223                       -seq=> 'attTcgcatgT' );
224 ok my $xseq = $seq->transcribe;
225 is $xseq->alphabet, 'rna';
226 ok !($xseq->seq =~ /[tT]/);
227 is $xseq->seq, 'auuUcgcaugU';
228 ok !$xseq->transcribe;
229 ok $seq = $xseq->rev_transcribe;
230 is $seq->seq, 'attTcgcatgT';
231 is $seq->alphabet, 'dna';