fix #105 - find the real "first orf"
[bioperl-live.git] / t / Seq / PrimarySeq.t
blob93d430e9ba9273123a7eb779739ca22c982f4ec4
1 # -*-Perl-*- Test Harness script for Bioperl
2 # $Id$
4 use strict;
5 use Data::Dumper;
7 BEGIN {
8     use lib '.';
9     use Bio::Root::Test;
10     test_begin( -tests => 312 );
12     use_ok('Bio::PrimarySeq');
13     use_ok('Bio::Location::Simple');
14     use_ok('Bio::Location::Fuzzy');
15     use_ok('Bio::Location::Split');
19 # Bare object
20 ok my $seq = Bio::PrimarySeq->new(), 'Bare object';
21 isa_ok $seq, 'Bio::PrimarySeqI';
22 is $seq->id, undef;
23 is $seq->seq, undef;
24 is $seq->length, 0;
25 is $seq->alphabet, undef;
26 is $seq->is_circular, undef;
29 # Empty sequence
30 ok $seq = Bio::PrimarySeq->new( -seq => '', -nowarnonempty => 1);
31 is $seq->seq, '';
32 is $seq->length, 0;
33 is $seq->alphabet, undef;
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'
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;
53 # check IdentifiableI and DescribableI interfaces
54 isa_ok $seq, 'Bio::IdentifiableI';
55 isa_ok $seq, 'Bio::DescribableI';
57 # make sure all methods are implemented
58 is $seq->authority("bioperl.org"), "bioperl.org";
59 is $seq->authority, "bioperl.org";
60 is $seq->namespace("t"), "t";
61 is $seq->namespace, "t";
62 is $seq->version(0), 0;
63 is $seq->version, 0;
64 is $seq->lsid_string(), "bioperl.org:t:X677667";
65 is $seq->namespace_string, "t:X677667.0";
66 is $seq->version(47), 47;
67 is $seq->version, 47;
68 is $seq->namespace_string, "t:X677667.47";
69 is $seq->description, 'Sample Bio::Seq object';
70 is $seq->display_name, "new-id";
73 # Test subseq
74 is $seq->subseq(2, 5), 'TGGT';
76 is $seq->subseq( -start => 1, -end => 15), 'TTGGTGGCGTCAACT';
78 my $location = Bio::Location::Simple->new(
79     '-start'  => 2,
80     '-end'    => 5,
81     '-strand' => -1
83 is $seq->subseq($location), 'ACCA';
85 my $splitlocation = Bio::Location::Split->new();
86 $splitlocation->add_sub_Location(
87     Bio::Location::Simple->new(
88         '-start'  => 1,
89         '-end'    => 4,
90         '-strand' => 1
91     )
94 $splitlocation->add_sub_Location(
95     Bio::Location::Simple->new(
96         '-start'  => 7,
97         '-end'    => 12,
98         '-strand' => -1
99     )
102 is $seq->subseq($splitlocation), 'TTGGTGACGC';
104 my $fuzzy = Bio::Location::Fuzzy->new(
105     -start  => '<3',
106     -end    => '8',
107     -strand => 1
110 is $seq->subseq($fuzzy), 'GGTGGC';
113     ok my $seq = Bio::PrimarySeq->new( -seq => 'TT-GTGGCGTCAACT' );
114     is $seq->subseq(2, 5, 'nogap'), 'TGT';
115     is $seq->subseq( -start => 2, -end => 5, -nogap => 1 ), 'TGT';
116     my $location = Bio::Location::Simple->new(
117        '-start'  => 2,
118        '-end'    => 5,
119        '-strand' => 1
120     );
121     is $seq->subseq( $location, -nogap => 1), 'TGT';
123     is $seq->subseq(-start=>2, -end=>5, -replace_with=>'aa'), 'T-GT';
124     is $seq->seq, 'TaaGGCGTCAACT';
126     throws_ok { $seq->subseq(-start=>2, -end=>5, -replace_with=>'?!'); } qr/.+/;
130     ok my $seq = Bio::PrimarySeq->new( -seq => 'AACCGGTT', -is_circular => 1 );
131     is $seq->subseq( -start => 7, -end => 10 ), 'TTAA';
134 ### Test for Bug #2936
135 # Without strand input argument (case: user don't think is necessary)
136 my $split_loc_obj1 = Bio::Location::Split->new();
137 $split_loc_obj1->add_sub_Location(
138     Bio::Location::Simple->new(
139         '-start'  => 1,
140         '-end'    => 10
141     )
143 $split_loc_obj1->add_sub_Location(
144     Bio::Location::Simple->new(
145         '-start'  => 20,
146         '-end'    => 30
147     )
149 # With strand input argument (case: user provides the argument)
150 my $split_loc_obj2 = Bio::Location::Split->new();
151 $split_loc_obj2->add_sub_Location(
152     Bio::Location::Simple->new(
153         '-start'  => 1,
154         '-end'    => 10,
155         '-strand' => 1
156     )
158 $split_loc_obj2->add_sub_Location(
159     Bio::Location::Simple->new(
160         '-start'  => 20,
161         '-end'    => 30,
162         '-strand' => 1
163     )
165 is $split_loc_obj1->to_FTstring, "join(1..10,20..30)";
166 is $split_loc_obj2->to_FTstring, "join(1..10,20..30)";
167 $split_loc_obj1->flip_strand;
168 $split_loc_obj2->flip_strand;
169 is $split_loc_obj1->to_FTstring, "complement(join(1..10,20..30))";
170 is $split_loc_obj2->to_FTstring, "complement(join(1..10,20..30))";
173 # Test trunc
174 my $trunc = $seq->trunc( 1, 4 );
175 isa_ok $trunc, 'Bio::PrimarySeqI';
176 is $trunc->seq(), 'TTGG' or diag( "Expecting TTGG. Got " . $trunc->seq() );
178 $trunc = $seq->trunc($splitlocation);
179 isa_ok $trunc, 'Bio::PrimarySeqI' ;
180 is $trunc->seq(), 'TTGGTGACGC';
182 $trunc = $seq->trunc($fuzzy);
183 isa_ok $trunc, 'Bio::PrimarySeqI';
184 is $trunc->seq(), 'GGTGGC';
186 my $rev = $seq->revcom();
187 isa_ok $rev, 'Bio::PrimarySeqI';
189 is $rev->seq(), 'AGTTGACGCCACCAA'
190   or diag( 'revcom() failed, was ' . $rev->seq() );
192 is $rev->display_id,         'new-id';
193 is $rev->display_name(),     'new-id';
194 is $rev->accession_number(), 'X677667';
195 is $rev->alphabet,           'dna';
196 is $rev->description,        'Sample Bio::Seq object';
197 is $rev->is_circular(),      0;
198 is $rev->version,            47;
199 is $rev->authority,          'bioperl.org';
200 is $rev->namespace,          't';
201 is $rev->namespace_string(), 't:X677667.47';
204 # Translate
207 my $aa = $seq->translate();    # TTG GTG GCG TCA ACT
208 is $aa->seq, 'LVAST', "Translation: " . $aa->seq;
210 # tests for non-standard initiator codon coding for
211 # M by making translate() look for an initiator codon and
212 # terminator codon ("complete", the 5th argument below)
213 $seq->seq('TTGGTGGCGTCAACTTAA');    # TTG GTG GCG TCA ACT TAA
214 $aa = $seq->translate( undef, undef, undef, undef, 1 );
215 is $aa->seq, 'MVAST', "Translation: " . $aa->seq;
217 # same test as previous, but using named parameter
218 $aa = $seq->translate( -complete => 1 );
219 is $aa->seq, 'MVAST', "Translation: " . $aa->seq;
221 # find ORF, ignore codons outside the ORF or CDS
222 $seq->seq('TTTTATGGTGGCGTCAACTTAATTT');    # ATG GTG GCG TCA ACT
223 $aa = $seq->translate( -orf => 1 );
224 is $aa->seq, 'MVAST*', "Translation: " . $aa->seq;
226 # smallest possible ORF
227 $seq->seq("ggggggatgtagcccc");             # atg tga
228 $aa = $seq->translate( -orf => 1 );
229 is $aa->seq, 'M*', "Translation: " . $aa->seq;
231 # same as previous but complete, so * is removed
232 $aa = $seq->translate(
233     -orf      => 1,
234     -complete => 1
236 is $aa->seq, 'M', "Translation: " . $aa->seq;
238 # ORF without termination codon
239 # should warn, let's change it into throw for testing
240 $seq->verbose(2);
241 $seq->seq("ggggggatgtggcccc");    # atg tgg ccc
242 eval { $seq->translate( -orf => 1 ); };
243 like( $@, qr/\batgtggccc\b/i );
244 $seq->verbose(-1);
245 $aa = $seq->translate( -orf => 1 );
246 is $aa->seq, 'MWP', "Translation: MWP";
247 $seq->verbose(0);
249 # use non-standard codon table where terminator is read as Q
250 $seq->seq('ATGGTGGCGTCAACTTAG');    # ATG GTG GCG TCA ACT TAG
251 $aa = $seq->translate( -codontable_id => 6 );
252 is $aa->seq, 'MVASTQ' or diag( "Translation: " . $aa->seq );
254 # insert an odd character instead of terminating with *
255 $aa = $seq->translate( -terminator => 'X' );
256 is $aa->seq, 'MVASTX' or diag( "Translation: " . $aa->seq );
258 # change frame from default
259 $aa = $seq->translate( -frame => 1 );    # TGG TGG CGT CAA CTT AG
260 is $aa->seq, 'WWRQL' or diag( "Translation: " . $aa->seq );
262 $aa = $seq->translate( -frame => 2 );    # GGT GGC GTC AAC TTA G
263 is $aa->seq, 'GGVNL' or diag( "Translation: " . $aa->seq );
265 # TTG is initiator in Standard codon table? Afraid so.
266 $seq->seq("ggggggttgtagcccc");           # ttg tag
267 $aa = $seq->translate( -orf => 1 );
268 is $aa->seq, 'L*' or diag( "Translation: " . $aa->seq );
270 # Replace L at 1st position with M by setting complete to 1
271 $seq->seq("ggggggttgtagcccc");           # ttg tag
272 $aa = $seq->translate(
273     -orf      => 1,
274     -complete => 1
276 is $aa->seq, 'M' or diag( "Translation: " . $aa->seq );
278 # Ignore non-ATG initiators (e.g. TTG) in codon table
279 $seq->seq("ggggggttgatgtagcccc");        # atg tag
280 $aa = $seq->translate(
281     -orf      => 1,
282     -start    => "atg",
283     -complete => 1
285 is $aa->seq, 'M' or diag( "Translation: " . $aa->seq );
287 # test for character '?' in the sequence string
288 is $seq->seq('TTGGTGGCG?CAACT'), 'TTGGTGGCG?CAACT';
290 # issue #105 - when there are starts and stops in both frame 0 and
291 # frame 1, frame 0 start < frame 1 start, then should return the frame
292 # 0 ORF per the pod ('the first orf') even if frame 1 stop < frame 0 stop
294 $seq->seq('ATGAATGTAAATAA');
295 $aa = $seq->translate( -orf => 1 );
296 my $aa0 = $seq->translate(-frame => 0);
297 is $aa->seq, $aa0->seq, "frame 0 start, frame 1 stop < frame 0 stop";
298 $seq->seq('AAATGAATGTAAATAA');
299 $aa = $seq->translate( -orf => 1, -frame=>1 );
300 my $aa2 = $seq->translate(-frame => 2);
301 is $aa->seq, $aa2->seq, "frame 1 start, frame 2 stop < frame 1 stop";
303 # test for some aliases
304 $seq = Bio::PrimarySeq->new(
305     -id          => 'aliasid',
306     -description => 'Alias desc'
308 is $seq->description, 'Alias desc';
309 is $seq->display_id,  'aliasid';
311 # Test alphabet
313 ok $seq->seq('actgx');
314 is $seq->alphabet, 'protein', 'Alphabet';
315 ok $seq->seq('actge');
316 is $seq->alphabet, 'protein';
317 ok $seq->seq('actgf');
318 is $seq->alphabet, 'protein';
319 ok $seq->seq('actgi');
320 is $seq->alphabet, 'protein';
321 ok $seq->seq('actgj');
322 is $seq->alphabet, 'protein';
323 ok $seq->seq('actgl');
324 is $seq->alphabet, 'protein';
325 ok $seq->seq('actgo');
326 is $seq->alphabet, 'protein';
327 ok $seq->seq('actgp');
328 is $seq->alphabet, 'protein';
329 ok $seq->seq('actgq');
330 is $seq->alphabet, 'protein';
331 ok $seq->seq('actgz');
332 is $seq->alphabet, 'protein';
333 ok $seq->seq('actgn');
334 is $seq->alphabet, 'dna';
335 ok $seq->seq('acugn');
336 is $seq->alphabet, 'rna';
337 ok $seq->seq('bdhkm');
338 is $seq->alphabet, 'protein';
339 ok $seq->seq('rsvwx');
340 is $seq->alphabet, 'protein';
341 ok $seq->seq('AAACTYAAAAGAATTGRCGG'); # valid degenerate DNA PCR primer sequence (90% ACGTN)
342 is $seq->alphabet, 'dna';
343 ok $seq->seq('AAACTYAAAKGAATTGRCGG'); # another primer previously detected as protein (85% ACGTN)
344 is $seq->alphabet, 'dna';
345 ok $seq->seq('YWACTYAAAKGARTTGRCGG'); # 70% ACGTNWSRM. Everything <= 70% is considered a protein
346 is $seq->alphabet, 'dna';
347 ok $seq->seq('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'); # Bug 2438
348 is $seq->alphabet, 'protein', 'Bug 2438';
349 ok $seq->seq('CAGTCXXXXXXXXXXXXXXXXXXXXXXXXXXXCAGCG');
350 is $seq->alphabet, 'protein';
351 ok $seq->seq('WTGGGGCTATGAAAAAAAAAWTTKMGMMAAAAAWTTWTKRWMRATC'); # showed up on MAKER list
352 is $seq->alphabet, 'dna';
354 ok $seq->seq('actgn', 'protein'); # accept specified alphabet, no matter what
355 is $seq->alphabet, 'protein';
356 ok $seq->seq('bdhkm', 'dna');
357 is $seq->alphabet, 'dna';
360 # Bug #2864:
362 $seq = Bio::PrimarySeq->new( -display_id => 0, -seq => 'GATC' );
364 is $seq->display_id, 0, "Bug #2864";
366 # Test that the check for terminators inside the translated protein
367 # works when the terminator isn't '*':
369 $seq = Bio::PrimarySeq->new(-seq=>'ATGCTCTAAGCAGGGTAA'); # ML*AG*
370 eval { $aa = $seq->translate(-complete=>1, -throw=>1, -terminator=>'#') };
371 my $error = $@;
372 ok $error =~ /\QTerminator codon inside CDS!\E/, 'Terminator + inside sequence';
374 $seq = Bio::PrimarySeq->new(-seq=>'ATGCTCGCAGGGTAA'); # MLAG*
375 $aa = $seq->translate(-complete=>1, -throw=>1, -terminator=>'#');
376 is $aa->seq, 'MLAG';
379 # Test length method
380 ok $seq = Bio::PrimarySeq->new(), 'Length method';
381 is $seq->length, 0;
382 ok $seq->length(123);
383 is $seq->length, 123;
385 ok $seq = Bio::PrimarySeq->new( -seq => 'ATGCTCTAAGCAGGGTAA' );
386 is $seq->length, 18;
387 ok $seq->seq('ATGCTCTAAG');
388 is $seq->length, 10;
389 is $seq->seq(undef), undef;
390 is $seq->length, 0;
392 ok $seq = Bio::PrimarySeq->new( -length => 123 );
393 is $seq->length, 123;
395 ok $seq = Bio::PrimarySeq->new( -seq => 'ATGCTCTAAGCAGGGTAA' );
396 is $seq->length, 18;
397 ok $seq->length( $seq->length ); # save memory by removing seq
398 is $seq->seq( undef ), undef;    # ... but keeping a record of length
399 is $seq->length, 18;
400 is $seq->seq, undef;
401 ok $seq->seq('ACGT');
402 is $seq->length, 4; # manually-specified length changed when sequence is changed
404 throws_ok { $seq->length(666); } qr/.+/; # Cannot lie about length
407 # Sequence validation method
408 is $seq->validate_seq( undef    ), 1;
409 is $seq->validate_seq( ''       ), 1;
410 is $seq->validate_seq( 'acgt'   ), 1;
411 is $seq->validate_seq( 'ACGT'   ), 1;
412 is $seq->validate_seq( 'XFRH'   ), 1;
413 is $seq->validate_seq( '-~'     ), 1; # gap symbols
414 is $seq->validate_seq( '-.*?=~' ), 1; # other valid symbols
415 is $seq->validate_seq( '0'      ), 0;
416 is $seq->validate_seq( '   '    ), 0;
417 is $seq->validate_seq( 'AAAA$'  ), 0;
418 is $seq->validate_seq( 'tt&t!'  ), 0;
420 throws_ok { $seq->validate_seq('tt&t!', 1); } qr/.+/;
423 # Test direct option (no sequence validation)
424 throws_ok { $seq = Bio::PrimarySeq->new(-seq => 'A\T$AGQ+T'); } qr/.+/, 'Validation';
425 ok $seq = Bio::PrimarySeq->new( -seq => 'A\T$AGQ+T', -direct => 1 );
426 is $seq->seq, 'A\T$AGQ+T';
427 throws_ok { $seq->seq('NT@/') } qr/.+/;
429 # Set a sequence by reference
430 my $string = 'AAAACCCCGGGGTTTT';
431 ok $seq = Bio::PrimarySeq->new( -ref_to_seq => \$string );
432 is $seq->seq, 'AAAACCCCGGGGTTTT';
435 # Test internal PrimarySeqI _find_orfs function and translate( -orf => 'longest' )
437     my @tests = (
438         #tiny test
439         ['TTTTATGGTGGCGTCAACTTAATTT',
440          [[4,22,18,1]],
441         ],
443         #bigger test (this is a tomato unigene)
444         ['GAAGGCTGGTTCTGAGTTGGATCTATGTTTGATGAAGGGAAGTAGACCGGAGGTCTTGCATCAGCAATATTAGTACCAAATCCAGGTGGAGGCGCATCCTGTCTCCGTTGCATTTCAACTTTCATTTCAGCAATCTGTTGCATCAGTTGCATGATCAATTCATTCTGTTCCACTACAGTGGGCTGAGCGACCACAACGTCAGTAAGACGCCCTTCGTCATTGTTGTCTCCCATAACTGTTTTTCCTTTATCTGAATTTGATCGAGGGAAGGAATCTGTAGGACCTTTCGATCTGGTGAAGTAAGGATGATCTGCCAGCTTTATTGACACAGATCAGTAAAAAGGTACCTGAAAGGTAAAAACAACTCAAAGGCAAATTTGTTAGTGCATATCCAGAGTACAAAATGCTTAATATCGCACATAAAACCGATAAACACACAAGTCGTTTTGTTTGAGGATATCTTAACCCACGAATAAGGACGGATATATATTTTGAACAAACAGGAATTTGTTTGTTTGGCGTTATCTTGGGAAATCTG',
445          [[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]],
446         ],
449        );
450     foreach my $test (@tests) {
451         my ($test_seq, $orfs) = @$test;
452         my @orfs = Bio::PrimarySeqI::_find_orfs_nucleotide(
453             undef,
454             $test_seq,
455             Bio::Tools::CodonTable->new,
456             undef,
457            ); # ATG GTG GCG TCA ACT
458         is_deeply( \@orfs, $orfs, '_find_orfs 1')
459             or diag "for $test_seq, _find_orfs returned:\n"
460                     .Dumper([map [@$_], @orfs]);
462         is_deeply( $orfs->[0],
463                    (sort {$b->[2] <=> $a->[2]} @$orfs)[0],
464                    'orfs are sorted by descending length'
465                   );
467         # make sure we get the same sequence by taking the longest orf
468         # nucleotide from the test data and translating it, as by
469         # calling translate with -orf => 'longest'
470         is(
471             Bio::PrimarySeq
472               ->new( -seq => $test_seq, -id => 'fake_id' )
473               ->translate( -orf => 'longest' )
474               ->seq,
476             Bio::PrimarySeq
477               ->new( -seq => substr( $test_seq, $orfs->[0][0], $orfs->[0][2] ),
478                      -id => 'foo'
479                     )
480               ->translate
481               ->seq,
482             'got correct -orf => "longest" seq',
483            );
484     }
487 #####
488 # Extensive location and subsequence tests
489 ok $seq = Bio::PrimarySeq->new('-seq' => 'AAAAACCCCCGGGGGTTTTT',);
490 ok $seq->is_circular(1);
492 # NOTE: "_no_strand" variables tests the possibility that the user didn't set
493 # Strand for positive coordinates (or the object comes from
494 # Bio::Factory::FTLocationFactory->from_string)
496 # Single location
497 # Coordinates: 1..5 => AAAAA
498 # Revcom: complement(1..5) => TTTTT
499 ok my $loc1_strand    = Bio::Location::Simple->new('-start' => 1, '-end' => 5,'-strand' => 1);
500 ok my $loc1_no_strand = Bio::Location::Simple->new('-start' => 1, '-end' => 5);
501 is $seq->subseq($loc1_strand),    'AAAAA';
502 is $seq->subseq($loc1_no_strand), 'AAAAA';
503 is $loc1_strand->to_FTstring,     '1..5';
504 is $loc1_no_strand->to_FTstring,  '1..5';
505 $loc1_strand->flip_strand;
506 $loc1_no_strand->flip_strand;
507 is $seq->subseq($loc1_strand),    'TTTTT';
508 is $seq->subseq($loc1_no_strand), 'TTTTT';
509 is $loc1_strand->to_FTstring,     'complement(1..5)';
510 is $loc1_no_strand->to_FTstring,  'complement(1..5)';
511 is $loc1_strand->length,    5;
512 is $loc1_no_strand->length, 5;
514 # Basic split, both locations in positive strand
515 # Coords: join(6..10,16..20) => CCCCCTTTTT
516 # Revcom: complement(join(6..10,16..20)) => AAAAAGGGGG
517 ok my $loc2_strand    = Bio::Location::Split->new();
518 ok my $loc2_no_strand = Bio::Location::Split->new();
519 ok $loc2_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 6,  '-end' => 10, '-strand' => 1) );
520 ok $loc2_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 16, '-end' => 20, '-strand' => 1) );
521 ok $loc2_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 6,  '-end' => 10) );
522 ok $loc2_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 16, '-end' => 20) );
523 is $seq->subseq($loc2_strand),    'CCCCCTTTTT';
524 is $seq->subseq($loc2_no_strand), 'CCCCCTTTTT';
525 is $loc2_strand->to_FTstring,     'join(6..10,16..20)';
526 is $loc2_no_strand->to_FTstring,  'join(6..10,16..20)';
527 $loc2_strand->flip_strand;
528 $loc2_no_strand->flip_strand;
529 is $seq->subseq($loc2_strand),    'AAAAAGGGGG';
530 is $seq->subseq($loc2_no_strand), 'AAAAAGGGGG';
531 is $loc2_strand->to_FTstring,     'complement(join(6..10,16..20))';
532 is $loc2_no_strand->to_FTstring,  'complement(join(6..10,16..20))';
533 is $loc2_strand->length,    15;
534 is $loc2_no_strand->length, 15;
536 # Basic split, both locations in negative strand
537 # Coords: complement(join(6..10,16..20)) => AAAAAGGGGG
538 # Revcom: join(6..10,16..20) => CCCCCTTTTT
539 my $loc3_strand    = Bio::Location::Split->new();
540 $loc3_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 6,  '-end' => 10, '-strand' => -1) );
541 $loc3_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 16, '-end' => 20, '-strand' => -1) );
542 is $seq->subseq($loc3_strand),    'AAAAAGGGGG';
543 is $loc3_strand->to_FTstring,     'complement(join(6..10,16..20))';
544 $loc3_strand->flip_strand;
545 is $seq->subseq($loc3_strand),    'CCCCCTTTTT';
546 is $loc3_strand->to_FTstring,     'join(6..10,16..20)';
547 is $loc3_strand->length, 15;
549 ## Cut by origin-split, same strand, single sequence that pass through origin
550 #Coords: join(16..20,1..2) => TTTTTAA
551 #Revcom: complement(join(16..20,1..2)) => TTAAAAA
552 my $loc4_strand    = Bio::Location::Split->new();
553 my $loc4_no_strand = Bio::Location::Split->new();
554 $loc4_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 16, '-end' => 20, '-strand' => 1) );
555 $loc4_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 1,  '-end' => 2,  '-strand' => 1) );
556 $loc4_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 16, '-end' => 20) );
557 $loc4_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 1,  '-end' => 2)  );
558 is $seq->subseq($loc4_strand),    'TTTTTAA';
559 is $seq->subseq($loc4_no_strand), 'TTTTTAA';
560 is $loc4_strand->to_FTstring,     'join(16..20,1..2)';
561 is $loc4_no_strand->to_FTstring,  'join(16..20,1..2)';
562 $loc4_strand->flip_strand;
563 $loc4_no_strand->flip_strand;
564 is $seq->subseq($loc4_strand),    'TTAAAAA';
565 is $seq->subseq($loc4_no_strand), 'TTAAAAA';
566 is $loc4_strand->to_FTstring,     'complement(join(16..20,1..2))';
567 is $loc4_no_strand->to_FTstring,  'complement(join(16..20,1..2))';
568 is $loc4_strand->length,    7;
569 is $loc4_no_strand->length, 7;
571 ## Cut by origin-combo split, same strand, 2 sequences with 1st passing through origin
572 #Coords: join(19..20,1..2,11..13) => TTAAGGG
573 #Revcom: complement(join(19..20,1..2,11..13)) => CCCTTAA
574 my $loc5_strand    = Bio::Location::Split->new();
575 my $loc5_no_strand = Bio::Location::Split->new();
576 $loc5_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 19, '-end' => 20, '-strand' => 1) );
577 $loc5_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 1,  '-end' => 2,  '-strand' => 1) );
578 $loc5_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 11, '-end' => 13, '-strand' => 1) );
579 $loc5_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 19, '-end' => 20) );
580 $loc5_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 1,  '-end' => 2)  );
581 $loc5_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 11, '-end' => 13) );
582 is $seq->subseq($loc5_strand),    'TTAAGGG';
583 is $seq->subseq($loc5_no_strand), 'TTAAGGG';
584 is $loc5_strand->to_FTstring,     'join(19..20,1..2,11..13)';
585 is $loc5_no_strand->to_FTstring,  'join(19..20,1..2,11..13)';
586 $loc5_strand->flip_strand;
587 $loc5_no_strand->flip_strand;
588 is $seq->subseq($loc5_strand),    'CCCTTAA';
589 is $seq->subseq($loc5_no_strand), 'CCCTTAA';
590 is $loc5_strand->to_FTstring,     'complement(join(19..20,1..2,11..13))';
591 is $loc5_no_strand->to_FTstring,  'complement(join(19..20,1..2,11..13))';
592 is $loc5_strand->length,    15;
593 is $loc5_no_strand->length, 15;
595 ## Cut by origin-combo split, same strand, 2 sequences with 2nd passing through origin
596 #Coords: join(6..10,19..20,1..4) => CCCCCTTAAAA
597 #Revcom: complement(join(6..10,19..20,1..4)) => TTTTAAGGGGG
598 my $loc6_strand    = Bio::Location::Split->new();
599 my $loc6_no_strand = Bio::Location::Split->new();
600 $loc6_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 6,  '-end' => 10, '-strand' => 1) );
601 $loc6_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 19, '-end' => 20, '-strand' => 1) );
602 $loc6_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 1,  '-end' => 4,  '-strand' => 1) );
603 $loc6_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 6,  '-end' => 10) );
604 $loc6_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 19, '-end' => 20) );
605 $loc6_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 1,  '-end' => 4)  );
606 is $seq->subseq($loc6_strand),    'CCCCCTTAAAA';
607 is $seq->subseq($loc6_no_strand), 'CCCCCTTAAAA';
608 is $loc6_strand->to_FTstring,     'join(6..10,19..20,1..4)';
609 is $loc6_no_strand->to_FTstring,  'join(6..10,19..20,1..4)';
610 $loc6_strand->flip_strand;
611 $loc6_no_strand->flip_strand;
612 is $seq->subseq($loc6_strand),    'TTTTAAGGGGG';
613 is $seq->subseq($loc6_no_strand), 'TTTTAAGGGGG';
614 is $loc6_strand->to_FTstring,     'complement(join(6..10,19..20,1..4))';
615 is $loc6_no_strand->to_FTstring,  'complement(join(6..10,19..20,1..4))';
616 is $loc6_strand->length,    19;
617 is $loc6_no_strand->length, 19;
619 ## Trans-splicing, 2 sequences in different strands, 2nd in complement
620 #Coords: join(6..10,complement(16..20)) => CCCCCAAAAA
621 #Revcom: join(16..20,complement(6..10)) => TTTTTGGGGG
622 my $loc7_strand    = Bio::Location::Split->new();
623 my $loc7_no_strand = Bio::Location::Split->new();
624 $loc7_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 6,  '-end' => 10, '-strand' =>  1) );
625 $loc7_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 16, '-end' => 20, '-strand' => -1) );
626 $loc7_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 6,  '-end' => 10) );
627 $loc7_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 16, '-end' => 20, '-strand' => -1) );
628 is $seq->subseq($loc7_strand),    'CCCCCAAAAA';
629 is $seq->subseq($loc7_no_strand), 'CCCCCAAAAA';
630 is $loc7_strand->to_FTstring,     'join(6..10,complement(16..20))';
631 is $loc7_no_strand->to_FTstring,  'join(6..10,complement(16..20))';
632 $loc7_strand->flip_strand;
633 $loc7_no_strand->flip_strand;
634 is $seq->subseq($loc7_strand),    'TTTTTGGGGG';
635 is $seq->subseq($loc7_no_strand), 'TTTTTGGGGG';
636 is $loc7_strand->to_FTstring,     'join(16..20,complement(6..10))';
637 is $loc7_no_strand->to_FTstring,  'join(16..20,complement(6..10))';
638 is $loc7_strand->length,    10;
639 is $loc7_no_strand->length, 10;
641 ## Trans-splicing, 2 sequences in different strands, 1st in complement
642 #Coords: join(complement(16..20),6..10) => AAAAACCCCC
643 #Revcom: join(complement(6..10),16..20) => GGGGGTTTTT
644 my $loc8_strand    = Bio::Location::Split->new();
645 my $loc8_no_strand = Bio::Location::Split->new();
646 $loc8_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 16, '-end' => 20, '-strand' => -1) );
647 $loc8_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 6,  '-end' => 10, '-strand' =>  1) );
648 $loc8_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 16, '-end' => 20, '-strand' => -1) );
649 $loc8_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 6,  '-end' => 10) );
650 is $seq->subseq($loc8_strand),    'AAAAACCCCC';
651 is $seq->subseq($loc8_no_strand), 'AAAAACCCCC';
652 is $loc8_strand->to_FTstring,     'join(complement(16..20),6..10)';
653 is $loc8_no_strand->to_FTstring,  'join(complement(16..20),6..10)';
654 $loc8_strand->flip_strand;
655 $loc8_no_strand->flip_strand;
656 is $seq->subseq($loc8_strand),    'GGGGGTTTTT';
657 is $seq->subseq($loc8_no_strand), 'GGGGGTTTTT';
658 is $loc8_strand->to_FTstring,     'join(complement(6..10),16..20)';
659 is $loc8_no_strand->to_FTstring,  'join(complement(6..10),16..20)';
660 is $loc8_strand->length,    10;
661 is $loc8_no_strand->length, 10;
663 ## Trans-splicing w/cut by origin, 2 sequences with 1st passing through origin, 2nd in complement
664 #Coords: join(19..20,1..3,complement(11..13)) => TTAAACCC
665 #Revcom: join(11..13,complement(1..3),complement(19..20)) => GGGTTTAA
666 my $loc9_strand    = Bio::Location::Split->new();
667 my $loc9_no_strand = Bio::Location::Split->new();
668 $loc9_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 19, '-end' => 20, '-strand' =>  1) );
669 $loc9_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 1,  '-end' => 3,  '-strand' =>  1) );
670 $loc9_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 11, '-end' => 13, '-strand' => -1) );
671 $loc9_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 19, '-end' => 20) );
672 $loc9_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 1,  '-end' => 3)  );
673 $loc9_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 11, '-end' => 13, '-strand' => -1) );
674 is $seq->subseq($loc9_strand),    'TTAAACCC';
675 is $seq->subseq($loc9_no_strand), 'TTAAACCC';
676 is $loc9_strand->to_FTstring,     'join(19..20,1..3,complement(11..13))';
677 is $loc9_no_strand->to_FTstring,  'join(19..20,1..3,complement(11..13))';
678 $loc9_strand->flip_strand;
679 $loc9_no_strand->flip_strand;
680 is $seq->subseq($loc9_strand),    'GGGTTTAA';
681 is $seq->subseq($loc9_no_strand), 'GGGTTTAA';
682 is $loc9_strand->to_FTstring,     'join(11..13,complement(1..3),complement(19..20))';
683 is $loc9_no_strand->to_FTstring,  'join(11..13,complement(1..3),complement(19..20))';
684 is $loc9_strand->length,    8;
685 is $loc9_no_strand->length, 8;
687 ## Trans-splicing w/cut by origin, 2 sequences with 1st passing through origin, 1st in complement
688 #Coords: join(complement(1..3),complement(19..20),11..13) => TTTAAGGG
689 #Revcom: join(complement(11..13),19..20,1..3) => CCCTTAAA
690 my $loc10_strand    = Bio::Location::Split->new();
691 my $loc10_no_strand = Bio::Location::Split->new();
692 $loc10_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 1,  '-end' => 3,  '-strand' => -1) );
693 $loc10_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 19, '-end' => 20, '-strand' => -1) );
694 $loc10_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 11, '-end' => 13, '-strand' =>  1) );
695 $loc10_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 1,  '-end' => 3,  '-strand' => -1) );
696 $loc10_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 19, '-end' => 20, '-strand' => -1) );
697 $loc10_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 11, '-end' => 13) );
698 is $seq->subseq($loc10_strand),    'TTTAAGGG';
699 is $seq->subseq($loc10_no_strand), 'TTTAAGGG';
700 is $loc10_strand->to_FTstring,     'join(complement(1..3),complement(19..20),11..13)';
701 is $loc10_no_strand->to_FTstring,  'join(complement(1..3),complement(19..20),11..13)';
702 $loc10_strand->flip_strand;
703 $loc10_no_strand->flip_strand;
704 is $seq->subseq($loc10_strand),    'CCCTTAAA';
705 is $seq->subseq($loc10_no_strand), 'CCCTTAAA';
706 is $loc10_strand->to_FTstring,     'join(complement(11..13),19..20,1..3)';
707 is $loc10_no_strand->to_FTstring,  'join(complement(11..13),19..20,1..3)';
708 is $loc10_strand->length,    8;
709 is $loc10_no_strand->length, 8;
711 ## Trans-splicing w/cut by origin, 2 sequences with 2nd passing through origin, 2nd in complement
712 #Coords: join(6..10,complement(1..2),complement(18..20)) => CCCCCTTAAA
713 #Revcom: join(18..20,1..2,complement(6..10)) => TTTAAGGGGG
714 my $loc11_strand    = Bio::Location::Split->new();
715 my $loc11_no_strand = Bio::Location::Split->new();
716 $loc11_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 6,  '-end' => 10, '-strand' =>  1) );
717 $loc11_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 1,  '-end' => 2,  '-strand' => -1) );
718 $loc11_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 18, '-end' => 20, '-strand' => -1) );
719 $loc11_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 6,  '-end' => 10) );
720 $loc11_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 1,  '-end' => 2,  '-strand' => -1) );
721 $loc11_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 18, '-end' => 20, '-strand' => -1) );
722 is $seq->subseq($loc11_strand),    'CCCCCTTAAA';
723 is $seq->subseq($loc11_no_strand), 'CCCCCTTAAA';
724 is $loc11_strand->to_FTstring,     'join(6..10,complement(1..2),complement(18..20))';
725 is $loc11_no_strand->to_FTstring,  'join(6..10,complement(1..2),complement(18..20))';
726 $loc11_strand->flip_strand;
727 $loc11_no_strand->flip_strand;
728 is $seq->subseq($loc11_strand),    'TTTAAGGGGG';
729 is $seq->subseq($loc11_no_strand), 'TTTAAGGGGG';
730 is $loc11_strand->to_FTstring,     'join(18..20,1..2,complement(6..10))';
731 is $loc11_no_strand->to_FTstring,  'join(18..20,1..2,complement(6..10))';
732 is $loc11_strand->length,    10;
733 is $loc11_no_strand->length, 10;
735 ## Trans-splicing w/cut by origin, 2 sequences with 2nd passing through origin, 1st in complement
736 #Coords: join(complement(6..10),18..20,1..2) => GGGGGTTTAA
737 #Revcom: join(complement(1..2),complement(18..20),6..10) => TTAAACCCCC
738 my $loc12_strand    = Bio::Location::Split->new();
739 my $loc12_no_strand = Bio::Location::Split->new();
740 $loc12_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 6,  '-end' => 10, '-strand' => -1) );
741 $loc12_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 18, '-end' => 20, '-strand' =>  1) );
742 $loc12_strand->add_sub_Location(    Bio::Location::Simple->new('-start'  => 1,  '-end' => 2,  '-strand' =>  1) );
743 $loc12_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 6,  '-end' => 10, '-strand' => -1) );
744 $loc12_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 18, '-end' => 20) );
745 $loc12_no_strand->add_sub_Location( Bio::Location::Simple->new('-start'  => 1,  '-end' => 2)  );
746 is $seq->subseq($loc12_strand),    'GGGGGTTTAA';
747 is $seq->subseq($loc12_no_strand), 'GGGGGTTTAA';
748 is $loc12_strand->to_FTstring,     'join(complement(6..10),18..20,1..2)';
749 is $loc12_no_strand->to_FTstring,  'join(complement(6..10),18..20,1..2)';
750 $loc12_strand->flip_strand;
751 $loc12_no_strand->flip_strand;
752 is $seq->subseq($loc12_strand),    'TTAAACCCCC';
753 is $seq->subseq($loc12_no_strand), 'TTAAACCCCC';
754 is $loc12_strand->to_FTstring,     'join(complement(1..2),complement(18..20),6..10)';
755 is $loc12_no_strand->to_FTstring,  'join(complement(1..2),complement(18..20),6..10)';
756 is $loc12_strand->length,    10;
757 is $loc12_no_strand->length, 10;