Bio::EnsEMBL::Variation::DBSQL LDFeatureContainerAdaptor
SummaryIncluded librariesPackage variablesSynopsisDescriptionGeneral documentationMethods
Toolbar
WebCvsRaw content
Summary
Bio::EnsEMBL::Variation::DBSQL::LDFeatureContainerAdaptor
Package variables
No package variables defined.
Included modules
Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor
Bio::EnsEMBL::Utils::Exception qw ( throw warning )
Bio::EnsEMBL::Variation::LDFeatureContainer
Data::Dumper
FileHandle
POSIX
Inherit
Bio::EnsEMBL::DBSQL::BaseAdaptor
Synopsis
  $vdb = Bio::EnsEMBL::Variation::DBSQL::DBAdaptor->new(...);
$db = Bio::EnsEMBL::DBSQL::DBAdaptor->new(...);
$sa = $db->get_SliceAdaptor(); $lda = $vdb->get_LDFeatureContainerAdaptor(); $vfa = $vdb->get_VariationFeatureAdaptor(); # Get a LDFeatureContainer in a region $slice = $sa->fetch_by_region('chromosome', 'X', 1e6, 2e6); $ldContainer = $lda->fetch_by_Slice($slice); print "Name of the ldContainer is: ", $ldContainer->name(); # fetch ld featureContainer for a particular variation feature $vf = $vfa->fetch_by_dbID(145); $ldContainer = $lda->fetch_by_VariationFeature($vf); print "Name of the ldContainer: ", $ldContainer->name();
Description
This adaptor provides database connectivity for LDFeature objects.
LD Features may be retrieved from the Ensembl variation database by
several means using this module.
Methods
_convert_genotype
No description
Code
_get_LD_populations
No description
Code
_get_siblings
No description
Code
_objs_from_sth
No description
Code
_objs_from_sth_temp_file
No description
Code
_store_genotype
No description
Code
executable
No description
Code
fetch_by_SliceDescriptionCode
fetch_by_VariationFeatureDescriptionCode
temp_path
No description
Code
Methods description
fetch_by_Slicecode    nextTop
  Arg [1]    : Bio::EnsEMBL::Slice $slice
The slice to fetch genes on. Assuming it is always correct (in the top level)
Arg [2] : (optional) Bio::EnsEMBL::Variation::Population $population. Population where
we want to select the LD information
Example : $ldFeatureContainer = $ldfeaturecontainer_adaptor->fetch_by_Slice($slice);
Description: Overwrites superclass method to add the name of the slice to the LDFeatureContainer.
Returntype : Bio::EnsEMBL::Variation::LDFeatureContainer
Exceptions : thrown on bad argument
Caller : general
Status : At Risk
fetch_by_VariationFeaturecodeprevnextTop
  Arg [1]    : Bio::EnsEMBL:Variation::VariationFeature $vf
