Archive Ensembl HomeArchive Ensembl Home
GenomicAlign.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::GenomicAlign - Defines one of the sequences involved in a genomic alignment
00022 
00023 =head1 SYNOPSIS
00024 
00025   use Bio::EnsEMBL::Compara::GenomicAlign; 
00026   my $genomic_align = new Bio::EnsEMBL::Compara::GenomicAlign(
00027           -adaptor => $genomic_align_adaptor,
00028           -genomic_align_block => $genomic_align_block,
00029           -method_link_species_set => $method_link_species_set,
00030           -dnafrag => $dnafrag,
00031           -dnafrag_start => 100001,
00032           -dnafrag_end => 100050,
00033           -dnafrag_strand => -1,
00034           -aligned_sequence => "TTGCAGGTAGGCCATCTGCAAGC----TGAGGAGCAAGGACTCCAGTCGGAGTC"
00035           -visible => 1,
00036         );
00037 
00038 
00039 SET VALUES
00040   $genomic_align->adaptor($adaptor);
00041   $genomic_align->dbID(12);
00042   $genomic_align->genomic_align_block($genomic_align_block);
00043   $genomic_align->genomic_align_block_id(1032);
00044   $genomic_align->method_link_species_set($method_link_species_set);
00045   $genomic_align->method_link_species_set_id(3);
00046   $genomic_align->dnafrag($dnafrag);
00047   $genomic_align->dnafrag_id(134);
00048   $genomic_align->dnafrag_start(100001);
00049   $genomic_align->dnafrag_end(100050);
00050   $genomic_align->dnafrag_strand(-1);
00051   $genomic_align->aligned_sequence("TTGCAGGTAGGCCATCTGCAAGC----TGAGGAGCAAGGACTCCAGTCGGAGTC");
00052   $genomic_align->original_sequence("TTGCAGGTAGGCCATCTGCAAGCTGAGGAGCAAGGACTCCAGTCGGAGTC");
00053   $genomic_align->cigar_line("23M4D27M");
00054   $genomic_align->visible(1);
00055 
00056 GET VALUES
00057   $adaptor = $genomic_align->adaptor;
00058   $dbID = $genomic_align->dbID;
00059   $genomic_align_block = $genomic_align->genomic_block;
00060   $genomic_align_block_id = $genomic_align->genomic_align_block_id;
00061   $method_link_species_set = $genomic_align->method_link_species_set;
00062   $method_link_species_set_id = $genomic_align->method_link_species_set_id;
00063   $dnafrag = $genomic_align->dnafrag;
00064   $dnafrag_id = $genomic_align->dnafrag_id;
00065   $dnafrag_start = $genomic_align->dnafrag_start;
00066   $dnafrag_end = $genomic_align->dnafrag_end;
00067   $dnafrag_strand = $genomic_align->dnafrag_strand;
00068   $aligned_sequence = $genomic_align->aligned_sequence;
00069   $original_sequence = $genomic_align->original_sequence;
00070   $cigar_line = $genomic_align->cigar_line;
00071   $visible = $genomic_align->visible;
00072   $slice = $genomic_align->get_Slice();
00073 
00074 =head1 DESCRIPTION
00075 
00076 The GenomicAlign object stores information about a single sequence within an alignment.
00077 
00078 =head1 OBJECT ATTRIBUTES
00079 
00080 =over
00081 
00082 =item dbID
00083 
00084 corresponds to genomic_align.genomic_align_id
00085 
00086 =item adaptor
00087 
00088 Bio::EnsEMBL::Compara::DBSQL::GenomicAlignAdaptor object to access DB
00089 
00090 =item genomic_align_block_id
00091 
00092 corresponds to genomic_align_block.genomic_align_block_id (ext. reference)
00093 
00094 =item genomic_align_block
00095 
00096 Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlock object corresponding to genomic_align_block_id
00097 
00098 =item method_link_species_set_id
00099 
00100 corresponds to method_link_species_set.method_link_species_set_id (external ref.)
00101 
00102 =item method_link_species_set
00103 
00104 Bio::EnsEMBL::Compara::DBSQL::MethodLinkSpeciesSet object corresponding to method_link_species_set_id
00105 
00106 =item dnafrag_id
00107 
00108 corresponds to dnafrag.dnafrag_id (external ref.)
00109 
00110 =item dnafrag
00111 
00112 Bio::EnsEMBL::Compara::DnaFrag object corresponding to dnafrag_id
00113 
00114 =item dnafrag_start
00115 
00116 corresponds to genomic_align.dnafrag_start
00117 
00118 =item dnafrag_end
00119 
00120 corresponds to genomic_align.dnafrag_end
00121 
00122 =item dnafrag_strand
00123 
00124 corresponds to genomic_align.dnafrag_strand
00125 
00126 =item cigar_line
00127 
00128 corresponds to genomic_align.cigar_line
00129 
00130 =item visible
00131 
00132 corresponds to genomic_align.visible
00133 
00134 =item aligned_sequence
00135 
00136 corresponds to the sequence rebuilt using dnafrag and cigar_line
00137 
00138 =item original_sequence
00139 
00140 corresponds to the original sequence. It can be rebuilt from the aligned_sequence, the dnafrag object or can be used
00141 in conjuction with cigar_line to get the aligned_sequence.
00142 
00143 =back
00144 
00145 =head1 APPENDIX
00146 
00147 The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
00148 
00149 =cut
00150 
00151 
00152 # Let the code begin...
00153 
00154 
00155 package Bio::EnsEMBL::Compara::GenomicAlign;
00156 use strict;
00157 
00158 use Bio::EnsEMBL::Utils::Exception qw(throw warning deprecate verbose);
00159 use Bio::EnsEMBL::Utils::Argument qw(rearrange);
00160 use Scalar::Util qw(weaken);
00161 use Bio::EnsEMBL::Compara::BaseGenomicAlignSet;
00162 use Bio::EnsEMBL::Compara::MethodLinkSpeciesSet;
00163 use Bio::EnsEMBL::Mapper;
00164 
00165 
00166 =head2 new (CONSTRUCTOR)
00167 
00168   Arg [-DBID] : (opt.) int $dbID (the database internal ID for this object)
00169   Arg [-ADAPTOR]
00170               : (opt.) Bio::EnsEMBL::Compara::DBSQL::GenomicAlignAdaptor $adaptor
00171                 (the adaptor for connecting to the database)
00172   Arg [-GENOMIC_ALIGN_BLOCK]
00173               : (opt.) Bio::EnsEMBL::Compara::GenomicAlignBlock $genomic_align_block
00174                 (the block to which this Bio::EnsEMBL::Compara::GenomicAlign object
00175                 belongs to)
00176   Arg [-GENOMIC_ALIGN_BLOCK_ID]
00177               : (opt.) int $genomic_align_block_id (the database internal ID for the
00178                 $genomic_align_block)
00179   Arg [-METHOD_LINK_SPECIES_SET]
00180               : (opt.) Bio::EnsEMBL::Compara::MethodLinkSpeciesSet $mlss
00181                 (this defines the type of alignment and the set of species used
00182                 to get this GenomicAlignBlock)
00183   Arg [-METHOD_LINK_SPECIES_SET_ID]
00184               : (opt.) int $mlss_id (the database internal ID for the $mlss)
00185   Arg [-DNAFRAG]
00186               : (opt.) Bio::EnsEMBL::Compara::DnaFrag $dnafrag (the genomic
00187                 sequence object to which this object refers to)
00188   Arg [-DNAFRAG_ID]
00189               : (opt.) int $dnafrag_id (the database internal ID for the $dnafrag)
00190   Arg [-DNAFRAG_START]
00191               : (opt.) int $dnafrag_start (the starting position of this
00192                 Bio::EnsEMBL::Compara::GenomicAling within its corresponding $dnafrag)
00193   Arg [-DNAFRAG_END]
00194               : (opt.) int $dnafrag_end (the ending position of this
00195                 Bio::EnsEMBL::Compara::GenomicAling within its corresponding $dnafrag)
00196   Arg [-DNAFRAG_STRAND]
00197               : (opt.) int $dnafrag_strand (1 or -1; defines in which strand of its
00198                 corresponding $dnafrag this Bio::EnsEMBL::Compara::GenomicAlign is)
00199   Arg [-ALIGNED_SEQUENCE]
00200               : (opt.) string $aligned_sequence (the sequence of this object, including
00201                 gaps and all)
00202   Arg [-CIGAR_LINE]
00203               : (opt.) string $cigar_line (a compressed way of representing the indels in
00204                 the $aligned_sequence of this object)
00205   Arg [-VISIBLE]
00206               : (opt.) int $visible. Used in self alignments to ensure only one Bio::EnsEMBL::Compara::GenomicAlignBlock is visible when you have more than 1 block covering the same region.
00207   Example     : my $genomic_align = new Bio::EnsEMBL::Compara::GenomicAlign(
00208                         -adaptor => $genomic_align_adaptor,
00209                         -genomic_align_block => $genomic_align_block,
00210                         -method_link_species_set => $method_link_species_set,
00211                         -dnafrag => $dnafrag,
00212                         -dnafrag_start => 100001,
00213                         -dnafrag_end => 100050,
00214                         -dnafrag_strand => -1,
00215                         -aligned_sequence => "TTGCAGGTAGGCCATCTGCAAGC----TGAGGAGCAAGGACTCCAGTCGGAGTC"
00216                         -visible => 1,
00217                       );
00218   Description : Creates a new Bio::EnsEMBL::Compara::GenomicAlign object
00219   Returntype  : Bio::EnsEMBL::Compara::GenomicAlign object
00220   Exceptions  : none
00221   Caller      : general
00222   Status      : Stable
00223 
00224 =cut
00225 
00226 sub new {
00227     my($class, @args) = @_;
00228 
00229     my $self = {};
00230     bless $self,$class;
00231 
00232     my ($cigar_line, $adaptor,
00233         $dbID, $genomic_align_block, $genomic_align_block_id, $method_link_species_set,
00234         $method_link_species_set_id, $dnafrag, $dnafrag_id,
00235         $dnafrag_start, $dnafrag_end, $dnafrag_strand,
00236         $aligned_sequence, $visible, $node_id ) = 
00237       
00238       rearrange([qw(
00239           CIGAR_LINE ADAPTOR
00240           DBID GENOMIC_ALIGN_BLOCK GENOMIC_ALIGN_BLOCK_ID METHOD_LINK_SPECIES_SET
00241           METHOD_LINK_SPECIES_SET_ID DNAFRAG DNAFRAG_ID
00242           DNAFRAG_START DNAFRAG_END DNAFRAG_STRAND
00243           ALIGNED_SEQUENCE VISIBLE NODE_ID)], @args);
00244 
00245     $self->adaptor( $adaptor ) if defined $adaptor;
00246     $self->cigar_line( $cigar_line ) if defined $cigar_line;
00247     
00248     $self->dbID($dbID) if (defined($dbID));
00249     $self->genomic_align_block($genomic_align_block) if (defined($genomic_align_block));
00250     $self->genomic_align_block_id($genomic_align_block_id) if (defined($genomic_align_block_id));
00251     $self->method_link_species_set($method_link_species_set) if (defined($method_link_species_set));
00252     $self->method_link_species_set_id($method_link_species_set_id) if (defined($method_link_species_set_id));
00253     $self->dnafrag($dnafrag) if (defined($dnafrag));
00254     $self->dnafrag_id($dnafrag_id) if (defined($dnafrag_id));
00255     $self->dnafrag_start($dnafrag_start) if (defined($dnafrag_start));
00256     $self->dnafrag_end($dnafrag_end) if (defined($dnafrag_end));
00257     $self->dnafrag_strand($dnafrag_strand) if (defined($dnafrag_strand));
00258     $self->aligned_sequence($aligned_sequence) if (defined($aligned_sequence));
00259     $self->visible($visible) if (defined($visible));
00260     $self->node_id($node_id) if (defined($node_id));
00261 
00262     return $self;
00263 }
00264 
00265 =head2 new_fast
00266 
00267   Arg [1]    : hash reference $hashref
00268   Example    : none
00269   Description: This is an ultra fast constructor which requires knowledge of
00270                the objects internals to be used.
00271   Returntype :
00272   Exceptions : none
00273   Caller     :
00274   Status     : Stable
00275 
00276 =cut
00277 
00278 sub new_fast {
00279   my $class = shift;
00280   my $hashref = shift;
00281 
00282   return bless $hashref, $class;
00283 }
00284 
00285 
00286 =head2 copy (CONSTRUCTOR)
00287 
00288   Arg         : -none-
00289   Example     : my $new_genomic_align = $genomic_align->copy();
00290   Description : Create a new object with the same attributes
00291                 as this one.
00292   Returntype  : Bio::EnsEMBL::Compara::GenomicAlign (or subclassed) object
00293   Exceptions  :
00294   Status      : Stable
00295 
00296 =cut
00297 
00298 sub copy {
00299   my ($self) = @_;
00300   my $new_copy = {};
00301   bless $new_copy, ref($self);
00302 
00303   while (my ($key, $value) = each %$self) {
00304     $new_copy->{$key} = $value;
00305   }
00306 
00307   return $new_copy;
00308 }
00309 
00310 
00311 =head2 adaptor
00312 
00313   Arg [1]    : Bio::EnsEMBL::Compara::DBSQL::GenomicAlignAdaptor
00314   Example    : $adaptor = $genomic_align->adaptor;
00315   Example    : $genomic_align->adaptor($adaptor);
00316   Description: Getter/Setter for the adaptor this object uses for database
00317                interaction.
00318   Returntype : Bio::EnsEMBL::Compara::DBSQL::GenomicAlignAdaptor
00319   Exceptions : thrown if $adaptor is not a
00320                Bio::EnsEMBL::Compara::DBSQL::GenomicAlignAdaptor object
00321   Caller     : object->methodname
00322   Status     : Stable
00323 
00324 =cut
00325 
00326 sub adaptor {
00327   my $self = shift;
00328 
00329   if (@_) {
00330     $self->{'adaptor'} = shift;
00331     throw($self->{'adaptor'}." is not a Bio::EnsEMBL::Compara::DBSQL::GenomicAlignAdaptor object")
00332         if ($self->{'adaptor'} and !UNIVERSAL::isa($self->{'adaptor'}, "Bio::EnsEMBL::Compara::DBSQL::GenomicAlignAdaptor"));
00333   }
00334 
00335   return $self->{'adaptor'};
00336 }
00337 
00338 
00339 =head2 dbID
00340 
00341   Arg [1]    : integer $dbID
00342   Example    : $dbID = $genomic_align->dbID;
00343   Example    : $genomic_align->dbID(12);
00344   Description: Getter/Setter for the attribute dbID
00345   Returntype : integer
00346   Exceptions : none
00347   Caller     : object->methodname
00348   Status     : Stable
00349 
00350 =cut
00351 
00352 sub dbID {
00353   my ($self, $dbID) = @_;
00354 
00355   if (defined($dbID)) {
00356      $self->{'dbID'} = $dbID;
00357   }
00358 
00359   return $self->{'dbID'};
00360 }
00361 
00362 
00363 =head2 genomic_align_block
00364 
00365   Arg [1]    : Bio::EnsEMBL::Compara::GenomicAlignBlock $genomic_align_block
00366   Example    : $genomic_align_block = $genomic_align->genomic_align_block;
00367   Example    : $genomic_align->genomic_align_block($genomic_align_block);
00368   Description: Getter/Setter for the attribute genomic_align_block
00369   Returntype : Bio::EnsEMBL::Compara::GenomicAlignBlock object. If no
00370                argument is given, the genomic_align_block is not defined but
00371                both the genomic_align_block_id and the adaptor are, it tries
00372                to fetch the data using the genomic_align_block_id.
00373   Exception  : throws if $genomic_align_block is not a
00374                Bio::EnsEMBL::Compara::GenomicAlignBlock object or if 
00375                $genomic_align_block does not match a previously defined
00376                genomic_align_block_id
00377   Warning    : warns if getting data from other sources fails.
00378   Caller     : object->methodname
00379   Status     : Stable
00380 
00381 =cut
00382 
00383 sub genomic_align_block {
00384   my ($self, $genomic_align_block) = @_;
00385 
00386   if (defined($genomic_align_block)) {
00387     throw("$genomic_align_block is not a Bio::EnsEMBL::Compara::BaseGenomicAlignSet object")
00388         if (!$genomic_align_block->isa("Bio::EnsEMBL::Compara::BaseGenomicAlignSet"));
00389     weaken($self->{'genomic_align_block'} = $genomic_align_block);
00390 
00391     ## Add adaptor to genomic_align_block object if possible and needed
00392     if (!defined($genomic_align_block->{'adaptor'}) and !defined($genomic_align_block->{'_adaptor'}) and defined($self->{'adaptor'})) {
00393       $genomic_align_block->adaptor($self->adaptor->db->get_GenomicAlignBlockAdaptor);
00394     }
00395 
00396     if ($genomic_align_block->isa("Bio::EnsEMBL::Compara::GenomicAlignBlock")) {
00397     if ($self->{'genomic_align_block_id'}) {
00398         if (!$self->{'genomic_align_block'}->{'dbID'}) {
00399         $self->{'genomic_align_block'}->dbID($self->{'genomic_align_block_id'});
00400         }
00401         #       warning("Defining both genomic_align_block_id and genomic_align_block");
00402         throw("dbID of genomic_align_block object does not match previously defined".
00403           " genomic_align_block_id. If you want to override a".
00404           " Bio::EnsEMBL::Compara::GenomicAlign object, you can reset the ".
00405           "genomic_align_block_id using \$genomic_align->genomic_align_block_id(0)")
00406           if ($self->{'genomic_align_block'}->{'dbID'} != $self->{'genomic_align_block_id'});
00407     } else {
00408         $self->{'genomic_align_block_id'} = $genomic_align_block->{'dbID'};
00409     }
00410     }
00411 
00412   } elsif (!defined($self->{'genomic_align_block'})) {
00413     # Try to get the genomic_align_block from other sources...
00414     if (defined($self->genomic_align_block_id) and defined($self->{'adaptor'})) {
00415       # ...from the genomic_align_block_id. Uses genomic_align_block_id function
00416       # and not the attribute in the <if> clause because the attribute can be retrieved from other
00417       # sources if it has not been set before.
00418       my $genomic_align_block_adaptor = $self->{'adaptor'}->db->get_GenomicAlignBlockAdaptor;
00419       $self->{'genomic_align_block'} = $genomic_align_block_adaptor->fetch_by_dbID(
00420               $self->{'genomic_align_block_id'});
00421     } else {
00422 #      warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->genomic_align_block".
00423 #          " You either have to specify more information (see perldoc for".
00424 #          " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00425     }
00426   }
00427 
00428   return $self->{'genomic_align_block'};
00429 }
00430 
00431 
00432 =head2 genomic_align_block_id
00433 
00434   Arg [1]    : integer $genomic_align_block_id
00435   Example    : $genomic_align_block_id = $genomic_align->genomic_align_block_id;
00436   Example    : $genomic_align->genomic_align_block_id(1032);
00437   Description: Getter/Setter for the attribute genomic_align_block_id. If no
00438                argument is given and the genomic_align_block_id is not defined, it
00439                tries to get the data from other sources like the corresponding
00440                Bio::EnsEMBL::Compara::GenomicAlignBlock object or the database using
00441                the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object.
00442                Use 0 as argument to clear this attribute.
00443   Returntype : integer
00444   Exceptions : thrown if $genomic_align_block_id does not match a previously defined
00445                genomic_align_block
00446   Warning    : warns if getting data from other sources fails.
00447   Caller     : object->methodname
00448   Status     : Stable
00449 
00450 =cut
00451 
00452 sub genomic_align_block_id {
00453   my ($self, $genomic_align_block_id) = @_;
00454 
00455   if (defined($genomic_align_block_id)) {
00456 
00457     $self->{'genomic_align_block_id'} = ($genomic_align_block_id or undef);
00458     if (defined($self->{'genomic_align_block'}) and $self->{'genomic_align_block_id'}) {
00459 #       warning("Defining both genomic_align_block_id and genomic_align_block");
00460       throw("genomic_align_block_id does not match previously defined genomic_align_block object")
00461           if ($self->{'genomic_align_block'} and
00462               $self->{'genomic_align_block'}->dbID != $self->{'genomic_align_block_id'});
00463     }
00464   } elsif (!($self->{'genomic_align_block_id'})) {
00465     # Try to get the ID from other sources...
00466     if (defined($self->{'genomic_align_block'})) {
00467     if ($self->{genomic_align_block}->isa("Bio::EnsEMBL::Compara::GenomicAlignBlock")and defined($self->{'genomic_align_block'}->dbID)) {
00468         # ...from the corresponding Bio::EnsEMBL::Compara::GenomicAlignBlock object
00469         $self->{'genomic_align_block_id'} = $self->{'genomic_align_block'}->dbID;
00470     }
00471     } elsif (defined($self->{'adaptor'}) and defined($self->{'dbID'})) {
00472       # ...from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object
00473       $self->adaptor->retrieve_all_direct_attributes($self);
00474     } else {
00475 #      warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->genomic_align_block_id".
00476 #          " You either have to specify more information (see perldoc for".
00477 #          " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00478     }
00479   }
00480 
00481   return $self->{'genomic_align_block_id'};
00482 }
00483 
00484 
00485 =head2 method_link_species_set
00486 
00487   Arg [1]    : Bio::EnsEMBL::Compara::MethodLinkSpeciesSet $method_link_species_set
00488   Example    : $method_link_species_set = $genomic_align->method_link_species_set;
00489   Example    : $genomic_align->method_link_species_set($method_link_species_set);
00490   Description: Getter/Setter for the attribute method_link_species_set. If no
00491                argument is given and the method_link_species_set is not defined, it
00492                tries to get the data from other sources like the corresponding
00493                Bio::EnsEMBL::Compara::GenomicAlignBlock object or from
00494                the method_link_species_set_id.
00495   Returntype : Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object
00496   Exceptions : thrown if $method_link_species_set is not a
00497                Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object or if 
00498                $method_link_species_set does not match a previously defined
00499                method_link_species_set_id
00500   Warning    : warns if getting data from other sources fails.
00501   Caller     : object->methodname
00502   Status     : Stable
00503 
00504 =cut
00505 
00506 sub method_link_species_set {
00507   my ($self, $method_link_species_set) = @_;
00508 
00509   if (defined($method_link_species_set)) {
00510     throw("$method_link_species_set is not a Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object")
00511         if (!$method_link_species_set->isa("Bio::EnsEMBL::Compara::MethodLinkSpeciesSet"));
00512     $self->{'method_link_species_set'} = $method_link_species_set;
00513     if ($self->{'method_link_species_set_id'}) {
00514       if (!$self->{'method_link_species_set'}->dbID) {
00515         $self->{'method_link_species_set'}->dbID($self->{'method_link_species_set_id'});
00516       } else {
00517         $self->{'method_link_species_set_id'} = $self->{'method_link_species_set'}->dbID();
00518       }
00519     } else {
00520       $self->{'method_link_species_set_id'} = $self->{'method_link_species_set'}->dbID;
00521     }
00522   
00523   } elsif (!defined($self->{'method_link_species_set'})) {
00524     # Try to get the object from other sources...
00525     if (defined($self->genomic_align_block) and ($self->{'genomic_align_block'}->method_link_species_set)) {
00526       # ...from the corresponding Bio::EnsEMBL::Compara::GenomicAlignBlock object. Uses genomic_align_block
00527       # function and not the attribute in the <if> clause because the attribute can be retrieved from other
00528       # sources if it has not been already defined.
00529       $self->{'method_link_species_set'} = $self->genomic_align_block->method_link_species_set;
00530     } elsif (defined($self->method_link_species_set_id) and defined($self->{'adaptor'})) {
00531       # ...from the method_link_species_set_id. Uses method_link_species_set_id function and not the attribute
00532       # in the <if> clause because the attribute can be retrieved from other sources if it has not been
00533       # already defined.
00534       my $method_link_species_set_adaptor = $self->adaptor->db->get_MethodLinkSpeciesSetAdaptor;
00535       $self->{'method_link_species_set'} = $method_link_species_set_adaptor->fetch_by_dbID(
00536               $self->{'method_link_species_set_id'});
00537     } else {
00538       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->method_link_species_set".
00539           " You either have to specify more information (see perldoc for".
00540           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00541     }
00542   }
00543 
00544   return $self->{'method_link_species_set'};
00545 }
00546 
00547 
00548 =head2 method_link_species_set_id
00549 
00550   Arg [1]    : integer $method_link_species_set_id
00551   Example    : $method_link_species_set_id = $genomic_align->method_link_species_set_id;
00552   Example    : $genomic_align->method_link_species_set_id(3);
00553   Description: Getter/Setter for the attribute method_link_species_set_id. If no
00554                argument is given and the method_link_species_set_id is not defined, it
00555                tries to get the data from other sources like the corresponding
00556                Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object or the database
00557                using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object.
00558                Use 0 as argument to clear this attribute.
00559   Returntype : integer
00560   Exceptions : thrown if $method_link_species_set_id does not match a previously defined
00561                method_link_species_set
00562   Warning    : warns if getting data from other sources fails.
00563   Caller     : object->methodname
00564   Status     : Stable
00565 
00566 =cut
00567 
00568 sub method_link_species_set_id {
00569   my ($self, $method_link_species_set_id) = @_;
00570 
00571   if (defined($method_link_species_set_id)) {
00572     $self->{'method_link_species_set_id'} = $method_link_species_set_id;
00573     if (defined($self->{'method_link_species_set'}) and $self->{'method_link_species_set_id'}) {
00574       $self->{'method_link_species_set'} = undef;
00575     }
00576   } elsif (!$self->{'method_link_species_set_id'}) {
00577     # Try to get the ID from other sources...
00578     if (defined($self->{'method_link_species_set'}) and $self->{'method_link_species_set'}->dbID) {
00579       # ...from the corresponding Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object
00580       $self->{'method_link_species_set_id'} = $self->{'method_link_species_set'}->dbID;
00581     } elsif (defined($self->{'dbID'}) and defined($self->{'adaptor'})) {
00582       # ...from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object
00583       $self->adaptor->retrieve_all_direct_attributes($self);
00584     } else {
00585       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->method_link_species_set_id".
00586           " You either have to specify more information (see perldoc for".
00587           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00588     }
00589   }
00590 
00591   return $self->{'method_link_species_set_id'};
00592 }
00593 
00594 
00595 =head2 genome_db
00596 
00597   Arg [1]    : Bio::EnsEMBL::Compara::GenomeDB $genome_db
00598   Example    : $genome_db = $genomic_align->genome_db;
00599   Example    : $genomic_align->genome_db($genome_db);
00600   Description: Getter/Setter for the attribute genome_db of
00601                the dnafrag. This method is a short cut for
00602                $genomic_align->dnafrag->genome_db()
00603   Returntype : Bio::EnsEMBL::Compara::DnaFrag object
00604   Exceptions : thrown if $genomic_align->dnafrag is not
00605                defined and cannot be fetched from other
00606                sources.
00607   Caller     : object->methodname
00608   Status     : Stable
00609 
00610 =cut
00611 
00612 sub genome_db {
00613   my ($self, $genome_db) = @_;
00614 
00615   if (defined($genome_db)) {
00616     throw("$genome_db is not a Bio::EnsEMBL::Compara::GenomeDB object")
00617         if (!UNIVERSAL::isa($genome_db, "Bio::EnsEMBL::Compara::GenomeDB"));
00618     my $dnafrag = $self->dnafrag();
00619     if (!$dnafrag) {
00620       throw("Cannot set genome_db if dnafrag does not exist");
00621     } else {
00622       $self->dnafrag->genome_db($genome_db);
00623     }
00624   }
00625   return $self->dnafrag->genome_db;
00626 }
00627 
00628 
00629 =head2 dnafrag
00630 
00631   Arg [1]    : Bio::EnsEMBL::Compara::DnaFrag $dnafrag
00632   Example    : $dnafrag = $genomic_align->dnafrag;
00633   Example    : $genomic_align->dnafrag($dnafrag);
00634   Description: Getter/Setter for the attribute dnafrag. If no
00635                argument is given, the dnafrag is not defined but
00636                both the dnafrag_id and the adaptor are, it tries
00637                to fetch the data using the dnafrag_id
00638   Returntype : Bio::EnsEMBL::Compara::DnaFrag object
00639   Exceptions : thrown if $dnafrag is not a Bio::EnsEMBL::Compara::DnaFrag
00640                object or if $dnafrag does not match a previously defined
00641                dnafrag_id
00642   Warning    : warns if getting data from other sources fails.
00643   Caller     : object->methodname
00644   Status     : Stable
00645 
00646 =cut
00647 
00648 sub dnafrag {
00649   my ($self, $dnafrag) = @_;
00650 
00651   if (defined($dnafrag)) {
00652     throw("$dnafrag is not a Bio::EnsEMBL::Compara::DnaFrag object")
00653         if (!$dnafrag->isa("Bio::EnsEMBL::Compara::DnaFrag"));
00654     $self->{'dnafrag'} = $dnafrag;
00655     if ($self->{'dnafrag_id'}) {
00656       if (!$self->{'dnafrag'}->dbID) {
00657         $self->{'dnafrag'}->dbID($self->{'dnafrag_id'});
00658       }
00659 #       warning("Defining both dnafrag_id and dnafrag");
00660       throw("dnafrag object does not match previously defined dnafrag_id")
00661           if ($self->{'dnafrag'}->dbID != $self->{'dnafrag_id'});
00662     } else {
00663       $self->{'dnafrag_id'} = $self->{'dnafrag'}->dbID;
00664     }
00665   
00666   } elsif (!defined($self->{'dnafrag'})) {
00667 
00668     # Try to get data from other sources...
00669     if (defined($self->dnafrag_id) and defined($self->{'adaptor'})) {
00670       # ...from the dnafrag_id. Use dnafrag_id function and not the attribute in the <if>
00671       # clause because the attribute can be retrieved from other sources if it has not been already defined.
00672       my $dnafrag_adaptor = $self->adaptor->db->get_DnaFragAdaptor;
00673       $self->{'dnafrag'} = $dnafrag_adaptor->fetch_by_dbID($self->{'dnafrag_id'});
00674     } else {
00675       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->dnafrag".
00676           " You either have to specify more information (see perldoc for".
00677           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00678     }
00679   }
00680   return $self->{'dnafrag'};
00681 }
00682 
00683 
00684 =head2 dnafrag_id
00685 
00686   Arg [1]    : integer $dnafrag_id
00687   Example    : $dnafrag_id = $genomic_align->dnafrag_id;
00688   Example    : $genomic_align->dnafrag_id(134);
00689   Description: Getter/Setter for the attribute dnafrag_id. If no
00690                argument is given and the dnafrag_id is not defined, it tries to
00691                get the ID from other sources like the corresponding
00692                Bio::EnsEMBL::Compara::DnaFrag object or the database using the dbID
00693                of the Bio::EnsEMBL::Compara::GenomicAlign object.
00694                Use 0 as argument to clear this attribute.
00695   Returntype : integer
00696   Exceptions : thrown if $dnafrag_id does not match a previously defined
00697                dnafrag
00698   Warning    : warns if getting data from other sources fails.
00699   Caller     : object->methodname
00700   Status     : Stable
00701 
00702 =cut
00703 
00704 sub dnafrag_id {
00705   my ($self, $dnafrag_id) = @_;
00706 
00707   if (defined($dnafrag_id)) {
00708     $self->{'dnafrag_id'} = $dnafrag_id;
00709     if (defined($self->{'dnafrag'}) and $self->{'dnafrag_id'}) {
00710 #       warning("Defining both dnafrag_id and dnafrag");
00711       throw("dnafrag_id does not match previously defined dnafrag object")
00712           if ($self->{'dnafrag'} and $self->{'dnafrag'}->dbID != $self->{'dnafrag_id'});
00713     }
00714 
00715   } elsif (!($self->{'dnafrag_id'})) {
00716     # Try to get the ID from other sources...
00717     if (defined($self->{'dnafrag'}) and defined($self->{'dnafrag'}->dbID)) {
00718       # ...from the corresponding Bio::EnsEMBL::Compara::DnaFrag object
00719       $self->{'dnafrag_id'} = $self->{'dnafrag'}->dbID;
00720     } elsif (defined($self->{'adaptor'}) and defined($self->{'dbID'})) {
00721       # ...from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object
00722       $self->adaptor->retrieve_all_direct_attributes($self);
00723     } else {
00724       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->dnafrag_id".
00725           " You either have to specify more information (see perldoc for".
00726           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00727     }
00728   }
00729 
00730   return $self->{'dnafrag_id'};
00731 }
00732 
00733 
00734 =head2 dnafrag_start
00735 
00736   Arg [1]    : integer $dnafrag_start
00737   Example    : $dnafrag_start = $genomic_align->dnafrag_start;
00738   Example    : $genomic_align->dnafrag_start(1233354);
00739   Description: Getter/Setter for the attribute dnafrag_start. If no argument is given, the
00740                dnafrag_start is not defined but both the dbID and the adaptor are, it tries
00741                to fetch and set all the direct attributes from the database using the
00742                dbID of the Bio::EnsEMBL::Compara::GenomicAlign object.
00743   Returntype : integer
00744   Exceptions : none
00745   Warning    : warns if getting data from other sources fails.
00746   Caller     : object->methodname
00747   Status     : Stable
00748 
00749 =cut
00750 
00751 sub dnafrag_start {
00752   my ($self, $dnafrag_start) = @_;
00753 
00754   if (defined($dnafrag_start)) {
00755      $self->{'dnafrag_start'} = $dnafrag_start;
00756 
00757    } elsif (!defined($self->{'dnafrag_start'})) {
00758     if (defined($self->{'dbID'}) and defined($self->{'adaptor'})) {
00759       # Try to get the values from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object
00760       $self->adaptor->retrieve_all_direct_attributes($self);
00761     } else {
00762       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->dnafrag_start".
00763           " You either have to specify more information (see perldoc for".
00764           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00765     }
00766   }
00767 
00768   return $self->{'dnafrag_start'};
00769 }
00770 
00771 
00772 =head2 dnafrag_end
00773 
00774   Arg [1]    : integer $dnafrag_end
00775   Example    : $dnafrag_end = $genomic_align->dnafrag_end;
00776   Example    : $genomic_align->dnafrag_end(1235320);
00777   Description: Getter/Setter for the attribute dnafrag_end. If no argument is given, the
00778                dnafrag_end is not defined but both the dbID and the adaptor are, it tries
00779                to fetch and set all the direct attributes from the database using the
00780                dbID of the Bio::EnsEMBL::Compara::GenomicAlign object.
00781   Returntype : integer
00782   Exceptions : none
00783   Warning    : warns if getting data from other sources fails.
00784   Caller     : object->methodname
00785   Status     : Stable
00786 
00787 =cut
00788 
00789 sub dnafrag_end {
00790   my ($self, $dnafrag_end) = @_;
00791 
00792   if (defined($dnafrag_end)) {
00793      $self->{'dnafrag_end'} = $dnafrag_end;
00794 
00795   } elsif (!defined($self->{'dnafrag_end'})) {
00796     if (defined($self->{'dbID'}) and defined($self->{'adaptor'})) {
00797       # Try to get the values from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object
00798       $self->adaptor->retrieve_all_direct_attributes($self);
00799     } else {
00800       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->dnafrag_end".
00801           " You either have to specify more information (see perldoc for".
00802           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00803     }
00804   }
00805 
00806   return $self->{'dnafrag_end'};
00807 }
00808 
00809 
00810 =head2 dnafrag_strand
00811 
00812   Arg [1]    : integer $dnafrag_strand (1 or -1)
00813   Example    : $dnafrag_strand = $genomic_align->dnafrag_strand;
00814   Example    : $genomic_align->dnafrag_strand(1);
00815   Description: Getter/Setter for the attribute dnafrag_strand. If no argument is given, the
00816                dnafrag_strand is not defined but both the dbID and the adaptor are, it tries
00817                to fetch and set all the direct attributes from the database using the
00818                dbID of the Bio::EnsEMBL::Compara::GenomicAlign object.
00819   Returntype : integer
00820   Exceptions : none
00821   Warning    : warns if getting data from other sources fails.
00822   Caller     : object->methodname
00823   Status     : Stable
00824 
00825 =cut
00826 
00827 sub dnafrag_strand {
00828   my ($self, $dnafrag_strand) = @_;
00829 
00830   if (defined($dnafrag_strand)) {
00831      $self->{'dnafrag_strand'} = $dnafrag_strand;
00832 
00833   } elsif (!defined($self->{'dnafrag_strand'})) {
00834     if (defined($self->{'dbID'}) and defined($self->{'adaptor'})) {
00835       # Try to get the values from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object
00836       $self->adaptor->retrieve_all_direct_attributes($self);
00837     } else {
00838       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->dnafrag_strand".
00839           " You either have to specify more information (see perldoc for".
00840           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00841     }
00842   }
00843 
00844   return $self->{'dnafrag_strand'};
00845 }
00846 
00847 
00848 =head2 aligned_sequence
00849 
00850   Arg [1...] : string $aligned_sequence or string @flags
00851   Example    : $aligned_sequence = $genomic_align->aligned_sequence
00852   Example    : $aligned_sequence = $genomic_align->aligned_sequence("+FIX_SEQ");
00853   Example    : $genomic_align->aligned_sequence("ACTAGTTAGCT---TATCT--TTAAA")
00854   Description: With no arguments, rebuilds the alignment string for this sequence
00855                using the cigar_line information and the original sequence if needed.
00856                This sequence depends on the strand defined by the dnafrag_strand attribute.
00857   Flags      : +FIX_SEQ
00858                    With this flag, the method will return a sequence that could be
00859                    directly aligned with the original_sequence of the reference
00860                    genomic_align.
00861   Returntype : string $aligned_sequence
00862   Exceptions : thrown if sequence contains unknown symbols
00863   Warning    : warns if getting data from other sources fails.
00864   Caller     : object->methodname
00865   Status     : Stable
00866 
00867 =cut
00868 
00869 sub aligned_sequence {
00870   my ($self, @aligned_sequence_or_flags) = @_;
00871   my $aligned_sequence;
00872 
00873   my $fix_seq = 0;
00874   my $fake_seq = 0;
00875   foreach my $flag (@aligned_sequence_or_flags) {
00876     if ($flag =~ /^\+/) {
00877       if ($flag eq "+FIX_SEQ") {
00878         $fix_seq = 1;
00879       } elsif ($flag eq "+FAKE_SEQ") {
00880         $fake_seq = 1;
00881       } else {
00882         warning("Unknow flag $flag when calling".
00883             " Bio::EnsEMBL::Compara::GenomicAlign::aligned_sequence()");
00884       }
00885     } else {
00886       $aligned_sequence = $flag;
00887     }
00888   }
00889 
00890   if (defined($aligned_sequence)) {
00891     $aligned_sequence =~ s/[\r\n]+$//;
00892     
00893     if ($aligned_sequence) {
00894       ## Check sequence
00895       throw("Unreadable sequence ($aligned_sequence)") if ($aligned_sequence !~ /^[\-\.A-Z]+$/i);
00896       $self->{'aligned_sequence'} = $aligned_sequence;
00897     } else {
00898       $self->{'aligned_sequence'} = undef;
00899     }
00900   } elsif (!defined($self->{'aligned_sequence'})) {
00901     # Try to get the aligned_sequence from other sources...
00902     if (defined($self->cigar_line) and $fake_seq) {
00903       # ...from the corresponding cigar_line (using a fake seq)
00904       $aligned_sequence = _get_fake_aligned_sequence_from_cigar_line(
00905           $self->{'cigar_line'});
00906     
00907     } elsif (defined($self->cigar_line) and defined($self->original_sequence)) {
00908       my $original_sequence = $self->original_sequence;
00909       # ...from the corresponding orginial_sequence and cigar_line
00910       $aligned_sequence = _get_aligned_sequence_from_original_sequence_and_cigar_line(
00911           $original_sequence, $self->{'cigar_line'});
00912       $self->{'aligned_sequence'} = $aligned_sequence;
00913 
00914     } else {
00915       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->aligned_sequence".
00916           " You either have to specify more information (see perldoc for".
00917           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
00918     }
00919   }
00920 
00921   $aligned_sequence = $self->{'aligned_sequence'} if (defined($self->{'aligned_sequence'}));
00922   if ($aligned_sequence and $fix_seq) {
00923     $aligned_sequence = _get_aligned_sequence_from_original_sequence_and_cigar_line(
00924         $aligned_sequence, $self->genomic_align_block->reference_genomic_align->cigar_line, $fix_seq);
00925   } 
00926 
00927   return $aligned_sequence;
00928 }
00929 
00930 
00931 =head2 length
00932 
00933   Arg [1]    : -none-
00934   Example    : $length = $genomic_align->length;
00935   Description: get the length of the aligned sequence. This method will try to
00936                get the length from the aligned_sequence if already set or by
00937                parsing the cigar_line otherwise
00938   Returntype : int
00939   Exceptions : none
00940   Warning    : 
00941   Caller     : object->methodname
00942   Status     : Stable
00943 
00944 =cut
00945 
00946 sub length {
00947   my $self = shift;
00948 
00949   if ($self->{aligned_sequence}) {
00950     return length($self->{aligned_sequence});
00951   } elsif ($self->{cigar_line}) {
00952     my $length = 0;
00953     my $cigar_line = $self->{cigar_line};
00954     my @cig = ( $cigar_line =~ /(\d*[GMDXI])/g );
00955     for my $cigElem ( @cig ) {
00956       my $cigType = substr( $cigElem, -1, 1);
00957       my $cigCount = substr( $cigElem, 0 , -1);
00958       $cigCount = 1 if ($cigCount eq "");
00959       $length += $cigCount unless ($cigType eq "I");
00960     }
00961     return $length;
00962   }
00963 
00964   return undef;
00965 }
00966 
00967 =head2 cigar_line
00968 
00969   Arg [1]    : string $cigar_line
00970   Example    : $cigar_line = $genomic_align->cigar_line;
00971   Example    : $genomic_align->cigar_line("35M2D233M7D23MD100M");
00972   Description: get/set for attribute cigar_line.
00973                If no argument is given, the cigar line has not been
00974                defined yet but the aligned sequence was, it calculates
00975                the cigar line based on the aligned (gapped) sequence.
00976                If no argument is given, the cigar_line is not defined but both
00977                the dbID and the adaptor are, it tries to fetch and set all
00978                the direct attributes from the database using the dbID of the
00979                Bio::EnsEMBL::Compara::GenomicAlign object. You can reset this
00980                attribute using an empty string as argument.
00981                The cigar_line depends on the strand defined by the dnafrag_strand
00982                attribute.
00983   Returntype : string
00984   Exceptions : none
00985   Warning    : warns if getting data from other sources fails.
00986   Caller     : object->methodname
00987   Status     : Stable
00988 
00989 =cut
00990 
00991 sub cigar_line {
00992   my ($self, $arg) = @_;
00993 
00994   if (defined($arg)) {
00995     if ($arg) {
00996       $self->{'cigar_line'} = $arg;
00997     } else {
00998       $self->{'cigar_line'} = undef;
00999     }
01000 
01001   } elsif (!defined($self->{'cigar_line'})) {
01002     # Try to get the cigar_line from other sources...
01003     if (defined($self->{'aligned_sequence'})) {
01004       # ...from the aligned sequence
01005 
01006 
01007       my $cigar_line = _get_cigar_line_from_aligned_sequence($self->{'aligned_sequence'});
01008       $self->cigar_line($cigar_line);
01009     
01010     } elsif (defined($self->{'dbID'}) and defined($self->{'adaptor'})) {
01011       # ...from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object
01012       $self->adaptor->retrieve_all_direct_attributes($self);
01013     } else {
01014       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->cigar_line".
01015           " You either have to specify more information (see perldoc for".
01016           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
01017     }
01018   }
01019 
01020   return $self->{'cigar_line'};
01021 }
01022 
01023 
01024 =head2 visible
01025 
01026   Arg [1]    : int $visible
01027   Example    : $visible = $genomic_align->visible
01028   Example    : $genomic_align->visible(1);
01029   Description: get/set for attribute visible. If no argument is given, visible
01030                is not defined but both the dbID and the adaptor are, it tries to
01031                fetch and set all the direct attributes from the database using the
01032                dbID of the Bio::EnsEMBL::Compara::GenomicAlign object.
01033   Returntype : int
01034   Exceptions : none
01035   Warning    : warns if getting data from other sources fails.
01036   Caller     : object->methodname
01037   Status     : Stable
01038 
01039 =cut
01040 
01041 sub visible {
01042   my ($self, $visible) = @_;
01043 
01044   if (defined($visible)) {
01045     $self->{'visible'} = $visible;
01046 
01047   } elsif (!defined($self->{'visible'})) {
01048     if (defined($self->{'dbID'}) and defined($self->{'adaptor'})) {
01049       # Try to get the values from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object
01050       $self->adaptor->retrieve_all_direct_attributes($self);
01051     } else {
01052       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->visible".
01053           " You either have to specify more information (see perldoc for".
01054           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
01055     }
01056   }
01057 
01058   return $self->{'visible'};
01059 }
01060 
01061 =head2 node_id
01062 
01063   Arg [1]    : [optional] int $node_id
01064   Example    : $node_id = $genomic_align->node_id;
01065   Example    : $genomic_align->node_id(5530000000004);
01066   Description: get/set for the node_id.This links the Bio::EnsEMBL::Compara::GenomicAlign to the 
01067                Bio::EnsEMBL::Compara::GenomicAlignTree. The default value is NULL. If no argument is given, the node_id
01068                is not defined but both the dbID and the adaptor are, it tries to
01069                fetch and set all the direct attributes from the database using the
01070                dbID of the Bio::EnsEMBL::Compara::GenomicAlign object.
01071   Returntype : int
01072   Exceptions : none
01073   Warning    : warns if getting data from other sources fails.
01074   Caller     : object->methodname
01075   Status     : At risk 
01076 
01077 =cut
01078 
01079 sub node_id {
01080   my ($self, $node_id) = @_;
01081 
01082   if (defined($node_id)) {
01083     $self->{'node_id'} = $node_id;
01084   } elsif (!defined($self->{'node_id'})) {
01085     if (defined($self->{'dbID'}) and defined($self->{'adaptor'})) {
01086       # Try to get the values from the database using the dbID of the Bio::EnsEMBL::Compara::GenomicAlign object
01087       $self->adaptor->retrieve_all_direct_attributes($self);
01088     }
01089   }
01090   return $self->{'node_id'};
01091 }
01092 
01093 =head2 genomic_align_group
01094 
01095   Arg [2]    : [optional] Bio::EnsEMBL::Compara::GenomicAlignGroup $genomic_align_group
01096   Example    : $genomic_align_group = $genomic_align->genomic_align_group();
01097   Example    : $genomic_align->genomic_align_group($genomic_align_group);
01098   Description: get/set for the Bio::EnsEMBL::Compara::GenomicAlginGroup object
01099                corresponding to this Bio::EnsEMBL::Compara::GenomicAlign object 
01100   Returntype : int
01101   Exceptions : none
01102   Warning    : warns if getting data from other sources fails.
01103   Caller     : object->methodname
01104   Status     : Stable
01105 
01106 =cut
01107 
01108 sub genomic_align_group {
01109   my ($self, $genomic_align_group) = @_;
01110 
01111   deprecate("Removed genomic_align_group table");
01112 
01113 }
01114 
01115 
01116 =head2 genomic_align_group_id
01117 
01118   Arg [2]    : [optional] int $genomic_align_group_id
01119   Example    : $genomic_align_group_id = $genomic_align->genomic_align_group_id();
01120   Example    : $genomic_align->genomic_align_group_id(18);
01121   Description: get/set for the genomic_align_group_id corresponding to this
01122                Bio::EnsEMBL::Compara::GenomicAlign object
01123   Returntype : int
01124   Exceptions : none
01125   Warning    : warns if getting data from other sources fails.
01126   Caller     : object->methodname
01127   Status     : Stable
01128 
01129 =cut
01130 
01131 sub genomic_align_group_id {
01132   my ($self, $genomic_align_group_id) = @_;
01133 
01134   deprecate("Removed genomic_align_group table");
01135 }
01136 
01137 
01138 =head2 original_sequence
01139 
01140   Arg [1]    : none
01141   Example    : $original_sequence = $genomic_align->original_sequence
01142   Description: get/set original sequence. If no argument is given and the original_sequence
01143                is not defined, it tries to fetch the data from other sources like the
01144                aligned sequence or the the Bio::EnsEMBL::Compara:DnaFrag object. You can
01145                reset this attribute using an empty string as argument.
01146                This sequence depends on the strand defined by the dnafrag_strand attribute.
01147   Returntype : string $original_sequence
01148   Exceptions : 
01149   Caller     : object->methodname
01150   Status     : Stable
01151 
01152 =cut
01153 
01154 sub original_sequence {
01155   my ($self, $original_sequence) = @_;
01156 
01157   if (defined($original_sequence)) {
01158     if ($original_sequence) {
01159       $self->{'original_sequence'} = $original_sequence;
01160     } else {
01161       $self->{'original_sequence'} = undef;
01162     }
01163 
01164   } elsif (!defined($self->{'original_sequence'})) {
01165     # Try to get the data from other sources...
01166     if ($self->{'aligned_sequence'} and $self->{'cigar_line'} !~ /I/) {
01167       # ...from the aligned sequence
01168       $self->{'original_sequence'} = $self->{'aligned_sequence'};
01169       $self->{'original_sequence'} =~ s/\-//g;
01170 
01171     } elsif (!defined($self->{'original_sequence'}) and defined($self->dnafrag)
01172           and defined($self->dnafrag_start) and defined($self->dnafrag_end)
01173           and defined($self->dnafrag_strand) and defined($self->dnafrag->slice)) {
01174       # ...from the dnafrag object. Uses dnafrag, dnafrag_start and dnafrag_methods instead of the attibutes
01175       # in the <if> clause because the attributes can be retrieved from other sources if they have not been
01176       # already defined.
01177       $self->{'original_sequence'} = $self->dnafrag->slice->subseq(
01178               $self->dnafrag_start,
01179               $self->dnafrag_end,
01180               $self->dnafrag_strand
01181           );
01182     } else {
01183       warning("Fail to get data from other sources in Bio::EnsEMBL::Compara::GenomicAlign->genomic_align_groups".
01184           " You either have to specify more information (see perldoc for".
01185           " Bio::EnsEMBL::Compara::GenomicAlign) or to set it up directly");
01186     }
01187   }
01188 
01189   return $self->{'original_sequence'};
01190 }
01191 
01192 =head2 _get_cigar_line_from_aligned_sequence
01193 
01194   Arg [1]    : string $aligned_sequence
01195   Example    : $cigar_line = _get_cigar_line_from_aligned_sequence("CGT-AACTGATG--TTA")
01196   Description: get cigar line from gapped sequence
01197   Returntype : string $cigar_line
01198   Exceptions : 
01199   Caller     : methodname
01200   Status     : Stable
01201 
01202 =cut
01203 
01204 sub _get_cigar_line_from_aligned_sequence {
01205   my ($aligned_sequence) = @_;
01206   my $cigar_line = "";
01207   
01208   my @pieces = grep {$_} split(/(\-+)|(\.+)/, $aligned_sequence);
01209   foreach my $piece (@pieces) {
01210     my $mode;
01211     if ($piece =~ /\-/) {
01212       $mode = "D"; # D for gaps (deletions)
01213     } elsif ($piece =~ /\./) {
01214       $mode = "X"; # X for pads (in 2X genomes)
01215     } else {
01216       $mode = "M"; # M for matches/mismatches
01217     }
01218     if (CORE::length($piece) == 1) {
01219       $cigar_line .= $mode;
01220     } elsif (CORE::length($piece) > 1) { #length can be 0 if the sequence starts with a gap
01221       $cigar_line .= CORE::length($piece).$mode;
01222     }
01223   }
01224 
01225   return $cigar_line;
01226 }
01227 
01228 
01229 =head2 _get_aligned_sequence_from_original_sequence_and_cigar_line
01230 
01231   Arg [1]    : string $original_sequence
01232   Arg [1]    : string $cigar_line
01233   Example    : $aligned_sequence = _get_aligned_sequence_from_original_sequence_and_cigar_line(
01234                    "CGTAACTGATGTTA", "3MD8M2D3M")
01235   Description: get gapped sequence from original one and cigar line
01236   Returntype : string $aligned_sequence
01237   Exceptions : thrown if cigar_line does not match sequence length
01238   Caller     : methodname
01239   Status     : Stable
01240 
01241 =cut
01242 
01243 sub _get_aligned_sequence_from_original_sequence_and_cigar_line {
01244   my ($original_sequence, $cigar_line, $fix_seq) = @_;
01245   my $aligned_sequence = "";
01246 
01247   return undef if (!defined($original_sequence) or !$cigar_line);
01248 
01249   my $seq_pos = 0;
01250   my @cig = ( $cigar_line =~ /(\d*[GMDXI])/g );
01251 
01252   for my $cigElem ( @cig ) {
01253     my $cigType = substr( $cigElem, -1, 1 );
01254     my $cigCount = substr( $cigElem, 0 ,-1 );
01255     $cigCount = 1 unless ($cigCount =~ /^\d+$/);
01256 
01257     if( $cigType eq "M" ) {
01258       $aligned_sequence .= substr($original_sequence, $seq_pos, $cigCount);
01259       $seq_pos += $cigCount;
01260     } elsif( $cigType eq "I") {
01261       $seq_pos += $cigCount;
01262     } elsif( $cigType eq "X") {
01263       $aligned_sequence .=  "." x $cigCount;
01264     } elsif( $cigType eq "G" || $cigType eq "D") {
01265       if ($fix_seq) {
01266         $seq_pos += $cigCount;
01267       } else {
01268         $aligned_sequence .=  "-" x $cigCount;
01269       }
01270     }
01271   }
01272   throw("Cigar line ($seq_pos) does not match sequence length (".CORE::length($original_sequence).")")
01273       if ($seq_pos != CORE::length($original_sequence));
01274 
01275   return $aligned_sequence;
01276 }
01277 
01278 
01279 =head2 _get_fake_aligned_sequence_from_cigar_line
01280 
01281   Arg [1]    : string $cigar_line
01282   Example    : $aligned_sequence = _get_fake_aligned_sequence_from_cigar_line(
01283                    "3MD8M2D3M")
01284   Description: get gapped sequence of N\'s from the cigar line
01285   Returntype : string $fake_aligned_sequence or undef if no $cigar_line
01286   Exceptions : 
01287   Caller     : methodname
01288   Status     : Stable
01289 
01290 =cut
01291 
01292 sub _get_fake_aligned_sequence_from_cigar_line {
01293   my ($cigar_line, $fix_seq) = @_;
01294   my $fake_aligned_sequence = "";
01295 
01296   return undef if (!$cigar_line);
01297 
01298   my $seq_pos = 0;
01299 
01300   my @cig = ( $cigar_line =~ /(\d*[GMDXI])/g );
01301   for my $cigElem ( @cig ) {
01302     my $cigType = substr( $cigElem, -1, 1 );
01303     my $cigCount = substr( $cigElem, 0 ,-1 );
01304     $cigCount = 1 if ($cigCount eq "");
01305 
01306     if( $cigType eq "M" ) {
01307       $fake_aligned_sequence .= "N" x $cigCount;
01308       $seq_pos += $cigCount;
01309     } elsif( $cigType eq "I") {
01310       $seq_pos += $cigCount;
01311     } elsif( $cigType eq "X") {
01312       $fake_aligned_sequence .=  "." x $cigCount;
01313     } elsif( $cigType eq "G" || $cigType eq "D") {
01314       if ($fix_seq) {
01315         $seq_pos += $cigCount;
01316       } else {
01317         $fake_aligned_sequence .=  "-" x $cigCount;
01318       }
01319     }
01320   }
01321 
01322   return $fake_aligned_sequence;
01323 }
01324 
01325 
01326 =head2 _print
01327 
01328   Arg [1]    : ref to a FILEHANDLE
01329   Example    : $genomic_align->_print
01330   Description: print attributes of the object to the STDOUT or to the FILEHANDLE.
01331                Used for debuging purposes.
01332   Returntype : none
01333   Exceptions : 
01334   Caller     : object->methodname
01335   Status     : At risk
01336 
01337 =cut
01338 
01339 sub _print {
01340   my ($self, $FILEH) = @_;
01341 
01342   my $verbose = verbose;
01343   verbose(0);
01344   
01345   $FILEH ||= \*STDOUT;
01346 
01347   print $FILEH
01348 "Bio::EnsEMBL::Compara::GenomicAlign object ($self)
01349   dbID = ".($self->dbID or "-undef-")."
01350   adaptor = ".($self->adaptor or "-undef-")."
01351   genomic_align_block = ".($self->genomic_align_block or "-undef-")."
01352   genomic_align_block_id = ".($self->genomic_align_block_id or "-undef-")."
01353   method_link_species_set = ".($self->method_link_species_set or "-undef-")."
01354   method_link_species_set_id = ".($self->method_link_species_set_id or "-undef-")."
01355   dnafrag = ".($self->dnafrag or "-undef-")."
01356   dnafrag_id = ".($self->dnafrag_id or "-undef-")."
01357   dnafrag_start = ".($self->dnafrag_start or "-undef-")."
01358   dnafrag_end = ".($self->dnafrag_end or "-undef-")."
01359   dnafrag_strand = ".($self->dnafrag_strand or "-undef-")."
01360   cigar_line = ".($self->cigar_line or "-undef-")."
01361   visible = ".($self->visible or "-undef-")."
01362   original_sequence = ".($self->original_sequence or "-undef-")."
01363   aligned_sequence = ".($self->aligned_sequence or "-undef-")."
01364   
01365 ";
01366   verbose($verbose);
01367 
01368 }
01369 
01370 =head2 display_id
01371 
01372   Args       : none
01373   Example    : my $id = $genomic_align->display_id;
01374   Description: returns string describing this genomic_align which can be used
01375                as display_id of a Bio::Seq object or in a fasta file. The actual form is
01376                taxon_id:genome_db_id:coord_system_name:dnafrag_name:dnafrag_start:dnafrag_end:dnafrag_strand
01377                e.g.
01378                9606:1:chromosome:14:50000000:51000000:-1
01379 
01380                Uses dnafrag information in addition to start and end.
01381   Returntype : string
01382   Exceptions : none
01383   Caller     : general
01384   Status     : Stable
01385 
01386 =cut
01387 
01388 sub display_id {
01389   my $self = shift;
01390 
01391   my $dnafrag = $self->dnafrag;
01392   return "" unless($dnafrag);
01393   my $id = join(':',
01394                 $dnafrag->genome_db->taxon_id,
01395                 $dnafrag->genome_db->dbID,
01396                 $dnafrag->coord_system_name,
01397                 $dnafrag->name,
01398                 $self->dnafrag_start,
01399                 $self->dnafrag_end,
01400                 $self->dnafrag_strand);
01401   return $id;
01402 }
01403 
01404 =head2 reverse_complement
01405 
01406   Args       : none
01407   Example    : none
01408   Description: reverse complement the object modifing dnafrag_strand and cigar_line
01409   Returntype : none
01410   Exceptions : none
01411   Caller     : general
01412   Status     : Stable
01413 
01414 =cut
01415 
01416 sub reverse_complement {
01417   my ($self) = @_;
01418 
01419   # reverse strand
01420   #$self->dnafrag_strand($self->dnafrag_strand * -1);
01421   $self->dnafrag_strand($self->{'dnafrag_strand'} * -1);
01422 
01423   # reverse original and aligned sequences if cached
01424   my $original_sequence = $self->{'original_sequence'};
01425   if ($original_sequence) {
01426     $original_sequence = reverse $original_sequence;
01427     $original_sequence =~ tr/ATCGatcg/TAGCtagc/;
01428     $self->original_sequence($original_sequence);
01429   }
01430   my $aligned_sequence = $self->{'aligned_sequence'};
01431   if ($aligned_sequence) {
01432     $aligned_sequence = reverse $aligned_sequence;
01433     $aligned_sequence =~ tr/ATCGatcg/TAGCtagc/;
01434     $self->aligned_sequence($aligned_sequence);
01435   }
01436   
01437   # reverse cigar_string as consequence
01438   my $cigar_line = $self->{'cigar_line'};
01439   
01440   #$cigar_line = join("", reverse grep {$_} split(/(\d*[GDMIX])/, $cigar_line));
01441   $cigar_line = join("", reverse ($cigar_line=~(/(\d*[GDMIX])/g)));
01442   $self->cigar_line($cigar_line);
01443 }
01444 
01445 
01446 =head2 get_Mapper
01447 
01448   Arg[1]     : [optional] integer $cache (default = FALSE)
01449   Arg[2]     : [optional] boolean $condensed (default = FALSE)
01450   Example    : $this_mapper = $genomic_align->get_Mapper();
01451   Example    : $mapper1 = $genomic_align1->get_Mapper();
01452                $mapper2 = $genomic_align2->get_Mapper();
01453   Description: creates and returns a Bio::EnsEMBL::Mapper to map coordinates from
01454                the original sequence of this Bio::EnsEMBL::Compara::GenomicAlign
01455                to the aligned sequence, i.e. the alignment. In order to map a sequence
01456                from this Bio::EnsEMBL::Compara::GenomicAlign object to another
01457                Bio::EnsEMBL::Compara::GenomicAlign of the same
01458                Bio::EnsEMBL::Compara::GenomicAlignBlock object, you may use this mapper
01459                to transform coordinates into the "alignment" coordinates and then to
01460                the other Bio::EnsEMBL::Compara::GenomicAlign coordinates using the
01461                corresponding Bio::EnsEMBL::Mapper.
01462                The coordinates of the "alignment" starts with the start
01463                position of the GenomicAlignBlock if available or 1 otherwise.
01464                With the $cache argument you can decide whether you want to cache the
01465                result or not. Result is *not* cached by default.
01466   Returntype : Bio::EnsEMBL::Mapper object
01467   Exceptions : throw if no cigar_line can be found
01468   Status     : Stable
01469 
01470 =cut
01471 
01472 sub get_Mapper {
01473   my ($self, $cache, $condensed) = @_;
01474   my $mapper;
01475   $cache = 0 if (!defined($cache));
01476   my $mode = "expanded";
01477   if (defined($condensed) and $condensed) {
01478     $mode = "condensed";
01479   }
01480 
01481   if (!defined($self->{$mode.'_mapper'})) {
01482     if ($mode eq "condensed") {
01483 
01484       $mapper = Bio::EnsEMBL::Mapper->new("sequence", "alignment");
01485 
01486       my $rel_strand = $self->dnafrag_strand; # This call ensures all direct attribs have been fetched
01487       my $ref_cigar_line = $self->genomic_align_block->reference_genomic_align->cigar_line;
01488 
01489       my $aln_pos = (eval{$self->genomic_align_block->start} or 1);
01490 
01491       #if the reference genomic_align, I only need a simple 1 to 1 mapping
01492       if ($self eq $self->genomic_align_block->reference_genomic_align) {
01493       $mapper->add_map_coordinates(
01494               'sequence',
01495               $self->{'dnafrag_start'},
01496               $self->{'dnafrag_end'},
01497               $self->{'dnafrag_strand'},
01498               'alignment',
01499           $self->genomic_align_block->start,
01500           $self->genomic_align_block->end,
01501           );
01502       return $mapper if (!$cache);
01503 
01504       $self->{$mode.'_mapper'} = $mapper;
01505       return $self->{$mode.'_mapper'};
01506       }
01507 
01508       my $aln_seq_pos = 0;
01509       my $seq_pos = 0;
01510 
01511       my $insertions = 0;
01512       my $target_cigar_pieces;
01513       @$target_cigar_pieces = $self->{'cigar_line'} =~ /(\d*[GMDXI])/g;
01514       my $ref_cigar_pieces;
01515       @$ref_cigar_pieces = $ref_cigar_line =~ /(\d*[GMDXI])/g;
01516       my $i = 0;
01517       my $j = 0;
01518       my ($ref_num, $ref_type) = $ref_cigar_pieces->[$i] =~ /(\d*)([GMDXI])/;
01519       $ref_num = 1 if (!defined($ref_num) or $ref_num eq "");
01520       my ($target_num, $target_type) = $target_cigar_pieces->[$j] =~ /(\d*)([GMDXI])/;
01521       $target_num = 1 if (!defined($target_num) or $target_num eq "");
01522 
01523       while ($i < @$ref_cigar_pieces and $j<@$target_cigar_pieces) {
01524       while ($ref_type eq "I") {
01525           $aln_pos += $ref_num;
01526           $i++;
01527           last if ($i >= @$ref_cigar_pieces);
01528           ($ref_num, $ref_type) = $ref_cigar_pieces->[$i] =~ /(\d*)([GMDXI])/;
01529           $ref_num = 1 if (!defined($ref_num) or $ref_num eq "");
01530       }
01531       while ($target_type eq "I") {
01532           $seq_pos += $target_num;
01533           $j++;
01534           last if ($j >= @$target_cigar_pieces);
01535           ($target_num, $target_type) = $target_cigar_pieces->[$j] =~ /(\d*)([GMDXI])/;
01536           $target_num = 1 if (!defined($target_num) or $target_num eq "");
01537       }
01538 
01539         my $length;
01540 
01541     if ($ref_num == $target_num) {
01542       $length = $ref_num;
01543     } elsif ($ref_num > $target_num) {
01544       $length = $target_num;
01545     } elsif ($ref_num < $target_num) {
01546       $length = $ref_num;
01547         }
01548     my $this_piece_of_cigar_line = $length.$target_type;
01549 
01550     if ($ref_type eq "M") {
01551           my $this_mapper;
01552           if ($rel_strand == 1) {
01553             _add_cigar_line_to_Mapper($this_piece_of_cigar_line, $aln_pos,
01554                 $seq_pos + $self->dnafrag_start, 1, $mapper);
01555           } else {
01556             _add_cigar_line_to_Mapper($this_piece_of_cigar_line, $aln_pos, $self->dnafrag_end - $seq_pos, -1, $mapper);
01557           }
01558       $aln_pos += $length;
01559         }
01560     my $gaps = 0;
01561     if ($target_type eq "D" || $target_type eq "X") {
01562         $gaps += $length;
01563     }
01564 
01565         $seq_pos -= $gaps;
01566     $seq_pos += $length;
01567 
01568     if ($ref_num == $target_num) {
01569       $i++;
01570       $j++;
01571       last if ($i >= @$ref_cigar_pieces);
01572       last if ($j >= @$target_cigar_pieces);
01573       ($ref_num, $ref_type) = $ref_cigar_pieces->[$i] =~ /(\d*)([GMDXI])/;
01574       $ref_num = 1 if (!defined($ref_num) or $ref_num eq "");
01575       ($target_num, $target_type) = $target_cigar_pieces->[$j] =~ /(\d*)([GMDXI])/;
01576       $target_num = 1 if (!defined($target_num) or $target_num eq "");
01577     } elsif ($ref_num > $target_num) {
01578       $j++;
01579       $ref_num -= $target_num;
01580       last if ($j >= @$target_cigar_pieces);
01581       ($target_num, $target_type) = $target_cigar_pieces->[$j] =~ /(\d*)([GMDXI])/;
01582       $target_num = 1 if (!defined($target_num) or $target_num eq "");
01583     } elsif ($ref_num < $target_num) {
01584       $i++;
01585       $target_num -= $ref_num;
01586       last if ($i >= @$ref_cigar_pieces);
01587       ($ref_num, $ref_type) = $ref_cigar_pieces->[$i] =~ /(\d*)([GMDXI])/;
01588       $ref_num = 1 if (!defined($ref_num) or $ref_num eq "");
01589         }
01590       }
01591     } else {
01592       my $cigar_line = $self->cigar_line;
01593       if (!$cigar_line) {
01594         throw("[$self] has no cigar_line and cannot be retrieved by any means");
01595       }
01596       my $alignment_position = (eval{$self->genomic_align_block->start} or 1);
01597       my $sequence_position = $self->dnafrag_start;
01598       my $rel_strand = $self->dnafrag_strand;
01599       if ($rel_strand == 1) {
01600         $sequence_position = $self->dnafrag_start;
01601       } else {
01602         $sequence_position = $self->dnafrag_end;
01603       }
01604       $mapper = _get_Mapper_from_cigar_line($cigar_line, $alignment_position, $sequence_position, $rel_strand);
01605     }
01606 
01607     return $mapper if (!$cache);
01608 
01609     $self->{$mode.'_mapper'} = $mapper;
01610   }
01611 
01612   return $self->{$mode.'_mapper'};
01613 }
01614 
01615 sub get_MapperOLD {
01616   my ($self, $cache, $condensed) = @_;
01617   my $mapper;
01618   $cache = 0 if (!defined($cache));
01619   my $mode = "expanded";
01620   if (defined($condensed) and $condensed) {
01621     $mode = "condensed";
01622   }
01623 
01624   if (!defined($self->{$mode.'_mapper'})) {
01625     if ($mode eq "condensed") {
01626       $mapper = Bio::EnsEMBL::Mapper->new("sequence", "alignment");
01627       my $rel_strand = $self->dnafrag_strand;
01628       my $ref_cigar_line = $self->genomic_align_block->reference_genomic_align->cigar_line;
01629       my $this_aligned_seq = $self->aligned_sequence("+FAKE_SEQ");
01630 
01631       my $aln_pos = (eval{$self->genomic_align_block->start} or 1);
01632       my $aln_seq_pos = 0;
01633       my $seq_pos = 0;
01634 
01635       my $target_cigar_pieces;
01636       @$target_cigar_pieces = $self->cigar_line =~ /(\d*[GMDXI])/g;
01637 
01638       my $insertions = 0;
01639       my $array_index = 0;
01640       my $this_target_pos = 0;
01641       foreach my $cigar_piece ($ref_cigar_line =~ /(\d*[GMDX])/g) {
01642         my ($cig_count, $cig_mode) = $cigar_piece =~ /(\d*)([GMDX])/;
01643         $cig_count = 1 if (!defined($cig_count) or $cig_count eq "");
01644 
01645     #because of 2X genomes, need different method for extracting the
01646     #cigar_line
01647     #my $this_piece_of_cigar_line = _get_sub_cigar_line_slow($target_cigar_pieces, $aln_seq_pos, $cig_count);
01648 
01649     #quicker method which keeps track of how far through the 
01650     #target_cigar_pieces array we are
01651     my $this_piece_of_cigar_line;
01652     ($this_piece_of_cigar_line,$array_index, $this_target_pos) = _get_sub_cigar_line($target_cigar_pieces, $aln_seq_pos, $cig_count, $array_index, $this_target_pos);
01653 
01654     #find number of each cigar_line mode in cigar_line
01655     my $num_cigar_elements = _count_cigar_elements($this_piece_of_cigar_line);
01656         if ($cig_mode eq "M") {
01657 
01658           my $this_mapper;
01659           if ($rel_strand == 1) {
01660             $this_mapper = _get_Mapper_from_cigar_line($this_piece_of_cigar_line, $aln_pos,
01661                 $seq_pos + $self->dnafrag_start, 1);
01662           } else {
01663             $this_mapper = _get_Mapper_from_cigar_line($this_piece_of_cigar_line, $aln_pos,
01664                 $self->dnafrag_end - $seq_pos, -1);
01665           }
01666           $mapper->add_Mapper($this_mapper);
01667           $aln_pos += $cig_count;
01668 
01669       $insertions = $num_cigar_elements->{"I"};
01670 
01671       $seq_pos += $insertions;
01672         }
01673     my $gaps = $num_cigar_elements->{"D"};
01674     $gaps += $num_cigar_elements->{"X"};
01675 
01676         $seq_pos -= $gaps;
01677         $seq_pos += $cig_count;
01678         $aln_seq_pos += $cig_count;
01679       }
01680 
01681     } else {
01682       my $cigar_line = $self->cigar_line;
01683       if (!$cigar_line) {
01684         throw("[$self] has no cigar_line and cannot be retrieved by any means");
01685       }
01686       my $alignment_position = (eval{$self->genomic_align_block->start} or 1);
01687       my $sequence_position = $self->dnafrag_start;
01688       my $rel_strand = $self->dnafrag_strand;
01689       if ($rel_strand == 1) {
01690         $sequence_position = $self->dnafrag_start;
01691       } else {
01692         $sequence_position = $self->dnafrag_end;
01693       }
01694       $mapper = _get_Mapper_from_cigar_line($cigar_line, $alignment_position, $sequence_position, $rel_strand);
01695     }
01696 
01697     return $mapper if (!$cache);
01698 
01699     $self->{$mode.'_mapper'} = $mapper;
01700   }
01701 
01702   return $self->{$mode.'_mapper'};
01703 }
01704 
01705 =head2 _count_cigar_elements
01706 
01707   Arg[1]     : string $cigar_line 
01708   Example    : $num_elements = _count_cigar_elements("5M3D2M5D")
01709   Description: Counts the number of each cigar_line mode in a cigar_line
01710                and stores them in a hash reference. In the above example
01711                $num_elements->{"M"} is 7, $num_elements->{"D"} is 8
01712   Returntype : hash reference
01713   Exceptions : None
01714   Status     : At risk
01715 
01716 =cut
01717 sub _count_cigar_elements {
01718     my ($cigar_line) = @_;
01719 
01720     my $this_count = 0;
01721     my $num_elements;
01722 
01723     #initialise each element to 0
01724     foreach my $mode (qw(G M D X I)) {
01725     $num_elements->{$mode} = 0;
01726     }
01727 
01728     foreach my $cigar_piece ($cigar_line =~ /(\d*[GMDXI])/g) {
01729         my ($cig_count, $cig_mode) = $cigar_piece =~ /(\d*)([GMDXI])/;
01730         $cig_count = 1 if (!defined($cig_count) or $cig_count eq "");
01731     $num_elements->{$cig_mode} += $cig_count;
01732     }
01733     return $num_elements;
01734 }
01735 
01736 =head2 _get_sub_cigar_line
01737 
01738   Arg[1]     : ref to array of target cigar_line elements
01739   Arg[2]     : int $offset start position
01740   Arg[3]     : int $length amount to extract
01741   Arg[4]     : int $start_array_index current element in target array
01742   Arg[5]     : int $start_target_pos current position in target coords
01743   Example    : my $new_cigar_line = _get_sub_cigar_line($target_cigar_pieces, $pos, $count);
01744   Description: Extracts a cigar_line of size $length starting at $offset
01745   Returntype : string
01746   Exceptions : None
01747   Status     : At risk
01748 
01749 =cut
01750 sub _get_sub_cigar_line {
01751     my ($target_cigar_pieces, $offset, $length, $start_array_index, $start_target_pos) = @_;
01752     my $ref_pos = $offset + $length;
01753 
01754     my $i = $start_array_index;
01755     my $target_pos = $start_target_pos;
01756 
01757     #current target element
01758     my ($target_cig_count, $target_cig_mode) = $target_cigar_pieces->[$i] =~ /(\d*)([GMDXI])/;
01759     $target_cig_count = 1 if (!defined($target_cig_count) or $target_cig_count eq "");
01760 
01761     my $new_cigar_line = "";
01762     #check to see if previous target overlaps this ref_pos
01763     if ($offset) {
01764     if ($target_pos > $offset) {
01765 
01766         #need to only add on cig_count amount
01767         my $new_count;
01768         if ($target_pos - $offset < $length) {
01769         $new_count = ($target_pos - $offset);
01770         } else {
01771         $new_count = $length;
01772         }
01773         #$new_cigar_line .= $new_count . $target_cig_mode;
01774         $new_cigar_line .= _cigar_element($target_cig_mode,$new_count);
01775         #print "here1 $target_cig_mode $new_count\n";
01776     }
01777     #increment to next target element
01778     $i++;
01779     }
01780     while ($target_pos < $ref_pos && $i < @$target_cigar_pieces) {
01781     ($target_cig_count, $target_cig_mode) = $target_cigar_pieces->[$i++] =~ /(\d*)([GMDXI])/;
01782     $target_cig_count = 1 if (!defined($target_cig_count) or $target_cig_count eq "");
01783 
01784     #first piece
01785     if (!$target_pos) {
01786         if ($target_cig_count >= $length) {
01787         $new_cigar_line .= _cigar_element($target_cig_mode,$length);
01788         } else {
01789         $new_cigar_line .= _cigar_element($target_cig_mode,$target_cig_count);
01790         }
01791     } else {
01792         if ($target_cig_mode ne "I" && 
01793         $target_cig_count + $target_pos > $ref_pos) {
01794         #if new target piece extends beyond ref_piece but is not I 
01795         #(since this doesn't count to target_pos) need to shorten it
01796         my $count = $ref_pos - $target_pos;
01797         $new_cigar_line .= _cigar_element($target_cig_mode,$count);
01798         } else {
01799         $new_cigar_line .= _cigar_element($target_cig_mode,$target_cig_count);
01800         }
01801     }
01802     $target_pos += $target_cig_count unless $target_cig_mode eq "I";
01803     }
01804     #need to check if the next element is an I which doesn't count to 
01805     #target_pos but need to add it to cigar_line
01806     if ($i < @$target_cigar_pieces) {
01807     ($target_cig_count, $target_cig_mode) = $target_cigar_pieces->[$i] =~ /(\d*)([GMDXI])/;
01808     $target_cig_count = 1 if (!defined($target_cig_count) or $target_cig_count eq "");
01809     if ($target_cig_mode eq "I") {
01810         $new_cigar_line .= _cigar_element($target_cig_mode,$target_cig_count);
01811     }
01812     }
01813 
01814     #decrement to return current target element 
01815     if ( $i > 0) {
01816     $i--;
01817     }
01818     return ($new_cigar_line, $i, $target_pos);
01819 }
01820 
01821 sub _get_sub_cigar_line_slow {
01822     my ($target_cigar_pieces, $offset, $length) = @_;
01823     my $i = 0;
01824     my $ref_pos = $offset + $length;
01825     my ($target_cig_count, $target_cig_mode);
01826     my $target_pos = 0;
01827 
01828     #skip through target_cigar_line until get to correct position
01829     while ($target_pos < $offset && $i < @$target_cigar_pieces) {
01830     ($target_cig_count, $target_cig_mode) = $target_cigar_pieces->[$i++] =~ /(\d*)([GMDXI])/;
01831     $target_cig_count = 1 if (!defined($target_cig_count) or $target_cig_count eq "");
01832 
01833     $target_pos += $target_cig_count unless $target_cig_mode eq "I";
01834     }
01835 
01836     my $new_cigar_line = "";
01837     #check to see if previous target overlaps this ref_pos
01838     if ($offset) {
01839     if ($target_pos > $offset) {
01840 
01841         #need to only add on cig_count amount
01842         my $new_count;
01843         if ($target_pos - $offset < $length) {
01844         $new_count = ($target_pos - $offset);
01845         } else {
01846         $new_count = $length;
01847         }
01848         #$new_cigar_line .= $new_count . $target_cig_mode;
01849         $new_cigar_line .= _cigar_element($target_cig_mode,$new_count);
01850     }
01851     }
01852 
01853     while ($target_pos < $ref_pos && $i < @$target_cigar_pieces) {
01854     ($target_cig_count, $target_cig_mode) = $target_cigar_pieces->[$i++] =~ /(\d*)([GMDXI])/;
01855     $target_cig_count = 1 if (!defined($target_cig_count) or $target_cig_count eq "");
01856 
01857     #first piece
01858     if (!$target_pos) {
01859         if ($target_cig_count >= $length) {
01860         $new_cigar_line .= _cigar_element($target_cig_mode,$length);
01861         } else {
01862         $new_cigar_line .= _cigar_element($target_cig_mode,$target_cig_count);
01863         }
01864     } else {
01865         if ($target_cig_mode ne "I" && 
01866         $target_cig_count + $target_pos > $ref_pos) {
01867         #if new target piece extends beyond ref_piece but is not I 
01868         #(since this doesn't count to target_pos) need to shorten it
01869         my $count = $ref_pos - $target_pos;
01870         $new_cigar_line .= _cigar_element($target_cig_mode,$count);
01871         } else {
01872         $new_cigar_line .= _cigar_element($target_cig_mode,$target_cig_count);
01873         }
01874     }
01875     $target_pos += $target_cig_count unless $target_cig_mode eq "I";
01876     }
01877     #need to check if the next element is an I which doesn't count to 
01878     #target_pos but need to add it to cigar_line
01879     if ($i < @$target_cigar_pieces) {
01880     ($target_cig_count, $target_cig_mode) = $target_cigar_pieces->[$i++] =~ /(\d*)([GMDXI])/;
01881     $target_cig_count = 1 if (!defined($target_cig_count) or $target_cig_count eq "");
01882     if ($target_cig_mode eq "I") {
01883         $new_cigar_line .= _cigar_element($target_cig_mode,$target_cig_count);
01884     }
01885     }
01886     return $new_cigar_line;
01887 }
01888 
01889 =head2 _cigar_element
01890 
01891   Arg[1]     : char $mode valid cigar_line mode
01892   Arg[2]     : int $length size of element
01893   Example    : $elem = _cigar_element("M", 5);
01894   Description: Creates a valid cigar element
01895   Returntype : integer
01896   Exceptions : None
01897   Status     : At risk
01898 
01899 =cut
01900 sub _cigar_element {
01901     my ($mode, $len) = @_;
01902     my $elem;
01903     if ($len == 1) {
01904     $elem = $mode;
01905     #} elsif ($len > 1) { #length can be 0 if the sequence starts with a gap
01906     } else { #length can be 0 if the sequence starts with a gap
01907     $elem = $len.$mode;
01908     }
01909     return $elem;
01910 }
01911 
01912 =head2 _get_Mapper_from_cigar_line
01913 
01914   Arg[1]     : $cigar_line
01915   Arg[2]     : $alignment_position
01916   Arg[3]     : $sequence_position
01917   Arg[4]     : $relative_strand
01918   Example    : $this_mapper = _get_Mapper_from_cigar_line($cigar_line, 
01919                 $aln_pos, $seq_pos, 1);
01920   Description: creates a new Bio::EnsEMBL::Mapper object for mapping between
01921                sequence and alignment coordinate systems using the cigar_line
01922                and starting from the $alignment_position and sequence_position.
01923   Returntype : Bio::EnsEMBL::Mapper object
01924   Exceptions : None
01925   Status     : Stable
01926 
01927 =cut
01928 
01929 sub _get_Mapper_from_cigar_line {
01930   my ($cigar_line, $alignment_position, $sequence_position, $rel_strand) = @_;
01931 
01932   my $mapper = Bio::EnsEMBL::Mapper->new("sequence", "alignment");
01933 
01934   my @cigar_pieces = ($cigar_line =~ /(\d*[GMDXI])/g);
01935   if ($rel_strand == 1) {
01936     foreach my $cigar_piece (@cigar_pieces) {
01937       my $cigar_type = substr($cigar_piece, -1, 1);
01938       my $cigar_count = substr($cigar_piece, 0, -1);
01939       $cigar_count = 1 if ($cigar_count eq "");
01940       next if ($cigar_count < 1);
01941   
01942       if( $cigar_type eq "M" ) {
01943          $mapper->add_map_coordinates(
01944                 "sequence", #$self->dbID,
01945                 $sequence_position,
01946                 $sequence_position + $cigar_count - 1,
01947                 $rel_strand,
01948                 "alignment", #$self->genomic_align_block->dbID,
01949                 $alignment_position,
01950                 $alignment_position + $cigar_count - 1
01951             );
01952         $sequence_position += $cigar_count;
01953         $alignment_position += $cigar_count;
01954       } elsif( $cigar_type eq "I") {
01955     #add to sequence_position but not alignment_position
01956     $sequence_position += $cigar_count;
01957       } elsif( $cigar_type eq "G" || $cigar_type eq "D" || $cigar_type eq "X") {
01958         $alignment_position += $cigar_count;
01959       }
01960     }
01961   } else {
01962     foreach my $cigar_piece (@cigar_pieces) {
01963       my $cigar_type = substr($cigar_piece, -1, 1);
01964       my $cigar_count = substr($cigar_piece, 0 ,-1);
01965       $cigar_count = 1 if ($cigar_count eq "");
01966       next if ($cigar_count < 1);
01967   
01968       if( $cigar_type eq "M" ) {
01969         $mapper->add_map_coordinates(
01970                 "sequence", #$self->dbID,
01971                 $sequence_position - $cigar_count + 1,
01972                 $sequence_position,
01973                 $rel_strand,
01974                 "alignment", #$self->genomic_align_block->dbID,
01975                 $alignment_position,
01976                 $alignment_position + $cigar_count - 1
01977             );
01978         $sequence_position -= $cigar_count;
01979         $alignment_position += $cigar_count;
01980       } elsif( $cigar_type eq "I") {
01981     #add to sequence_position but not alignment_position
01982     $sequence_position -= $cigar_count;
01983       } elsif( $cigar_type eq "G" || $cigar_type eq "D" || $cigar_type eq "X") {
01984         $alignment_position += $cigar_count;
01985       }
01986     }
01987   }
01988 
01989   return $mapper;
01990 }
01991 
01992 sub _add_cigar_line_to_Mapper {
01993   my ($cigar_line, $alignment_position, $sequence_position, $rel_strand, $mapper) = @_;
01994 
01995   my @cigar_pieces = ($cigar_line =~ /(\d*[GMDXI])/g);
01996   if ($rel_strand == 1) {
01997     foreach my $cigar_piece (@cigar_pieces) {
01998       my $cigar_type = substr($cigar_piece, -1, 1 );
01999       my $cigar_count = substr($cigar_piece, 0 ,-1 );
02000       $cigar_count = 1 unless ($cigar_count =~ /^\d+$/);
02001       next if ($cigar_count < 1);
02002   
02003       if( $cigar_type eq "M" ) {
02004         $mapper->add_map_coordinates(
02005                 "sequence", #$self->dbID,
02006                 $sequence_position,
02007                 $sequence_position + $cigar_count - 1,
02008                 $rel_strand,
02009                 "alignment", #$self->genomic_align_block->dbID,
02010                 $alignment_position,
02011                 $alignment_position + $cigar_count - 1
02012             );
02013         $sequence_position += $cigar_count;
02014         $alignment_position += $cigar_count;
02015       } elsif( $cigar_type eq "I") {
02016     #add to sequence_position but not alignment_position
02017     $sequence_position += $cigar_count;
02018       } elsif( $cigar_type eq "G" || $cigar_type eq "D" || $cigar_type eq "X") {
02019         $alignment_position += $cigar_count;
02020       }
02021     }
02022   } else {
02023     foreach my $cigar_piece (@cigar_pieces) {
02024       my $cigar_type = substr($cigar_piece, -1, 1 );
02025       my $cigar_count = substr($cigar_piece, 0 ,-1 );
02026       $cigar_count = 1 unless ($cigar_count =~ /^\d+$/);
02027       next if ($cigar_count < 1);
02028   
02029       if( $cigar_type eq "M" ) {
02030         $mapper->add_map_coordinates(
02031                 "sequence", #$self->dbID,
02032                 $sequence_position - $cigar_count + 1,
02033                 $sequence_position,
02034                 $rel_strand,
02035                 "alignment", #$self->genomic_align_block->dbID,
02036                 $alignment_position,
02037                 $alignment_position + $cigar_count - 1
02038             );
02039         $sequence_position -= $cigar_count;
02040         $alignment_position += $cigar_count;
02041       } elsif( $cigar_type eq "I") {
02042     #add to sequence_position but not alignment_position
02043     $sequence_position -= $cigar_count;
02044       } elsif( $cigar_type eq "G" || $cigar_type eq "D" || $cigar_type eq "X") {
02045         $alignment_position += $cigar_count;
02046       }
02047     }
02048   }
02049 
02050   return $mapper;
02051 }
02052 
02053 
02054 =head2 get_Slice
02055 
02056   Arg[1]     : -none-
02057   Example    : $slice = $genomic_align->get_Slice();
02058   Description: creates and returns a Bio::EnsEMBL::Slice which corresponds to
02059                this Bio::EnsEMBL::Compara::GenomicAlign
02060   Returntype : Bio::EnsEMBL::Slice object
02061   Exceptions : return -undef- if slice cannot be created (this is likely to
02062                happen if the Registry is misconfigured)
02063   Status     : Stable
02064 
02065 =cut
02066 
02067 sub get_Slice {
02068   my ($self) = @_;
02069 
02070   my $slice = $self->dnafrag->slice;
02071   return undef if (!defined($slice));
02072 
02073   $slice = $slice->sub_Slice(
02074               $self->dnafrag_start,
02075               $self->dnafrag_end,
02076               $self->dnafrag_strand
02077           );
02078 
02079   return $slice;
02080 }
02081 
02082 
02083 =head2 restrict
02084 
02085   Arg[1]     : int start
02086   Arg[1]     : int end
02087   Example    : my $genomic_align = $genomic_align->restrict(10, 20);
02088   Description: restrict (trim) this GenomicAlign to the start and end
02089                positions (in alignment coordinates). If no trimming is
02090                required, the original object is returned instead.
02091   Returntype : Bio::EnsEMBL::Compara::GenomicAlign object
02092   Exceptions :
02093   Status     : At risk
02094 
02095 =cut
02096 
02097 sub restrict {
02098   my ($self, $start, $end, $aligned_seq_length) = @_;
02099   throw("Wrong arguments") if (!$start or !$end);
02100 
02101   my $restricted_genomic_align = $self->copy();
02102   delete($restricted_genomic_align->{dbID});
02103   delete($restricted_genomic_align->{genomic_align_block_id});
02104   delete($restricted_genomic_align->{original_sequence});
02105   delete($restricted_genomic_align->{aligned_sequence});
02106   delete($restricted_genomic_align->{cigar_line});
02107   $restricted_genomic_align->{original_dbID} = $self->dbID if ($self->dbID);
02108 
02109   # Need to calculate the original aligned sequence length myself
02110   if (!$aligned_seq_length) {
02111     my @cigar = grep {$_} split(/(\d*[GDMXI])/, $self->cigar_line);
02112     foreach my $num_type (@cigar) {
02113       my $type = substr($num_type, -1, 1, "");
02114       $num_type = 1 if ($num_type eq "");
02115       $aligned_seq_length += $num_type unless ($type eq "I");
02116     }
02117   }
02118 
02119   my $final_aligned_length = $end - $start + 1;
02120   my $number_of_columns_to_trim_from_the_start = $start - 1;
02121   my $number_of_columns_to_trim_from_the_end = $aligned_seq_length - $end;
02122 
02123   my @cigar = grep {$_} split(/(\d*[GDMXI])/, $self->cigar_line);
02124 
02125   ## Trim start of cigar_line if needed
02126   if ($number_of_columns_to_trim_from_the_start >= 0) {
02127     my $counter_of_trimmed_columns_from_the_start = 0;
02128     my $counter_of_trimmed_base_pairs = 0; # num of bp we trim (from the start)
02129     ## Loop through the cigar pieces
02130     while (my $cigar = shift(@cigar)) {
02131       # Parse each cigar piece
02132       my $type = substr( $cigar, -1, 1 );
02133       my $num = substr( $cigar, 0 ,-1 );
02134       $num = 1 if ($num eq "");
02135 
02136       # Insertions are not part of the alignment, don't count them
02137       if ($type ne "I") {
02138         $counter_of_trimmed_columns_from_the_start += $num;
02139       }
02140 
02141       # Matches and insertions are actual base pairs in the sequence
02142       if ($type eq "M" || $type eq "I") {
02143         $counter_of_trimmed_base_pairs += $num;
02144       }
02145 
02146       # If this cigar piece is too long and we overshoot the number of columns we want to trim,
02147       # we substitute this cigar piece by a shorter one
02148       if ($counter_of_trimmed_columns_from_the_start >= $number_of_columns_to_trim_from_the_start) {
02149         my $new_cigar_piece;
02150         # length of the new cigar piece
02151         my $length = $counter_of_trimmed_columns_from_the_start - $number_of_columns_to_trim_from_the_start;
02152         if ($length > 1) {
02153           $new_cigar_piece = $length.$type;
02154         } elsif ($length == 1) {
02155           $new_cigar_piece = $type;
02156         }
02157         unshift(@cigar, $new_cigar_piece) if ($new_cigar_piece);
02158         if ($type eq "M") {
02159           $counter_of_trimmed_base_pairs -= $length;
02160         }
02161 
02162         ## We don't want to start with an insertion. Trim it!
02163         while (@cigar and $cigar[0] =~ /[I]/) {
02164       my ($num, $type) = ($cigar[0] =~ /^(\d*)([DIGMX])/);  
02165           #my $type = substr( $cigar, -1, 1 );
02166           #my $num = substr( $cigar, 0 ,-1 );
02167           $num = 1 if ($num eq "");
02168           $counter_of_trimmed_base_pairs += $num;
02169           shift(@cigar);
02170         }
02171         last;
02172       }
02173     }
02174     if ($self->dnafrag_strand == 1) {
02175       $restricted_genomic_align->dnafrag_start($self->dnafrag_start + $counter_of_trimmed_base_pairs);
02176     } else {
02177       $restricted_genomic_align->dnafrag_end($self->dnafrag_end - $counter_of_trimmed_base_pairs);
02178     }
02179   }
02180 
02181   ## Trim end of cigar_line if needed
02182   if ($number_of_columns_to_trim_from_the_end >= 0) {
02183     my $counter_of_trimmed_columns_from_the_end = 0;
02184     my $counter_of_trimmed_base_pairs = 0; # num of bp we trim (from the start)
02185     ## Loop through the cigar pieces
02186     while (my $cigar = pop(@cigar)) {
02187       # Parse each cigar piece
02188       my $type = substr( $cigar, -1, 1 );
02189       my $num = substr( $cigar, 0 ,-1 );
02190       $num = 1 if ($num eq "");
02191 
02192       # Insertions are not part of the alignment, don't count them
02193       if ($type ne "I") {
02194         $counter_of_trimmed_columns_from_the_end += $num;
02195       }
02196 
02197       # Matches and insertions are actual base pairs in the sequence
02198       if ($type eq "M" || $type eq "I") {
02199         $counter_of_trimmed_base_pairs += $num;
02200       }
02201       # If this cigar piece is too long and we overshoot the number of columns we want to trim,
02202       # we substitute this cigar piece by a shorter one
02203       if ($counter_of_trimmed_columns_from_the_end >= $number_of_columns_to_trim_from_the_end) {
02204         my $new_cigar_piece;
02205         # length of the new cigar piece
02206         my $length = $counter_of_trimmed_columns_from_the_end - $number_of_columns_to_trim_from_the_end;
02207         if ($length > 1) {
02208           $new_cigar_piece = $length.$type;
02209         } elsif ($length == 1) {
02210           $new_cigar_piece = $type;
02211         }
02212         push(@cigar, $new_cigar_piece) if ($new_cigar_piece);
02213         if ($type eq "M") {
02214           $counter_of_trimmed_base_pairs -= $length;
02215         }
02216 
02217         ## We don't want to end with an insertion. Trim it!
02218         while (@cigar and $cigar[-1] =~ /[I]/) {
02219       my ($num, $type) = ($cigar[-1] =~ /^(\d*)([DIGMX])/);
02220           #my $type = substr( $cigar, -1, 1 );
02221           #my $num = substr( $cigar, 0 ,-1 );
02222           $num = 1 if ($num eq "");
02223           $counter_of_trimmed_base_pairs += $num;
02224           pop(@cigar);
02225         }
02226         last;
02227       }
02228     }
02229     if ($self->dnafrag_strand == 1) {
02230       $restricted_genomic_align->dnafrag_end($restricted_genomic_align->dnafrag_end - $counter_of_trimmed_base_pairs);
02231     } else {
02232       $restricted_genomic_align->dnafrag_start($restricted_genomic_align->dnafrag_start + $counter_of_trimmed_base_pairs);
02233     }
02234   }
02235 
02236   ## Save genomic_align's cigar_line
02237   $restricted_genomic_align->aligned_sequence(0);
02238   $restricted_genomic_align->cigar_line(join("", @cigar));
02239 
02240   return $restricted_genomic_align;
02241 }
02242 
02243 1;