Bio::EnsEMBL::Variation::DBSQL ReadCoverageAdaptor
SummaryIncluded librariesPackage variablesSynopsisDescriptionGeneral documentationMethods
Toolbar
WebCvsRaw content
Summary
Bio::EnsEMBL::Variation::DBSQL::ReadCoverageAdaptor
Package variables
No package variables defined.
Included modules
Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor
Bio::EnsEMBL::Mapper::RangeRegistry
Bio::EnsEMBL::Utils::Exception qw ( throw warning )
Bio::EnsEMBL::Variation::ReadCoverage
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
_columns
No description
Code
_initialize_range_registry
No description
Code
_invert_pair
No description
Code
_objs_from_sth
No description
Code
_register_range_level
No description
Code
_tables
No description
Code
fetch_all_by_Slice_Sample_depthDescriptionCode
fetch_all_regions_covered
No description
Code
get_coverage_levelsDescriptionCode
Methods description
fetch_all_by_Slice_Sample_depthcode    nextTop
    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
get_coverage_levelscodeprevnextTop
    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
_columnsdescriptionprevnextTop
sub _columns {
    return qw (rc.seq_region_id rc.seq_region_start rc.seq_region_end rc.level rc.sample_id);
}
_initialize_range_registrydescriptionprevnextTop
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;
}
_invert_pairdescriptionprevnextTop
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
}
_objs_from_sthdescriptionprevnextTop
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;
}
_register_range_leveldescriptionprevnextTop
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
}
_tablesdescriptionprevnextTop
sub _tables {
return (['read_coverage','rc']		     )
}
fetch_all_by_Slice_Sample_depthdescriptionprevnextTop
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
}
fetch_all_regions_covereddescriptionprevnextTop
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);
}
get_coverage_levelsdescriptionprevnextTop
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
AUTHOR - Daniel RiosTop
CONTACTTop
Post questions to the Ensembl development list ensembl-dev@ebi.ac.uk