[ensembl-dev] unable to recover high-r2 pair using fetch_by_VariationFeatures

[andrew126 at mac.com]

Thanks!

> On Apr 24, 2017, at 5:18 AM, Anja Thormann <anja at ebi.ac.uk> wrote:
> 
> Hi Andrew,
> 
> both is now fixed on the release/88 branch. Please refresh the code.
> 
> Thank you,
> Anja
> 
>> On 20 Apr 2017, at 17:30, andrew126 at mac.com wrote:
>> 
>> Hi Anja,
>> 
>> Thanks for the response.
>> 
>> I’ve installed API 88 and things seem to work so far (the expected r2 value is returned .. smile), with the following observations:
>> 
>> 1)
>> Running my code now generates looooots of these messages (regardless of if I use remote or local vcfs):
>> 
>> 	Use of uninitialized value in hash element at /home/andrew/gte/src_88/ensembl-variation/modules/Bio/EnsEMBL/Variation/VCFCollection.pm line 488.
>> 
>> This is the offending line:
>> 	$synonyms{$syn->[0]} = $sample->name;
>> 
>> I presume these are benign messages?
>> 
>> 2)
>> ensemble-variation/C_code/ld_vcf.c wanted:
>> 	return 0;
>> at the end of main().
>> 
>> Thanks for the help.
>> 
>> Best,
>> 
>> Andrew
>> 
>> 
>>> On Apr 20, 2017, at 1:51 AM, Anja Thormann <anja at ebi.ac.uk> wrote:
>>> 
>>> 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/
>>> 
>>> _______________________________________________
>>> 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/
>> 
>> _______________________________________________
>> 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/
> 
> _______________________________________________
> 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/