Toolbar

Summary

GlovarHaplotypeAdaptor -
Database adaptor for Glovar haplotypes

Package variables

Included modules

Inherit

Bio::EnsEMBL::ExternalData::Glovar::GlovarAdaptor

Synopsis

$glodb = Bio::EnsEMBL::ExternalData::Glovar::DBAdaptor->new(
-user => 'ensro',
-pass => 'secret',
-dbname => 'snp',
-host => 'go_host',
-driver => 'Oracle');
my $glovar_adaptor = $glodb->get_GlovarHaplotypeAdaptor;
$var_listref = $glovar_adaptor->fetch_all_by_clone_accession(
'AL100005', 'AL100005', 1, 10000);

Description

This module is an entry point into a Glovar database. It allows you to retrieve
haplotype data from the database.

Methods

Methods description

  Arg[1]      : none
  Example     : my @coord_systems = $glovar_adaptor->coordinate_systems;
  Description : This method returns a list of coordinate systems which are
                implemented by this class. A minimum of one valid coordinate
                system must be implemented. Valid coordinate systems are:
                'SLICE', 'ASSEMBLY', 'CONTIG', and 'CLONE'.
  Return type : list of strings
  Exceptions  : none
  Caller      : internal

  Arg[1]      : clone internal ID
  Arg[2]      : clone embl accession
  Arg[3]      : clone start coordinate
  Arg[4]      : clone end coordinate
  Example     : @list = @{$glovar_adaptor->fetch_all_by_clone_accession(
                    'AL100005', 'AL100005', 1, 10000)};
  Description: Retrieves haplotypes (blocks of SNPs) on a clone in clone
               coordinates. Since Glovar stores now haplotype start/end data,
               the haplotypes have to be retrieved in a two-step process:
               1. get IDs of all haplotypes in the requested region
               2. for all these IDs, get tagSNPs assigned to them
               Haplotype start/end have to be inferred from the tagSNP
               coordinates.
  Returntype : Listref of Bio::EnsEMBL::ExternalData::Glovar::Haplotype objects
  Exceptions : none
  Caller     : $self->fetch_all_by_Clone

  Arg[1]      : none
  Example     : my $track_name = $haplotype_adaptor->track_name;
  Description : returns the track name
  Return type : String - track name
  Exceptions  : none
  Caller      : Bio::EnsEMBL::Slice,
                Bio::EnsEMBL::ExternalData::ExternalFeatureAdaptor

Methods code

sub coordinate_systems {

    return ('CLONE');

}

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([]);
    }
    
    # get info on clone
    my @cloneinfo = $self->fetch_clone_by_accession($embl_acc);
    return([]) unless (@cloneinfo);
    my ($nt_name, $id_seq, $clone_start, $clone_end, $clone_strand) = @cloneinfo;

    # get only features in the desired region of the clone
    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;
    
    # get all haplotype blocks
    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);

    # now get all tagSNPs for the haplotypes found
    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([]);
    }

    # get a clone slice and projectc it to chromosomal coordinates
    # this is needed for tag SNP and haplotype coordinate calculations
    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'});

        ## filter SNPs without strand (this is gruft in mapped_snp)
        next unless $row->{'SNP_STRAND'};

        ## calculate clone coords depending on clone orientation
        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;

        # create a haplotype object if it doesn't already exist
        $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'}};

        # record haplotype start/end, ie max(snp_start)/max(snp_end)
        $hap->start($start) if ($hap->start > $start);
        $hap->end($end) if ($hap->end < $end);

        # calculate SNP chromosomal coords
        my $snp = Bio::EnsEMBL::Feature->new(
                    -SLICE  => $clone,
                    -START  => $start,
                    -END    => $end,
                    -STRAND => $strand,
        );
        $snp = $snp->transform('chromosome');
        # skip SNPs that don't transform
        next unless $snp;
        
        # push tag SNP coords on haplotype
        $hap->add_tagSNP($snp->start, $snp->end);
    }
    
    # return listref of haplotype objects
    my @haps = ();
    foreach my $key (keys %$haplotypes) {
        # trim haplotype to golden portion of clone (otherwise it will be 
        # discarded during the mapping to chromosome)
        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;

}

sub track_name {

    my ($self) = @_;    
    return("GlovarHaplotype");
}

1;

}

General documentation

Patrick Meidl <pm2@sanger.ac.uk>

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