VariationFeature coordinate test

#!/usr/bin/perl

use strict;
use warnings;

use Bio::EnsEMBL::Registry;

# example variant rs1800460 is protein coding and corresponds to
# NG_012137.1:g.26260C>G
# NM_000367.2:c.500C>G,
# NP_000358.1:p.Ala167Gly.
# ENST00000309983.4:c.500C>G
# ENSP00000312304.4:p.Ala167Gly
# http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=rs74423290;source=dbSNP
# http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=74423290

my $rs = 'rs74423290';
my $id = 'ENST00000309983';
my $cds_pos = 500;
my $allele_string = 'C/G';

my %ens = ensembl_init();

my $rs_v = $ens{va}->fetch_by_name($rs);
my ($rs_vf) = map {@{$ens{vfa}->fetch_all_by_Variation($_)}} $rs_v;
print_vf($rs_vf);

my ($tx) = $ens{ta}->fetch_by_stable_id($id);
my ($gstart) = $tx->cdna2genomic($cds_pos + $tx->cdna_coding_start,
                                 $cds_pos + $tx->cdna_coding_start);
my $vf = Bio::EnsEMBL::Variation::VariationFeature->new(
  -variation_name => 'test',
  -strand => 1,
  -start => $gstart->start,
  -end => $gstart->start,
  -slice => $tx->slice,
  -allele_string => $allele_string,
  -adapter => $ens{vfa}
 );
print_vf($vf);


sub print_vf {
  my ($vf) = @_;
  printf("* %s @%s:%s\n",
         $vf->variation_name,
         $vf->seq_region_name,
         $vf->seq_region_start);
}


sub ensembl_init {
  my %ens;
  $ens{registry} = 'Bio::EnsEMBL::Registry';
  $ens{registry}->load_registry_from_db(
    -host => 'outcast.locusdev.net',
    -user => 'reece',
   );
  $ens{sa} = $ens{registry}->get_adaptor( 'Human', 'Core', 'Slice' );
  $ens{ta} = $ens{registry}->get_adaptor( 'Human', 'Core', 'Transcript' );
  $ens{va} = $ens{registry}->get_adaptor( 'Human', 'Variation', 'Variation' );
  $ens{vfa} = $ens{registry}->get_adaptor( 'Human', 'Variation', 'VariationFeature' );

  return %ens;
}