transform() and transfer() of Ensembl Perl API

Dear friends,

When converting coordinates between assemblies such as from grch38 to grch37 or vice-versa using ensembl perl API using transfer() and transform(), how do we know if the coordinates are converted between the assemblies? In the tutorial (http://asia.ensembl.org/info/docs/api/core/core_tutorial.html#coordinates) there is no mention of grch38 or grch37 or any other assembly.

I am using ensemble release 95 and want to convert the coordinates of a specific region from grch38 to 37. Any comments will be very helpful.

Thanks,

ensembl api coordinates conversion

• 338 views • link

updated 7.3 years ago by Emily

• written 7.3 years ago by DanielC

1 answer

You want to use the project method on the slice to convert coordinates. Something like:
my $37_slice = $38_slice->project( 'chromosome', 'GRCh37' );

• 0 views • link

written 7.3 years ago by Emily

Thanks Emily! I wrote this code to convert the coordinates from 20000 to 30000 on chromosome 10 in grch37 using ensembl release 95:

use warnings;
use strict;
use Bio::EnsEMBL::Registry;

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

$registry->load_registry_from_db(
    -host => 'ensembldb.ensembl.org', # alternatively 'useastdb.ensembl.org'
    -user => 'anonymous'
);


my $slice_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Slice' );
my $slice = $slice_adaptor->fetch_by_region( 'chromosome', '10', 20000, 30000 );
my $coord_sys1  = $slice->coord_system()->name();
my $seq_region1 = $slice->seq_region_name();
my $start1      = $slice->start();
my $end1        = $slice->end();
my $strand1     = $slice->strand();
print "Slice: $coord_sys1 $seq_region1 $start1-$end1 ($strand1)\n";
my $cs_adaptor = $registry->get_adaptor( 'Human', 'Core', 'CoordSystem' );
my $cs = $cs_adaptor->fetch_by_name('chromosome');

printf "Coordinate system: %s %s\n", $cs->name(), $cs->version();


my $slice2 = $slice->project( 'chromosome', 'GRCh37' );
my $coord_sys  = $slice2->coord_system()->name();
my $seq_region = $slice2->seq_region_name();
my $start      = $slice2->start();
my $end        = $slice2->end();
my $strand     = $slice2->strand();

print "New_Slice: $coord_sys $seq_region $start-$end ($strand)\n";

And I get his output:

Slice: chromosome 10 25000-30000 (1)
Coordinate system: chromosome GRCh38
Can't call method "coord_system" on unblessed reference at test.pl line 41.

Could you please provide your comments on resolving this error? Thanks much!

• 0 views • link

written 7.3 years ago by DanielC

The project method doesn't return a slice, but a list of arrays, which are each [$start,$end,$slice]. This is because a region does not necessarily project to only one region on the other assembly, it might in fact be two or more.

You'll need to dereference the array, get the items one-by-one and get the slice. Something like:

my @projection = @{  $slice->project( 'chromosome', 'GRCh37' ) };
foreach my $projection_segment (@projection) {   
    my $slice2 = $projection_segment[2]
    my $coord_sys  = $slice2->coord_system()->name();
    # Do the rest of your stuff with $slice2
}

• 0 views • link

written 7.3 years ago by Emily

I figured out (thanks for your reply anyway), the script is this:

use warnings;
use strict;
use Bio::EnsEMBL::Registry;
my $registry = 'Bio::EnsEMBL::Registry';

$registry->load_registry_from_db(
    -host => 'ensembldb.ensembl.org', # alternatively 'useastdb.ensembl.org'
    -user => 'anonymous'
);

my $slice_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Slice' );
my $slice = $slice_adaptor->fetch_by_region( 'chromosome', '10', 25000, 30000 );
my $coord_sys1  = $slice->coord_system()->name();
my $seq_region1 = $slice->seq_region_name();
my $start1      = $slice->start();
my $end1        = $slice->end();
my $strand1     = $slice->strand();
print "Slice: $coord_sys1 $seq_region1 $start1-$end1 ($strand1)\n";
my $cs_adaptor = $registry->get_adaptor( 'Human', 'Core', 'CoordSystem' );
my $cs = $cs_adaptor->fetch_by_name('chromosome');
printf "Coordinate system: %s %s\n", $cs->name(), $cs->version();
my $slice2 = $slice->project( 'chromosome', 'GRCh37' );
foreach my $segment (@$slice2) {
      my $chromosome = $segment->to_Slice();
      print $slice->seq_region_name(), ':', $segment->from_start(), '-',
            $segment->from_end(), ' -> ',
            $chromosome->seq_region_name(), ':', $chromosome->start(), '-',$chromosome->end(),
            ':', $chromosome->strand(), "\n";
    }

The output is:

Slice: chromosome 10 25000-30000 (1)
Coordinate system: chromosome GRCh38
10:1-846 -> 10:70936-71781:1
10:847-1247 -> 10:71784-72184:1
10:1250-2609 -> 10:72185-73544:1
10:2610-5001 -> 10:73546-75937:1

One question: Why in the output I am seeing is different region for "slice", where the region I am converting is from 25000 to 30000:

10:1-846 
    10:847-1247 
    10:1250-2609 
    10:2610-5001

Please let me know your comments. Thanks!

• 0 views • link

written 7.3 years ago by DanielC

That's 1-846 of your region of interest. It's just saying that the first 846 bases of your original region correspond to 10:70936-71781:1 in GRCh37.

• 0 views • link

written 7.3 years ago by Emily

I see, OK, perfect! Thanks!

• 0 views • link

written 7.3 years ago by DanielC