Raw content of Bio::EnsEMBL::Variation::DBSQL::LDFeatureContainerAdaptor
<![CDATA[#
# Ensembl module for Bio::EnsEMBL::Variation::DBSQL::LDFeatureContainerAdaptor
#
# Copyright (c) 2004 Ensembl
#
# You may distribute this module under the same terms as perl itself
#
#
=head1 NAME
Bio::EnsEMBL::Variation::DBSQL::LDFeatureContainerAdaptor
=head1 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();
=head1 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.
=head1 AUTHOR - Daniel Rios
=head1 CONTACT
Post questions to the Ensembl development list ensembl-dev@ebi.ac.uk
=head1 METHODS
=cut
use strict;
use warnings;
package Bio::EnsEMBL::Variation::DBSQL::LDFeatureContainerAdaptor;
use Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor;
use Bio::EnsEMBL::Variation::LDFeatureContainer;
use vars qw(@ISA);
use Data::Dumper;
use POSIX;
use FileHandle;
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use constant MAX_SNP_DISTANCE => 100_000;
use base qw(Bio::EnsEMBL::DBSQL::BaseAdaptor);
our $MAX_SNP_DISTANCE = 100000;
our $BINARY_FILE = '';
our $TMP_PATH = '';
=head2 fetch_by_Slice
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
=cut
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;
}
sub temp_path {
my $self = shift;
$TMP_PATH = shift if @_;
$TMP_PATH ||= '/tmp';
return $TMP_PATH;
}
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;
}
=head2 fetch_by_VariationFeature
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
=cut
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
#
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)
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
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
#
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';
}
}
}
}
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]);
}
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;
]]>