view release on metacpan or  search on metacpan

lib/Bio/Polloc/LociGroup.pm  view on Meta::CPAN

   return if not $force and defined $self->{'_fixed_strands'} and $self->{'_fixed_strands'} == $#{$self->loci};
   $self->{'_fixed_strands'} = $#{$self->loci};
   $self->_load_module('Bio::Polloc::GroupCriteria');
   $self->_load_module('Bio::Tools::Run::Alignment::Muscle');
   return unless $#{$self->loci}>0; # No need to check
   $size ||= 500;
   
   my $factory = Bio::Tools::Run::Alignment::Muscle->new();
   $factory->quiet(1);
   
   # Find a suitable reference
   my $ref = [undef, undef];
   LOCUS: for my $lk (1 .. $#{$self->loci}){
      my $ref_test = [
      		Bio::Polloc::GroupCriteria->_build_subseq(
				$self->loci->[$lk]->seq,
				$self->loci->[$lk]->from - $size,
				$self->loci->[$lk]->from),
      		Bio::Polloc::GroupCriteria->_build_subseq(
				$self->loci->[$lk]->seq,
				$self->loci->[$lk]->to,
				$self->loci->[$lk]->to + $size)
		];
      if(defined $ref->[0] and defined $ref->[1]){
         # Longer pair:
	 $ref = $ref_test
	 	if  defined $ref_test->[0] and defined $ref_test->[1]
		and $ref_test->[0]->length >= $ref->[0]->length
		and $ref_test->[1]->length >= $ref->[1]->length;
      }elsif(defined $ref->[0] or defined $ref->[1]){
         # Both sequences defined:
	 $ref = $ref_test if defined $ref_test->[0] and defined $ref_test->[1];
      }else{
         # At least one sequence defined:
	 $ref = $ref_test if defined $ref_test->[0] or defined $ref_test->[1];
      }
   }
   unless(defined $ref->[0] or defined $ref->[1]){
      $self->debug('Impossible to find a suitable reference');
      return;
   }
   $ref = defined $ref->[0] ?
   		( defined $ref->[1] ?
			Bio::Seq->new(-seq=>$ref->[0]->seq . ("N"x20) . $ref->[1]->seq)
			: $ref->[0]
		) : $ref->[1];
   
   $ref->id('ref');
   $self->loci->[0]->strand('+');
   
   # Compare
   LOCUS: for my $k (0 .. $#{$self->loci}){
      my $tgt = Bio::Polloc::GroupCriteria->_build_subseq(
      		$self->loci->[$k]->seq,
		$self->loci->[$k]->from-$size,
		$self->loci->[$k]->to+$size);
      next LOCUS unless $tgt; # <- This may be way too paranoic!
      $tgt->id('tgt');
      my $tgtrc = $tgt->revcom;
      $self->debug("Setting strand for ".$self->loci->[$k]->id) if defined $self->loci->[$k]->id;
      my $eval_fun = 'average_percentage_identity';
      #$eval_fun = 'overall_percentage_identity';
      if($factory->align([$ref, $tgt])->$eval_fun
      		< $factory->align([$ref,$tgtrc])->$eval_fun){
         $self->debug("Assuming negative strand, setting locus orientation");
	 $self->loci->[$k]->strand('-');
      }else{
         $self->debug("Assuming positive strand, setting locus orientation");
         $self->loci->[$k]->strand('+');
      }
   } # LOCUS
}

=head1 INTERNAL METHODS

Methods intended to be used only within the scope of Bio::Polloc::*

=head2 _initialize

=cut

sub _initialize {
   my ($self, @args) = @_;
   my($name, $featurename, $genomes) = $self->_rearrange([qw(NAME FEATURENAME GENOMES)], @args);
   $self->name($name);
   $self->featurename($featurename);
   $self->genomes($genomes);
}

1;