[ensembl-dev] unable to recover high-r2 pair using fetch_by_VariationFeatures
Anja Thormann
anja at ebi.ac.uk
Thu Apr 20 09:51:02 BST 2017
Hi Andrew,
we had problems with fetch_by_VariationFeatures for that release. Would it be possible for you to use the most recent version 88 of our API?
Thank you,
Anja
> On 17 Apr 2017, at 05:56, andrew126 at mac.com wrote:
>
> Hi,
>
> I am using API version 86.
>
> There are two snps within 1kb of each other that have a high r2 value:
>
> rs79796976 (chr1, pos 247,390,677)
> rs56285508 (chr1, pos 247,391,186)
> r^2 0.983802 in 1000GENOMES:phase_3:EUR
>
> I can recover rs56285508 as a high-LD pair when I use fetch_by_VariationFeature (no “s” at the end), querying with rs79796976.
>
> However, I get no result when I use fetch_by_VariationFeatures (with an “s” at the end), querying with both snps.
>
> Below is the perl script using fetch_by_VariationFeatures.
>
> Can you please let me know if I am using fetch_by_VarationFeatures incorrectly?
>
> Thanks for any guidance.
>
> Best regards,
>
> Andrew
>
>
>
> use strict;
> $|=1;
> use IPC::Run;
> use Bio::EnsEMBL::Registry;
> use Bio::EnsEMBL::ApiVersion;
> printf( "The API version used is %s\n", software_version() );
>
> my $registry = 'Bio::EnsEMBL::Registry';
> $registry->load_all();
> $registry->set_reconnect_when_lost(1);
> $registry->load_registry_from_db(
> -host => 'ensembldb.ensembl.org', # alternatively ensembldb 'useastdb.ensembl.org'
> -user => 'anonymous'
> );
>
> my $chrom = "1";
> my $var1 = "rs79796976";
> my $region_start1 = "247390677";
> my $region_end1 = "247390677";
> my $var2 = "rs56285508";
> my $region_start2 = "247391186";
> my $region_end2 = "247391186";
>
> my $pop_adaptor = $registry->get_adaptor("human","variation","population");
> $pop_adaptor->db->use_vcf(1);
> my $ldFeatureContainerAdaptor = $registry->get_adaptor('human', 'variation', 'ldfeaturecontainer');
> $ldFeatureContainerAdaptor->db->use_vcf(1);
> my $variation_adaptor = $registry->get_adaptor( 'human', 'variation', 'variation' );
> $variation_adaptor->db->use_vcf(1);
>
> my $population_name = "1000GENOMES:phase_3:EUR";
> my $population = $pop_adaptor->fetch_by_name($population_name);
>
> my $variation1 = $variation_adaptor->fetch_by_name("$var1");
> my $vfref1 = $variation1->get_all_VariationFeatures();
> my $query_variation_feature1;
> my $found=0;
> foreach my $vf1 (@{$vfref1}) {
> if ($chrom eq $vf1->seq_region_name) {
> if ($region_start1 eq $vf1->seq_region_start && $region_end1 eq $vf1->seq_region_end) {
> $found=1;
> $query_variation_feature1 = $vf1;
> }
> }
> }
> if ($found==0) {die;}
>
> my $variation2 = $variation_adaptor->fetch_by_name("$var2");
> my $vfref2 = $variation2->get_all_VariationFeatures();
> my $query_variation_feature2;
> $found=0;
> foreach my $vf2 (@{$vfref2}) {
> if ($chrom eq $vf2->seq_region_name) {
> if ($region_start2 eq $vf2->seq_region_start && $region_end2 eq $vf2->seq_region_end) {
> $found=1;
> $query_variation_feature2 = $vf2;
> }
> }
> }
> if ($found==0) {die;}
>
> my $ldFeatureContainer = $ldFeatureContainerAdaptor->fetch_by_VariationFeatures([$query_variation_feature1, $query_variation_feature2],$population);
>
> foreach my $r_square (@{$ldFeatureContainer->get_all_r_square_values}){
> if ($r_square->{r2} >= 0.1){
> print $r_square->{variation1}->variation_name."\n";
> print "\t".$r_square->{variation1}->seq_region_start."\n";
> print "\t".$r_square->{variation1}->seq_region_end."\n";
> print $r_square->{variation2}->variation_name."\n";
> print "\t".$r_square->{variation2}->seq_region_start."\n";
> print "\t".$r_square->{variation2}->seq_region_end."\n";
> print "r^2 ".$r_square->{r2}."\n";
> }
> }
>
>
>
> _______________________________________________
> Dev mailing list Dev at ensembl.org
> Posting guidelines and subscribe/unsubscribe info: http://lists.ensembl.org/mailman/listinfo/dev
> Ensembl Blog: http://www.ensembl.info/
More information about the Dev
mailing list