Toolbar

Summary

Bio::EnsEMBL::Variation::DBSQL::ReadCoverageAdaptor

Package variables

Included modules

Inherit

Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor

Synopsis

  $vdb = Bio::EnsEMBL::Variation::DBSQL::DBAdaptor->new(...);
  $db  = Bio::EnsEMBL::DBSQL::DBAdaptor->new(...);

  # tell the variation database where core database information can be
  # be found
  $vdb->dnadb($db);

  $rca = $vdb->get_ReadCoverageAdaptor();
  $sa  = $db->get_SliceAdaptor();
  $pa  = $vdb->get_PopulationAdaptor();

  # get read coverage in a region for a certain population in in a certain level
  $slice = $sa->fetch_by_region('chromosome', 'X', 1e6, 2e6);
  $population = $pa->fetch_by_name("UNKNOWN");
  $level = 1;
  foreach $rc (@{$vfa->fetch_all_by_Slice_Sample_depth($slice,$population,$level)}) {
    print $rc->start(), '-', $rc->end(), "\n";
  }

Description

This adaptor provides database connectivity for ReadCoverage objects.
Coverage information for reads can be obtained from the database using this
adaptor. See the base class BaseFeatureAdaptor for more information.

Methods

Methods description

    Arg[0]      : Bio::EnsEMBL::Slice $slice
    Arg[1]      : (optional) Bio::EnsEMBL::Variation::Sample $sample
    Arg[2]      : (optional) int $level
    Example     : my $features = $rca->fetch_all_by_Slice_Sample_depth($slice,$sample,$level); or
                  my $features = $rca->fetch_all_by_Slice_Sample_depth($slice, $sample); or
                  my $features = $rca->fetch_all_by_Slice_Sample_depth($slice, $level); or
                  my $features = $rca->fetch_all_by_Slice_Sample_depth($slice); 
    Description : Gets all the read coverage features for a given sample(strain or individual) in a certain level
                  in the provided slice
    ReturnType  : listref of Bio::EnsEMBL::Variation::ReadCoverage
    Exceptions  : thrown on bad arguments
    Caller      : general
    Status      : At Risk

    Args        : none
    Example     : my @coverage_levels = @{$rca->fetch_coverage_levels()};
    Description : Gets the read coverage depths calculated in the database
    ReturnType  : listref of integer
    Exceptions  : none
    Caller      : general
    Status      : At Risk

Methods code

sub _columns {

    return qw (rc.seq_region_id rc.seq_region_start rc.seq_region_end rc.level rc.sample_id);

}

sub _initialize_range_registry {

    my $range_registry = shift;
    my $max_level = shift;

    foreach my $level (1..$max_level){
	$range_registry->[$level] = Bio::EnsEMBL::Mapper::RangeRegistry->new();
    }

    return;

}

sub _invert_pair {

    my $range = shift; #initial range of the array
    my $pairs = shift; #listref with the pairs that have been added to the range

    my @inverted_pairs;
    my $inverted;

    my $rr = Bio::EnsEMBL::Mapper::RangeRegistry->new();

    foreach my $pair (@{$pairs}){
	$rr->check_and_register(1,$pair->[0],$pair->[1]);
    }
    return $rr->check_and_register(1,$range->[0],$range->[1]); #register again the range

}

