This module is an entry point into a Glovar database. It allows you to retrieve
haplotype data from the database.
sub fetch_all_by_clone_accession
{ my ($self, $embl_acc, $embl_version, $cl_start, $cl_end) = @_;
my $dnadb = Bio::EnsEMBL::Registry->get_DNAAdaptor($ENV{'ENSEMBL_SPECIES'}, 'glovar');
unless ($dnadb) {
warning("ERROR: No dnadb attached to Glovar.\n");
return([]);
}
my @cloneinfo = $self->fetch_clone_by_accession($embl_acc);
return([]) unless (@cloneinfo);
my ($nt_name, $id_seq, $clone_start, $clone_end, $clone_strand) = @cloneinfo;
my ($q_start, $q_end);
if ($clone_strand == 1) {
$q_start = $clone_start + $cl_start - 1;
$q_end = $clone_start + $cl_end + 1;
} else{
$q_start = $clone_end - $cl_end - 1;
$q_end = $clone_end - $cl_start + 1;
}
my $clone_length = $clone_end - $clone_start + 1;
my $q2 = qq(
SELECT
sb.id_block as block,
ss.is_private as private_snp,
bs.is_private as private_block
FROM
mapped_snp ms,
snp_summary ss,
snp_block sb,
block b,
block_set bs
WHERE ms.id_sequence = ?
AND ms.id_snp = ss.id_snp
AND ss.id_snp = sb.id_snp
AND sb.id_block = b.id_block
AND b.id_block_set = bs.id_block_set
AND ms.position BETWEEN $q_start AND $q_end
);
my $sth;
eval {
$sth = $self->prepare($q2);
$sth->execute($id_seq);
};
if ($@){
warning("ERROR: SQL failed in " . (caller(0))[3] . "\n$@");
return([]);
}
my @blocks;
while (my $b = $sth->fetchrow_hashref()) {
return([]) unless keys %{$b};
warning("WARNING: private haplotype/SNP!") if ($b->{'PRIVATE_SNP'} ||
$b->{'PRIVATE_BLOCK'});
push @blocks, $b->{'BLOCK'};
}
return([]) unless (@blocks);
my $blocklist = join(",", @blocks);
my $q3 = qq(
SELECT
ms.position as snp_start,
ms.end_position as snp_end,
ms.is_revcomp as snp_strand,
sb.id_block as internal_id,
b.name as block_name,
b.length as block_length,
b.num_snps as num_snps,
p.description as population,
ss.is_private as private_snp,
bs.is_private as private_block
FROM
mapped_snp ms,
snp_summary ss,
snp_block sb,
block b,
block_set bs,
population p
WHERE ms.id_sequence = ?
AND sb.id_block IN ($blocklist)
AND ms.id_snp = ss.id_snp
AND ss.id_snp = sb.id_snp
AND sb.id_block = b.id_block
AND b.id_block_set = bs.id_block_set
AND bs.id_pop = p.id_pop
);
eval {
$sth = $self->prepare($q3);
$sth->execute($id_seq);
};
if ($@){
warning("ERROR: SQL failed in " . (caller(0))[3] . "\n$@");
return([]);
}
my $clone = $dnadb->get_SliceAdaptor->fetch_by_region('clone', $embl_acc);
my $projected_clone = $clone->project('chromosome');
unless ($projected_clone) {
warning("Clone $embl_acc doesn't project to chromosome.");
return([]);
}
my $haplotypes;
while (my $row = $sth->fetchrow_hashref()) {
return([]) unless keys %{$row};
warn "WARNING: private haplotype/SNP!" if ($row->{'PRIVATE_SNP'} ||
$row->{'PRIVATE_BLOCK'});
next unless $row->{'SNP_STRAND'};
my ($start, $end);
$row->{'SNP_END'} ||= $row->{'SNP_START'};
if ($clone_strand == 1) {
$start = $row->{'SNP_START'} - $clone_start + 1;
$end = $row->{'SNP_END'} - $clone_start + 1;
} else {
$start = $clone_end - $row->{'SNP_END'} + 1;
$end = $clone_end - $row->{'SNP_START'} + 1;
}
my $strand = $row->{'SNP_STRAND'}*$clone_strand;
$haplotypes->{$row->{'INTERNAL_ID'}} ||=
Bio::EnsEMBL::ExternalData::Glovar::Haplotype->new_fast({
'analysis' => 'glovar_haplotype',
'display_id' => $row->{'BLOCK_NAME'},
'dbID' => $row->{'INTERNAL_ID'},
'start' => $start,
'end' => $end,
'strand' => $clone_strand,
'seqname' => $embl_acc,
'population' => $row->{'POPULATION'},
'num_snps' => $row->{'NUM_SNPS'},
});
my $hap = $haplotypes->{$row->{'INTERNAL_ID'}};
$hap->start($start) if ($hap->start > $start);
$hap->end($end) if ($hap->end < $end);
my $snp = Bio::EnsEMBL::Feature->new(
-SLICE => $clone,
-START => $start,
-END => $end,
-STRAND => $strand,
);
$snp = $snp->transform('chromosome');
next unless $snp;
$hap->add_tagSNP($snp->start, $snp->end);
}
my @haps = ();
foreach my $key (keys %$haplotypes) {
if ($haplotypes->{$key}->start < $projected_clone->[0]->from_start) {
$haplotypes->{$key}->start($projected_clone->[0]->from_start);
}
if ($haplotypes->{$key}->end > $projected_clone->[0]->from_end) {
$haplotypes->{$key}->end($projected_clone->[0]->from_end);
}
push @haps, $haplotypes->{$key};
}
return\@ haps; } |