Arg [2] : (optional) int $population_id. Population where we want to select the LD information
Example : my $ldFeatureContainer = $ldFetureContainerAdaptor->fetch_by_VariationFeature($vf); Description: Retrieves LDFeatureContainer for a given variation feature. Most
variations should only hit the genome once and only a return
a single variation feature.
Returntype : reference to Bio::EnsEMBL::Variation::LDFeatureContainer
Exceptions : throw on bad argument
Caller : general
Status : At Risk
Methods code
_convert_genotypedescriptionprevnextTop
sub _convert_genotype {
    my $self = shift;
    my $alleles_variation = shift; #reference to the hash containing the alleles for the variation present in the genotypes
my $individual_information = shift; #reference to a hash containing the values to be written to the file
my @alleles_ordered; #the array will contain the alleles ordered by apparitions in the genotypes (only 2 values possible)
@alleles_ordered = sort({$alleles_variation->{$b} <=> $alleles_variation->{$a}} keys %{$alleles_variation}); #let's convert the allele_1 allele_2 to a genotype in the AA, Aa or aa format, where A corresponds to the major allele and a to the minor
foreach my $individual_id (keys %{$individual_information}){ #if both alleles are different, this is the Aa genotype
if ($individual_information->{$individual_id}{allele_1} ne $individual_information->{$individual_id}{allele_2}){ $individual_information->{$individual_id}{genotype} = 'Aa'; } #when they are the same, must find out which is the major
else{ if ($alleles_ordered[0] eq $individual_information->{$individual_id}{allele_1}){ #it is the major allele
$individual_information->{$individual_id}{genotype} = 'AA'; } else{ $individual_information->{$individual_id}{genotype} = 'aa'; } } }
}
_get_LD_populationsdescriptionprevnextTop
sub _get_LD_populations {
    my $self = shift;
    my $siblings = shift;
    my ($pop_id,$population_name);
    my $sth = $self->db->dbc->prepare(qq{SELECT s.sample_id, s.name
				     FROM population p, sample s
				     WHERE (s.name like 'PERLEGEN:AFD%'
				     OR s.name like 'CSHL-HAPMAP%')
				     AND s.sample_id = p.sample_id});

    $sth->execute();
    $sth->bind_columns(\$pop_id,\$population_name);
    #get all the children that we do not want in the genotypes
my @pops; while($sth->fetch){ if($population_name =~ /CEU|YRI/){ $self->_get_siblings($pop_id,$siblings); } push @pops, $pop_id; } my $in_str = " IN (" . join(',', @pops). ")"; return $in_str if (defined $pops[0]); return '' if (!defined $pops[0]);
}
_get_siblingsdescriptionprevnextTop
sub _get_siblings {
    my $self = shift;
    my $population_id = shift;
    my $siblings = shift;

    my $sth_individual = $self->db->dbc->prepare(qq{SELECT i.sample_id
							     FROM individual i, individual_population ip
							     WHERE ip.individual_sample_id = i.sample_id
							     AND ip.population_sample_id = ? 
							     AND i.father_individual_sample_id IS NOT NULL
							     AND i.mother_individual_sample_id IS NOT NULL
							 });
    my ($individual_id);
    $sth_individual->execute($population_id);
    $sth_individual->bind_columns(\$individual_id);
    while ($sth_individual->fetch){
	$siblings->{$population_id.'-'.$individual_id}++; #necessary to have in the key the population, since some individuals are shared between
#populations
} } #reads one line from the compress_genotypes table, uncompress the data, and writes it to the different hashes: one containing the number of bases for the variation and the other with the actual genotype information we need to print in the file
}
_objs_from_sthdescriptionprevnextTop
sub _objs_from_sth {
  my $self = shift;
  return $self->_objs_from_sth_temp_file( @_ );
  my $sth  = shift;
  my $slice = shift;
  my $siblings = shift;

  my ($sample_id,$ld_region_id,$ld_region_start,$ld_region_end,$d_prime,$r2,$sample_count);
  my ($vf_id1,$vf_id2);

  my %feature_container = ();
  my %vf_objects = ();

  #get all Variation Features in Slice
my $vfa = $self->db->get_VariationFeatureAdaptor(); my $variations = $vfa->fetch_all_by_Slice($slice); #retrieve all variation features
#create a hash that maps the position->vf_id
my %pos_vf = (); my $region_Slice = $slice->seq_region_Slice(); map {$pos_vf{$_->seq_region_start} = $_->transfer($region_Slice)} @{$variations}; my %alleles_variation = (); #will contain a record of the alleles in the variation. A will be the major, and a the minor. When more than 2 alleles
#, the genotypes for that variation will be discarded
my %individual_information = (); #to store the relation of snps->individuals
my %regions; #will contain all the regions in the population and the number of genotypes in each one
my $previous_seq_region_id = 0; my %_pop_ids; my ($individual_id, $seq_region_id, $seq_region_start,$seq_region_end,$genotypes, $population_id); my $bin = $self->executable; #open the pipe between processes if the binary file exists in the PATH
if( ! $bin ) { warning("Binary file calc_genotypes not found. Please, read the ensembl-variation/C_code/README.txt file if you want to use LD calculation\n"); goto OUT; } my $pid; eval "require IPC::Run"; #check wether the IPC::Run module it is installed
if( $@ ){ warning("IPC::Run it is not installed in you system. Please, read ensembl-variation/C_code/README.txt if you want to use the LD calculation"); goto OUT; } else{ use IPC::Run qw(start finish); $pid = start [$bin], '<pipe',\* IN, '>pipe',\* OUT, '2>pipe',\* ERR || die "returned $?" ; } #set autoflush
my $piid = fork; if (!defined $piid){ throw("Not possible to fork: $!\n"); } elsif ($piid !=0){ close IN || die "Could not close writer filehandle: $!\n"; #you are the father, read from the pipe
while(<OUT>){ my %ld_values = (); # 936 965891 164284 166818 0.628094 0.999996 120
#get the ouput into the hashes
chomp; ($sample_id,$ld_region_id,$ld_region_start,$ld_region_end,$r2,$d_prime,$sample_count) = split /\s/; $ld_values{'d_prime'} = $d_prime; $ld_values{'r2'} = $r2; $ld_values{'sample_count'} = $sample_count; $vf_id1 = $pos_vf{$ld_region_start}->dbID(); $vf_id2 = $pos_vf{$ld_region_end}->dbID(); $feature_container{$vf_id1 . '-' . $vf_id2}->{$sample_id} =\% ld_values; $vf_objects{$vf_id1} = $pos_vf{$ld_region_start}; $vf_objects{$vf_id2} = $pos_vf{$ld_region_end}; $_pop_ids{$sample_id} = 1; } waitpid($piid,0); close OUT || die "Could not close filehandle: $!\n"; finish $pid || die "Could not finish fork..: $!\n"; } else{ close OUT || die "Could not close filehandle: $!\n"; #the parent dumps the data from the database and writes in the pipe
$sth->bind_columns(\$individual_id,\$ seq_region_id,\$ seq_region_start,\$ seq_region_end,\$ genotypes,\$ population_id); while($sth->fetch()) { #only print genotypes without parents genotyped
if (!exists $siblings->{$population_id . '-' . $individual_id}){ #necessary to use the population_id
$self->_store_genotype(\%individual_information,\%alleles_variation, $individual_id, $seq_region_start, $genotypes, $population_id, $slice); $previous_seq_region_id = $seq_region_id; } } $sth->finish(); select(IN); #necessary to flush the pipe
$| = 1; #we have to print the variations
foreach my $snp_start (sort{$a<=>$b} keys %alleles_variation){ foreach my $population (keys %{$alleles_variation{$snp_start}}){ #if the variation has 2 alleles, print all the genotypes to the file
if (keys %{$alleles_variation{$snp_start}{$population}} == 2){ $self->_convert_genotype($alleles_variation{$snp_start}{$population},$individual_information{$population}{$snp_start}); foreach my $individual_id (keys %{$individual_information{$population}{$snp_start}}){ print IN join("\t",$previous_seq_region_id,$snp_start, $snp_start, $population, $individual_id, $individual_information{$population}{$snp_start}{$individual_id}{genotype})."\n" || warn $!; } } } } close IN || die "Could not close filehandle: $! and $?\n"; POSIX:_exit(0); } OUT: my $t = Bio::EnsEMBL::Variation::LDFeatureContainer->new( '-ldContainer'=>\% feature_container, '-name' => '', '-variationFeatures' =>\% vf_objects ); $t->{'_pop_ids'} =\%_pop_ids; return $t; } #for a given population, gets all individuals that are children (have father or mother)
}
_objs_from_sth_temp_filedescriptionprevnextTop
sub _objs_from_sth_temp_file {
  my $self = shift;
  my $sth  = shift;
  my $slice = shift;
  my $siblings = shift;

  my ($sample_id,$ld_region_id,$ld_region_start,$ld_region_end,$d_prime,$r2,$sample_count);
  my ($vf_id1,$vf_id2);

  my %feature_container = ();
  my %vf_objects = ();

  #get all Variation Features in Slice
my $vfa = $self->db->get_VariationFeatureAdaptor(); my $variations = $vfa->fetch_all_by_Slice($slice); #retrieve all variation features
#create a hash that maps the position->vf_id
my %pos_vf = (); my $region_Slice = $slice->seq_region_Slice(); map {$pos_vf{$_->seq_region_start} = $_->transfer($region_Slice)} @{$variations}; my %alleles_variation = (); #will contain a record of the alleles in the variation. A will be the major, and a the minor. When more than 2 alleles
#, the genotypes for that variation will be discarded
my %individual_information = (); #to store the relation of snps->individuals
my %regions; #will contain all the regions in the population and the number of genotypes in each one
my $previous_seq_region_id = 0; my %_pop_ids; my ($individual_id, $seq_region_id, $seq_region_start,$seq_region_end,$genotypes, $population_id); my @cmd = qw(calc_genotypes); #open the pipe between processes if the binary file exists in the PATH
my $bin = $self->executable; #open the pipe between processes if the binary file exists in the PATH
if( ! $bin ) { warning("Binary file calc_genotypes not found. Please, read the ensembl-variation/C_code/README.txt file if you want to use LD calculation\n"); goto OUT; } my $pid; my $file= $self->temp_path."/".sprintf( "ld%08x%08x%08x", $$, time, rand( 0x7fffffff) ); open IN, ">$file.in"; $sth->bind_columns(\$individual_id,\$ seq_region_id,\$ seq_region_start,\$ seq_region_end,\$ genotypes,\$ population_id); while($sth->fetch()) { #only print genotypes without parents genotyped
if (!exists $siblings->{$population_id . '-' . $individual_id}){ #necessary to use the population_id
$self->_store_genotype(\%individual_information,\%alleles_variation, $individual_id, $seq_region_start, $genotypes, $population_id, $slice); $previous_seq_region_id = $seq_region_id; } } $sth->finish(); #we have to print the variations
foreach my $snp_start (sort{$a<=>$b} keys %alleles_variation){ foreach my $population (keys %{$alleles_variation{$snp_start}}){ #if the variation has 2 alleles, print all the genotypes to the file
if (keys %{$alleles_variation{$snp_start}{$population}} == 2){ $self->_convert_genotype($alleles_variation{$snp_start}{$population},$individual_information{$population}{$snp_start}); foreach my $individual_id (keys %{$individual_information{$population}{$snp_start}}){ print IN join("\t",$previous_seq_region_id,$snp_start, $snp_start, $population, $individual_id, $individual_information{$population}{$snp_start}{$individual_id}{genotype})."\n" || warn $!; } } } } close IN; # `calc_genotypes <$IN >$OUT`;
`$bin <$file.in >$file.out`; open OUT, "$file.out"; while(<OUT>){ my %ld_values = (); # 936 965891 164284 166818 0.628094 0.999996 120
#get the ouput into the hashes
chomp; ($sample_id,$ld_region_id,$ld_region_start,$ld_region_end,$r2,$d_prime,$sample_count) = split /\s/; $ld_values{'d_prime'} = $d_prime; $ld_values{'r2'} = $r2; $ld_values{'sample_count'} = $sample_count; if (!defined $pos_vf{$ld_region_start} || !defined $pos_vf{$ld_region_end}){ next; #problem to fix in the compressed genotype table: some of the positions seem to be wrong
} $vf_id1 = $pos_vf{$ld_region_start}->dbID(); $vf_id2 = $pos_vf{$ld_region_end}->dbID(); $feature_container{$vf_id1 . '-' . $vf_id2}->{$sample_id} =\% ld_values; $vf_objects{$vf_id1} = $pos_vf{$ld_region_start}; $vf_objects{$vf_id2} = $pos_vf{$ld_region_end}; $_pop_ids{$sample_id} = 1; } close OUT || die "Could not close filehandle: $!\n"; unlink( "$file.in" ); unlink( "$file.out" ); OUT: my $t = Bio::EnsEMBL::Variation::LDFeatureContainer->new( '-ldContainer'=>\% feature_container, '-name' => '', '-variationFeatures' =>\% vf_objects ); $t->{'_pop_ids'} =\%_pop_ids; return $t; } 1;
}
_store_genotypedescriptionprevnextTop
sub _store_genotype {
    my $self = shift;
    my $individual_information = shift;
    my $alleles_variation = shift;
    my $individual_id = shift;
    my $seq_region_start = shift;
    my $genotype = shift;
    my $population_id = shift;
    my $slice = shift;

    #get the first byte of the string, and unpack it (the genotype, without the gaps)
my $blob = substr($genotype,2); #the array contains the uncompressed value of genotype, always in the format number_gaps . genotype
my @genotypes = unpack("naa" x (length($blob)/4),$blob);
unshift @genotypes, substr($genotype,1,1); #add the second allele of the first genotype
unshift @genotypes, substr($genotype,0,1); #add the first allele of the first genotype
unshift @genotypes, 0; #the first SNP is in the position indicated by the seq_region1
my $snp_start; my $allele_1; my $allele_2; for (my $i=0; $i < @genotypes -1;$i+=3){ #number of gaps
if ($i == 0){ $snp_start = $seq_region_start; #first SNP is in the beginning of the region
} else{ #ignore when there is more than 1 genotype in the same position
if ($genotypes[$i] == 0){ $snp_start += 1; next; } $snp_start += $genotypes[$i] +1; } #genotype
$allele_1 = $genotypes[$i+1]; $allele_2 = $genotypes[$i+2]; #only get genotypes in the range
if (($snp_start >= $slice->start) && ($snp_start <= $slice->end)){ #store in structure
if ($allele_1 ne 'N' and $allele_2 ne 'N'){ $alleles_variation->{$snp_start}->{$population_id}->{$allele_1}++; $alleles_variation->{$snp_start}->{$population_id}->{$allele_2}++; $individual_information->{$population_id}->{$snp_start}->{$individual_id}->{allele_1} = $allele_1; $individual_information->{$population_id}->{$snp_start}->{$individual_id}->{allele_2} = $allele_2; } } } } #
# Converts the genotype into the required format for the calculation of the pairwise_ld value: AA, Aa or aa
# From the Allele table, will select the alleles and compare to the alleles in the genotype
#
}
executabledescriptionprevnextTop
sub executable {
  my $self = shift;
  $BINARY_FILE = shift if @_;
  unless( $BINARY_FILE ) {
    my $binary_name = 'calc_genotypes';
    ($BINARY_FILE) = grep {-e $_} map {"$_/calc_genotypes"} split /:/,$ENV{'PATH'};
  }
  return $BINARY_FILE;
}
fetch_by_SlicedescriptionprevnextTop
sub fetch_by_Slice {
    my $self = shift;
    my $slice = shift;
    my $population = shift;

    if(!ref($slice) || !$slice->isa('Bio::EnsEMBL::Slice')) {
	throw('Bio::EnsEMBL::Slice arg expected');
    }
    my $sth;
    my $in_str;
    my $siblings = {};
    #when there is no population selected, return LD in the HapMap and PerlEgen populations
$in_str = $self->_get_LD_populations($siblings); #if a population is passed as an argument, select the LD in the region with the population
if ($population){ if(!ref($population) || !$population->isa('Bio::EnsEMBL::Variation::Population')) { throw('Bio::EnsEMBL::Variation::Population arg expected'); } my $population_id = $population->dbID; $in_str = " = $population_id"; # if ($in_str =~ /$population_id/){
# $in_str = "IN ($population_id)";
#' }
# else{
# warning("Not possible to calculate LD for a non HapMap or PerlEgen population: $population_id");
# return {};
# }
} if ($in_str eq ''){ #there is no population, not a human specie or not passed as an argument, return the empy container
my $t = Bio::EnsEMBL::Variation::LDFeatureContainer->new( '-ldContainer'=> {}, '-name' => $slice->name, '-variationFeatures' => {} ); return $t } $sth = $self->prepare(qq{SELECT c.sample_id,c.seq_region_id,c.seq_region_start,c.seq_region_end,c.genotypes,ip.population_sample_id FROM compressed_genotype_single_bp c, individual_population ip WHERE ip.individual_sample_id = c.sample_id AND ip.population_sample_id $in_str AND c.seq_region_id = ? AND c.seq_region_start >= ? and c.seq_region_start <= ? AND c.seq_region_end >= ? ORDER BY c.seq_region_id, c.seq_region_start},{mysql_use_result => 1}); $sth->bind_param(1,$slice->get_seq_region_id,SQL_INTEGER); $sth->bind_param(2,$slice->start - MAX_SNP_DISTANCE,SQL_INTEGER) if ($slice->start - MAX_SNP_DISTANCE >= 1); $sth->bind_param(2,1,SQL_INTEGER) if ($slice->start - MAX_SNP_DISTANCE < 1); $sth->bind_param(3,$slice->end,SQL_INTEGER); $sth->bind_param(4,$slice->start,SQL_INTEGER); $sth->execute(); my $ldFeatureContainer = $self->_objs_from_sth($sth,$slice,$siblings); $sth->finish(); #and store the name of the slice in the Container
$ldFeatureContainer->name($slice->name()); return $ldFeatureContainer;
}
fetch_by_VariationFeaturedescriptionprevnextTop
sub fetch_by_VariationFeature {
  my $self = shift;
  my $vf  = shift;
  my $pop = shift;

  if(!ref($vf) || !$vf->isa('Bio::EnsEMBL::Variation::VariationFeature')) {
    throw('Bio::EnsEMBL::Variation::VariationFeature arg expected');
  }

  if(!defined($vf->dbID())) {
    throw("VariationFeature arg must have defined dbID");
  }
  my $ldFeatureContainer = $self->fetch_by_Slice($vf->feature_Slice->expand(MAX_SNP_DISTANCE,MAX_SNP_DISTANCE),$pop);
  #we need to filter and only return LD with the vf associated
my %feature_container = (); my %vf_objects = (); my %pop_ids = (); my $vf_id = $vf->dbID; #iterate all the LD container looking for values related to this variation Feature
foreach my $ld_key (keys %{$ldFeatureContainer->{'ldContainer'}}){ if ($ld_key =~ /$vf_id/){ #add the value to the new hash
my ($vf_id1,$vf_id2) = split /-/,$ld_key; $feature_container{$vf_id1 . '-' . $vf_id2} = $ldFeatureContainer->{'ldContainer'}->{$ld_key}; $vf_objects{$vf_id1} = $ldFeatureContainer->{'variationFeatures'}->{$vf_id1}; $vf_objects{$vf_id2} = $ldFeatureContainer->{'variationFeatures'}->{$vf_id2}; map {$pop_ids{$_}++} keys %{$ldFeatureContainer->{'ldContainer'}->{$ld_key}}; } } my $new_ldFeatureContainer = Bio::EnsEMBL::Variation::LDFeatureContainer->new( '-ldContainer'=>\% feature_container, '-name' => $vf->dbID, '-variationFeatures' =>\% vf_objects ); $new_ldFeatureContainer->{'_pop_ids'} =\% pop_ids; return $new_ldFeatureContainer; } #
# private method, creates ldfeatureContainer objects from an executed statement handle
# ordering of columns must be consistant
#
}
temp_pathdescriptionprevnextTop
sub temp_path {
  my $self = shift;
  $TMP_PATH = shift if @_;
  $TMP_PATH ||= '/tmp'; 
  return $TMP_PATH;
}
General documentation
AUTHOR - Daniel RiosTop
CONTACTTop
Post questions to the Ensembl development list ensembl-dev@ebi.ac.uk