Bio/Assembly/Scaffold.pm

   1 #
   2 # BioPerl module for Bio::Assembly::Scaffold
   3 #
   4 # Copyright by Robson F. de Souza
   5 #
   6 # You may distribute this module under the same terms as perl itself
   7 #
   8 # POD documentation - main docs before the code
   9
  10 =head1 NAME
  11
  12 Bio::Assembly::Scaffold - Perl module to hold and manipulate sequence assembly
  13 data.
  14
  15 =head1 SYNOPSIS
  16 #
  17     # Module loading
  18     use Bio::Assembly::IO;
  19
  20     # Assembly loading methods
  21     my $aio = Bio::Assembly::IO->new(-file=>"test.ace.1", -format=>'phrap');
  22     my $assembly = $aio->next_assembly;
  23
  24     foreach my $contig ($assembly->all_contigs) {
  25         # do something... (see Bio::Assembly::Contig)
  26     }
  27
  28 =head1 DESCRIPTION
  29
  30 Bio::Assembly::Scaffold was developed to store and manipulate data
  31 from sequence assembly programs like Phrap. It implements the
  32 ScaffoldI interface and intends to be generic enough to be used by
  33 Bio::Assembly::IO drivers written to programs other than Phrap.
  34
  35 =head1 FEEDBACK
  36
  37 =head2 Mailing Lists
  38
  39 User feedback is an integral part of the evolution of this and other
  40 Bioperl modules. Send your comments and suggestions preferably to the
  41 Bioperl mailing lists  Your participation is much appreciated.
  42
  43   bioperl-l@bioperl.org                  - General discussion
  44   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  45
  46 =head2 Support
  47
  48 Please direct usage questions or support issues to the mailing list:
  49
  50 I<bioperl-l@bioperl.org>
  51
  52 rather than to the module maintainer directly. Many experienced and
  53 reponsive experts will be able look at the problem and quickly
  54 address it. Please include a thorough description of the problem
  55 with code and data examples if at all possible.
  56
  57 =head2 Reporting Bugs
  58
  59 Report bugs to the Bioperl bug tracking system to help us keep track
  60 the bugs and their resolution.  Bug reports can be submitted via the
  61 web:
  62
  63   https://github.com/bioperl/bioperl-live/issues
  64
  65 =head1 AUTHOR - Robson Francisco de Souza
  66
  67 rfsouza@citri.iq.usp.br
  68
  69 =head1 APPENDIX
  70
  71 The rest of the documentation details each of the object
  72 methods. Internal methods are usually preceded with a _
  73
  74 =cut
  75
  76 package Bio::Assembly::Scaffold;
  77
  78 use strict;
  79 use Bio::Annotation::Collection;
  80
  81 use base qw(Bio::Root::Root Bio::Assembly::ScaffoldI);
  82
  83 =head2 new ()
  84
  85     Title   : new
  86     Usage   : $scaffold = new ( -id       => "assembly 1",
  87                                 -source   => 'program_name',
  88                                 -contigs  => \@contigs,
  89                                 -singlets => \@singlets );
  90     Function: creates a new scaffold object
  91     Returns : Bio::Assembly::Scaffold
  92     Args    : -id       : [string] scaffold name
  93               -source   : [string] sequence assembly program
  94               -contigs  : reference to array of Bio::Assembly::Contig objects
  95               -singlets : reference to array of Bio::Assembly::Singlet objects
  96
  97
  98 =cut
  99
 100 sub new {
 101   my($class, @args) = @_;
 102   my $self = $class->SUPER::new(@args);
 103   my ($id, $src, $contigs, $singlets) = $self->_rearrange(
 104       [qw(ID SOURCE CONTIGS SINGLETS)], @args);
 105
 106   # Scaffold defaults
 107   $self->{'_id'} = 'NoName';
 108   $self->{'_source'} = 'Unknown';
 109   $self->{'_contigs'} = {};
 110   $self->{'_singlets'} = {};
 111   $self->{'_seqs'} = {};
 112   $self->{'_annotation'} = Bio::Annotation::Collection->new();
 113
 114   # Import manual info
 115   $self->{'_id'} = $id if (defined $id);
 116   $self->{'_source'} = $src if (defined $src);
 117
 118   # Add contigs and singlets to scaffold
 119   if (defined $contigs && ref($contigs = 'ARRAY')) {
 120     for my $contig (@{$contigs}) {
 121       if( ! ref $contig || ! $contig->isa('Bio::Assembly::Contig') ) {
 122         $self->throw("Bio::Assembly::Scaffold::new is unable to process non".
 123           "Bio::Assembly::Contig object [", ref($contig), "]");
 124       }
 125       $self->add_contig($contig);
 126     }
 127   }
 128   if (defined $singlets && ref($singlets = 'ARRAY')) {
 129     for my $singlet (@{$singlets}) {
 130       if( ! ref $singlet || ! $singlet->isa('Bio::Assembly::Singlet') ) {
 131         $self->throw("Bio::Assembly::Scaffold::new is unable to process non".
 132           "Bio::Assembly::Singlet object [", ref($singlet), "]");
 133       }
 134       $self->add_singlet($singlet);
 135     }
 136   }
 137
 138   return $self;
 139 }
 140
 141 =head1 Accessing general assembly data
 142
 143 =cut
 144
 145 =head2 id
 146
 147     Title   : id
 148     Usage   : $assembly->id()
 149     Function: Get/Set assembly ID
 150     Returns : string or undef
 151     Args    : string
 152
 153 =cut
 154
 155 sub id {
 156     my ($self, $id) = @_;
 157     return $self->{'_id'} = $id if (defined $id);
 158     return $self->{'_id'};
 159 }
 160
 161 =head2 annotation
 162
 163     Title   : annotation
 164     Usage   : $assembly->annotation()
 165     Function: Get/Set assembly annotation object
 166     Returns : Bio::Annotation::Collection
 167     Args    : none
 168
 169 =cut
 170
 171 sub annotation {
 172     my ($self, $ref) = shift;
 173     $self->{'_annotation'} = $ref if (defined $ref);
 174     return $self->{'_annotation'};
 175 }
 176
 177 =head2 get_nof_contigs
 178
 179     Title   : get_nof_contigs
 180     Usage   : $assembly->get_nof_contigs()
 181     Function: Get the number of contigs included in the scaffold
 182     Returns : integer
 183     Args    : none
 184
 185 =cut
 186
 187 sub get_nof_contigs {
 188     my $self = shift;
 189     return scalar( $self->get_contig_ids() );
 190 }
 191
 192 =head2 get_nof_contig_seqs
 193
 194     Title   : get_nof_contig_seqs
 195     Usage   : $assembly->get_nof_contig_seqs()
 196     Function: Get the number of sequences included in contigs of the
 197               scaffold (no consensus sequences or singlets)
 198     Returns : integer
 199     Args    : none
 200
 201 =cut
 202
 203 sub get_nof_contig_seqs {
 204     my $self = shift;
 205
 206     my $nof_seqs = 0;
 207     foreach my $contig ($self->all_contigs) {
 208         $nof_seqs += scalar( $contig->get_seq_ids() );
 209     }
 210
 211     return $nof_seqs;
 212 }
 213 # function alias for backward compatibility
 214 *get_nof_sequences_in_contigs = \&get_nof_contig_seqs;
 215
 216
 217 =head2 get_nof_singlets (get_nof_singlet_seqs)
 218
 219     Title   : nof_singlets
 220     Usage   : $assembly->nof_singlets()
 221     Function: Get the number of singlets included in the assembly
 222     Returns : integer
 223     Args    : none
 224
 225 =cut
 226
 227 sub get_nof_singlets {
 228     my $self = shift;
 229     return scalar( $self->get_singlet_ids() );
 230 }
 231 *get_nof_singlet_seqs = \&get_nof_singlets;
 232
 233 =head2 get_all_seq_ids
 234
 235     Title   : get_all_seq_ids
 236     Usage   : $assembly->get_all_seq_ids()
 237     Function: Get the ID of all sequences making up the scaffold
 238               (sequences from contigs and singlets, not consensus).
 239     Returns : array of strings
 240     Args    : none
 241
 242 =cut
 243
 244 sub get_all_seq_ids {
 245     my $self = shift;
 246     return keys %{ $self->{'_seqs'} };
 247 }
 248
 249 =head2 get_nof_seqs
 250
 251     Title   : get_nof_seqs
 252     Usage   : $assembly->get_nof_seqs()
 253     Function: Get total number of sequences making up the scaffold
 254               (sequences from contigs and singlets, not consensus).
 255     Returns : integer
 256     Args    : none
 257
 258 =cut
 259
 260 sub get_nof_seqs {
 261     my $self = shift;
 262     return scalar $self->get_all_seq_ids;
 263 }
 264
 265 =head2 get_contig_seq_ids
 266
 267     Title   : get_contig_seq_ids
 268     Usage   : $assembly->get_contig_seq_ids()
 269     Function: Get the ID of all sequences in contigs
 270     Returns : array of strings
 271     Args    : none
 272
 273 =cut
 274
 275 sub get_contig_seq_ids {
 276     my $self = shift;
 277     my @ids;
 278     for my $contig ( $self->all_contigs ) {
 279         push @ids, $contig->get_seq_ids;
 280     }
 281     return @ids;
 282 }
 283 # function alias for backward compatibility
 284 *get_seq_ids = \&get_contig_seq_ids;
 285
 286 =head2 get_contig_ids
 287
 288     Title   : get_contig_ids
 289     Usage   : $assembly->get_contig_ids()
 290     Function: Access list of contig IDs from assembly
 291     Returns : an array, if there are any contigs in the
 292               assembly. An empty array otherwise
 293     Args    : none
 294
 295 =cut
 296
 297 sub get_contig_ids {
 298     my $self = shift;
 299
 300     return wantarray
 301         ? sort keys %{$self->{'_contigs'}}
 302         : scalar keys %{$self->{'_contigs'}};
 303 }
 304
 305 =head2 get_singlet_ids (get_singlet_seq_ids)
 306
 307     Title   : get_singlet_ids
 308     Usage   : $assembly->get_singlet_ids()
 309     Function: Access list of singlet IDs from assembly
 310     Returns : array of strings if there are any singlets
 311               otherwise an empty array
 312     Args    : none
 313
 314 =cut
 315
 316 sub get_singlet_ids {
 317     my $self = shift;
 318
 319     return wantarray
 320         ? sort keys %{$self->{'_singlets'}}
 321         : scalar keys %{$self->{'_singlets'}};
 322 }
 323 *get_singlet_seq_ids = \&get_singlet_ids;
 324
 325 =head2 get_seq_by_id
 326
 327     Title   : get_seq_by_id
 328     Usage   : $assembly->get_seq_by_id($id)
 329     Function: Get a reference for an sequence making up the scaffold
 330               (from a contig or singlet, not consensus)
 331     Returns : a Bio::LocatableSeq object
 332               undef if sequence $id is not found in the scaffold
 333     Args    : [string] sequence identifier (id)
 334
 335 =cut
 336
 337 sub get_seq_by_id {
 338     my $self = shift;
 339     my $seqID = shift;
 340
 341     return unless (exists $self->{'_seqs'}{$seqID});
 342
 343     return $self->{'_seqs'}{$seqID}->get_seq_by_name($seqID);
 344 }
 345
 346 =head2 get_contig_by_id
 347
 348     Title   : get_contig_by_id
 349     Usage   : $assembly->get_contig_by_id($id)
 350     Function: Get a reference for a contig
 351     Returns : a Bio::Assembly::Contig object or undef
 352     Args    : [string] contig unique identifier (ID)
 353
 354 =cut
 355
 356 sub get_contig_by_id {
 357     my $self = shift;
 358     my $contigID = shift;
 359
 360     return unless (exists $self->{'_contigs'}{$contigID});
 361
 362     return $self->{'_contigs'}{$contigID};
 363 }
 364
 365 =head2 get_singlet_by_id
 366
 367     Title   : get_singlet_by_id
 368     Usage   : $assembly->get_singlet_by_id()
 369     Function: Get a reference for a singlet
 370     Returns : Bio::Assembly::Singlet object or undef
 371     Args    : [string] a singlet ID
 372
 373 =cut
 374
 375 sub get_singlet_by_id {
 376     my $self = shift;
 377
 378     my $singletID = shift;
 379
 380     return unless (exists $self->{'_singlets'}{$singletID});
 381
 382     return $self->{'_singlets'}{$singletID};
 383 }
 384
 385 =head1 Modifier methods
 386
 387 =cut
 388
 389 =head2 add_contig
 390
 391     Title   : add_contig
 392     Usage   : $assembly->add_contig($contig)
 393     Function: Add a contig to the assembly
 394     Returns : 1 on success
 395     Args    : a Bio::Assembly::Contig object
 396           order (optional)
 397
 398 =cut
 399
 400 sub add_contig {
 401     my ($self, $contig) = @_;
 402
 403     # Input check
 404     if( !ref $contig || ! $contig->isa('Bio::Assembly::Contig') ) {
 405         $self->throw("Bio::Assembly::Scaffold::add_contig is unable to process".
 406             " non Bio::Assembly::Contig object [", ref($contig), "]");
 407     }
 408
 409     # Create and attribute contig ID
 410     my $contigID  = $contig->id();
 411     if( !defined $contigID ) {
 412         $contigID = 'Unknown_' . ($self->get_nof_contigs() + 1);
 413         $contig->id($contigID);
 414         $self->warn("Attributing ID $contigID to unnamed Bio::Assembly::Contig".
 415             " object.");
 416     }
 417
 418     # Adding contig to scaffold
 419     $self->warn("Replacing contig $contigID with a new contig object")
 420         if (exists $self->{'_contigs'}{$contigID});
 421     $self->{'_contigs'}{$contigID} = $contig;
 422     $contig->assembly($self); # weak circular reference
 423
 424     # Put contig sequences in the list of sequences belonging to the scaffold
 425     foreach my $seqID ($contig->get_seq_ids()) {
 426         if (exists $self->{'_seqs'}{$seqID} &&
 427             not($self->{'_seqs'}{$seqID} eq $contig) ) {
 428             $self->warn( "Sequence $seqID already assigned to object ".
 429                 $self->{'_seqs'}{$seqID}->id().". Moving to contig $contigID");
 430         }
 431         $self->{'_seqs'}{$seqID} = $contig;
 432     }
 433
 434     return 1;
 435 }
 436
 437 =head2 add_singlet
 438
 439     Title   : add_singlet
 440     Usage   : $assembly->add_singlet($seq)
 441     Function: Add a singlet to the assembly
 442     Returns : 1 on success
 443     Args    : a Bio::Assembly::Singlet object
 444               order (optional)
 445
 446 =cut
 447
 448 sub add_singlet {
 449     my ($self, $singlet) = @_;
 450
 451     # Input check
 452     if ( !ref $singlet || ! $singlet->isa('Bio::Assembly::Singlet') ) {
 453         $self->throw("Bio::Assembly::Scaffold::add_singlet is unable to process".
 454             " non Bio::Assembly::Singlet object [", ref($singlet), "]");
 455     }
 456
 457     # Create and attribute singlet ID
 458     my $singletID = $singlet->id();
 459     if( !defined $singletID ) {
 460         $singletID = 'Unknown_' . ($self->get_nof_singlets() + 1);
 461         $singlet->id($singletID);
 462         $self->warn("Attributing ID $singletID to unnamed Bio::Assembly::".
 463             "Singlet object.");
 464     }
 465
 466     # Adding singlet to scaffold
 467     $self->warn("Replacing singlet $singletID with a new singlet object")
 468         if (exists $self->{'_singlets'}{$singletID});
 469     $self->{'_singlets'}{$singletID} = $singlet;
 470     $singlet->assembly($self); # weak circular reference
 471
 472     # Put singlet sequence in the list of sequences belonging to the scaffold
 473     my $seqID = $singlet->seqref->id();
 474     if (exists $self->{'_seqs'}{$seqID} &&
 475         not($self->{'_seqs'}{$seqID} eq $singlet) ) {
 476         $self->warn( "Sequence $seqID already assigned to object ".
 477             $self->{'_seqs'}{$seqID}->id().". Moving to singlet $singletID");
 478     }
 479     $self->{'_seqs'}{$seqID} = $singlet;
 480
 481     return 1;
 482 }
 483
 484 =head2 update_seq_list
 485
 486     Title   : update_seq_list
 487     Usage   : $assembly->update_seq_list()
 488     Function:
 489
 490               Synchronizes the assembly registry for sequences in
 491               contigs and contig actual aligned sequences content. You
 492               probably want to run this after you remove/add a
 493               sequence from/to a contig in the assembly.
 494
 495     Returns : 1 for success
 496     Args    : none
 497
 498 =cut
 499
 500 sub update_seq_list {
 501     my $self = shift;
 502
 503     $self->{'_seqs'} = {};
 504
 505     # Put sequences in contigs in list of sequences belonging to the scaffold
 506     foreach my $contig ($self->all_contigs) {
 507         my $contigID = $contig->id();
 508         foreach my $seqID ($contig->get_seq_ids) {
 509             if (exists $self->{'_seqs'}{$seqID} &&
 510                 not($self->{'_seqs'}{$seqID} eq $contig) ) {
 511                 $self->warn( "Sequence $seqID already assigned to object ".
 512                     $self->{'_seqs'}{$seqID}->id().". Moving to contig $contigID");
 513             }
 514             $self->{'_seqs'}{$seqID} = $contig;
 515         }
 516     }
 517
 518     # Put singlet sequences in the list of sequences belonging to the scaffold
 519     foreach my $singlet ($self->all_singlets) {
 520         my $singletID = $singlet->id();
 521         my $seqID     = $singlet->seqref->id();
 522         if (exists $self->{'_seqs'}{$seqID} &&
 523             not($self->{'_seqs'}{$seqID} eq $singlet) ) {
 524             $self->warn( "Sequence $seqID already assigned to object ".
 525                 $self->{'_seqs'}{$seqID}->id().". Moving to singlet $singletID");
 526         }
 527         $self->{'_seqs'}{$seqID} = $singlet;
 528     }
 529
 530     return 1;
 531 }
 532
 533 =head2 remove_contigs
 534
 535     Title   : remove_contigs
 536     Usage   : $assembly->remove_contigs(1..4)
 537     Function: Remove contig from assembly object
 538     Returns : an array of removed Bio::Assembly::Contig
 539               objects
 540     Args    : an array of contig IDs
 541
 542     See function get_contig_ids() above
 543
 544 =cut
 545
 546 sub remove_contigs {
 547     my ($self, @args) = @_;
 548
 549     my @ret = ();
 550     foreach my $contigID (@args) {
 551         foreach my $seqID ($self->get_contig_by_id($contigID)->get_seq_ids()) {
 552             delete $self->{'_seqs'}{$seqID};
 553         }
 554         push(@ret, $self->{'_contigs'}{$contigID});
 555         delete $self->{'_contigs'}{$contigID};
 556     }
 557
 558     return @ret;
 559 }
 560
 561 =head2 remove_singlets
 562
 563     Title   : remove_singlets
 564     Usage   : $assembly->remove_singlets(@singlet_ids)
 565     Function: Remove singlet from assembly object
 566     Returns : the Bio::Assembly::Singlet objects removed
 567     Args    : a list of singlet IDs
 568
 569     See function get_singlet_ids() above
 570
 571 =cut
 572
 573 sub remove_singlets {
 574     my ($self,@args) = @_;
 575
 576     my @ret = ();
 577     foreach my $singletID (@args) {
 578         push(@ret,$self->{'_singlets'}{$singletID});
 579         delete $self->{'_singlets'}{$singletID};
 580     }
 581
 582     return @ret;
 583 }
 584
 585 =head2 remove_features_collection
 586
 587     Title   : remove_features_collection
 588     Usage   : $assembly->remove_features_collection()
 589     Function: Removes the collection of features associated to every
 590               contig and singlet of the scaffold. This can be useful
 591               to save some memory (when contig and singlet features are
 592               not needed).
 593     Returns   : none
 594     Argument  : none
 595
 596 =cut
 597
 598 sub remove_features_collection {
 599     my ($self) = @_;
 600     for my $obj ( $self->all_contigs, $self->all_singlets ) {
 601         $obj->remove_features_collection;
 602     }
 603     return;
 604 }
 605
 606
 607 =head1 Contig and singlet selection methods
 608
 609 =cut
 610
 611 =head2 select_contigs
 612
 613     Title   : select_contigs
 614     Usage   : $assembly->select_contigs(@list)
 615     Function: Select an array of contigs from the assembly
 616     Returns : an array of Bio::Assembly::Contig objects
 617     Args    : an array of contig ids
 618
 619     See function get_contig_ids() above
 620
 621 =cut
 622
 623 sub select_contigs {
 624     my ($self,@args) = @_;
 625
 626     my @contigs = ();
 627     foreach my $contig (@args) {
 628     unless (exists $self->{'_contigs'}{$contig}) {
 629         $self->warn("$contig contig not found. Ignoring...");
 630         next;
 631     }
 632     push(@contigs, $self->{'_contigs'}{$contig});
 633     }
 634
 635     return @contigs;
 636 }
 637
 638 =head2 select_singlets
 639
 640     Title   : select_singlets
 641     Usage   : $assembly->select_singlets(@list)
 642     Function: Selects an array of singlets from the assembly
 643     Returns : an array of Bio::Assembly::Singlet objects
 644     Args    : an array of singlet ids
 645
 646     See function get_singlet_ids() above
 647
 648 =cut
 649
 650 sub select_singlets {
 651     my ($self,@args) = @_;
 652
 653     my @singlets = ();
 654     foreach my $singlet (@args) {
 655     unless (exists $self->{'_singlets'}{$singlet}) {
 656         $self->warn("$singlet singlet not found. Ignoring...");
 657         next;
 658     }
 659     push(@singlets, $self->{'_singlets'}{$singlet});
 660     }
 661
 662     return @singlets;
 663 }
 664
 665 =head2 all_contigs
 666
 667     Title   : all_contigs
 668     Usage   : my @contigs = $assembly->all_contigs
 669     Function:
 670
 671               Returns a list of all contigs in this assembly.  Contigs
 672               are both clusters and alignments of one or more reads,
 673               with an associated consensus sequence.
 674
 675     Returns : array of Bio::Assembly::Contig (in lexical id order)
 676     Args    : none
 677
 678 =cut
 679
 680 sub all_contigs {
 681     my ($self) = @_;
 682
 683     my @contigs = ();
 684     foreach my $contig (sort { $a cmp $b } keys %{ $self->{'_contigs'} }) {
 685         push(@contigs, $self->{'_contigs'}{$contig});
 686     }
 687
 688     return @contigs;
 689 }
 690
 691 =head2 all_singlets
 692
 693     Title   : all_singlets
 694     Usage   : my @singlets = $assembly->all_singlets
 695     Function:
 696
 697               Returns a list of all singlets in this assembly.
 698           Singlets are isolated reads, without non-vector
 699           matches to any other read in the assembly.
 700
 701     Returns : array of Bio::Assembly::Singlet objects (in lexical order by id)
 702     Args    : none
 703
 704 =cut
 705
 706 sub all_singlets {
 707     my ($self) = @_;
 708
 709     my @singlets = ();
 710     foreach my $singlet (sort { $a cmp $b } keys %{ $self->{'_singlets'} }) {
 711     push(@singlets, $self->{'_singlets'}{$singlet});
 712     }
 713
 714     return @singlets;
 715 }
 716
 717
 718 1;