Error while performing "lift-over" of genomic coordinates between two species using Ensembl Compara API

Question

I'm trying to convert coordinates between the reference mouse genome and CAST/EiJ mice strain. Since Ensebml already has alignment available in the browser, I decided to use Ensembl perl API to lift over the alignment for the regions I'm interested in. The problem is that it doesn't seem to be working. I'm getting following error message:

-------------------- EXCEPTION --------------------
MSG: The database do not contain any CACTUS_HAL data for Bio::EnsEMBL::Compara::GenomeDB=HASH(0x557db0d661e8)!
STACK toplevel coord.pl:42
Date (localtime)    = Mon Oct 15 07:42:17 2018
Ensembl API version = 94
---------------------------------------------------

which I assume means that I'm either searching for the wrong alignment record or the alignment itself is not available. Following is the perl script which I modified from the Ensembl Compara API tutorial:

use strict;
use Bio::EnsEMBL::Registry;
use Bio::EnsEMBL::Utils::Exception qw(throw);
use Bio::AlignIO;
use Getopt::Long;

my $species = "murinae";
my $targetSpecies = "mus_musculus_casteij";
my $alignment_type = "CACTUS_HAL";
my $output_file = undef;
my $coord_system = "chromosome";
my $seq_region = "14";
my $seq_region_start = 5000000;
my $seq_region_end = 5010000;
my $output_format = "clustalw";
my $genome_dbs;

use Bio::EnsEMBL::Registry;
Bio::EnsEMBL::Registry->load_registry_from_db(
    -host => 'ensembldb.ensembl.org',
    -user => 'anonymous',
    -port => 5306);

my $genome_db_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'GenomeDB');

throw("Cannot connect to Compara") if (!$genome_db_adaptor);

my $this_genome_db = $genome_db_adaptor->fetch_by_name_assembly($species);
my $target_genome_db = $genome_db_adaptor->fetch_by_name_assembly($targetSpecies);
push(@$genome_dbs, $this_genome_db);
push(@$genome_dbs, $target_genome_db);

my $method_link_species_set_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'MethodLinkSpeciesSet');

my $method_link_species_set =
$method_link_species_set_adaptor->fetch_by_method_link_type_GenomeDBs(
      $alignment_type, 
  $genome_dbs);

throw("The database do not contain any $alignment_type data for $target_genome_db!")
    if (!$method_link_species_set);


my $slice_adaptor = Bio::EnsEMBL::Registry->get_adaptor($species, 'core', 'Slice');

throw("Registry configuration file has no data for connecting to <$species>")
if (!$slice_adaptor);

my $query_slice = $slice_adaptor->fetch_by_region('toplevel', $seq_region, $seq_region_start, $seq_region_end);

throw("No Slice can be created with coordinates $seq_region:$seq_region_start-".
    "$seq_region_end") if (!$query_slice);

# Fetching all the GenomicAlignBlock corresponding to this Slice:
my $genomic_align_block_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'GenomicAlignBlock');

my $genomic_align_blocks =
$genomic_align_block_adaptor->fetch_all_by_MethodLinkSpeciesSet_Slice(
      $method_link_species_set, 
  $query_slice);

my $all_aligns;

# Get a Bio::SimpleAlign object from every GenomicAlignBlock
foreach my $this_genomic_align_block (@$genomic_align_blocks) {
    my $simple_align = $this_genomic_align_block->get_SimpleAlign;
    push(@$all_aligns, $simple_align);
}

# print all the genomic alignments using a Bio::AlignIO object
my $alignIO = Bio::AlignIO->newFh(
-interleaved => 0,
-fh => \*STDOUT,
-format => $output_format,
    -idlength => 10
);

foreach my $this_align (@$all_aligns) {
    print $alignIO $this_align;
}

exit; 

Any pointers on what I'm doing wrong ?

Answer 1

While I wouldn't consider this as the best answer or solution - I decided to use Ensembl REST api to extract the coordinates of alignment blocks.

The quick fix - in python looks like this:

#This is the main function which uses Ensembl REST API to access the Enesembl Compara database 
#And retrive the alignment blocks which correspond to the MM10 standard. 
#The function accepts chromosome, start and array and returns the alignment blocks
#reported by the API. Since the blocks can be on different chromosomes/and/or/contigs, it 
#returns the whole output in form of dataframe - containing chromosome, start and stop of the
#alignment blocks corresponding to input and all the other alignment blocks corresponding the 
#16 mice species, if available. 
#!!!!!!!----------Caution start----------!!!!!!!
#The API is a bit buggy, so it is advised not to request for specific species.
#!!!!!!!-----------Caution end-----------!!!!!!!
def liftover(chrom, start, end):
    mmStart, mmEnd, mmchr = [], [], []
    caStart, caEnd, cachr = [], [], []
    #ryStart, ryEnd, rychr = [], [], []
    server = "https://rest.ensembl.org"
    ext = "/alignment/region/mus_musculus/" + str(chrom) + ":" + str(start) + "-" + str(end) + ":1?species_set_group=murinae;method=CACTUS_HAL;content-type=application/json"
    headers={ "Content-Type" : "application/json"}
    r = requests.get(server+ext, headers={ "Content-Type" : "application/json"})
    if(r.ok):
        json_string = json.dumps(r.json())
        jData = json.loads(r.content)
        results = r.json()
        print("The response contains {0} properties".format(len(jData)))
        print("\n")
    else:
        r.raise_for_status()

    for element in jData:
        element.pop('tree', None)

    for element in jData:
        for key, value in element.items():
            for z in value:
                if (z["species"]=="mus_musculus"):
                    mmStart.append(z["start"])
                    mmEnd.append(z["end"])
                    mmchr.append(z["seq_region"])
                if (z["species"]=="mus_musculus_casteij"):
                    caStart.append(z["start"])
                    caEnd.append(z["end"])
                    cachr.append(z["seq_region"])
    if (len(mmStart) != len(caStart)):
        print("Length mismatch between species blocks - proceed with caution!")
    return [min(mmStart), max(mmEnd), min(caStart), max(caEnd), len(mmStart)]

Again, it depends on how alignment blocks are reported - and requires further digging to resolve which mm10 alignment block resembles CAST alignment block.