t/SeqTools/SeqUtils.t

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