%%bash
#Create an installation directory
mkdir src
cd src
#Take a clone of BioPerl from GitHub and change to the correct version:
git clone https://github.com/bioperl/bioperl-live.git
cd bioperl-live
git checkout bioperl-release-1-6-1
cd -
#Install the Ensembl git tools to allow easy management of the Ensembl repos.
git clone https://github.com/Ensembl/ensembl-git-tools.git
export PATH=$PWD/ensembl-git-tools/bin:$PATH
#Install the APIs that you need. You can install all the APIs using the git ensembl command:
git ensembl --clone api
/home/dhimmels/Documents/serg/rephetio/construct/gwas-catalog/ensembl/src * Processing 'api' * Working with module 'ensembl' * Cloning from remote 'https://github.com/Ensembl/ensembl.git' * Enabling git hooks * Working with module 'ensembl-compara' * Cloning from remote 'https://github.com/Ensembl/ensembl-compara.git' * Enabling git hooks * Working with module 'ensembl-funcgen' * Cloning from remote 'https://github.com/Ensembl/ensembl-funcgen.git' * Enabling git hooks * Working with module 'ensembl-io' * Cloning from remote 'https://github.com/Ensembl/ensembl-io.git' * Enabling git hooks * Working with module 'ensembl-variation' * Cloning from remote 'https://github.com/Ensembl/ensembl-variation.git' * Enabling git hooks
Cloning into 'bioperl-live'... Note: checking out 'bioperl-release-1-6-1'. You are in 'detached HEAD' state. You can look around, make experimental changes and commit them, and you can discard any commits you make in this state without impacting any branches by performing another checkout. If you want to create a new branch to retain commits you create, you may do so (now or later) by using -b with the checkout command again. Example: git checkout -b new_branch_name HEAD is now at 83dcde3... * bump point version * sync minor doc changes Cloning into 'ensembl-git-tools'... Cloning into 'ensembl'... Cloning into 'ensembl-compara'... Cloning into 'ensembl-funcgen'... Cloning into 'ensembl-io'... Cloning into 'ensembl-variation'...
snps = ['rs1333049', 'rs2248462', 'rs225190', 'rs225212', 'rs2252865', 'rs2252996']
with open('SNP_list_sorted.txt', 'w') as write_file:
write_file.write('\n'.join(sorted(snps)))
%%perl
use lib "src/bioperl-live";
use lib "src/ensembl/modules";
use lib "src/ensembl-compara/modules";
use lib "src/ensembl-variation/modules";
use lib "src/ensembl-funcgen/modules";
use strict;
use warnings;
use Bio::EnsEMBL::Registry;
open (IN, "<SNP_list_sorted.txt");
open (OUT, ">SNPs_in_LD.tsv");
my $reg = 'Bio::EnsEMBL::Registry';
$reg->load_registry_from_db(-host => 'ensembldb.ensembl.org',-user => 'anonymous');
my $va = $reg->get_adaptor( 'human', 'variation', 'variation' );
while (<IN>) {
print $_;
my $v = $va->fetch_by_name($_);
if (!$v) {
print 'unfetched ' . $_;
next;
}
my @vfs = @{ $v->get_all_VariationFeatures };
foreach my $vf (@vfs){
my $ld = $vf->get_all_LD_values;
my @pops = @{ $vf->get_all_LD_Populations };
my @ldvs = @{ $ld->get_variations };
foreach my $pop (@pops) {
if ($pop->name =~ /1000GENOMES/) {
foreach my $ldv (@ldvs) {
if ($ldv->stable_id ne $_) {
my @ldvfs = @{ $ldv->get_all_VariationFeatures };
foreach my $ldvf (@ldvfs) {
my @tvs = @{ $ldvf->get_all_TranscriptVariations };
my $r2 = $ld->get_r_square($vf, $ldvf, $pop);
foreach my $tv (@tvs) {
my $gene = $tv->transcript->get_Gene;
if ($r2) {
print OUT $v->stable_id, "\t", $ldv->stable_id, "\t", $gene->external_name, "\t", $r2, "\t", $pop->name, "\n";
}
}
}
}
}
}
}
}
}
rs1333049 unfetched rs1333049 rs2248462 unfetched rs2248462 rs225190 unfetched rs225190 rs225212 unfetched rs225212 rs2252865 unfetched rs2252865 rs2252996
UNIVERSAL->import is deprecated and will be removed in a future perl at src/bioperl-live/Bio/Tree/TreeFunctionsI.pm line 94. -------------------- WARNING ---------------------- MSG: Binary file calc_genotypes not found. Please, read the ensembl-variation/C_code/README.txt file if you want to use LD calculation FILE: Variation/DBSQL/LDFeatureContainerAdaptor.pm LINE: 674 CALLED BY: Variation/DBSQL/LDFeatureContainerAdaptor.pm LINE: 607 Date (localtime) = Mon Jun 22 18:03:46 2015 Ensembl API version = 80 ---------------------------------------------------