File Coverage

blib/lib/CracTools/GenomeMask.pm

Criterion	Covered	Total	%
statement	94	99	94.9
branch	18	30	60.0
condition	2	3	66.6
subroutine	15	15	100.0
pod	10	10	100.0
total	139	157	88.5

package CracTools::GenomeMask;

{
  $CracTools::GenomeMask::DIST = 'CracTools';
}
# ABSTRACT: A bit vector mask over the whole genome
$CracTools::GenomeMask::VERSION = '1.251';
use strict;
use warnings;

use CracTools::BitVector;
use CracTools::Utils;
use Carp;

#You can also define is the genome should be consider as uniquely stranded or as double stranded
#with the C argument.


sub new {
  my $class = shift;

  my %args = @_;

  # the genome mask is not stranded by default
  #my $is_stranded = defined $args{is_stranded}? $args{is_stranded} : 0;
  my $verbose = defined $args{verbose}? $args{verbose} : 0;

  my %bit_vectors;

  if(defined $args{genome}) {

    foreach my $chr (keys %{$args{genome}}) {
      if(!defined $bit_vectors{$chr}) {
        $bit_vectors{$chr} = CracTools::BitVector->new($args{genome}->{$chr});
      } else {
        croak "Multiple definition of sequence $chr inf the genome lengths";
      }
    }
  } elsif(defined $args{crac_index_conf}) {
    print STDERR "Creating GenomeMask from crac index conf file : $args{crac_index_conf}\n" if $verbose;
    my $conf_fh = CracTools::Utils::getReadingFileHandle($args{crac_index_conf});
    my $nb_chr = <$conf_fh>;
    my $nb_chr_found = 0;
    while(<$conf_fh>) {
      my $chr = $_;
      chomp $chr;
      my $chr_length = <$conf_fh>;
      chomp $chr_length;
      if(defined $chr_length) {
        print STDERR "\tCreating bitvecor for chr $chr of length $chr_length\n" if $verbose;
        $bit_vectors{$chr} = CracTools::BitVector->new($chr_length);
        $nb_chr_found++;
      } else {
        croak "Missing genome length for chromosome $chr";
      }
    }
    croak "There is less chromosome found ($nb_chr_found) in $args{crac_index_conf} than expected ($nb_chr)" if $nb_chr_found < $nb_chr;
  } elsif(defined $args{sam_reader}) {
    my $refseq_lengths = $args{sam_reader}->allRefSeqLengths();
    foreach my $chr (keys %{$refseq_lengths}) {
      $bit_vectors{$chr} = CracTools::BitVector->new($refseq_lengths->{$chr});
    }
  } else {
    croak "There is no valid argument to extract the chromosomes names and length";
  }

  my $self = bless {
    bit_vectors => \%bit_vectors,
    #is_stranded => $is_stranded,
  }, $class;

  return $self;
}

#=head2 isStranded
#
#  Description : Return true is genomeMask is double stranded
#
#=cut
#
#sub isStranded {
#  my $self = shift;
#  return $self->{is_stranded};
#}

#Arg [2] : Integer (1,-1) - strand

sub getBitvector {
  my $self = shift;
  my $chr = shift;
  #my $strand = shift;
  if(defined $self->{bit_vectors}->{$chr}) {
    return $self->{bit_vectors}->{$chr};
  } else {
    carp "There is no bitvector for sequence $chr in the genome mask";
    return undef;
  }
}


sub getChrLength {
  my $self = shift;
  my $chr = shift;
  my $bv = $self->getBitvector($chr);
  return defined $bv? $bv->length : undef;
}


sub setPos {
  my ($self,$chr,$pos) = @_;
  my $bv = $self->getBitvector($chr);
  $bv->set($pos) if defined $bv;
}


sub setRegion {
  my ($self,$chr,$start,$end) = @_;
  for(my $i = $start; $i <= $end; $i++) {
    $self->setPos($chr,$i);
  }
}


sub getPos {
  my ($self,$chr,$pos) = @_;
  my $bv = $self->getBitvector($chr);
  return $bv->get($pos) if defined $bv;
}


sub getPosSetInRegion {
  my ($self,$chr,$start,$end) = @_;
  my $bv = $self->getBitvector($chr);
  my @pos;
  if(defined $bv) {
    for(my $i = $start; $i <= $end; $i++) {
      push(@pos,$i) if $bv->get($i) == 1;
    }
  }
  return \@pos;
}


sub getNbBitsSetInRegion {
  my $self = shift;
  return scalar @{$self->getPosSetInRegion(@_)};
}


sub rank {
  my ($self,$chr,$pos) = @_;
  my $cumulated_bits = 0;
  my $i = 0;
  my @chr_sorted = sort keys %{$self->{bit_vectors}};
  while($chr_sorted[$i] ne $chr) {
    $cumulated_bits += $self->getBitvector($chr_sorted[$i])->nb_set;
    $i++;
  }
  return $cumulated_bits + $self->getBitvector($chr_sorted[$i])->rank($pos);
}


sub select {
  my $self = shift;
  my $i = shift;
  my $cumulated_bits = 0;
  my @chr_sorted = sort keys %{$self->{bit_vectors}};
  my $j = 0;
  while($j < @chr_sorted && $cumulated_bits + $self->getBitvector($chr_sorted[$j])->nb_set < $i) {
    my $chr = $chr_sorted[$j];
    my $bv = $self->getBitvector($chr);
    $cumulated_bits += $self->getBitvector($chr)->nb_set;
    $j++;
  }
  my $chr = $chr_sorted[$j-1];
  my $pos = $self->getBitvector($chr_sorted[$j-1])->select($i - $cumulated_bits + 1);
  return ($chr,$pos);
}

1;

__END__