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 $seq_obj->add_SeqFeature( $composite_feat1, $feature1, $feature2, $feature3, $feature4, $feature5, $feature6);
 391
 392 my $coll = Bio::Annotation::Collection->new;
 393 $coll->add_Annotation(
 394   'comment', Bio::Annotation::Comment->new( '-text' => 'a comment on the whole sequence')
 395 );
 396 $seq_obj->annotation($coll);
 397
 398
 399 my $fragment_obj = Bio::Seq->new(
 400   -seq =>'atatatatat',
 401   -display_id => 'fragment1',
 402   -desc       => 'some fragment to insert'
 403 );
 404 my $frag_feature1 = Bio::SeqFeature::Generic->new(
 405   -primary_tag => 'frag_feat1',
 406   -seq_id      => 'fragment1',
 407   -start       => 2,
 408   -end         => 4,
 409   -strand      => -1,
 410 );
 411 $fragment_obj->add_SeqFeature( $frag_feature1 );
 412 my $frag_coll = Bio::Annotation::Collection->new;
 413 $frag_coll->add_Annotation(
 414   'comment', Bio::Annotation::Comment->new( '-text' => 'a comment on the fragment')
 415 );
 416 $fragment_obj->annotation($frag_coll);
 417
 418 # delete
 419 my $product;
 420 lives_ok(
 421   sub {
 422     $product = Bio::SeqUtils->delete( $seq_obj, 11, 20 );
 423   },
 424   "No error thrown when deleting a segment of the sequence"
 425 );
 426
 427 my ($seq_obj_comment) = $seq_obj->annotation->get_Annotations('comment');
 428 my ($product_comment) = $product->annotation->get_Annotations('comment');
 429 is( $seq_obj_comment, $product_comment, 'annotation of whole sequence has been moved to new molecule');
 430
 431 ok(
 432   grep ($_ eq 'deletion of 10bp',
 433     map ($_->get_tag_values('note'),
 434       grep ($_->primary_tag eq 'misc_feature', $product->get_SeqFeatures)
 435     )
 436   ),
 437   "the product has an additional 'misc_feature' and the note specifies the lengths of the deletion'"
 438 );
 439
 440 my ($composite_feat1_del) = grep ($_->primary_tag eq 'comp_feat1', $product->get_SeqFeatures);
 441 ok ($composite_feat1_del, "The composite feature is still present");
 442 isa_ok( $composite_feat1_del, 'Bio::SeqFeature::Generic');
 443 isa_ok( $composite_feat1_del->location, 'Bio::Location::Split', "a composite feature that spanned the deletion site has been split up, Location");
 444
 445 is( $composite_feat1_del->get_SeqFeatures, 2, 'one of the sub-eatures of the composite feature has been deleted completely');
 446 my ($subfeat1_del) = grep ($_->primary_tag eq 'sf1', $composite_feat1_del->get_SeqFeatures);
 447 ok ($subfeat1_del, "sub-feature 1 of the composite feature is still present");
 448 is ($subfeat1->end, 12, "the original end of sf1 is 12");
 449 is ($subfeat1_del->end, 10, "after deletion, the end of sf1 is 1nt before the deletion site");
 450 is ($subfeat1->location->end_pos_type, 'EXACT', 'the original end location of sf1 EXACT');
 451
 452 my ($subfeat1_comment) = $subfeat1->annotation->get_Annotations('comment');
 453 my ($subfeat1_del_comment) = $subfeat1_del->annotation->get_Annotations('comment');
 454 is( $subfeat1_comment, $subfeat1_del_comment, 'annotation of subeature 1 has been moved to new molecule');
 455
 456 my ($feature1_del) = grep ($_->primary_tag eq 'feat1', $product->get_SeqFeatures);
 457 ok ($feature1_del, "feature1 is till present");
 458 isa_ok( $feature1_del->location, 'Bio::Location::Split', 'feature1 location has now been split by the deletion and location object');
 459 is( my @feature1_del_sublocs = $feature1_del->location->each_Location, 2, 'feature1 has two locations after the deletion');
 460 is( $feature1_del_sublocs[0]->start, 2, 'feature1 start is unaffected by the deletion');
 461 is( $feature1_del_sublocs[0]->end, 10, 'feature1 end of first split is 1nt before deletion site');
 462 is( $feature1_del_sublocs[1]->start, 11, 'feature1 start of second split is 1nt after deletion site');
 463 is( $feature1_del_sublocs[1]->end, 15, 'feature1 end of second split has been adjusted correctly');
 464 my @fd1_notes = $feature1_del->get_tag_values('note');
 465 is( @fd1_notes,1, 'split feature now has a note');
 466 is (shift @fd1_notes, '10bp internal deletion between pos 10 and 11', 'got the expected note about length and position of deletion');
 467
 468 my ($feature3_del) = grep ($_->primary_tag eq 'feat3', $product->get_SeqFeatures);
 469 ok ($feature3_del, "feature3 is till present");
 470 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');
 471
 472 my ($feature4_del) = grep ($_->primary_tag eq 'feat4', $product->get_SeqFeatures);
 473 ok ($feature4_del, "feature4 is till present");
 474 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');
 475
 476 my ($feature2_del) = grep ($_->primary_tag eq 'feat2', $product->get_SeqFeatures);
 477 ok ($feature2_del, "feature2 is till present");
 478 is ( $feature2_del->start, 11, 'start pos of a feature that started in the deletion site has been altered accordingly');
 479 my @fd2_notes = $feature2_del->get_tag_values('note');
 480 is( @fd2_notes,1, 'feature 2 now has a note');
 481 is (shift @fd2_notes, "6bp deleted from feature 3' end", "note added to feature2 about deletion at 3' end");
 482
 483 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');
 484
 485 my ($feature6_del) = grep ($_->primary_tag eq 'feat6', $product->get_SeqFeatures);
 486 ok ($feature6_del, "feature6 is till present");
 487 is ( $feature6_del->start, 11, 'start pos of a feature that started in the deletion site has been altered accordingly');
 488 is ( $feature6_del->end, 15, 'end pos of a feature that started in the deletion site has been altered accordingly');
 489
 490
 491 # insert
 492 lives_ok(
 493   sub {
 494     $product = Bio::SeqUtils->insert( $seq_obj, $fragment_obj, 10 );
 495   },
 496   "No error thrown when inserting a fragment into recipient sequence"
 497 );
 498 ($seq_obj_comment) = $seq_obj->annotation->get_Annotations('comment');
 499 ($product_comment) = $product->annotation->get_Annotations('comment');
 500 is( $seq_obj_comment, $product_comment, 'annotation of whole sequence has been moved to new molecule');
 501
 502 my ($composite_feat1_ins) = grep ($_->primary_tag eq 'comp_feat1', $product->get_SeqFeatures);
 503 ok ($composite_feat1_ins, "The composite feature is still present");
 504 isa_ok( $composite_feat1_ins, 'Bio::SeqFeature::Generic');
 505 isa_ok( $composite_feat1_ins->location, 'Bio::Location::Split', "a composite feature that spanned the insertion site has been split up, Location");
 506 is( $composite_feat1_ins->get_SeqFeatures, 3, 'all of the parts of the composite feature are still present');
 507
 508 my ($subfeat1_ins) = grep ($_->primary_tag eq 'sf1', $composite_feat1_ins->get_SeqFeatures);
 509 ok ($subfeat1_ins, "sub-feature 1 of the composite feature is still present");
 510 is ($subfeat1->end, 12, "the original end of sf1 is 12");
 511 is ($subfeat1_ins->end, $subfeat1->end + $fragment_obj->length, "after insertion, the end of sf1 has been shifted by the length of the insertion");
 512 isa_ok( $subfeat1_ins->location, 'Bio::Location::Split', 'sub-feature 1 (spans insertion site) is now split up and');
 513 is_deeply (
 514   [$subfeat1->location->end_pos_type, $subfeat1->location->start_pos_type],
 515   [$subfeat1_ins->location->end_pos_type, $subfeat1_ins->location->start_pos_type],
 516   'the start and end position types of sub-feature1 have not changed'
 517 );
 518 ($subfeat1_comment) = $subfeat1->annotation->get_Annotations('comment');
 519 my ($subfeat1_ins_comment) = $subfeat1_ins->annotation->get_Annotations('comment');
 520 is( $subfeat1_comment, $subfeat1_ins_comment, 'annotation of subeature 1 has been moved to new molecule');
 521 my @sf1ins_notes = $subfeat1_ins->get_tag_values('note');
 522 is( @sf1ins_notes,1, 'split feature now has a note');
 523 is (shift @sf1ins_notes, '10bp internal insertion between pos 10 and 21', 'got the expected note about length and position of insertion');
 524
 525 my ($feature3_ins) = grep ($_->primary_tag eq 'feat3', $product->get_SeqFeatures);
 526 ok ($feature3_ins, "feature3 is till present");
 527 is_deeply (
 528   [$feature3_ins->start, $feature3_ins->end],
 529   [$feature3->start + $fragment_obj->length, $feature3->end + $fragment_obj->length],
 530   'a feature downstream of the insertion site is shifted entirely to the left by the length of the insertion');
 531
 532 my ($feature4_ins) = grep ($_->primary_tag eq 'feat4', $product->get_SeqFeatures);
 533 ok ($feature4_ins, "feature4 is till present");
 534 is_deeply ( [$feature4_ins->start, $feature4_ins->end], [$feature4->start, $feature4->end], 'a feature upstream of the insertion site is not repositioned');
 535
 536 my ($frag_feature1_ins) = grep ($_->primary_tag eq 'frag_feat1', $product->get_SeqFeatures);
 537 ok( $frag_feature1_ins, 'a feature on the inserted fragment is present on the product molecule');
 538 is_deeply (
 539   [$frag_feature1_ins->start, $frag_feature1_ins->end],
 540   [12, 14],
 541   'position of the feature on the insert has been adjusted to product coordinates'
 542 );
 543 is( $frag_feature1_ins->strand, $frag_feature1->strand, 'strand of the feature on insert has not changed');
 544 like( $product->desc, qr/some fragment to insert/, 'desctription of the product contains description of the fragment');
 545 like( $product->desc, qr/some sequence for testing/, 'desctription of the product contains description of the recipient');
 546
 547 ok(
 548   grep ($_ eq 'inserted fragment',
 549     map ($_->get_tag_values('note'),
 550       grep ($_->primary_tag eq 'misc_feature', $product->get_SeqFeatures)
 551     )
 552   ),
 553   "the product has an additional 'misc_feature' with note='inserted fragment'"
 554 );
 555
 556 # ligate
 557 lives_ok(
 558   sub {
 559     $product = Bio::SeqUtils->ligate(
 560       -recipient => $seq_obj,
 561       -fragment  => $fragment_obj,
 562       -left      => 10,
 563       -right     => 31,
 564       -flip      => 1
 565     );
 566   },
 567   "No error thrown using 'ligate' of fragment into recipient"
 568 );
 569
 570 is ($product->length, 30, 'product has the expected length');
 571 is ($product->subseq(11,20), 'atatatatat', 'the sequence of the fragment is inserted into the product');
 572
 573 my ($inserted_fragment_feature) = grep(
 574   grep($_ eq 'inserted fragment', $_->get_tag_values('note')),
 575   grep( $_->has_tag('note'), $product->get_SeqFeatures)
 576 );
 577
 578 ok($inserted_fragment_feature, 'we have a feature annotating the ligated fragment');
 579 is_deeply (
 580   [$inserted_fragment_feature->start, $inserted_fragment_feature->end],
 581   [11, 20],
 582   'coordinates of the feature annotating the ligated feature are correct'
 583 );
 584
 585 my ($fragment_feat_lig) = grep ($_->primary_tag eq 'frag_feat1', $product->get_SeqFeatures);
 586 ok( $fragment_feat_lig, 'the fragment feature1 is now a feature of the product');
 587 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');
 588
 589 # test clone_obj option (create new objects via clone not 'new')
 590 my $foo_seq_obj = Bio::Seq::Foo->new(
 591   -seq =>'aaaaaaaaaaccccccccccggggggggggtttttttttt',
 592   -display_id => 'seq1',
 593   -desc       => 'some sequence for testing'
 594 );
 595 for ($composite_feat1, $feature1, $feature2, $feature3, $feature4, $feature5) {
 596     $foo_seq_obj->add_SeqFeature( $_ );
 597 }
 598 $foo_seq_obj->annotation($coll);
 599
 600 dies_ok(
 601   sub {
 602     $product = Bio::SeqUtils->delete( $foo_seq_obj, 11, 20, { clone_obj => 0} );
 603   },
 604   "Trying to delete from an object of a custom Bio::Seq subclass that doesn't allow calling 'new' throws an error"
 605 );
 606
 607 lives_ok(
 608   sub {
 609     $product = Bio::SeqUtils->delete( $foo_seq_obj, 11, 20, { clone_obj => 1} );
 610   },
 611   "Deleting from Bio::Seq::Foo does not throw an error when using the 'clone_obj' option to clone instead of calling 'new'"
 612 );
 613
 614 isa_ok( $product, 'Bio::Seq::Foo');
 615
 616 # just repeat some of the tests for the cloned feature
 617 ok(
 618   grep ($_ eq 'deletion of 10bp',
 619     map ($_->get_tag_values('note'),
 620       grep ($_->primary_tag eq 'misc_feature', $product->get_SeqFeatures)
 621     )
 622   ),
 623   "the product has an additional 'misc_feature' and the note specifies the lengths of the deletion'"
 624 );
 625 ($composite_feat1_del) = grep ($_->primary_tag eq 'comp_feat1', $product->get_SeqFeatures);
 626 ok ($composite_feat1_del, "The composite feature is still present");
 627 isa_ok( $composite_feat1_del, 'Bio::SeqFeature::Generic');
 628 isa_ok( $composite_feat1_del->location, 'Bio::Location::Split', "a composite feature that spanned the deletion site has been split up, Location");
 629
 630 # ligate with clone_obj
 631 dies_ok(
 632   sub {
 633     $product = Bio::SeqUtils->ligate(
 634       -recipient => $foo_seq_obj,
 635       -fragment  => $fragment_obj,
 636       -left      => 10,
 637       -right     => 31,
 638       -flip      => 1
 639     );
 640   },
 641   "'ligate' without clone_obj option dies with a Bio::Seq::Foo object that can't call new"
 642 );
 643
 644 lives_ok(
 645   sub {
 646     $product = Bio::SeqUtils->ligate(
 647       -recipient => $foo_seq_obj,
 648       -fragment  => $fragment_obj,
 649       -left      => 10,
 650       -right     => 31,
 651       -flip      => 1,
 652       -clone_obj => 1,
 653     );
 654   },
 655   "'ligate' with clone_obj option works with a Bio::Seq::Foo object that can't call new"
 656 );
 657
 658 sub uniq_sort {
 659     my @args = @_;
 660     my %uniq;
 661     @args = sort @args;
 662     @uniq{@args} = (0..$#args);
 663     return sort {$uniq{$a} <=> $uniq{$b}} keys %uniq;
 664 }
 665
 666 package Bio::Seq::Foo;
 667
 668 use base 'Bio::Seq';
 669
 670 sub can_call_new { 0 }