sub _objs_from_sth {

  my ($self, $sth, $mapper, $dest_slice) = @_;
  
  my $sa = $self->db()->dnadb()->get_SliceAdaptor();

  #
  # This code is ugly because an attempt has been made to remove as many
  # function calls as possible for speed purposes.  Thus many caches and
  # a fair bit of gymnastics is used.
  #

  my ($seq_region_id, $seq_region_start, $seq_region_end, $level, $sample_id);
  my $seq_region_strand = 1; #assume all reads are in the + strand

  $sth->bind_columns(\$seq_region_id,\$ seq_region_start,\$ seq_region_end,\$ level,\$ sample_id);

  my @features;
  my %seen_pops;  #hash containing the populations seen
  my %slice_hash;
  my %sr_name_hash;
  my %sr_cs_hash;
  my $asm_cs;
  my $cmp_cs;
  my $asm_cs_vers;
  my $asm_cs_name;
  my $cmp_cs_vers;
  my $cmp_cs_name;

  if($mapper) {
      $asm_cs = $mapper->assembled_CoordSystem();
      $cmp_cs = $mapper->component_CoordSystem();
      $asm_cs_name = $asm_cs->name();
      $asm_cs_vers = $asm_cs->version();
      $cmp_cs_name = $cmp_cs->name();
      $cmp_cs_vers = $cmp_cs->version();
  }
  
  my $dest_slice_start;
  my $dest_slice_end;
  my $dest_slice_strand;
  my $dest_slice_length;

  if($dest_slice) {
      $dest_slice_start  = $dest_slice->start();
      $dest_slice_end    = $dest_slice->end();
      $dest_slice_strand = $dest_slice->strand();
      $dest_slice_length = $dest_slice->length();
  }

  my $read_coverage;
  FEATURE: while ($sth->fetch()){

      #get the population_adaptor object
      my $pa = $self->db()->get_PopulationAdaptor();
      my $ia = $self->db()->get_IndividualAdaptor();
      #get the slice object
      my $slice = $slice_hash{"ID:".$seq_region_id};
      if(!$slice) {
	  $slice = $sa->fetch_by_seq_region_id($seq_region_id);
	  $slice_hash{"ID:".$seq_region_id} = $slice;
	  $sr_name_hash{$seq_region_id} = $slice->seq_region_name();
	  $sr_cs_hash{$seq_region_id} = $slice->coord_system();
      }
      #
      # remap the feature coordinates to another coord system
      # if a mapper was provided
      #
      if($mapper) {
	  my $sr_name = $sr_name_hash{$seq_region_id};
	  my $sr_cs   = $sr_cs_hash{$seq_region_id};
	  
	  ($sr_name,$seq_region_start,$seq_region_end,$seq_region_strand) =
	      $mapper->fastmap($sr_name, $seq_region_start, $seq_region_end,
			       $seq_region_strand, $sr_cs);
	  
	  #skip features that map to gaps or coord system boundaries
	  next FEATURE if(!defined($sr_name));
	  
	  #get a slice in the coord system we just mapped to
	  if($asm_cs == $sr_cs || ($cmp_cs != $sr_cs && $asm_cs->equals($sr_cs))) {
	      $slice = $slice_hash{"NAME:$sr_name:$cmp_cs_name:$cmp_cs_vers"} ||=
		  $sa->fetch_by_region($cmp_cs_name, $sr_name,undef, undef, undef,
				       $cmp_cs_vers);
	  } else {
	      $slice = $slice_hash{"NAME:$sr_name:$asm_cs_name:$asm_cs_vers"} ||=
		  $sa->fetch_by_region($asm_cs_name, $sr_name, undef, undef, undef,
				       $asm_cs_vers);
	  }
      }
      
      #
      # If a destination slice was provided convert the coords
      # If the dest_slice starts at 1 and is foward strand, nothing needs doing
      #
      if($dest_slice) {
	  if($dest_slice_start != 1 || $dest_slice_strand != 1) {
	      if($dest_slice_strand == 1) {
		  $seq_region_start = $seq_region_start - $dest_slice_start + 1;
		  $seq_region_end   = $seq_region_end   - $dest_slice_start + 1;
	      } else {
		  my $tmp_seq_region_start = $seq_region_start;
		  $seq_region_start = $dest_slice_end - $seq_region_end + 1;
		  $seq_region_end   = $dest_slice_end - $tmp_seq_region_start + 1;
		  $seq_region_strand *= -1;
	      }
	      
	      #throw away features off the end of the requested slice
	      if($seq_region_end < 1 || $seq_region_start > $dest_slice_length) {
		  next FEATURE;
	      }
	  }
	  $slice = $dest_slice;
      }

      my $sample;
      if($sample_id){
	  $sample = $seen_pops{$sample_id} ||= $pa->fetch_by_dbID($sample_id);
	  $sample = $seen_pops{$sample_id} ||= $ia->fetch_by_dbID($sample_id);
      }
      $read_coverage = Bio::EnsEMBL::Variation::ReadCoverage->new(-start => $seq_region_start,
								  -end   => $seq_region_end,
								  -slice => $slice,
								  -adaptor  => $self,
								  -level => $level,
								  -sample => $sample,
								  -strand => 1
								  );
      push @features, $read_coverage;
  }
  $sth->finish();
  return\@ features;
}

1;

}

