t/SeqTools/SeqUtils.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7 #    use List::MoreUtils qw(uniq);
   8     use Bio::Root::Test;
   9
  10     test_begin(-tests => 133);
  11
  12     use_ok('Bio::PrimarySeq');
  13     use_ok('Bio::SeqUtils');
  14     use_ok('Bio::SeqFeature::Generic');
  15     use_ok('Bio::Annotation::SimpleValue');
  16     use_ok('Bio::Annotation::Collection');
  17     use_ok('Bio::Annotation::Comment');
  18 }
  19
  20 my ($seq, $util, $ascii, $ascii_aa, $ascii3);
  21
  22 # Entire alphabet now IUPAC-endorsed and used in GenBank (Oct 2006)
  23 $ascii =    'ABCDEFGHIJKLMNOPQRSTUVWXYZ';
  24 $ascii_aa = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ';
  25
  26 $ascii3 =
  27     'AlaAsxCysAspGluPheGlyHisIleXleLysLeuMetAsnPylProGlnArgSerThrSecValTrpXaaTyrGlx';
  28
  29 $seq = Bio::PrimarySeq->new('-seq'=> $ascii,
  30                             '-alphabet'=>'protein',
  31                             '-id'=>'test');
  32
  33 # one letter amino acid code to three letter code
  34 ok $util = Bio::SeqUtils->new();
  35 is $util->seq3($seq), $ascii3;
  36
  37 #using anonymous hash
  38 is (Bio::SeqUtils->seq3($seq), $ascii3);
  39 is (Bio::SeqUtils->seq3($seq, undef, ','),
  40     'Ala,Asx,Cys,Asp,Glu,Phe,Gly,His,Ile,Xle,Lys,'.
  41     'Leu,Met,Asn,Pyl,Pro,Gln,Arg,Ser,Thr,Sec,Val,Trp,Xaa,Tyr,Glx');
  42
  43 $seq->seq('asd-KJJK-');
  44 is (Bio::SeqUtils->seq3($seq, '-', ':'),
  45     'Ala:Ser:Asp:Ter:Lys:Xle:Xle:Lys:Ter');
  46
  47 # three letter amino acid code to one letter code
  48 ok (Bio::SeqUtils->seq3in($seq, 'AlaPYHCysAspGlu'));
  49 is $seq->seq, 'AXCDE';
  50 is (Bio::SeqUtils->seq3in($seq, $ascii3)->seq, $ascii_aa);
  51
  52 #
  53 # Tests for multiframe translations
  54 #
  55
  56 $seq = Bio::PrimarySeq->new('-seq'=> 'agctgctgatcggattgtgatggctggatggcttgggatgctgg',
  57                             '-alphabet'=>'dna',
  58                             '-id'=>'test2');
  59
  60 my @a = $util->translate_3frames($seq);
  61 is scalar @a, 3;
  62 #foreach $a (@a) {
  63 #    print 'ID: ', $a->id, ' ', $a->seq, "\n";
  64 #}
  65
  66 @a = $util->translate_6frames($seq);
  67 is scalar @a, 6;
  68 #foreach $a (@a) {
  69 #    print 'ID: ', $a->id, ' ', $a->seq, "\n";
  70 #}
  71
  72 #
  73 # test for valid AA return
  74 #
  75
  76 my @valid_aa = sort Bio::SeqUtils->valid_aa;
  77 is(@valid_aa, 27);
  78 is($valid_aa[1], 'A');
  79
  80 @valid_aa = sort Bio::SeqUtils->valid_aa(1);
  81 is(@valid_aa, 27);
  82 is ($valid_aa[1], 'Arg');
  83
  84 my %valid_aa = Bio::SeqUtils->valid_aa(2);
  85 is keys %valid_aa, 54;
  86 is($valid_aa{'C'}, 'Cys');
  87 is( $valid_aa{'Cys'}, 'C');
  88
  89
  90 #
  91 # Mutate
  92 #
  93
  94 SKIP: {
  95     test_skip(-tests => 4,
  96               -requires_modules => ['Bio::LiveSeq::Mutation']);
  97     use_ok('Bio::LiveSeq::Mutation');
  98     my $string1 = 'aggt';
  99     $seq = Bio::PrimarySeq->new('-seq'=> 'aggt',
 100                                 '-alphabet'=>'dna',
 101                                 '-id'=>'test3');
 102
 103     # point
 104     Bio::SeqUtils->mutate($seq,
 105                           Bio::LiveSeq::Mutation->new(-seq => 'c',
 106                                                       -pos => 3
 107                           )
 108         );
 109     is $seq->seq, 'agct';
 110
 111     # insertion and deletion
 112     my @mutations = (
 113         Bio::LiveSeq::Mutation->new(-seq => 'tt',
 114                                     -pos => 2,
 115                                     -len => 0
 116         ),
 117         Bio::LiveSeq::Mutation->new(-pos => 2,
 118                                     -len => 2
 119         )
 120         );
 121
 122     Bio::SeqUtils->mutate($seq, @mutations);
 123     is $seq->seq, 'agct';
 124
 125     # insertion to the end of the sequence
 126     Bio::SeqUtils->mutate($seq,
 127                           Bio::LiveSeq::Mutation->new(-seq => 'aa',
 128                                                       -pos => 5,
 129                                                       -len => 0
 130                           )
 131         );
 132     is $seq->seq, 'agctaa';
 133 }
 134
 135
 136 #
 137 # testing Bio::SeqUtils->cat
 138 #
 139
 140 # PrimarySeqs
 141
 142 my $primseq1 = Bio::PrimarySeq->new(-id => 1, -seq => 'acgt', -description => 'master');
 143 my $primseq2 = Bio::PrimarySeq->new(-id => 2, -seq => 'tgca');
 144
 145 Bio::SeqUtils->cat($primseq1, $primseq2);
 146 is $primseq1->seq, 'acgttgca';
 147 is $primseq1->description, 'master';
 148
 149 #should work for Bio::LocatableSeq
 150 #should work for Bio::Seq::MetaI Seqs?
 151
 152
 153 # Bio::SeqI
 154
 155 my $seq1 = Bio::Seq->new(-id => 1, -seq => 'aaaa', -description => 'first');
 156 my $seq2 = Bio::Seq->new(-id => 2, -seq => 'tttt', -description => 'second');
 157 my $seq3 = Bio::Seq->new(-id => 3, -seq => 'cccc', -description => 'third');
 158
 159
 160 #  annotations
 161 my $ac2 = Bio::Annotation::Collection->new();
 162 my $simple1 = Bio::Annotation::SimpleValue->new(
 163                                                 -tagname => 'colour',
 164                                                 -value   => 'blue'
 165                                                ), ;
 166 my $simple2 = Bio::Annotation::SimpleValue->new(
 167                                                 -tagname => 'colour',
 168                                                 -value   => 'black'
 169                                                ), ;
 170 $ac2->add_Annotation('simple',$simple1);
 171 $ac2->add_Annotation('simple',$simple2);
 172 $seq2->annotation($ac2);
 173
 174 my $ac3 = Bio::Annotation::Collection->new();
 175 my $simple3 = Bio::Annotation::SimpleValue->new(
 176                                                 -tagname => 'colour',
 177                                                 -value   => 'red'
 178                                                );
 179 $ac3->add_Annotation('simple',$simple3);
 180 $seq3->annotation($ac3);
 181
 182
 183 ok (Bio::SeqUtils->cat($seq1, $seq2, $seq3));
 184 is $seq1->seq, 'aaaattttcccc';
 185 is scalar $seq1->annotation->get_Annotations, 3;
 186
 187
 188 # seq features
 189 my $ft2 = Bio::SeqFeature::Generic->new( -start => 1,
 190                                          -end => 4,
 191                                          -strand => 1,
 192                                          -primary => 'source',
 193                                          -tag     => {note => 'note2'},
 194                                        );
 195
 196
 197 my $ft3 = Bio::SeqFeature::Generic->new( -start => 3,
 198                                          -end => 3,
 199                                          -strand => 1,
 200                                          -primary => 'hotspot',
 201                                          -tag     => {note => ['note3a','note3b'],
 202                                                       comment => 'c1'},
 203                                        );
 204
 205 $seq2->add_SeqFeature($ft2);
 206 $seq2->add_SeqFeature($ft3);
 207
 208 my $seq1_length = $seq1->length;
 209 ok (Bio::SeqUtils->cat($seq1, $seq2));
 210 is $seq1->seq, 'aaaattttcccctttt';
 211 is scalar $seq1->annotation->get_Annotations, 5;
 212 is_deeply([uniq_sort(map{$_->get_all_tags}$seq1->get_SeqFeatures)], [sort qw(note comment)], 'cat - has expected tags');
 213 is_deeply([sort map{$_->get_tagset_values('note')}$seq1->get_SeqFeatures], [sort qw(note2 note3a note3b)], 'cat - has expected tag values');
 214 my @tags;
 215 lives_ok {
 216   @tags = map{$_->get_tag_values(q(note))}$seq1->get_SeqFeatures ;
 217 } 'cat - note tag transfered (no throw)';
 218 cmp_ok(scalar(@tags),'==',3, 'cat - note tag values transfered (correct count)') ;
 219 my ($ft3_precat) = grep ($_->primary_tag eq 'hotspot', $seq2->get_SeqFeatures);
 220 is ($ft3_precat->start, 3, "get correct start of feature before 'cat'");
 221 my ($ft3_cat) = grep ($_->primary_tag eq 'hotspot', $seq1->get_SeqFeatures);
 222 is ($ft3_cat->start, 3+$seq1_length, "get correct start of feature after 'cat'");
 223
 224
 225 my $protseq = Bio::PrimarySeq->new(-id => 2, -seq => 'MVTF'); # protein seq
 226
 227 throws_ok {
 228     Bio::SeqUtils->cat($seq1, $protseq);
 229 } qr/different alphabets:/, 'different alphabets' ;
 230
 231
 232 #
 233 # evolve()
 234 #
 235
 236 $seq = Bio::PrimarySeq->new('-seq'=> 'aaaaaaaaaa',
 237                             '-id'=>'test');
 238
 239
 240
 241 $util = Bio::SeqUtils->new(-verbose => 0);
 242 ok my $newseq = $util->evolve($seq, 60, 4);
 243
 244 #  annotations
 245
 246 $seq2 = Bio::Seq->new(-id => 2, -seq => 'ggttaaaa', -description => 'second');
 247 $ac3 = Bio::Annotation::Collection->new();
 248 $simple3 = Bio::Annotation::SimpleValue->new(
 249                                                 -tagname => 'colour',
 250                                                 -value   => 'red'
 251                                             );
 252 $ac3->add_Annotation('simple',$simple3);
 253 $seq2->annotation($ac3);
 254 $ft2 = Bio::SeqFeature::Generic->new( -start => 1,
 255                                       -end => 4,
 256                                       -strand => 1,
 257                                       -primary => 'source',
 258                                       -tag     => {note => 'note2'},
 259                                     );
 260
 261
 262 $ft3 = Bio::SeqFeature::Generic->new( -start => 5,
 263                                       -end => 8,
 264                                       -strand => -1,
 265                                       -primary => 'hotspot',
 266                                       -tag     => {note => ['note3a','note3b'],
 267                                                    comment => 'c1'},
 268                                     );
 269
 270 my $ft4 = Bio::SeqFeature::Generic->new(-primary => 'CDS');
 271 $ft4->location(Bio::Location::Fuzzy->new(-start=>'<1',
 272                                          -end=>5,
 273                                          -strand=>-1));
 274
 275 $seq2->add_SeqFeature($ft2);
 276 $seq2->add_SeqFeature($ft3);
 277 $seq2->add_SeqFeature($ft4);
 278
 279 my $trunc=Bio::SeqUtils->trunc_with_features($seq2, 2, 7);
 280 is $trunc->seq, 'gttaaa';
 281 my @feat=$trunc->get_SeqFeatures;
 282 is $feat[0]->location->to_FTstring, '<1..3';
 283 is $feat[1]->location->to_FTstring, 'complement(4..>6)';
 284 is $feat[2]->location->to_FTstring, 'complement(<1..4)';
 285 is_deeply([uniq_sort(map{$_->get_all_tags}$trunc->get_SeqFeatures)], [sort qw(note comment)], 'trunc_with_features - has expected tags');
 286 is_deeply([sort map{$_->get_tagset_values('note')}$trunc->get_SeqFeatures], [sort qw(note2 note3a note3b)], 'trunc_with_features - has expected tag values');
 287
 288 my $revcom=Bio::SeqUtils->revcom_with_features($seq2);
 289 is $revcom->seq, 'ttttaacc';
 290 my ($rf1) = $revcom->get_SeqFeatures('hotspot');
 291 is $rf1->primary_tag, $ft3->primary_tag, 'primary_tag matches original feature...';
 292 is $rf1->location->to_FTstring, '1..4', 'but tagged sf is now revcom';
 293
 294 my ($rf2) = $revcom->get_SeqFeatures('source');
 295 is $rf2->primary_tag, $ft2->primary_tag, 'primary_tag matches original feature...';
 296 is $rf2->location->to_FTstring, 'complement(5..8)', 'but tagged sf is now revcom';
 297
 298 my ($rf3) = $revcom->get_SeqFeatures('CDS');
 299 is $rf3->primary_tag, $ft4->primary_tag, 'primary_tag matches original feature...';
 300 is $rf3->location->to_FTstring, '4..>8', 'but tagged sf is now revcom';
 301
 302 is_deeply([uniq_sort(map{$_->get_all_tags}$revcom->get_SeqFeatures)], [sort qw(note comment)], 'revcom_with_features - has expected tags');
 303 is_deeply([sort map{$_->get_tagset_values('note')}$revcom->get_SeqFeatures], [sort qw(note2 note3a note3b)], 'revcom_with_features - has expected tag values');
 304 # check circularity
 305 isnt($revcom->is_circular, 1, 'still not circular');
 306 $seq3 = Bio::Seq->new(-id => 3, -seq => 'ggttaaaa', -description => 'third', -is_circular => 1);
 307 is(Bio::SeqUtils->revcom_with_features($seq3)->is_circular, 1, 'still circular');
 308
 309
 310 # delete, insert and ligate
 311 # prepare some sequence objects
 312 my $seq_obj = Bio::Seq->new(
 313   -seq =>'aaaaaaaaaaccccccccccggggggggggtttttttttt',
 314   -display_id => 'seq1',
 315   -desc       => 'some sequence for testing'
 316 );
 317 my $subfeat1 = Bio::SeqFeature::Generic->new(
 318   -primary_tag => 'sf1',
 319   -seq_id      => 'seq1',
 320   -start       => 2,
 321   -end         => 12
 322 );
 323
 324 my $subfeat2 = Bio::SeqFeature::Generic->new(
 325   -primary_tag => 'sf2',
 326   -seq_id      => 'seq1',
 327   -start       => 14,
 328   -end         => 16
 329 );
 330 my $subfeat3 = Bio::SeqFeature::Generic->new(
 331   -primary_tag => 'sf3',
 332   -seq_id      => 'seq1',
 333   -start       => 21,
 334   -end         => 25
 335 );
 336
 337 my $composite_feat1 = Bio::SeqFeature::Generic->new(
 338   -primary_tag => 'comp_feat1',
 339   -seq_id      => 'seq1',
 340   -start       => 2,
 341   -end         => 30
 342 );
 343 my $coll_sf = Bio::Annotation::Collection->new;
 344 $coll_sf->add_Annotation(
 345   'comment', Bio::Annotation::Comment->new( '-text' => 'a comment on sf1')
 346 );
 347 $subfeat1->annotation($coll_sf);
 348
 349 $composite_feat1->add_SeqFeature( $subfeat1);
 350 $composite_feat1->add_SeqFeature( $subfeat2);
 351 $composite_feat1->add_SeqFeature( $subfeat3);
 352 my $feature1 = Bio::SeqFeature::Generic->new(
 353   -primary_tag => 'feat1',
 354   -seq_id      => 'seq1',
 355   -start       => 2,
 356   -end         => 25
 357 );
 358 my $feature2 = Bio::SeqFeature::Generic->new(
 359   -primary_tag => 'feat2',
 360   -seq_id      => 'seq1',
 361   -start       => 15,
 362   -end         => 25,
 363   -strand      => -1,
 364 );
 365 my $feature3 = Bio::SeqFeature::Generic->new(
 366   -primary_tag => 'feat3',
 367   -seq_id      => 'seq1',
 368   -start       => 30,
 369   -end         => 40
 370 );
 371 my $feature4 = Bio::SeqFeature::Generic->new(
 372   -primary_tag => 'feat4',
 373   -seq_id      => 'seq1',
 374   -start       => 1,
 375   -end         => 10
 376 );
 377 my $feature5 = Bio::SeqFeature::Generic->new(
 378   -primary_tag => 'feat5',
 379   -seq_id      => 'seq1',
 380   -start       => 11,
 381   -end         => 20
 382 );
 383
 384 my $feature6 = Bio::SeqFeature::Generic->new(
 385   -primary_tag => 'feat6',
 386   -seq_id      => 'seq1',
 387   -start       => 11,
 388   -end         => 25
 389 );
 390
 391 $seq_obj->add_SeqFeature($_)
 392     foreach ($composite_feat1, $feature1, $feature2, $feature3,
 393              $feature4, $feature5, $feature6);
 394
 395 my $coll = Bio::Annotation::Collection->new;
 396 $coll->add_Annotation(
 397   'comment', Bio::Annotation::Comment->new( '-text' => 'a comment on the whole sequence')
 398 );
 399 $seq_obj->annotation($coll);
 400
 401
 402 my $fragment_obj = Bio::Seq->new(
 403   -seq =>'atatatatat',
 404   -display_id => 'fragment1',
 405   -desc       => 'some fragment to insert'
 406 );
 407 my $frag_feature1 = Bio::SeqFeature::Generic->new(
 408   -primary_tag => 'frag_feat1',
 409   -seq_id      => 'fragment1',
 410   -start       => 2,
 411   -end         => 4,
 412   -strand      => -1,
 413 );
 414 $fragment_obj->add_SeqFeature( $frag_feature1 );
 415 my $frag_coll = Bio::Annotation::Collection->new;
 416 $frag_coll->add_Annotation(
 417   'comment', Bio::Annotation::Comment->new( '-text' => 'a comment on the fragment')
 418 );
 419 $fragment_obj->annotation($frag_coll);
 420
 421 # delete
 422 my $product;
 423 lives_ok(
 424   sub {
 425     $product = Bio::SeqUtils->delete( $seq_obj, 11, 20 );
 426   },
 427   "No error thrown when deleting a segment of the sequence"
 428 );
 429
 430 my ($seq_obj_comment) = $seq_obj->annotation->get_Annotations('comment');
 431 my ($product_comment) = $product->annotation->get_Annotations('comment');
 432 is( $seq_obj_comment, $product_comment, 'annotation of whole sequence has been moved to new molecule');
 433
 434 ok(
 435   grep ($_ eq 'deletion of 10bp',
 436     map ($_->get_tag_values('note'),
 437       grep ($_->primary_tag eq 'misc_feature', $product->get_SeqFeatures)
 438     )
 439   ),
 440   "the product has an additional 'misc_feature' and the note specifies the lengths of the deletion'"
 441 );
 442
 443 my ($composite_feat1_del) = grep ($_->primary_tag eq 'comp_feat1', $product->get_SeqFeatures);
 444 ok ($composite_feat1_del, "The composite feature is still present");
 445 isa_ok( $composite_feat1_del, 'Bio::SeqFeature::Generic');
 446 isa_ok( $composite_feat1_del->location, 'Bio::Location::Split', "a composite feature that spanned the deletion site has been split up, Location");
 447
 448 is( $composite_feat1_del->get_SeqFeatures, 2, 'one of the sub-eatures of the composite feature has been deleted completely');
 449 my ($subfeat1_del) = grep ($_->primary_tag eq 'sf1', $composite_feat1_del->get_SeqFeatures);
 450 ok ($subfeat1_del, "sub-feature 1 of the composite feature is still present");
 451 is ($subfeat1->end, 12, "the original end of sf1 is 12");
 452 is ($subfeat1_del->end, 10, "after deletion, the end of sf1 is 1nt before the deletion site");
 453 is ($subfeat1->location->end_pos_type, 'EXACT', 'the original end location of sf1 EXACT');
 454
 455 my ($subfeat1_comment) = $subfeat1->annotation->get_Annotations('comment');
 456 my ($subfeat1_del_comment) = $subfeat1_del->annotation->get_Annotations('comment');
 457 is( $subfeat1_comment, $subfeat1_del_comment, 'annotation of subeature 1 has been moved to new molecule');
 458
 459 my ($feature1_del) = grep ($_->primary_tag eq 'feat1', $product->get_SeqFeatures);
 460 ok ($feature1_del, "feature1 is till present");
 461 isa_ok( $feature1_del->location, 'Bio::Location::Split', 'feature1 location has now been split by the deletion and location object');
 462 is( my @feature1_del_sublocs = $feature1_del->location->each_Location, 2, 'feature1 has two locations after the deletion');
 463 is( $feature1_del_sublocs[0]->start, 2, 'feature1 start is unaffected by the deletion');
 464 is( $feature1_del_sublocs[0]->end, 10, 'feature1 end of first split is 1nt before deletion site');
 465 is( $feature1_del_sublocs[1]->start, 11, 'feature1 start of second split is 1nt after deletion site');
 466 is( $feature1_del_sublocs[1]->end, 15, 'feature1 end of second split has been adjusted correctly');
 467 my @fd1_notes = $feature1_del->get_tag_values('note');
 468 is( @fd1_notes,1, 'split feature now has a note');
 469 is (shift @fd1_notes, '10bp internal deletion between pos 10 and 11', 'got the expected note about length and position of deletion');
 470
 471 my ($feature3_del) = grep ($_->primary_tag eq 'feat3', $product->get_SeqFeatures);
 472 ok ($feature3_del, "feature3 is till present");
 473 is_deeply ( [$feature3_del->start, $feature3_del->end], [$feature3->start - 10, $feature3->end - 10], 'a feature downstream of the deletion site is shifted entirely by 10nt to the left');
 474
 475 my ($feature4_del) = grep ($_->primary_tag eq 'feat4', $product->get_SeqFeatures);
 476 ok ($feature4_del, "feature4 is till present");
 477 is_deeply ( [$feature4_del->start, $feature4_del->end], [$feature4->start, $feature4->end], 'a feature upstream of the deletion site is not repositioned by the deletion');
 478
 479 my ($feature2_del) = grep ($_->primary_tag eq 'feat2', $product->get_SeqFeatures);
 480 ok ($feature2_del, "feature2 is till present");
 481 is ( $feature2_del->start, 11, 'start pos of a feature that started in the deletion site has been altered accordingly');
 482 my @fd2_notes = $feature2_del->get_tag_values('note');
 483 is( @fd2_notes,1, 'feature 2 now has a note');
 484 is (shift @fd2_notes, "6bp deleted from feature 3' end", "note added to feature2 about deletion at 3' end");
 485
 486 ok (!grep ($_->primary_tag eq 'feat5', $product->get_SeqFeatures), 'a feature that was completely positioned inside the deletion site is not present on the new molecule');
 487
 488 my ($feature6_del) = grep ($_->primary_tag eq 'feat6', $product->get_SeqFeatures);
 489 ok ($feature6_del, "feature6 is till present");
 490 is ( $feature6_del->start, 11, 'start pos of a feature that started in the deletion site has been altered accordingly');
 491 is ( $feature6_del->end, 15, 'end pos of a feature that started in the deletion site has been altered accordingly');
 492
 493
 494 # insert
 495 lives_ok(
 496   sub {
 497     $product = Bio::SeqUtils->insert( $seq_obj, $fragment_obj, 10 );
 498   },
 499   "No error thrown when inserting a fragment into recipient sequence"
 500 );
 501 ($seq_obj_comment) = $seq_obj->annotation->get_Annotations('comment');
 502 ($product_comment) = $product->annotation->get_Annotations('comment');
 503 is( $seq_obj_comment, $product_comment, 'annotation of whole sequence has been moved to new molecule');
 504
 505 my ($composite_feat1_ins) = grep ($_->primary_tag eq 'comp_feat1', $product->get_SeqFeatures);
 506 ok ($composite_feat1_ins, "The composite feature is still present");
 507 isa_ok( $composite_feat1_ins, 'Bio::SeqFeature::Generic');
 508 isa_ok( $composite_feat1_ins->location, 'Bio::Location::Split', "a composite feature that spanned the insertion site has been split up, Location");
 509 is( $composite_feat1_ins->get_SeqFeatures, 3, 'all of the parts of the composite feature are still present');
 510
 511 my ($subfeat1_ins) = grep ($_->primary_tag eq 'sf1', $composite_feat1_ins->get_SeqFeatures);
 512 ok ($subfeat1_ins, "sub-feature 1 of the composite feature is still present");
 513 is ($subfeat1->end, 12, "the original end of sf1 is 12");
 514 is ($subfeat1_ins->end, $subfeat1->end + $fragment_obj->length, "after insertion, the end of sf1 has been shifted by the length of the insertion");
 515 isa_ok( $subfeat1_ins->location, 'Bio::Location::Split', 'sub-feature 1 (spans insertion site) is now split up and');
 516 is_deeply (
 517   [$subfeat1->location->end_pos_type, $subfeat1->location->start_pos_type],
 518   [$subfeat1_ins->location->end_pos_type, $subfeat1_ins->location->start_pos_type],
 519   'the start and end position types of sub-feature1 have not changed'
 520 );
 521 ($subfeat1_comment) = $subfeat1->annotation->get_Annotations('comment');
 522 my ($subfeat1_ins_comment) = $subfeat1_ins->annotation->get_Annotations('comment');
 523 is( $subfeat1_comment, $subfeat1_ins_comment, 'annotation of subeature 1 has been moved to new molecule');
 524 my @sf1ins_notes = $subfeat1_ins->get_tag_values('note');
 525 is( @sf1ins_notes,1, 'split feature now has a note');
 526 is (shift @sf1ins_notes, '10bp internal insertion between pos 10 and 21', 'got the expected note about length and position of insertion');
 527
 528 my ($feature3_ins) = grep ($_->primary_tag eq 'feat3', $product->get_SeqFeatures);
 529 ok ($feature3_ins, "feature3 is till present");
 530 is_deeply (
 531   [$feature3_ins->start, $feature3_ins->end],
 532   [$feature3->start + $fragment_obj->length, $feature3->end + $fragment_obj->length],
 533   'a feature downstream of the insertion site is shifted entirely to the left by the length of the insertion');
 534
 535 my ($feature4_ins) = grep ($_->primary_tag eq 'feat4', $product->get_SeqFeatures);
 536 ok ($feature4_ins, "feature4 is till present");
 537 is_deeply ( [$feature4_ins->start, $feature4_ins->end], [$feature4->start, $feature4->end], 'a feature upstream of the insertion site is not repositioned');
 538
 539 my ($frag_feature1_ins) = grep ($_->primary_tag eq 'frag_feat1', $product->get_SeqFeatures);
 540 ok( $frag_feature1_ins, 'a feature on the inserted fragment is present on the product molecule');
 541 is_deeply (
 542   [$frag_feature1_ins->start, $frag_feature1_ins->end],
 543   [12, 14],
 544   'position of the feature on the insert has been adjusted to product coordinates'
 545 );
 546 is( $frag_feature1_ins->strand, $frag_feature1->strand, 'strand of the feature on insert has not changed');
 547 like( $product->desc, qr/some fragment to insert/, 'desctription of the product contains description of the fragment');
 548 like( $product->desc, qr/some sequence for testing/, 'desctription of the product contains description of the recipient');
 549
 550 ok(
 551   grep ($_ eq 'inserted fragment',
 552     map ($_->get_tag_values('note'),
 553       grep ($_->primary_tag eq 'misc_feature', $product->get_SeqFeatures)
 554     )
 555   ),
 556   "the product has an additional 'misc_feature' with note='inserted fragment'"
 557 );
 558
 559 # ligate
 560 lives_ok(
 561   sub {
 562     $product = Bio::SeqUtils->ligate(
 563       -recipient => $seq_obj,
 564       -fragment  => $fragment_obj,
 565       -left      => 10,
 566       -right     => 31,
 567       -flip      => 1
 568     );
 569   },
 570   "No error thrown using 'ligate' of fragment into recipient"
 571 );
 572
 573 is ($product->length, 30, 'product has the expected length');
 574 is ($product->subseq(11,20), 'atatatatat', 'the sequence of the fragment is inserted into the product');
 575
 576 my ($inserted_fragment_feature) = grep(
 577   grep($_ eq 'inserted fragment', $_->get_tag_values('note')),
 578   grep( $_->has_tag('note'), $product->get_SeqFeatures)
 579 );
 580
 581 ok($inserted_fragment_feature, 'we have a feature annotating the ligated fragment');
 582 is_deeply (
 583   [$inserted_fragment_feature->start, $inserted_fragment_feature->end],
 584   [11, 20],
 585   'coordinates of the feature annotating the ligated feature are correct'
 586 );
 587
 588 my ($fragment_feat_lig) = grep ($_->primary_tag eq 'frag_feat1', $product->get_SeqFeatures);
 589 ok( $fragment_feat_lig, 'the fragment feature1 is now a feature of the product');
 590 is_deeply( [$fragment_feat_lig->start, $fragment_feat_lig->end], [17,19], 'start and end of a feature on the fragment are correct after insertion with "flip" option');
 591
 592 # test clone_obj option (create new objects via clone not 'new')
 593 my $foo_seq_obj = Bio::Seq::Foo->new(
 594   -seq =>'aaaaaaaaaaccccccccccggggggggggtttttttttt',
 595   -display_id => 'seq1',
 596   -desc       => 'some sequence for testing'
 597 );
 598 for ($composite_feat1, $feature1, $feature2, $feature3, $feature4, $feature5) {
 599     $foo_seq_obj->add_SeqFeature( $_ );
 600 }
 601 $foo_seq_obj->annotation($coll);
 602
 603 dies_ok(
 604   sub {
 605     $product = Bio::SeqUtils->delete( $foo_seq_obj, 11, 20, { clone_obj => 0} );
 606   },
 607   "Trying to delete from an object of a custom Bio::Seq subclass that doesn't allow calling 'new' throws an error"
 608 );
 609
 610 lives_ok(
 611   sub {
 612     $product = Bio::SeqUtils->delete( $foo_seq_obj, 11, 20, { clone_obj => 1} );
 613   },
 614   "Deleting from Bio::Seq::Foo does not throw an error when using the 'clone_obj' option to clone instead of calling 'new'"
 615 );
 616
 617 isa_ok( $product, 'Bio::Seq::Foo');
 618
 619 # just repeat some of the tests for the cloned feature
 620 ok(
 621   grep ($_ eq 'deletion of 10bp',
 622     map ($_->get_tag_values('note'),
 623       grep ($_->primary_tag eq 'misc_feature', $product->get_SeqFeatures)
 624     )
 625   ),
 626   "the product has an additional 'misc_feature' and the note specifies the lengths of the deletion'"
 627 );
 628 ($composite_feat1_del) = grep ($_->primary_tag eq 'comp_feat1', $product->get_SeqFeatures);
 629 ok ($composite_feat1_del, "The composite feature is still present");
 630 isa_ok( $composite_feat1_del, 'Bio::SeqFeature::Generic');
 631 isa_ok( $composite_feat1_del->location, 'Bio::Location::Split', "a composite feature that spanned the deletion site has been split up, Location");
 632
 633 # ligate with clone_obj
 634 dies_ok(
 635   sub {
 636     $product = Bio::SeqUtils->ligate(
 637       -recipient => $foo_seq_obj,
 638       -fragment  => $fragment_obj,
 639       -left      => 10,
 640       -right     => 31,
 641       -flip      => 1
 642     );
 643   },
 644   "'ligate' without clone_obj option dies with a Bio::Seq::Foo object that can't call new"
 645 );
 646
 647 lives_ok(
 648   sub {
 649     $product = Bio::SeqUtils->ligate(
 650       -recipient => $foo_seq_obj,
 651       -fragment  => $fragment_obj,
 652       -left      => 10,
 653       -right     => 31,
 654       -flip      => 1,
 655       -clone_obj => 1,
 656     );
 657   },
 658   "'ligate' with clone_obj option works with a Bio::Seq::Foo object that can't call new"
 659 );
 660
 661 sub uniq_sort {
 662     my @args = @_;
 663     my %uniq;
 664     @args = sort @args;
 665     @uniq{@args} = (0..$#args);
 666     return sort {$uniq{$a} <=> $uniq{$b}} keys %uniq;
 667 }
 668
 669 package Bio::Seq::Foo;
 670
 671 use base 'Bio::Seq';
 672
 673 sub can_call_new { 0 }