Archive Ensembl HomeArchive Ensembl Home
GenomicAlignBlock.pm
Go to the documentation of this file.
00001 =head1 LICENSE
00002 
00003   Copyright (c) 1999-2012 The European Bioinformatics Institute and
00004   Genome Research Limited.  All rights reserved.
00005 
00006   This software is distributed under a modified Apache license.
00007   For license details, please see
00008 
00009     http://www.ensembl.org/info/about/code_licence.html
00010 
00011 =head1 CONTACT
00012 
00013   Please email comments or questions to the public Ensembl
00014   developers list at <dev@ensembl.org>.
00015 
00016   Questions may also be sent to the Ensembl help desk at
00017   <helpdesk@ensembl.org>.
00018 
00019 =head1 NAME
00020 
00021 Bio::EnsEMBL::Compara::GenomicAlignBlock - Alignment of two or more pieces of genomic DNA
00022 
00023 =head1 SYNOPSIS
00024   
00025   use Bio::EnsEMBL::Compara::GenomicAlignBlock;
00026   
00027   my $genomic_align_block = new Bio::EnsEMBL::Compara::GenomicAlignBlock(
00028           -adaptor => $genomic_align_block_adaptor,
00029           -method_link_species_set => $method_link_species_set,
00030           -score => 56.2,
00031           -length => 1203,
00032           -genomic_align_array => [$genomic_align1, $genomic_align2...]
00033       );
00034 
00035 SET VALUES
00036   $genomic_align_block->adaptor($gen_ali_blk_adaptor);
00037   $genomic_align_block->dbID(12);
00038   $genomic_align_block->method_link_species_set($method_link_species_set);
00039   $genomic_align_block->reference_genomic_align_id(35123);
00040   $genomic_align_block->genomic_align_array([$genomic_align1, $genomic_align2]);
00041   $genomic_align_block->reference_slice($reference_slice);
00042   $genomic_align_block->reference_slice_start(1035);
00043   $genomic_align_block->reference_slice_end(1283);
00044   $genomic_align_block->score(56.2);
00045   $genomic_align_block->length(562);
00046 
00047 GET VALUES
00048   my $genomic_align_block_adaptor = $genomic_align_block->adaptor();
00049   my $dbID = $genomic_align_block->dbID();
00050   my $method_link_species_set = $genomic_align_block->method_link_species_set;
00051   my $genomic_aligns = $genomic_align_block->genomic_align_array();
00052   my $reference_genomic_align = $genomic_align_block->reference_genomic_align();
00053   my $non_reference_genomic_aligns = $genomic_align_block->get_all_non_reference_genomic_aligns();
00054   my $reference_slice = $genomic_align_block->reference_slice();
00055   my $reference_slice_start = $genomic_align_block->reference_slice_start();
00056   my $reference_slice_end = $genomic_align_block->reference_slice_end();
00057   my $score = $genomic_align_block->score();
00058   my $length = $genomic_align_block->length;
00059   my alignment_strings = $genomic_align_block->alignment_strings;
00060   my $genomic_align_block_is_on_the_original_strand =
00061       $genomic_align_block->get_original_strand;
00062 
00063 =head1 DESCRIPTION
00064 
00065 The GenomicAlignBlock object stores information about an alignment comprising of two or more pieces of genomic DNA.
00066 
00067 
00068 =head1 OBJECT ATTRIBUTES
00069 
00070 =over
00071 
00072 =item dbID
00073 
00074 corresponds to genomic_align_block.genomic_align_block_id
00075 
00076 =item adaptor
00077 
00078 Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlockAdaptor object to access DB
00079 
00080 =item method_link_species_set_id
00081 
00082 corresponds to method_link_species_set.method_link_species_set_id (external ref.)
00083 
00084 =item method_link_species_set
00085 
00086 Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object corresponding to method_link_species_set_id
00087 
00088 =item score
00089 
00090 corresponds to genomic_align_block.score
00091 
00092 =item perc_id
00093 
00094 corresponds to genomic_align_block.perc_id
00095 
00096 =item length
00097 
00098 corresponds to genomic_align_block.length
00099 
00100 =item group_id
00101 
00102 corresponds to the genomic_align_block.group_id
00103 
00104 =item reference_genomic_align_id
00105 
00106 When looking for genomic alignments in a given slice or dnafrag, the
00107 reference_genomic_align corresponds to the Bio::EnsEMBL::Compara::GenomicAlign
00108 included in the starting slice or dnafrag. The reference_genomic_align_id is
00109 the dbID corresponding to the reference_genomic_align. All remaining
00110 Bio::EnsEMBL::Compara::GenomicAlign objects included in the
00111 Bio::EnsEMBL::Compara::GenomicAlignBlock are the non_reference_genomic_aligns.
00112 
00113 =item reference_genomic_align
00114 
00115 Bio::EnsEMBL::Compara::GenomicAling object corresponding to reference_genomic_align_id
00116 
00117 =item genomic_align_array
00118 
00119 listref of Bio::EnsEMBL::Compara::GenomicAlign objects corresponding to this
00120 Bio::EnsEMBL::Compara::GenomicAlignBlock object
00121 
00122 =item reference_slice
00123 
00124 This is the Bio::EnsEMBL::Slice object used as argument to the
00125 Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlockAdaptor->fetch_all_by_Slice method.
00126 
00127 =item reference_slice_start
00128 
00129 starting position in the coordinates system defined by the reference_slice
00130 
00131 =item reference_slice_end
00132 
00133 ending position in the coordinates system defined by the reference_slice
00134 
00135 =back
00136 
00137 =head1 APPENDIX
00138 
00139 The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
00140 
00141 =cut
00142 
00143 
00144 # Let the code begin...
00145 
00146 
00147 package Bio::EnsEMBL::Compara::GenomicAlignBlock;
00148 use strict;
00149 
00150 # Object preamble
00151 use Bio::EnsEMBL::Utils::Argument qw(rearrange);
00152 use Bio::EnsEMBL::Utils::Exception qw(throw warning info deprecate verbose);
00153 use Bio::EnsEMBL::Compara::GenomicAlign;
00154 use Bio::SimpleAlign;
00155 use Bio::EnsEMBL::Compara::BaseGenomicAlignSet;
00156 
00157 our @ISA = qw(Bio::EnsEMBL::Compara::BaseGenomicAlignSet);
00158 
00159 =head2 new (CONSTRUCTOR)
00160 
00161   Arg [-DBID] : (opt.) int $dbID (the database internal ID for this object)
00162   Arg [-ADAPTOR]
00163               : (opt.) Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlockAdaptor $adaptor
00164                 (the adaptor for connecting to the database)
00165   Arg [-METHOD_LINK_SPECIES_SET]
00166               : (opt.) Bio::EnsEMBL::Compara::MethodLinkSpeciesSet $mlss
00167                 (this defines the type of alignment and the set of species used
00168                 to get this GenomicAlignBlock)
00169   Arg [-METHOD_LINK_SPECIES_SET_ID]
00170               : (opt.) int $mlss_id (the database internal ID for the $mlss)
00171   Arg [-SCORE]: (opt.) float $score (the score of this alignment)
00172   Arg [-PERC_ID]
00173               : (opt.) int $perc_id (the percentage of identity, only used for pairwise)
00174   Arg [-LENGTH]
00175               : (opt.) int $length (the length of this alignment, taking into account
00176                 gaps and all)
00177   Arg [-GROUP_ID]
00178               : (opt.) int $group)id (the group ID for this alignment)
00179   Arg [-REFERENCE_GENOMIC_ALIGN]
00180               : (opt.) Bio::EnsEMBL::Compara::GenomicAlign $reference_genomic_align (the
00181                 Bio::EnsEMBL::Compara::GenomicAlign corresponding to the requesting
00182                 Bio::EnsEMBL::Compara::DnaFrag or Bio::EnsEMBL::Slice when this
00183                 Bio::EnsEMBL::Compara::GenomicAlignBlock has been fetched from a
00184                 Bio::EnsEMBL::Compara::DnaFrag or a Bio::EnsEMBL::Slice)
00185   Arg [-REFERENCE_GENOMIC_ALIGN_ID]
00186               : (opt.) int $reference_genomic_align (the database internal ID of the
00187                 $reference_genomic_align)
00188   Arg [-GENOMIC_ALIGN_ARRAY]
00189               : (opt.) array_ref $genomic_aligns (a reference to the array of
00190                 Bio::EnsEMBL::Compara::GenomicAlign objects corresponding to this
00191                 Bio::EnsEMBL::Compara::GenomicAlignBlock object)
00192   Example    : my $genomic_align_block =
00193                    new Bio::EnsEMBL::Compara::GenomicAlignBlock(
00194                        -adaptor => $gaba,
00195                        -method_link_species_set => $method_link_species_set,
00196                        -score => 56.2,
00197                        -length => 1203,
00198                        -group_id => 1234,
00199                        -genomic_align_array => [$genomic_align1, $genomic_align2...]
00200                    );
00201   Description: Creates a new GenomicAlignBlock object
00202   Returntype : Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlock
00203   Exceptions : none
00204   Caller     : general
00205   Status     : Stable
00206 
00207 =cut
00208 
00209 sub new {
00210   my($class, @args) = @_;
00211   
00212   my $self = {};
00213   bless $self,$class;
00214     
00215   my ($adaptor, $dbID, $method_link_species_set, $method_link_species_set_id,
00216           $score, $perc_id, $length, $group_id, $level_id, $reference_genomic_align, $reference_genomic_align_id,
00217           $genomic_align_array, $starting_genomic_align_id, $ungapped_genomic_align_blocks) = 
00218     rearrange([qw(
00219         ADAPTOR DBID METHOD_LINK_SPECIES_SET METHOD_LINK_SPECIES_SET_ID
00220         SCORE PERC_ID LENGTH GROUP_ID LEVEL_ID REFERENCE_GENOMIC_ALIGN REFERENCE_GENOMIC_ALIGN_ID
00221         GENOMIC_ALIGN_ARRAY STARTING_GENOMIC_ALIGN_ID UNGAPPED_GENOMIC_ALIGN_BLOCKS)],
00222             @args);
00223 
00224   $self->{_original_strand} = 1;
00225 
00226   if (defined($ungapped_genomic_align_blocks)) {
00227     return $self->_create_from_a_list_of_ungapped_genomic_align_blocks($ungapped_genomic_align_blocks);
00228   }
00229   $self->adaptor($adaptor) if (defined ($adaptor));
00230   $self->dbID($dbID) if (defined ($dbID));
00231   $self->method_link_species_set($method_link_species_set) if (defined ($method_link_species_set));
00232   $self->method_link_species_set_id($method_link_species_set_id)
00233       if (defined ($method_link_species_set_id));
00234   $self->score($score) if (defined ($score));
00235   $self->perc_id($perc_id) if (defined ($perc_id));
00236   $self->length($length) if (defined ($length));
00237   $self->group_id($group_id) if (defined ($group_id));
00238   $self->level_id($level_id) if (defined ($level_id));
00239   $self->reference_genomic_align($reference_genomic_align)
00240       if (defined($reference_genomic_align));
00241   $self->reference_genomic_align_id($reference_genomic_align_id)
00242       if (defined($reference_genomic_align_id));
00243   $self->genomic_align_array($genomic_align_array) if (defined($genomic_align_array));
00244 
00245   $self->starting_genomic_align_id($starting_genomic_align_id) if (defined($starting_genomic_align_id));
00246 
00247   return $self;
00248 }
00249 
00250 =head2 new_fast
00251 
00252   Arg [1]    : hash reference $hashref
00253   Example    : none
00254   Description: This is an ultra fast constructor which requires knowledge of
00255                the objects internals to be used.
00256   Returntype :
00257   Exceptions : none
00258   Caller     :
00259   Status     : Stable
00260 
00261 =cut
00262 
00263 sub new_fast {
00264   my $class = shift;
00265   my $hashref = shift;
00266 
00267   return bless $hashref, $class;
00268 }
00269 
00270 
00271 =head2 adaptor
00272 
00273   Arg [1]    : Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlockAdaptor $adaptor
00274   Example    : my $gen_ali_blk_adaptor = $genomic_align_block->adaptor();
00275   Example    : $genomic_align_block->adaptor($gen_ali_blk_adaptor);
00276   Description: Getter/Setter for the adaptor this object uses for database
00277                interaction.
00278   Returntype : Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlockAdaptor object
00279   Exceptions : thrown if $adaptor is not a
00280                Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlockAdaptor object
00281   Caller     : general
00282   Status     : Stable
00283 
00284 =cut
00285 
00286 sub adaptor {
00287   my ($self, $adaptor) = @_;
00288 
00289   if (defined($adaptor)) {
00290     throw("$adaptor is not a Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlockAdaptor object")
00291         unless ($adaptor->isa("Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlockAdaptor"));
00292     $self->{'adaptor'} = $adaptor;
00293   }
00294 
00295   return $self->{'adaptor'};
00296 }
00297 
00298 
00299 =head2 dbID
00300 
00301   Arg [1]    : integer $dbID
00302   Example    : my $dbID = $genomic_align_block->dbID();
00303   Example    : $genomic_align_block->dbID(12);
00304   Description: Getter/Setter for the attribute dbID
00305   Returntype : integer
00306   Exceptions : none
00307   Caller     : general
00308   Status     : Stable
00309 
00310 =cut
00311 
00312 sub dbID {
00313   my ($self, $dbID) = @_;
00314 
00315   if (defined($dbID)) {
00316     $self->{'dbID'} = $dbID;
00317   }
00318 
00319   return $self->{'dbID'};
00320 }
00321 
00322 
00323 =head2 method_link_species_set
00324 
00325   Arg [1]    : Bio::EnsEMBL::Compara::MethodLinkSpeciesSet $method_link_species_set
00326   Example    : $method_link_species_set = $genomic_align_block->method_link_species_set;
00327   Example    : $genomic_align_block->method_link_species_set($method_link_species_set);
00328   Description: get/set for attribute method_link_species_set. If no
00329                argument is given, the method_link_species_set is not defined but
00330                both the method_link_species_set_id and the adaptor are, it tries
00331                to fetch the data using the method_link_species_set_id
00332   Returntype : Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object
00333   Exceptions : thrown if $method_link_species_set is not a
00334                Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object
00335   Caller     : general
00336   Status     : Stable
00337 
00338 =cut
00339 
00340 sub method_link_species_set {
00341   my ($self, $method_link_species_set) = @_;
00342 
00343   if (defined($method_link_species_set)) {
00344     throw("$method_link_species_set is not a Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object")
00345         unless ($method_link_species_set->isa("Bio::EnsEMBL::Compara::MethodLinkSpeciesSet"));
00346     $self->{'method_link_species_set'} = $method_link_species_set;
00347     if ($self->{'method_link_species_set_id'}) {
00348       $self->{'method_link_species_set_id'} = $self->{'method_link_species_set'}->dbID;
00349     }
00350     ## Update the MethodLinkSpeciesSet for the GenomicAligns included in this GenomicAlignBlock
00351     if (defined($self->{genomic_align_array})) {
00352       foreach my $this_genomic_align (@{$self->{genomic_align_array}}) {
00353         $this_genomic_align->method_link_species_set($method_link_species_set);
00354       }
00355     }
00356 
00357   } elsif (!defined($self->{'method_link_species_set'}) and defined($self->{'adaptor'})
00358           and defined($self->method_link_species_set_id)) {
00359     # Try to get object from ID. Use method_link_species_set_id function and not the attribute in the <if>
00360     # clause because the attribute can be retrieved from other sources if it has not been already defined.
00361     my $mlssa = $self->adaptor->db->get_MethodLinkSpeciesSetAdaptor;
00362     $self->{'method_link_species_set'} = $mlssa->fetch_by_dbID($self->{'method_link_species_set_id'})
00363   }
00364 
00365   return $self->{'method_link_species_set'};
00366 }
00367 
00368 
00369 =head2 method_link_species_set_id
00370 
00371   Arg [1]    : integer $method_link_species_set_id
00372   Example    : $method_link_species_set_id = $genomic_align_block->method_link_species_set_id;
00373   Example    : $genomic_align_block->method_link_species_set_id(3);
00374   Description: Getter/Setter for the attribute method_link_species_set_id. If no
00375                argument is given, the method_link_species_set_id is not defined but
00376                the method_link_species_set is, it tries to get the data from the
00377                method_link_species_set object. If this fails, it tries to get and set
00378                all the direct attributes from the database using the dbID of the
00379                Bio::Ensembl::Compara::GenomicAlignBlock object.
00380   Returntype : integer
00381   Exceptions : thrown if $method_link_species_set_id does not match a previously defined
00382                method_link_species_set
00383   Caller     : object::methodname
00384   Status     : Stable
00385 
00386 =cut
00387 
00388 sub method_link_species_set_id {
00389   my ($self, $method_link_species_set_id) = @_;
00390 
00391   if (defined($method_link_species_set_id)) {
00392     $self->{'method_link_species_set_id'} = $method_link_species_set_id;
00393     if (defined($self->{'method_link_species_set'}) and $self->{'method_link_species_set_id'}) {
00394       $self->{'method_link_species_set'} = undef;
00395     }
00396     ## Update the MethodLinkSpeciesSet for the GenomicAligns included in this GenomicAlignBlock
00397     if (defined($self->{genomic_align_array})) {
00398       foreach my $this_genomic_align (@{$self->{genomic_align_array}}) {
00399         $this_genomic_align->method_link_species_set_id($method_link_species_set_id);
00400       }
00401     }
00402 
00403   } elsif (!($self->{'method_link_species_set_id'})) {
00404     # Try to get the ID from other sources...
00405     if (defined($self->{'method_link_species_set'})) {
00406       # ...from the object
00407       $self->{'method_link_species_set_id'} = $self->{'method_link_species_set'}->dbID;
00408     } elsif (defined($self->{'adaptor'}) and defined($self->dbID)) {
00409       # ...from the database, using the dbID of the Bio::Ensembl::Compara::GenomicAlignBlock object
00410       $self->adaptor->retrieve_all_direct_attributes($self);
00411     }
00412   }
00413 
00414   return $self->{'method_link_species_set_id'};
00415 }
00416 
00417 =head2 starting_genomic_align_id (DEPRECATED)
00418 
00419   DEPRECATED! Use Bio::EnsEMBL::Compara::GenomicAlignBlock->reference_genomic_align_id() method instead
00420  
00421   Arg [1]    : integer $reference_genomic_align_id
00422   Example    : $genomic_align_block->starting_genomic_align_id(4321);
00423   Description: set for attribute reference_genomic_align_id. A value of 0 will set the
00424                reference_genomic_align_id attribute to undef. When looking for genomic
00425                alignments in a given slice or dnafrag, the reference_genomic_align
00426                corresponds to the Bio::EnsEMBL::Compara::GenomicAlign included in the
00427                starting slice or dnafrag. The reference_genomic_align_id is the dbID
00428                corresponding to the reference_genomic_align. All remaining
00429                Bio::EnsEMBL::Compara::GenomicAlign objects included in the
00430                Bio::EnsEMBL::Compara::GenomicAlignBlock are the non_reference_genomic_aligns.
00431   Returntype : none
00432   Exceptions : throw if $reference_genomic_align_id id not a postive number
00433   Caller     : $genomic_align_block->starting_genomic_align_id(int)
00434 
00435 =cut
00436 
00437 sub starting_genomic_align_id {
00438   my $self = shift;
00439   deprecate("Use Bio::EnsEMBL::Compara::GenomicAlignBlock->reference_genomic_align_id() method instead");
00440   return $self->reference_genomic_align_id(@_);
00441 }
00442 
00443 =head2 starting_genomic_align (DEPRECATED)
00444 
00445   DEPRECATED! Use Bio::EnsEMBL::Compara::GenomicAlignBlock->reference_genomic_align() method instead
00446  
00447   Arg [1]    : (none)
00448   Example    : $genomic_align = $genomic_align_block->starting_genomic_align();
00449   Description: get the reference_genomic_align. When looking for genomic alignments in
00450                a given slice or dnafrag, the reference_genomic_align corresponds to the
00451                Bio::EnsEMBL::Compara::GenomicAlign included in the starting slice or
00452                dnafrag. The reference_genomic_align_id is the dbID corresponding to the
00453                reference_genomic_align. All remaining Bio::EnsEMBL::Compara::GenomicAlign
00454                objects included in the Bio::EnsEMBL::Compara::GenomicAlignBlock are the
00455                non_reference_genomic_aligns.
00456   Returntype : Bio::EnsEMBL::Compara::GenomicAlign object
00457   Exceptions : warns if no reference_genomic_align_id has been set and returns a ref.
00458                to an empty array
00459   Exceptions : warns if no genomic_align_array has been set and returns a ref.
00460                to an empty array
00461   Exceptions : throw if reference_genomic_align_id does not match any of the
00462                Bio::EnsEMBL::Compara::GenomicAlign objects in the genomic_align_array
00463   Caller     : $genomic_align_block->starting_genomic_align()
00464 
00465 =cut
00466 
00467 sub starting_genomic_align {
00468   my $self = shift;
00469   deprecate("Use Bio::EnsEMBL::Compara::GenomicAlignBlock->reference_genomic_align() method instead");
00470   return $self->reference_genomic_align(@_);
00471 }
00472 
00473 
00474 =head2 resulting_genomic_aligns
00475 
00476   DEPRECATED! Use Bio::EnsEMBL::Compara::GenomicAlignBlock->get_all_non_reference_genomic_align()
00477   method instead
00478  
00479   Arg [1]    : (none)
00480   Example    : $genomic_aligns = $genomic_align_block->resulting_genomic_aligns();
00481   Description: get the all the non_reference_genomic_aligns. When looking for genomic
00482                alignments in a given slice or dnafrag, the reference_genomic_align
00483                corresponds to the Bio::EnsEMBL::Compara::GenomicAlign included in the
00484                reference slice or dnafrag. The reference_genomic_align_id is the dbID
00485                corresponding to the reference_genomic_align. All remaining
00486                Bio::EnsEMBL::Compara::GenomicAlign objects included in the
00487                Bio::EnsEMBL::Compara::GenomicAlignBlock are the
00488                non_reference_genomic_aligns.
00489   Returntype : a ref. to an array of Bio::EnsEMBL::Compara::GenomicAlign objects
00490   Exceptions : warns if no reference_genomic_align_id has been set and returns a ref.
00491                to an empty array
00492   Exceptions : warns if no genomic_align_array has been set and returns a ref.
00493                to an empty array
00494   Caller     : $genomic_align_block->resulting_genomic_aligns()
00495 
00496 =cut
00497 
00498 sub resulting_genomic_aligns {
00499   my ($self) = @_;
00500   deprecate("Use Bio::EnsEMBL::Compara::GenomicAlignBlock->get_all_non_reference_genomic_aligns() method instead");
00501   return $self->get_all_non_reference_genomic_aligns(@_);
00502 }
00503 
00504 
00505 =head2 reference_genomic_align
00506 
00507   Arg [1]    : (optional) Bio::EnsEMBL::Compara::GenomicAlign $reference_genomic_align
00508   Example    : $genomic_align_block->reference_genomic_align($reference_genomic_align);
00509   Example    : $genomic_align = $genomic_align_block->reference_genomic_align();
00510   Description: get/set the reference_genomic_align. When looking for genomic alignments in
00511                a given slice or dnafrag, the reference_genomic_align corresponds to the
00512                Bio::EnsEMBL::Compara::GenomicAlign included in the starting slice or
00513                dnafrag. The reference_genomic_align_id is the dbID corresponding to the
00514                reference_genomic_align. All remaining Bio::EnsEMBL::Compara::GenomicAlign
00515                objects included in the Bio::EnsEMBL::Compara::GenomicAlignBlock are the
00516                non_reference_genomic_aligns.
00517                Synchronises reference_genomic_align and reference_genomic_align_id
00518                attributes.
00519   Returntype : Bio::EnsEMBL::Compara::GenomicAlign object
00520   Exceptions : throw if reference_genomic_align is not a Bio::EnsEMBL::Compara::GenomicAlign
00521                object
00522   Exceptions : throw if reference_genomic_align_id does not match any of the
00523                Bio::EnsEMBL::Compara::GenomicAlign objects in the genomic_align_array
00524   Caller     : $genomic_align_block->reference_genomic_align()
00525   Status     : Stable
00526 
00527 =cut
00528 
00529 sub reference_genomic_align {
00530   my ($self, $reference_genomic_align) = @_;
00531 
00532   if (defined($reference_genomic_align)) {
00533     throw("[$reference_genomic_align] must be a Bio::EnsEMBL::Compara::GenomicAlign object")
00534         unless($reference_genomic_align and  ref($reference_genomic_align) and
00535             $reference_genomic_align->isa("Bio::EnsEMBL::Compara::GenomicAlign"));
00536     $self->{'reference_genomic_align'} = $reference_genomic_align;
00537 
00538     ## Synchronises reference_genomic_align and reference_genomic_align_id attributes
00539     if (defined($self->{'reference_genomic_align'}->dbID)) {
00540       $self->{'reference_genomic_align_id'} = $self->{'reference_genomic_align'}->dbID;
00541     }
00542 
00543   ## Try to get data from other sources...
00544   } elsif (!defined($self->{'reference_genomic_align'})) {
00545     
00546     ## ...from the reference_genomic_align_id attribute
00547     if (defined($self->{'reference_genomic_align_id'}) and @{$self->get_all_GenomicAligns}) {
00548       my $reference_genomic_align_id = $self->{'reference_genomic_align_id'};
00549       foreach my $this_genomic_align (@{$self->get_all_GenomicAligns}) {
00550         if ($this_genomic_align->dbID == $reference_genomic_align_id) {
00551           $self->{'reference_genomic_align'} = $this_genomic_align;
00552           return $this_genomic_align;
00553         }
00554       }
00555       throw("[$self] Cannot found Bio::EnsEMBL::Compara::GenomicAlign::reference_genomic_align_id".
00556           " ($reference_genomic_align_id) in the genomic_align_array");
00557     }
00558   }
00559 
00560   return $self->{'reference_genomic_align'};
00561 }
00562 
00563 =head2 reference_genomic_align_id
00564 
00565   Arg [1]    : integer $reference_genomic_align_id
00566   Example    : $genomic_align_block->reference_genomic_align_id(4321);
00567   Description: get/set for attribute reference_genomic_align_id. A value of 0 will set the
00568                reference_genomic_align_id attribute to undef. When looking for genomic
00569                alignments in a given slice or dnafrag, the reference_genomic_align
00570                corresponds to the Bio::EnsEMBL::Compara::GenomicAlign included in the
00571                starting slice or dnafrag. The reference_genomic_align_id is the dbID
00572                corresponding to the reference_genomic_align. All remaining
00573                Bio::EnsEMBL::Compara::GenomicAlign objects included in the
00574                Bio::EnsEMBL::Compara::GenomicAlignBlock are the
00575                non_reference_genomic_aligns.
00576                Synchronises reference_genomic_align and reference_genomic_align_id
00577                attributes.
00578   Returntype : integer
00579   Exceptions : throw if $reference_genomic_align_id id not a postive number
00580   Caller     : $genomic_align_block->reference_genomic_align_id(int)
00581   Status     : Stable
00582 
00583 =cut
00584 
00585 sub reference_genomic_align_id {
00586   my ($self, $reference_genomic_align_id) = @_;
00587  
00588   if (defined($reference_genomic_align_id)) {
00589     if ($reference_genomic_align_id !~ /^\d+$/) {
00590       throw "[$reference_genomic_align_id] should be a positive number.";
00591     }
00592     $self->{'reference_genomic_align_id'} = ($reference_genomic_align_id or undef);
00593 
00594     ## Synchronises reference_genomic_align and reference_genomic_align_id
00595     if (defined($self->{'reference_genomic_align'}) and
00596         defined($self->{'reference_genomic_align'}->dbID) and
00597         ($self->{'reference_genomic_align'}->dbID != ($self->{'reference_genomic_align_id'} or 0))) {
00598         $self->{'reference_genomic_align'} = undef; ## Attribute will be set on request
00599     }
00600 
00601   ## Try to get data from other sources...
00602   } elsif (!defined($self->{'reference_genomic_align_id'})) {
00603 
00604     ## ...from the reference_genomic_align attribute
00605     if (defined($self->{'reference_genomic_align'}) and
00606         defined($self->{'reference_genomic_align'}->dbID)) {
00607       $self->{'reference_genomic_align_id'} = $self->{'reference_genomic_align'}->dbID;
00608     }
00609   }
00610   
00611   return $self->{'reference_genomic_align_id'};
00612 }
00613 
00614 =head2 get_all_non_reference_genomic_aligns
00615 
00616   Arg [1]    : (none)
00617   Example    : $genomic_aligns = $genomic_align_block->get_all_non_reference_genomic_aligns();
00618   Description: get the non_reference_genomic_aligns. When looking for genomic alignments in
00619                a given slice or dnafrag, the reference_genomic_align corresponds to the
00620                Bio::EnsEMBL::Compara::GenomicAlign included in the starting slice or
00621                dnafrag. The reference_genomic_align_id is the dbID corresponding to the
00622                reference_genomic_align. All remaining Bio::EnsEMBL::Compara::GenomicAlign
00623                objects included in the Bio::EnsEMBL::Compara::GenomicAlignBlock are the
00624                non_reference_genomic_aligns.
00625   Returntype : a ref. to an array of Bio::EnsEMBL::Compara::GenomicAlign objects
00626   Exceptions : warns if no reference_genomic_align_id has been set and returns a ref.
00627                to an empty array
00628   Exceptions : warns if no genomic_align_array has been set and returns a ref.
00629                to an empty array
00630   Caller     : $genomic_align_block->non_reference_genomic_aligns()
00631   Status     : Stable
00632 
00633 =cut
00634 
00635 sub get_all_non_reference_genomic_aligns {
00636   my ($self) = @_;
00637   my $all_non_reference_genomic_aligns = [];
00638  
00639   my $reference_genomic_align_id = $self->reference_genomic_align_id;
00640   my $reference_genomic_align = $self->reference_genomic_align;
00641   if (!defined($reference_genomic_align_id) and !defined($reference_genomic_align)) {
00642     warning("Trying to get Bio::EnsEMBL::Compara::GenomicAlign::all_non_reference_genomic_aligns".
00643         " when no reference_genomic_align has been set before");
00644     return $all_non_reference_genomic_aligns;
00645   }
00646   my $genomic_aligns = $self->get_all_GenomicAligns; ## Lazy loading compliant
00647   if (!@$genomic_aligns) {
00648     warning("Trying to get Bio::EnsEMBL::Compara::GenomicAlign::all_non_reference_genomic_aligns".
00649         " when no genomic_align_array can be retrieved");
00650     return $all_non_reference_genomic_aligns;
00651   }
00652 
00653   foreach my $this_genomic_align (@$genomic_aligns) {
00654     if (($this_genomic_align->dbID or -1) != ($reference_genomic_align_id or -2) and
00655         $this_genomic_align != $reference_genomic_align) {
00656       push(@$all_non_reference_genomic_aligns, $this_genomic_align);
00657     }
00658   }
00659 
00660   return $all_non_reference_genomic_aligns;
00661 }
00662 
00663 
00664 =head2 genomic_align_array
00665 
00666   Arg [1]    : array reference containing Bio::EnsEMBL::Compara::GenomicAlign objects
00667   Example    : $genomic_aligns = $genomic_align_block->genomic_align_array();
00668                $genomic_align_block->genomic_align_array([$genomic_align1, $genomic_align2]);
00669   Description: get/set for attribute genomic_align_array. If no argument is given, the
00670                genomic_align_array is not defined but both the dbID and the adaptor are,
00671                it tries to fetch the data from the database using the dbID of the
00672                Bio::EnsEMBL::Compara::GenomicAlignBlock object.
00673                You can unset all cached GenomicAlign using 0 as argument. They will be
00674                loaded again from the database if needed.
00675   Returntype : array reference containing Bio::EnsEMBL::Compara::GenomicAlign objects
00676   Exceptions : none
00677   Caller     : general
00678   Status     : Stable
00679 
00680 =cut
00681 
00682 sub genomic_align_array {
00683   my ($self, $genomic_align_array) = @_;
00684  
00685   if (defined($genomic_align_array)) {
00686     if (!$genomic_align_array) {
00687       ## Clean cache.
00688       $self->{'genomic_align_array'} = undef;
00689       $self->{'reference_genomic_align'} = undef;
00690       return undef;
00691     }
00692     foreach my $genomic_align (@$genomic_align_array) {
00693       throw("[$genomic_align] is not a Bio::EnsEMBL::Compara::GenomicAlign object")
00694           unless (ref $genomic_align and $genomic_align->isa("Bio::EnsEMBL::Compara::GenomicAlign"));
00695       # Create weak circular reference to genomic_align_block from each genomic_align
00696       $genomic_align->genomic_align_block($self); 
00697     }
00698     $self->{'genomic_align_array'} = $genomic_align_array;
00699 
00700   } elsif (!defined($self->{'genomic_align_array'}) and defined($self->{'adaptor'})
00701         and defined($self->{'dbID'})) {
00702     # Fetch data from DB (allow lazy fetching of genomic_align_block objects)
00703     my $genomic_align_adaptor = $self->adaptor->db->get_GenomicAlignAdaptor();
00704     $self->{'genomic_align_array'} = 
00705         $genomic_align_adaptor->fetch_all_by_genomic_align_block_id($self->{'dbID'});
00706     foreach my $this_genomic_align (@{$self->{'genomic_align_array'}}) {
00707       $this_genomic_align->genomic_align_block($self);
00708     }
00709   }
00710   
00711   return $self->{'genomic_align_array'};
00712 }
00713 
00714 
00715 =head2 add_GenomicAlign
00716 
00717   Arg [1]    : Bio::EnsEMBL::Compara::GenomicAlign $genomic_align
00718   Example    : $genomic_align_block->add_GenomicAlign($genomic_align);
00719   Description: adds another Bio::EnsEMBL::Compara::GenomicAlign object to the set of
00720                Bio::EnsEMBL::Compara::GenomicAlign objects in the attribute
00721                genomic_align_array.
00722   Returntype : Bio::EnsEMBL::Compara::GenomicAlign object
00723   Exceptions : thrown if wrong argument
00724   Caller     : general
00725   Status     : Stable
00726 
00727 =cut
00728 
00729 sub add_GenomicAlign {
00730   my ($self, $genomic_align) = @_;
00731 
00732   throw("[$genomic_align] is not a Bio::EnsEMBL::Compara::GenomicAlign object")
00733       unless ($genomic_align and ref($genomic_align) and 
00734           $genomic_align->isa("Bio::EnsEMBL::Compara::GenomicAlign"));
00735   # Create weak circular reference to genomic_align_block from each genomic_align
00736   $genomic_align->genomic_align_block($self);
00737   push(@{$self->{'genomic_align_array'}}, $genomic_align);
00738 
00739   return $genomic_align;
00740 }
00741 
00742 
00743 =head2 get_all_GenomicAligns
00744 
00745   Arg [1]    : none
00746   Example    : $genomic_aligns = $genomic_align_block->get_all_GenomicAligns();
00747   Description: returns the set of Bio::EnsEMBL::Compara::GenomicAlign objects in
00748                the attribute genomic_align_array.
00749   Returntype : array reference containing Bio::EnsEMBL::Compara::GenomicAlign objects
00750   Exceptions : none
00751   Caller     : general
00752   Status     : Stable
00753 
00754 =cut
00755 
00756 sub get_all_GenomicAligns {
00757   my ($self) = @_;
00758  
00759   return ($self->genomic_align_array or []);
00760 }
00761 
00762 
00763 =head2 score
00764 
00765   Arg [1]    : double $score
00766   Example    : my $score = $genomic_align_block->score();
00767                $genomic_align_block->score(56.2);
00768   Description: get/set for attribute score. If no argument is given, the score
00769                is not defined but both the dbID and the adaptor are, it tries to
00770                fetch and set all the direct attributes from the database using the
00771                dbID of the Bio::EnsEMBL::Compara::GenomicAlignBlock object.
00772   Returntype : double
00773   Exceptions : none
00774   Caller     : general
00775   Status     : Stable
00776 
00777 =cut
00778 
00779 sub score {
00780   my ($self, $score) = @_;
00781 
00782   if (defined($score)) {
00783     $self->{'score'} = $score;
00784   } elsif (!defined($self->{'score'})) {
00785     # Try to get the ID from other sources...
00786     if (defined($self->{'adaptor'}) and defined($self->dbID)) {
00787       # ...from the database, using the dbID of the Bio::Ensembl::Compara::GenomicAlignBlock object
00788       $self->adaptor->retrieve_all_direct_attributes($self);
00789     }
00790   }
00791 
00792   return $self->{'score'};
00793 }
00794 
00795 
00796 =head2 perc_id
00797 
00798   Arg [1]    : double $perc_id
00799   Example    : my $perc_id = $genomic_align_block->perc_id;
00800   Example    : $genomic_align_block->perc_id(95.4);
00801   Description: get/set for attribute perc_id. If no argument is given, the perc_id
00802                is not defined but both the dbID and the adaptor are, it tries to
00803                fetch and set all the direct attributes from the database using the
00804                dbID of the Bio::EnsEMBL::Compara::GenomicAlignBlock object.
00805   Returntype : double
00806   Exceptions : none
00807   Caller     : general
00808   Status     : Stable
00809 
00810 =cut
00811 
00812 sub perc_id {
00813   my ($self, $perc_id) = @_;
00814  
00815   if (defined($perc_id)) {
00816     $self->{'perc_id'} = $perc_id;
00817   } elsif (!defined($self->{'perc_id'})) {
00818     # Try to get the ID from other sources...
00819     if (defined($self->{'adaptor'}) and defined($self->dbID)) {
00820       # ...from the database, using the dbID of the Bio::Ensembl::Compara::GenomicAlignBlock object
00821       $self->adaptor->retrieve_all_direct_attributes($self);
00822     }
00823   }
00824   
00825   return $self->{'perc_id'};
00826 }
00827 
00828 
00829 =head2 length
00830 
00831   Arg [1]    : integer $length
00832   Example    : my $length = $genomic_align_block->length;
00833   Example    : $genomic_align_block->length(562);
00834   Description: get/set for attribute length. If no argument is given, the length
00835                is not defined but both the dbID and the adaptor are, it tries to
00836                fetch and set all the direct attributes from the database using the
00837                dbID of the Bio::EnsEMBL::Compara::GenomicAlignBlock object.
00838   Returntype : integer
00839   Exceptions : none
00840   Caller     : general
00841   Status     : Stable
00842 
00843 =cut
00844 
00845 sub length {
00846   my ($self, $length) = @_;
00847  
00848   if (defined($length)) {
00849     $self->{'length'} = $length;
00850   } elsif (!defined($self->{'length'})) {
00851     # Try to get the ID from other sources...
00852     if (defined($self->{'adaptor'}) and defined($self->dbID)) {
00853       # ...from the database, using the dbID of the Bio::Ensembl::Compara::GenomicAlignBlock object
00854       $self->adaptor->retrieve_all_direct_attributes($self);
00855     } elsif (@{$self->get_all_GenomicAligns} and $self->get_all_GenomicAligns->[0]->aligned_sequence("+FAKE_SEQ")) {
00856       $self->{'length'} = CORE::length($self->get_all_GenomicAligns->[0]->aligned_sequence("+FAKE_SEQ"));
00857     }
00858   }
00859   
00860   return $self->{'length'};
00861 }
00862 
00863 =head2 group_id
00864 
00865   Arg [1]    : integer $group_id
00866   Example    : my $group_id = $genomic_align_block->group_id;
00867   Example    : $genomic_align_block->group_id(1234);
00868   Description: get/set for attribute group_id. 
00869   Returntype : integer
00870   Exceptions : none
00871   Caller     : general
00872   Status     : At risk
00873 
00874 =cut
00875 
00876 sub group_id {
00877     my ($self, $group_id) = @_;
00878 
00879     if (defined($group_id)) {
00880     $self->{'group_id'} = ($group_id);
00881     } elsif (!defined($self->{'group_id'})) {
00882     # Try to get the ID from other sources...
00883     if (defined($self->{'adaptor'}) and defined($self->dbID)) {
00884         # ...from the database, using the dbID of the Bio::Ensembl::Compara::GenomicAlignBlock object
00885         $self->adaptor->retrieve_all_direct_attributes($self);
00886     }
00887     }
00888     return $self->{'group_id'};
00889 }
00890 
00891 =head2 level_id
00892 
00893   Arg [1]    : int $level_id
00894   Example    : $level_id = $genomic_align->level_id;
00895   Example    : $genomic_align->level_id(1);
00896   Description: get/set for attribute level_id. If no argument is given, the level_id
00897                is not defined but both the dbID and the adaptor are, it tries to
00898                fetch and set all the direct attributes from the database using the
00899                dbID of the Bio::EnsEMBL::Compara::GenomicAlign object.
00900   Returntype : int
00901   Exceptions : none
00902   Warning    : warns if getting data from other sources fails.
00903   Caller     : object->methodname
00904   Status     : Stable
00905 
00906 =cut
00907 
00908 sub level_id {
00909   my ($self, $level_id) = @_;
00910 
00911   if (defined($level_id)) {
00912     $self->{'level_id'} = $level_id;
00913 
00914   } elsif (!defined($self->{'level_id'})) {
00915     if (defined($self->{'dbID'}) and defined($self->{'adaptor'})) {
00916       # Try to get the values from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlignBlock object
00917       $self->adaptor->retrieve_all_direct_attributes($self);
00918     } else {
00919 #      warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlignBlock->level_id".
00920 #          " You either have to specify more information (see perldoc for".
00921 #          " Bio::EnsEMBL::Compara::GenomicAlignBlock) or to set it up directly");
00922     }
00923   }
00924 
00925   return $self->{'level_id'};
00926 }
00927 
00928 =head2 requesting_slice (DEPRECATED)
00929 
00930   DEPRECATED! Use Bio::EnsEMBL::Compara::GenomicAlignBlock->reference_slice() method instead
00931  
00932   Arg [1]    : Bio::EnsEMBL::Slice $reference_slice
00933   Example    : my $reference_slice = $genomic_align_block->requesting_slice;
00934   Example    : $genomic_align_block->requesting_slice($reference_slice);
00935   Description: get/set for attribute reference_slice.
00936   Returntype : Bio::EnsEMBL::Slice object
00937   Exceptions : throw if $reference_slice is not a Bio::EnsEMBL::Slice
00938   Caller     : general
00939 
00940 =cut
00941 
00942 sub requesting_slice {
00943   my ($self) = shift;
00944   deprecate("Use Bio::EnsEMBL::Compara::GenomicAlignBlock->reference_slice() method instead");
00945   return $self->reference_slice(@_);
00946 }
00947 
00948 =head2 alignment_strings
00949 
00950   Arg [1]    : none
00951   Example    : $genomic_align_block->alignment_strings
00952   Description: Returns the alignment string of all the sequences in the
00953                alignment
00954   Returntype : array reference containing several strings
00955   Exceptions : none
00956   Caller     : general
00957   Status     : Stable
00958 
00959 =cut
00960 
00961 sub alignment_strings {
00962   my ($self) = @_;
00963   my $alignment_strings = [];
00964 
00965   foreach my $genomic_align (@{$self->get_all_GenomicAligns}) {
00966     push(@$alignment_strings, $genomic_align->aligned_sequence);
00967   }
00968 
00969   return $alignment_strings;
00970 }
00971 
00972 
00973 =head2 get_SimpleAlign
00974 
00975   Arg [1]    : list of string $flags
00976                "translated" = by default, the sequence alignment will be on nucleotide. With "translated" flag
00977                               the aligned sequences are translated.
00978                "uc" = by default aligned sequences are given in lower cases. With "uc" flag, the aligned
00979                       sequences are given in upper cases.
00980                "display_id" = by default the name of each sequence in the alignment is $dnafrag->name. With 
00981                               "dispaly_id" flag the name of each sequence is defined by the 
00982                               Bio::EnsEMBL::Compara::GenomicAlign display_id method.
00983   Example    : $daf->get_SimpleAlign or
00984                $daf->get_SimpleAlign("translated") or
00985                $daf->get_SimpleAlign("translated","uc")
00986   Description: Allows to rebuild the alignment string of all the genomic_align objects within
00987                this genomic_align_block using the cigar_line information
00988                and access to the core database slice objects
00989   Returntype : a Bio::SimpleAlign object
00990   Exceptions :
00991   Caller     : general
00992   Status     : Stable
00993 
00994 =cut
00995 
00996 sub get_SimpleAlign {
00997   my ( $self, @flags ) = @_;
00998 
00999   # setting the flags
01000   my $uc = 0;
01001   my $translated = 0;
01002   my $display_id = 0;
01003 
01004   for my $flag ( @flags ) {
01005     $uc = 1 if ($flag =~ /^uc$/i);
01006     $translated = 1 if ($flag =~ /^translated$/i);
01007     $display_id = 1 if ($flag =~ /^display_id$/i);
01008   }
01009 
01010   my $sa = Bio::SimpleAlign->new();
01011 
01012   #Hack to try to work with both bioperl 0.7 and 1.2:
01013   #Check to see if the method is called 'addSeq' or 'add_seq'
01014   my $bio07 = 0;
01015   if(!$sa->can('add_seq')) {
01016     $bio07 = 1;
01017   }
01018 
01019   my $all_genomic_aligns;
01020   if ($self->reference_genomic_align) {
01021     $all_genomic_aligns = [$self->reference_genomic_align,@{$self->get_all_non_reference_genomic_aligns}];
01022   } else {
01023     $all_genomic_aligns = $self->get_all_GenomicAligns();
01024   }
01025 
01026   foreach my $genomic_align (@$all_genomic_aligns) {
01027     my $alignSeq = $genomic_align->aligned_sequence;
01028     next if($alignSeq=~/^[\.\-]+$/);
01029     
01030     my $loc_seq = Bio::LocatableSeq->new(-SEQ    => $uc ? uc $alignSeq : lc $alignSeq,
01031                                          -START  => $genomic_align->dnafrag_start,
01032                                          -END    => $genomic_align->dnafrag_end,
01033                                          -ID     => $display_id ? $genomic_align->display_id : ($genomic_align->dnafrag->genome_db->name . "/" . $genomic_align->dnafrag->name),
01034                                          -STRAND => $genomic_align->dnafrag_strand);
01035 
01036     $loc_seq->seq($uc ? uc $loc_seq->translate->seq
01037                       : lc $loc_seq->translate->seq) if ($translated);
01038 
01039     if($bio07) { $sa->addSeq($loc_seq); }
01040     else       { $sa->add_seq($loc_seq); }
01041 
01042   }
01043   return $sa;
01044 }
01045 
01046 
01047 =head2 _create_from_a_list_of_ungapped_genomic_align_blocks (testing)
01048 
01049   Args       : listref of ungapped Bio::EnsEMBL::Compara::GenomicAlignBlocks
01050   Example    : $new_genomic_align_block =
01051                   $self->_create_from_a_list_of_ungapped_genomic_align_blocks(
01052                       $ungapped_genomic_align_blocks
01053                   );
01054   Description: Takes a list of ungapped Bio::EnsEMBL::Compara::GenomicAlignBlock
01055                objects and creates a new Bio::EnsEMBL::Compara::GenomicAlignBlock
01056   Returntype : Bio::EnsEMBL::Compara::GenomicAlignBlock object
01057   Exceptions : lots...
01058   Caller     : new()
01059   Status     : At risk
01060 
01061 =cut
01062 
01063 sub _create_from_a_list_of_ungapped_genomic_align_blocks {
01064   my ($self, $ungapped_genomic_align_blocks) = @_;
01065 
01066   ## Set adaptor
01067   my $adaptor = undef;
01068   foreach my $genomic_align_block (@$ungapped_genomic_align_blocks) {
01069     if ($genomic_align_block->adaptor) {
01070       $self->adaptor($adaptor);
01071       last;
01072     }
01073   }
01074   
01075   ## Set method_link_species_set
01076   my $method_link_species_set = undef;
01077   foreach my $genomic_align_block (@$ungapped_genomic_align_blocks) {
01078     if ($genomic_align_block->method_link_species_set) {
01079       if ($method_link_species_set) {
01080         if ($method_link_species_set->dbID != $genomic_align_block->method_link_species_set->dbID) {
01081           warning("Creating a GenomicAlignBlock from a list of ungapped GenomicAlignBlock with".
01082               " different method_link_species_set is not supported");
01083           return undef;
01084         }
01085       } else {
01086         $method_link_species_set = $genomic_align_block->method_link_species_set;
01087       }
01088     }
01089   }
01090   $self->method_link_species_set($method_link_species_set);
01091   
01092   ## Set method_link_species_set_id
01093   my $method_link_species_set_id = undef;
01094   foreach my $genomic_align_block (@$ungapped_genomic_align_blocks) {
01095     if ($genomic_align_block->method_link_species_set_id) {
01096       if ($method_link_species_set_id) {
01097         if ($method_link_species_set_id != $genomic_align_block->method_link_species_set_id) {
01098           warning("Creating a GenomicAlignBlock from a list of ungapped GenomicAlignBlock with".
01099               " different method_link_species_set_id is not supported");
01100           return undef;
01101         }
01102       } else {
01103         $method_link_species_set_id = $genomic_align_block->method_link_species_set_id;
01104       }
01105     }
01106   }
01107   $self->method_link_species_set_id($method_link_species_set_id);
01108 
01109   my $genomic_aligns;
01110   ## Check blocks and create new genomic_aligns
01111   foreach my $genomic_align_block (@$ungapped_genomic_align_blocks) {
01112     foreach my $this_genomic_align (@{$genomic_align_block->get_all_GenomicAligns}) {
01113       my $dnafrag_id = $this_genomic_align->dnafrag_id;
01114       if (!defined($genomic_aligns->{$dnafrag_id})) {
01115         $genomic_aligns->{$dnafrag_id} = new Bio::EnsEMBL::Compara::GenomicAlign(
01116                 -adaptor => $this_genomic_align->adaptor,
01117                 -method_link_species_set => $method_link_species_set,
01118                 -method_link_species_set_id => $method_link_species_set_id,
01119                 -dnafrag => $this_genomic_align->dnafrag,
01120                 -dnafrag_start => $this_genomic_align->dnafrag_start,
01121                 -dnafrag_end => $this_genomic_align->dnafrag_end,
01122                 -dnafrag_strand => $this_genomic_align->dnafrag_strand,
01123             );
01124       } else {
01125         ## Check strand
01126         if ($this_genomic_align->dnafrag_strand < $genomic_aligns->{$dnafrag_id}->dnafrag_strand) {
01127           warning("The list of ungapped GenomicAlignBlock is inconsistent in strand");
01128           return undef;
01129         }
01130 
01131         ## Check order and lengthen genomic_align
01132         if ($this_genomic_align->dnafrag_strand == -1) {
01133           if ($this_genomic_align->dnafrag_end >= $genomic_aligns->{$dnafrag_id}->dnafrag_start) {
01134             warning("The list of ungapped GenomicAlignBlock must be previously sorted");
01135             return undef;
01136           }
01137           $genomic_aligns->{$dnafrag_id}->dnafrag_start($this_genomic_align->dnafrag_start);
01138         } else {
01139           if ($this_genomic_align->dnafrag_start <= $genomic_aligns->{$dnafrag_id}->dnafrag_end) {
01140             warning("The list of ungapped GenomicAlignBlock must be previously sorted");
01141             return undef;
01142           }
01143           $genomic_aligns->{$dnafrag_id}->dnafrag_end($this_genomic_align->dnafrag_end);
01144         }
01145       }
01146     }
01147   }
01148   
01149   ## Create cigar lines
01150   my $cigar_lines;
01151   for (my $i=0; $i<@$ungapped_genomic_align_blocks; $i++) {
01152     my $genomic_align_block = $ungapped_genomic_align_blocks->[$i];
01153     my $block_length = 0;
01154     ## Calculate block length
01155     foreach my $this_genomic_align (@{$genomic_align_block->get_all_GenomicAligns}) {
01156       if ($block_length) {
01157         if ($block_length != CORE::length($this_genomic_align->aligned_sequence)) {
01158           warning("The list of ungapped GenomicAlignBlock is inconsistent in gaps");
01159           return undef;
01160         }
01161       } else {
01162         $block_length = CORE::length($this_genomic_align->aligned_sequence);
01163       }
01164     }
01165 
01166     next if ($block_length == 0); # Skip 0-length blocks (shouldn't happen)
01167     $block_length = "" if ($block_length == 1); # avoid a "1" in cigar_line
01168 
01169     ## Fix cigar line according to block length
01170     while (my ($id, $genomic_align) = each %{$genomic_aligns}) {
01171       my $is_included_in_this_block = 0;
01172       foreach my $this_genomic_align (@{$genomic_align_block->get_all_GenomicAligns}) {
01173         if ($this_genomic_align->dnafrag_id == $id) {
01174           $is_included_in_this_block = 1;
01175           $cigar_lines->{$id} .= $this_genomic_align->cigar_line;
01176           last;
01177         }
01178       }
01179       if (!$is_included_in_this_block) {
01180         $cigar_lines->{$id} .= $block_length."D";
01181       }
01182     }
01183 
01184     ## Add extra gaps between genomic_align_blocks
01185     if (defined($ungapped_genomic_align_blocks->[$i+1])) {
01186       foreach my $genomic_align1 (@{$genomic_align_block->get_all_GenomicAligns}) {
01187         foreach my $genomic_align2 (@{$ungapped_genomic_align_blocks->[$i+1]->get_all_GenomicAligns}) {
01188           next if ($genomic_align1->dnafrag_id != $genomic_align2->dnafrag_id);
01189           ## $gap is the piece of sequence of this dnafrag between this block and the next one
01190           my $gap;
01191           if ($genomic_align1->dnafrag_strand == 1) {
01192             $gap = $genomic_align2->dnafrag_start - $genomic_align1->dnafrag_end - 1;
01193           } else {
01194             $gap = $genomic_align1->dnafrag_start - $genomic_align2->dnafrag_end - 1;
01195           }
01196           if ($gap) {
01197             $gap = "" if ($gap == 1);
01198             foreach my $genomic_align3 (@{$genomic_align_block->get_all_GenomicAligns}) {
01199               if ($genomic_align1->dnafrag_id == $genomic_align3->dnafrag_id) {
01200                 ## Add (mis)matches for this sequence
01201                 $cigar_lines->{$genomic_align3->dnafrag_id} .= $gap."M";
01202               } else {
01203                 ## Add gaps for others
01204                 $cigar_lines->{$genomic_align3->dnafrag_id} .= $gap."D";
01205               }
01206             }
01207           }
01208         }
01209       }
01210     }
01211 
01212   }
01213 
01214   while (my ($id, $genomic_align) = each %{$genomic_aligns}) {
01215     $genomic_align->cigar_line($cigar_lines->{$id});
01216     $self->add_GenomicAlign($genomic_align);
01217   }
01218 
01219   return $self;
01220 }
01221 
01222 
01223 =head2 get_all_ungapped_GenomicAlignBlocks
01224 
01225   Args       : (optional) listref $genome_dbs
01226   Example    : my $ungapped_genomic_align_blocks =
01227                    $self->get_all_ungapped_GenomicAlignBlocks();
01228   Example    : my $ungapped_genomic_align_blocks =
01229                    $self->get_all_ungapped_GenomicAlignBlocks([$human_genome_db, $mouse_genome_db]);
01230   Description: split the GenomicAlignBlock object into a set of ungapped
01231                alignments. If a list of genome_dbs is provided, only those
01232                sequences will be taken into account. This can be used to extract
01233                ungapped pairwise alignments from multiple alignments.
01234   Returntype : listref of Bio::EnsEMBL::Compara::GenomicAlignBlocks objects
01235   Exceptions : none
01236   Caller     : general
01237   Status     : At risk
01238 
01239 =cut
01240 
01241 sub get_all_ungapped_GenomicAlignBlocks {
01242   my ($self, $genome_dbs) = @_;
01243   my $ungapped_genomic_align_blocks = [];
01244 
01245   my $genomic_aligns = $self->get_all_GenomicAligns;
01246   if ($genome_dbs and @$genome_dbs) {
01247     my $these_genomic_aligns = [];
01248     foreach my $this_genomic_align (@$genomic_aligns) {
01249       if (grep {$this_genomic_align->genome_db->name eq $_->name} @$genome_dbs) {
01250         push(@$these_genomic_aligns, $this_genomic_align);
01251       }
01252     }
01253     if (@$these_genomic_aligns > 1) {
01254       $genomic_aligns = $these_genomic_aligns;
01255     } else {
01256       return [];
01257     }
01258   }
01259   my $aln_length = CORE::length($genomic_aligns->[0]->aligned_sequence("+FAKE_SEQ"));
01260 #   foreach my $this_genomic_align (@$genomic_aligns) {
01261 #     print STDERR join(" - ", $this_genomic_align->dnafrag_start, $this_genomic_align->dnafrag_end,
01262 #         $this_genomic_align->dnafrag_strand, $this_genomic_align->aligned_sequence("+FAKE_SEQ")), "\n";
01263 #   }
01264 
01265   my $aln_pos = 0;
01266   my $gap;
01267   my $end_block_pos;
01268   do {
01269     $end_block_pos = undef;
01270     my $these_genomic_aligns_with_no_gaps;
01271 
01272     ## Get the (next) first gap from all the aligned sequences (sets: $gap_pos, $gap and $genomic_align_block_id)
01273     foreach my $this_genomic_align (@$genomic_aligns) {
01274       my $this_end_block_pos = index($this_genomic_align->aligned_sequence("+FAKE_SEQ"), "-", $aln_pos);
01275       if ($this_end_block_pos == $aln_pos) {
01276         ## try to find the end of the gaps
01277         my $gap_string = substr($this_genomic_align->aligned_sequence("+FAKE_SEQ"), $aln_pos);
01278         ($gap) = $gap_string =~ /^(\-+)/;
01279         my $gap_length = CORE::length($gap);
01280         $this_end_block_pos = $aln_pos+$gap_length;
01281       } else {
01282         $these_genomic_aligns_with_no_gaps->{$this_genomic_align} = $this_genomic_align;
01283       }
01284       $this_end_block_pos = CORE::length($this_genomic_align->aligned_sequence("+FAKE_SEQ")) if ($this_end_block_pos < 0); # no more gaps have been found in this sequence
01285 
01286       
01287       if (!defined($end_block_pos) or $this_end_block_pos < $end_block_pos) {
01288         $end_block_pos = $this_end_block_pos;
01289       }
01290     }
01291 
01292     if (scalar(keys(%$these_genomic_aligns_with_no_gaps)) > 1) {
01293       my $new_genomic_aligns;
01294       my $reference_genomic_align;
01295       foreach my $this_genomic_align (values %$these_genomic_aligns_with_no_gaps) {
01296         my $previous_seq = substr($this_genomic_align->aligned_sequence("+FAKE_SEQ"), 0, $aln_pos );
01297         $previous_seq =~ s/\-//g;
01298         my $dnafrag_start;
01299         my $dnafrag_end;
01300         my $cigar_line;
01301         my $cigar_length = ($end_block_pos - $aln_pos);
01302         $cigar_length = "" if ($cigar_length == 1);
01303         if ($this_genomic_align->dnafrag_strand == 1) {
01304           $dnafrag_start = $this_genomic_align->dnafrag_start + CORE::length($previous_seq);
01305           $dnafrag_end = $dnafrag_start + $end_block_pos - $aln_pos - 1;
01306           $cigar_line = $cigar_length."M";
01307         } else {
01308           $dnafrag_end = $this_genomic_align->dnafrag_end - CORE::length($previous_seq);
01309           $dnafrag_start = $dnafrag_end - $end_block_pos + $aln_pos + 1;
01310           $cigar_line = $cigar_length."M";
01311         }
01312         my $new_genomic_align = new Bio::EnsEMBL::Compara::GenomicAlign(
01313                 -adaptor => $this_genomic_align->adaptor,
01314                 -method_link_species_set => $this_genomic_align->method_link_species_set,
01315                 -dnafrag => $this_genomic_align->dnafrag,
01316                 -dnafrag_start => $dnafrag_start,
01317                 -dnafrag_end => $dnafrag_end,
01318                 -dnafrag_strand => $this_genomic_align->dnafrag_strand,
01319                 -cigar_line => $cigar_line,
01320             );
01321         $reference_genomic_align = $new_genomic_align
01322             if (defined($self->reference_genomic_align) and
01323                 $self->reference_genomic_align == $this_genomic_align);
01324         push(@$new_genomic_aligns, $new_genomic_align);
01325       }
01326       ## Create a new GenomicAlignBlock
01327       my $new_genomic_align_block = new Bio::EnsEMBL::Compara::GenomicAlignBlock(
01328               -method_link_species_set => $self->method_link_species_set,
01329               -length => $end_block_pos - $aln_pos,
01330               -genomic_align_array => $new_genomic_aligns,
01331           );
01332       $new_genomic_align_block->reference_genomic_align($reference_genomic_align) if (defined($reference_genomic_align));
01333       push(@$ungapped_genomic_align_blocks, $new_genomic_align_block);
01334     }
01335     $aln_pos = $end_block_pos;
01336 
01337   } while ($aln_pos < $aln_length); # exit loop if no gap has been found
01338 
01339   return $ungapped_genomic_align_blocks;
01340 }
01341 
01342 
01343 =head2 reverse_complement
01344 
01345   Args       : none
01346   Example    : none
01347   Description: reverse complement the ,
01348                modifying dnafrag_strand and cigar_line of each GenomicAlign in consequence
01349   Returntype : none
01350   Exceptions : none
01351   Caller     : general
01352   Status     : Stable
01353 
01354 =cut
01355 
01356 sub reverse_complement {
01357   my ($self) = @_;
01358 
01359   if (defined($self->{_original_strand})) {
01360     $self->{_original_strand} = 1 - $self->{_original_strand};
01361   } else {
01362     $self->{_original_strand} = 0;
01363   }
01364 
01365   my $gas = $self->get_all_GenomicAligns;
01366   foreach my $ga (@{$gas}) {
01367     $ga->reverse_complement;
01368   }
01369 }
01370 
01371 =head2 restrict_between_alignment_positions
01372 
01373   Arg[1]     : [optional] int $start, refers to the start of the alignment
01374   Arg[2]     : [optional] int $end, refers to the start of the alignment
01375   Arg[3]     : [optional] boolean $skip_empty_GenomicAligns
01376   Example    : none
01377   Description: restrict this GenomicAlignBlock. It returns a new object unless no
01378                restriction is needed. In that case, it returns the original unchanged
01379                object.
01380                This method uses coordinates relative to the alignment itself.
01381                For instance if you have an alignment like:
01382                             1    1    2    2    3
01383                    1   5    0    5    0    5    0
01384                    AAC--CTTGTGGTA-CTACTT-----ACTTT
01385                    AACCCCTT-TGGTATCTACTTACCTAACTTT
01386                and you restrict it between 5 and 25, you will get back a
01387                object containing the following alignment:
01388                             1    1
01389                    1   5    0    5
01390                    CTTGTGGTA-CTACTT----
01391                    CTT-TGGTATCTACTTACCT
01392 
01393                See restrict_between_reference_positions() elsewhere in this document
01394                for an alternative method using absolute genomic coordinates.
01395 
01396                NB: This method works only for GenomicAlignBlock which have been
01397                fetched from the DB as it is adjusting the dnafrag coordinates
01398                and the cigar_line only and not the actual sequences stored in the
01399                object if any. If you want to restrict an object with no coordinates
01400                a simple substr() will do!
01401 
01402   Returntype : Bio::EnsEMBL::Compara::GenomicAlignBlock object
01403   Exceptions : none
01404   Caller     : general
01405   Status     : At risk
01406 
01407 =cut
01408 
01409 sub restrict_between_alignment_positions {
01410   my ($self, $start, $end, $skip_empty_GenomicAligns) = @_;
01411   my $genomic_align_block;
01412   my $new_reference_genomic_align;
01413   my $new_genomic_aligns;
01414 
01415   $start = 1 if (!defined($start) or $start < 1);
01416   $end = $self->length if (!defined($end) or $end > $self->length);
01417 
01418   my $number_of_columns_to_trim_from_the_start = $start - 1;
01419   my $number_of_columns_to_trim_from_the_end = $self->length - $end;
01420 
01421   ## Skip if no restriction is needed. Return original object! We are still going on with the
01422   ## restriction when either excess_at_the_start or excess_at_the_end are 0 as a (multiple)
01423   ## alignment may start or end with gaps in the reference species. In that case, we want to
01424   ## trim these gaps from the alignment as they fall just outside of the region of interest
01425   return $self if ($number_of_columns_to_trim_from_the_start <= 0
01426       and $number_of_columns_to_trim_from_the_end <= 0);
01427 
01428   my $final_alignment_length = $end - $start + 1;
01429 
01430   ## Create a new Bio::EnsEMBL::Compara::GenomicAlignBlock object with restricted GenomicAligns
01431   my $length = $self->{length};
01432   foreach my $this_genomic_align (@{$self->get_all_GenomicAligns}) {
01433     my $new_genomic_align = $this_genomic_align->restrict($start, $end, $length);
01434     if ($self->reference_genomic_align and $this_genomic_align == $self->reference_genomic_align) {
01435       $new_reference_genomic_align = $new_genomic_align;
01436     }
01437     push(@$new_genomic_aligns, $new_genomic_align);
01438   }
01439   $genomic_align_block = new Bio::EnsEMBL::Compara::GenomicAlignBlock(
01440           -method_link_species_set => $self->method_link_species_set,
01441           -genomic_align_array => $new_genomic_aligns,
01442           -group_id => $self->group_id,
01443       -level_id => $self->level_id,
01444       );
01445   $genomic_align_block->{original_dbID} = ($self->dbID or $self->{original_dbID});
01446   $genomic_align_block->{_original_strand} = $self->{_original_strand};
01447   if ($new_reference_genomic_align) {
01448     $genomic_align_block->reference_genomic_align($new_reference_genomic_align);
01449   }
01450   $genomic_align_block->reference_slice($self->reference_slice);
01451 
01452   # The restriction might result in empty GenomicAligns. If the
01453   # skip_empty_GenomicAligns flag is set, remove them from the block.
01454   if ($skip_empty_GenomicAligns) {
01455     my $reference_genomic_align = $genomic_align_block->reference_genomic_align();
01456     my $genomic_align_array = $genomic_align_block->genomic_align_array;
01457     for (my $i=0; $i<@$genomic_align_array; $i++) {
01458       if ($genomic_align_array->[$i]->dnafrag_start > $genomic_align_array->[$i]->dnafrag_end) {
01459         splice(@$genomic_align_array, $i, 1);
01460         $i--;
01461       }
01462     }
01463     $genomic_align_block->reference_genomic_align($reference_genomic_align) if ($reference_genomic_align);
01464   }
01465   $genomic_align_block->length($final_alignment_length);
01466 
01467   return $genomic_align_block;
01468 }
01469 
01470 =head2 _print
01471 
01472   Arg [1]    : none
01473   Example    : $genomic_align->_print
01474   Description: print attributes of the object to the STDOUT. Used for debuging purposes.
01475   Returntype : none
01476   Exceptions : 
01477   Caller     : object::methodname
01478   Status     : At risk
01479 
01480 =cut
01481 
01482 sub _print {
01483   my ($self, $FILEH) = @_;
01484 
01485   $FILEH ||= \*STDOUT;
01486   print $FILEH
01487 "Bio::EnsEMBL::Compara::GenomicAlignBlock object ($self)
01488   dbID = ", ($self->dbID or "-undef-"), "
01489   adaptor = ", ($self->adaptor or "-undef-"), "
01490   method_link_species_set = ", ($self->method_link_species_set or "-undef-"), "
01491   method_link_species_set_id = ", ($self->method_link_species_set_id or "-undef-"), "
01492   genomic_aligns = ", (scalar(@{$self->genomic_align_array}) or "-undef-"), "
01493   score = ", ($self->score or "-undef-"), "
01494   length = ", ($self->length or "-undef-"), "
01495   alignments: \n";
01496   foreach my $this_genomic_align (@{$self->genomic_align_array()}) {
01497     my $slice = $this_genomic_align->get_Slice;
01498     if ($self->reference_genomic_align and $self->reference_genomic_align == $this_genomic_align) {
01499       print $FILEH "    * ", $this_genomic_align->genome_db->name, " ",
01500           ($slice?$slice->name:"--error--"), "\n";
01501     } else {
01502       print $FILEH "    - ", $this_genomic_align->genome_db->name, " ",
01503           ($slice?$slice->name:"--error--"), "\n";
01504     }
01505   }
01506 
01507 }
01508 
01509 
01510 #####################################################################
01511 #####################################################################
01512 
01513 =head1 METHODS FOR BACKWARDS COMPATIBILITY
01514 
01515 Consensus and Query DnaFrag are no longer used. DO NOT USE THOSE METHODS IN NEW SCRIPTS, THEY WILL DISSAPEAR!!
01516 
01517 For backwards compatibility, consensus_genomic_align correponds to the lower genome_db_id by convention.
01518 This convention works for pairwise alignment only! Trying to use the old API methods for multiple
01519 alignments will throw an exception.
01520 
01521 =cut
01522 
01523 #####################################################################
01524 #####################################################################
01525 
01526 
01527 =head2 get_old_consensus_genomic_align [FOR BACKWARDS COMPATIBILITY ONLY]
01528  
01529   Arg [1]    : none
01530   Example    : $old_consensus_genomic_aligns = $genomic_align_group->get_old_consensus_genomic_align();
01531   Description: get the Bio::EnsEMBL::Compara::GenomicAlign object following the convention for backwards
01532                compatibility
01533   Returntype : Bio::EnsEMBL::Compara::GenomicAlign object
01534   Exceptions : 
01535   Caller     : general
01536  
01537 =cut
01538 
01539 sub get_old_consensus_genomic_align {
01540   my ($self) = @_;
01541 
01542   my $genomic_aligns = $self->get_all_GenomicAligns;
01543   if (!@$genomic_aligns) {
01544     throw "Bio::EnsEMBL::Compara::GenomicAlignBlock ($self) does not have any associated".
01545         " Bio::EnsEMBL::Compara::GenomicAlign";
01546   }
01547 
01548   if (scalar(@{$genomic_aligns}) != 2) {
01549     throw "Trying to get old_consensus_genomic_align from Bio::EnsEMBL::Compara::GenomicAlignBlock".
01550         " ($self) holding a multiple alignment";
01551   }
01552 
01553   if ($genomic_aligns->[0]->dnafrag->genome_db->dbID > $genomic_aligns->[1]->dnafrag->genome_db->dbID) {
01554     return $genomic_aligns->[1];
01555 
01556   } elsif ($genomic_aligns->[0]->dnafrag->genome_db->dbID < $genomic_aligns->[1]->dnafrag->genome_db->dbID) {
01557     return $genomic_aligns->[0];
01558 
01559   ## If they belongs to the same genome_db, use the dnafrag_id instead
01560   } elsif ($genomic_aligns->[0]->dnafrag->dbID > $genomic_aligns->[1]->dnafrag->dbID) {
01561     return $genomic_aligns->[1];
01562 
01563   } elsif ($genomic_aligns->[0]->dnafrag->dbID < $genomic_aligns->[1]->dnafrag->dbID) {
01564     return $genomic_aligns->[0];
01565 
01566   ## If they belongs to the same genome_db and dnafrag, use the dnafrag_start instead
01567   } elsif ($genomic_aligns->[0]->dnafrag_start > $genomic_aligns->[1]->dnafrag_start) {
01568     return $genomic_aligns->[1];
01569 
01570   } elsif ($genomic_aligns->[0]->dnafrag_start < $genomic_aligns->[1]->dnafrag_start) {
01571     return $genomic_aligns->[0];
01572 
01573   ## If they belongs to the same genome_db and dnafrag and have the same danfrag_start, use the dnafrag_end instead
01574   } elsif ($genomic_aligns->[0]->dnafrag_end > $genomic_aligns->[1]->dnafrag_end) {
01575     return $genomic_aligns->[1];
01576 
01577   } elsif ($genomic_aligns->[0]->dnafrag_end < $genomic_aligns->[1]->dnafrag_end) {
01578     return $genomic_aligns->[0];
01579 
01580   ## If everithing else fails, use the dnafrag_strand
01581   } elsif ($genomic_aligns->[0]->dnafrag_strand > $genomic_aligns->[1]->dnafrag_strand) {
01582     return $genomic_aligns->[1];
01583 
01584   } elsif ($genomic_aligns->[0]->dnafrag_strand < $genomic_aligns->[1]->dnafrag_strand) {
01585     return $genomic_aligns->[0];
01586 
01587   # Whatever, they are the same. Use 0 for consensus and 1 for query
01588   } else {
01589     return $genomic_aligns->[0];
01590   }
01591 }
01592 
01593 
01594 =head2 get_old_query_genomic_align [FOR BACKWARDS COMPATIBILITY ONLY]
01595  
01596   Arg [1]    : none
01597   Example    : $old_query_genomic_aligns = $genomic_align_group->get_old_query_genomic_align();
01598   Description: get the Bio::EnsEMBL::Compara::GenomicAlign object following the convention for backwards
01599                compatibility
01600   Returntype : Bio::EnsEMBL::Compara::GenomicAlign object
01601   Exceptions : 
01602   Caller     : general
01603  
01604 =cut
01605 
01606 sub get_old_query_genomic_align {
01607   my ($self) = @_;
01608 
01609   my $genomic_aligns = $self->get_all_GenomicAligns;
01610   if (!@$genomic_aligns) {
01611     throw "Bio::EnsEMBL::Compara::GenomicAlignBlock ($self) does not have any associated".
01612         " Bio::EnsEMBL::Compara::GenomicAlign";
01613   }
01614   
01615   if (scalar(@{$genomic_aligns}) != 2) {
01616     throw "Trying to get old_consensus_genomic_align from Bio::EnsEMBL::Compara::GenomicAlignBlock".
01617         " ($self) holding a multiple alignment";
01618   }
01619 
01620   if ($genomic_aligns->[0]->dnafrag->genome_db->dbID > $genomic_aligns->[1]->dnafrag->genome_db->dbID) {
01621     return $genomic_aligns->[0];
01622 
01623   } elsif ($genomic_aligns->[0]->dnafrag->genome_db->dbID < $genomic_aligns->[1]->dnafrag->genome_db->dbID) {
01624     return $genomic_aligns->[1];
01625 
01626   ## If they belongs to the same genome_db, use the dnafrag_id instead
01627   } elsif ($genomic_aligns->[0]->dnafrag->dbID > $genomic_aligns->[1]->dnafrag->dbID) {
01628     return $genomic_aligns->[0];
01629 
01630   } elsif ($genomic_aligns->[0]->dnafrag->dbID < $genomic_aligns->[1]->dnafrag->dbID) {
01631     return $genomic_aligns->[1];
01632 
01633   ## If they belongs to the same genome_db and dnafrag, use the dnafrag_start instead
01634   } elsif ($genomic_aligns->[0]->dnafrag_start > $genomic_aligns->[1]->dnafrag_start) {
01635     return $genomic_aligns->[0];
01636 
01637   } elsif ($genomic_aligns->[0]->dnafrag_start < $genomic_aligns->[1]->dnafrag_start) {
01638     return $genomic_aligns->[1];
01639 
01640   ## If they belongs to the same genome_db and dnafrag and have the same danfrag_start, use the dnafrag_end instead
01641   } elsif ($genomic_aligns->[0]->dnafrag_end > $genomic_aligns->[1]->dnafrag_end) {
01642     return $genomic_aligns->[0];
01643 
01644   } elsif ($genomic_aligns->[0]->dnafrag_end < $genomic_aligns->[1]->dnafrag_end) {
01645     return $genomic_aligns->[1];
01646 
01647   ## If everithing else fails, use the dnafrag_strand
01648   } elsif ($genomic_aligns->[0]->dnafrag_strand > $genomic_aligns->[1]->dnafrag_strand) {
01649     return $genomic_aligns->[0];
01650 
01651   } elsif ($genomic_aligns->[0]->dnafrag_strand < $genomic_aligns->[1]->dnafrag_strand) {
01652     return $genomic_aligns->[1];
01653 
01654   # Whatever, they are the same. Use 0 for consensus and 1 for query
01655   } else {
01656     return $genomic_aligns->[1];
01657   }
01658 }
01659 
01660 1;