[ensembl-dev] Error/warning when using get_r_square of LDFeatureContainer

Johanne Håøy Horn johannhh at ifi.uio.no
Mon Feb 29 09:00:19 GMT 2016


Dear Ensembl dev team,

I try to get rsquare between two SNPs (on the same chromosome). When I use the ldFeatureContainer function get_r_square, however, I get a warning and no return rsquare value. Here is the function signature from the documentation:


$r_square = $obj->get_r_square($vf1,$vf2,$population_id);

The script I am using is the following:

# START SCRIPT
use strict;
use warnings;
use Bio::EnsEMBL::Registry;

my $start_run = time();

my $registry = 'Bio::EnsEMBL::Registry';

$registry->load_registry_from_db(
-host => 'ensembldb.ensembl.org<http://ensembldb.ensembl.org>',
-user => 'anonymous'
);

# Connect to the databases:
my $variation_adaptor = $registry->get_adaptor('homo_sapiens', 'variation', 'variation');
my $ldfc_adaptor = $registry->get_adaptor('homo_sapiens', 'variation', 'ldfeaturecontainer');
my $pop_adaptor = $registry->get_adaptor('homo_sapiens', 'variation', 'population');
$variation_adaptor->db->use_vcf(1);

my $ld_population = $pop_adaptor->fetch_by_name('1000GENOMES:phase_3:CEU');

# Build tagSNP list of SNPs on same chromosome:
my $tagSNP1 = 'rs6678275';
my $tagSNP2 = 'rs6656401';
my @variations;
push (@variations, $variation_adaptor->fetch_by_name($tagSNP1));
push (@variations, $variation_adaptor->fetch_by_name($tagSNP2));

my @all_features;
foreach my $variation (@variations) {
my @var_features = @{ $variation->get_all_VariationFeatures() };
foreach my $vf (@var_features) {
push (@all_features, $vf);
}
}

my $ldfc = $ldfc_adaptor->fetch_by_VariationFeatures(\@all_features, $ld_population);

# Print r2 between all SNPs:
foreach my $vf (@all_features) {
foreach my $vf2 (@all_features) {
my $rsid = $vf->name;
my $rsid2 = $vf2->name;
my $start = $vf->start;
my $start2 = $vf2->start;
my $region = $vf->seq_region_name;
my $region2 = $vf2->seq_region_name;

print "$region\t$rsid\t$start\n";
print "$region2\t$rsid2\t$start2\n";

my $r2 = $ldfc->get_r_square($vf, $vf2);
#print "\tr2=$r2\n";
}
}
# END SCRIPT

I get this printout:

$ perl rsquareTagSNPs.pl out_tagsnps.txt 0.9
1 rs6678275 193656103
1 rs6678275 193656103

-------------------- WARNING ----------------------
MSG: variation features have no pairwise ld information
FILE: EnsEMBL/Variation/LDFeatureContainer.pm LINE: 207
CALLED BY: rsquareTagSNPs.pl  LINE: 59
Date (localtime)    = Mon Feb 29 09:47:44 2016
Ensembl API version = 83
---------------------------------------------------
1 rs6678275 193656103
1 rs6656401 207518704

-------------------- WARNING ----------------------
MSG: variation features have no pairwise ld information
FILE: EnsEMBL/Variation/LDFeatureContainer.pm LINE: 207
CALLED BY: rsquareTagSNPs.pl  LINE: 59
Date (localtime)    = Mon Feb 29 09:47:44 2016
Ensembl API version = 83
---------------------------------------------------
1 rs6656401 207518704
1 rs6678275 193656103

-------------------- WARNING ----------------------
MSG: variation features have no pairwise ld information
FILE: EnsEMBL/Variation/LDFeatureContainer.pm LINE: 207
CALLED BY: rsquareTagSNPs.pl  LINE: 59
Date (localtime)    = Mon Feb 29 09:47:44 2016
Ensembl API version = 83
---------------------------------------------------
1 rs6656401 207518704
1 rs6656401 207518704

-------------------- WARNING ----------------------
MSG: variation features have no pairwise ld information
FILE: EnsEMBL/Variation/LDFeatureContainer.pm LINE: 207
CALLED BY: rsquareTagSNPs.pl  LINE: 59
Date (localtime)    = Mon Feb 29 09:47:44 2016
Ensembl API version = 83
---------------------------------------------------

Am I using the function on the wrong kind of ldfc? Or have I misunderstood the documentation? The ultimate goal of my script is to find r2 between all SNPs/rsIDs I pass to the script which are on the same chromosome.

Best,
Johanne Håøy Horn

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.ensembl.org/pipermail/dev_ensembl.org/attachments/20160229/b6fcb1b9/attachment.html>


More information about the Dev mailing list