sub _register_range_level {

    my $range_registry = shift;
    my $range = shift;
    my $level = shift;
    my $max_level = shift;
    
    return if ($level > $max_level);
    my $rr = $range_registry->[$level];
    my $pair = $rr->check_and_register(1,$range->[0],$range->[1]);
    my $pair_inverted = _invert_pair($range,$pair);
    return if (!defined $pair_inverted);
    foreach my $inverted_range (@{$pair_inverted}){
	_register_range_level($range_registry,$inverted_range,$level+1, $max_level);
    }
}

#for a given range and the one covered, returns the inverted

}

sub _tables {

return (['read_coverage','rc']		     )

}

sub fetch_all_by_Slice_Sample_depth {

    my $self = shift;
    my $slice = shift;
    my @args = @_; #can contain individual and/or level
    my $rcs;

    if(!ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
	throw('Bio::EnsEMBL::Slice arg expected');
    }
    if (defined $args[0]){ #contains, at least, 1 parameter, either a Individual or the level
	my $constraint;
	my $levels = $self->get_coverage_levels();
	if (defined $args[1]){ #contains both parameters, first the Individual and second the level
	    if(!ref($args[0]) || !$args[0]->isa('Bio::EnsEMBL::Variation::Sample')) {
		throw('Bio::EnsEMBL::Variation::Sample arg expected');
	    }
	    if(!defined($args[0]->dbID())) {
		throw("Sample arg must have defined dbID");
	    }
	    if ((grep {$args[1] == $_} @{$levels}) > 0){
	      $constraint = "rc.sample_id = " . $args[0]->dbID . " AND rc.level = " . $args[1];
	     # $constraint = "rc.sample_id = ? AND rc.level = ?";
	     # $self->bind_param_generic_fetch($args[0]->dbID,SQL_INTEGER);
	     # $self->bind_param_generic_fetch($args[1],SQL_INTEGER);
	    }
	    else{
		warning("Level must be a number of: " . join(",", @{$levels}));
		return [];
	    }
	}
	else{ #there is just 1 argument, can either be the Individual or the level
	    if (!ref($args[0])){
		#it should just contain the level
		if ((grep {$args[0] == $_} @{$levels}) > 0){
		  $constraint = "rc.level = " . $args[0];
		  #$constraint = "rc.level = ? ";
		  #$self->bind_param_generic_fetch($args[0],SQL_INTEGER);
		}
		else{
		    warning("Level must be a number of: " . join(",", @{$levels}));
		    return [];
		}									   
	    }
	    else{
		#it should contain the Individual
		if (!$args[0]->isa('Bio::EnsEMBL::Variation::Sample')){
		    throw('Bio::EnsEMBL::Variation::Sample arg expected');
		}
		$constraint = "rc.sample_id = " . $args[0]->dbID; 
		#$constraint = "rc.sample_id = ?";
		#$self->bind_param_generic_fetch($args[0]->dbID,SQL_INTEGER);
	    }
	}
	$rcs = $self->fetch_all_by_Slice_constraint($slice,$constraint);    
	return $rcs;
    }
    #call the method fetch_all_by_Slice
    $rcs = $self->fetch_all_by_Slice($slice);    
    return $rcs;
}

#returns a list of regions that are covered by all the sample given

}

sub fetch_all_regions_covered {

    my $self = shift;
    my $slice = shift;
    my $samples = shift;  #listref of sample names to get the coverage from

    my $ind_adaptor = $self->db->get_IndividualAdaptor;
    my $range_registry = [];
    my $max_level = scalar(@{$samples});
    _initialize_range_registry($range_registry,$max_level);
    foreach my $sample_name (@{$samples}){
	my $sample = shift@{$ind_adaptor->fetch_all_by_name($sample_name)};
	my $coverage = $self->fetch_all_by_Slice_Sample_depth($slice,$sample,1); #get coverage information
	foreach my $cv_feature (@{$coverage}){
	    my $range = [$cv_feature->seq_region_start,$cv_feature->seq_region_end]; #store toplevel coordinates
	    _register_range_level($range_registry,$range,1,$max_level);	   
	}
    }
    return $range_registry->[$max_level]->get_ranges(1);

}

sub get_coverage_levels {

    my $self = shift;
    my @levels;
    my $level_coverage;
    my $sth = $self->prepare(qq{SELECT meta_value from meta where meta_key = 'read_coverage.coverage_level'
				});
    $sth->execute();
    $sth->bind_columns(\$level_coverage);
    while ($sth->fetch()){
	push @levels, $level_coverage;
    }
    $sth->finish();

    return\@ levels;

}

General documentation

Post questions to the Ensembl development list ensembl-dev@ebi.ac.uk