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

andrew126 at mac.com andrew126 at mac.com
Thu Apr 20 17:30:18 BST 2017


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/




More information about the Dev mailing list