<div dir="ltr">Dear Stephen, <div><br></div><div>You were absolutely right. when retrieving just a single block (my ($ref_start, $ref_end) = (7247302, 7267366), for example) phyml works seamlessly.</div><div><br></div><div>
I however feel that this is just a workaround. As the web example also works with phyml and is composed by more than one block, I tend to think that the problem is somewhere else (maybe in the alignment you mentioned?). Anyway, finding this issue is out of the scope of my work right now: I am focused on testing this as a proof of concept of the technology, that has succeeded thanks to your help. I only hope that if anyone else tries to perform the same operation finds this useful as a starting point.</div>
<div><br></div><div>Thanks again for your fast and reliable support. Best regards, </div><div><br></div><div><br></div><div><br></div><div>Manuel</div><div><br></div><div class="gmail_extra"><br><br><div class="gmail_quote">
2014/1/14 Stephen Fitzgerald <span dir="ltr"><<a href="mailto:stephenf@ebi.ac.uk" target="_blank">stephenf@ebi.ac.uk</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi Manuel, yes, there is a difference between the data retrieved by the<br>
web and the data retrieved by the API code.<br>
Because the region (chr19:7112266-7294045) is quite large, it is actually<br>
composed of 6 primate alignment blocks:<br>
+-------------+-----------+<br>
| block_start | block_end |<br>
+-------------+-----------+<br>
| 7058566 | 7122884 |<br>
| 7122885 | 7192725 |<br>
| 7192726 | 7211509 |<br>
| 7211510 | 7226908 |<br>
| 7226909 | 7247301 |<br>
| 7247302 | 7267366 |<br>
+-------------+-----------+<br>
<br>
In this alignment:<br>
<a href="http://www.ensembl.org/Homo_sapiens/Gene/Compara_Alignments?align=548&db=core&g=ENSG00000171105&r=19%3A7112266-7294045" target="_blank">http://www.ensembl.org/Homo_<u></u>sapiens/Gene/Compara_<u></u>Alignments?align=548&db=core&<u></u>g=ENSG00000171105&r=19%<u></u>3A7112266-7294045</a><br>
The web-code retrieves the full alignment for the last 5 blocks and a<br>
truncated (restricted) region from the first block (7112266-7122884).<br>
However, the web-code will connect the blocks with unaligned<br>
reference sequence (human in this case) if there are any gaps<br>
between the blocks (in this region there are no gaps between the blocks)<br>
and pad out any upstream or downstream regions (where there are no<br>
alignment blocks) with reference sequence.<br>
In this case the downstream (on the +ve strand) 26679 bp in the<br>
7112266-7294045 region has no primate alignments (so, you get unaligned human sequence for this region).<br>
<br>
In contrast, the script will only retrieve the primate alignment blocks<br>
associated with the region (no extra reference padding sequence will be<br>
added) to connect or extend the blocks. In the above example, if you<br>
output the results of the script to a file, you will get six blocks in the<br>
one file. I have not run phyml/phylip, but it may expect only one block of<br>
alignment per file.<br>
<br>
Hope this helps.<span class="HOEnZb"><font color="#888888"><br>
Stephen.</font></span><div class="HOEnZb"><div class="h5"><br>
<br>
<br>
<br>
On Tue, 14 Jan 2014, Manuel Rodríguez Pascual wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi all, <br>
Thanks for your valuable suggestions. however, I still have some minor issues.<br>
<br>
I understand that in Stephen's code, in order to obtain results equivalent to the ones in the web, <br>
<br>
<a href="http://www.ensembl.org/Homo_sapiens/Gene/Compara_Alignments?align=548&db=core&g=ENSG00000171105&r=19%3A7112266-7294045" target="_blank">http://www.ensembl.org/Homo_<u></u>sapiens/Gene/Compara_<u></u>Alignments?align=548&db=core&<u></u>g=ENSG00000171105&r=19%<u></u>3A7112266-7294045</a><br>
<br>
variables ($ref_start, $ref_end) should be asigned to (7112266, 7294045), as my desired gene region is "Chromosome 19:<br>
7,112,266-7,294,045 reverse strand". I have called $ref_slice_adaptor->fetch_by_<u></u>region with strand 1 and -1.<br>
<br>
However, the results I get are not the same than in your web example. Also, when executing PHYML on the results, it suceeds on the web example<br>
but fails on the one obtained with Stephen's code. I don't know wheter it is a minor format issue or a deeper problem. Have you got any<br>
experience on running phyml/phylip on the fetched results?<br>
<br>
<br>
Thank you very much for your help. Best regards,<br>
<br>
<br>
Manuel<br>
<br>
<br>
<br>
<br>
<br>
<br>
2014/1/13 Stephen Fitzgerald <<a href="mailto:stephenf@ebi.ac.uk" target="_blank">stephenf@ebi.ac.uk</a>><br>
Hi Manuel, you can also try our REST service, which will give you the option of accessing the data (output in JSON format) in a<br>
number of ways:<br>
<a href="http://beta.rest.ensembl.org/documentation/info/genomic_alignment_block_region" target="_blank">http://beta.rest.ensembl.org/<u></u>documentation/info/genomic_<u></u>alignment_block_region</a><br>
for example (using human coords 19:7184302-7184344):<br>
bash> curl '<a href="http://beta.rest.ensembl.org/alignment/block/region/homo_sapiens/19:7184302-7184344:1?species_set_group=primates" target="_blank">http://beta.rest.ensembl.org/<u></u>alignment/block/region/homo_<u></u>sapiens/19:7184302-7184344:1?<u></u>species_set_group=primates</a>' -H<br>
'Content-type:application/<u></u>json'<br>
<br>
Also, if you wish to use the perl API, here is a script for retrieving the alignments in phylip format:<br>
bash> perl <a href="http://script.pl" target="_blank">script.pl</a> > out.phylip<br>
<br>
################ <a href="http://script.pl" target="_blank">script.pl</a><br>
<br>
use strict;<br>
use warnings;<br>
use Data::Dumper;<br>
use Bio::AlignIO;<br>
<br>
use Bio::EnsEMBL::Registry;<br>
<br>
#Auto-configure the registry<br>
Bio::EnsEMBL::Registry->load_<u></u>registry_from_db(<br>
-host=>"<a href="http://ensembldb.ensembl.org" target="_blank">ensembldb.ensembl.org</a>"<u></u>, -user=>"anonymous",<br>
-port=>'5306', db_version => 74);<br>
<br>
<br>
# Get the Compara Adaptor for MethodLinkSpeciesSets<br>
my $method_link_species_set_<u></u>adaptor =<br>
Bio::EnsEMBL::Registry->get_<u></u>adaptor(<br>
"Multi", "compara", "MethodLinkSpeciesSet");<br>
<br>
my $methodLinkSpeciesSet = $method_link_species_set_<u></u>adaptor-><br>
fetch_by_method_link_type_<u></u>species_set_name("EPO", "primates");<br>
<br>
# Define the start and end positions for the alignment<br>
my ($ref_start, $ref_end) = (7184302, 7184344);<br>
<br>
# Get the reference species *core* Adaptor for Slices<br>
my $ref_slice_adaptor =<br>
Bio::EnsEMBL::Registry->get_<u></u>adaptor(<br>
"homo_sapiens", "core", "Slice");<br>
<br>
# Get the slice corresponding to the region of interest<br>
my $ref_slice = $ref_slice_adaptor->fetch_by_<u></u>region(<br>
"chromosome", 19, $ref_start, $ref_end);<br>
<br>
# Get the Compara Adaptor for GenomicAlignBlocks<br>
my $genomic_align_block_adaptor =<br>
Bio::EnsEMBL::Registry->get_<u></u>adaptor(<br>
"Multi", "compara", "GenomicAlignBlock");<br>
<br>
# The fetch_all_by_<u></u>MethodLinkSpeciesSet_Slice() returns a ref.<br>
# to an array of GenomicAlingBlock objects (human is the reference species)<br>
my $all_genomic_align_blocks = $genomic_align_block_adaptor-><br>
fetch_all_by_<u></u>MethodLinkSpeciesSet_Slice(<br>
$methodLinkSpeciesSet, $ref_slice);<br>
<br>
# set up an AlignIO to format SimpleAlign output<br>
my $alignIO = Bio::AlignIO->newFh(-<u></u>interleaved => 0,<br>
-fh => \*STDOUT,<br>
-format => 'phylip',<br>
-idlength => 30);<br>
<br>
# print the restricted alignments<br>
foreach my $genomic_align_block( @{ $all_genomic_align_blocks }) {<br>
my $restricted_gab = $genomic_align_block-><u></u>restrict_between_reference_<u></u>positions($ref_start, $ref_end);<br>
print $alignIO $restricted_gab->get_<u></u>SimpleAlign;<br>
}<br>
<br>
################ out.phylip<br>
6 44<br>
homo_sapiens/19 cccc-<u></u>agacccacatccagaactcacttgctggaa<u></u>ttcatcgtg<br>
pongo_abelii/19 cccc-<u></u>agacctacatccagaactcacttgctggaa<u></u>ttcatcgtg<br>
pan_troglodytes/19 cccc-<u></u>agacccacaaccagaactcacttgctggaa<u></u>ttcatcgtg<br>
gorilla_gorilla/19 cccc-<u></u>agacccacatctagaactcacttgctggaa<u></u>ttcatcgtg<br>
callithrix_jacchus/22 ccccaagacccacatgcaggactcacttgc<u></u>tggaattcatcgtg<br>
macaca_mulatta/19 cccc-<u></u>ggacccacatacagaactcacttgctggaa<u></u>ttcatcgtg<br>
<br>
<br>
Cheers,<br>
Stephen.<br>
<br>
<br>
<br>
On Mon, 13 Jan 2014, Emily Pritchard wrote:<br>
<br>
Hi Manuel<br>
<br>
Have you seen our API course on EBI Train Online:<br>
<a href="http://www.ebi.ac.uk/training/online/course/ensembl-filmed-api-workshop" target="_blank">http://www.ebi.ac.uk/training/<u></u>online/course/ensembl-filmed-<u></u>api-workshop</a><br>
<br>
The Core section will introduce you to the API itself, and the Compara module for alignments.<br>
<br>
Hope this helps<br>
<br>
Emily<br>
<br>
On 13/01/2014 13:32, Manuel Rodríguez Pascual wrote:<br>
I am just starting working with Ensembl API. I am not really experienced neither with Ensembl or bioinformatics<br>
itself,so I am stuck<br>
in a problem that seems to be easy to solve.<br>
<br>
<br>
As a test and proof of concept, I am interested to retrieve information of a comparison in phylip format. In particular,<br>
the provided<br>
example<br>
<a href="http://www.ensembl.org/Homo_sapiens/Gene/Compara_Alignments?align=548&db=core&g=ENSG00000171105&r=19%3A7112266-7294045" target="_blank">http://www.ensembl.org/Homo_<u></u>sapiens/Gene/Compara_<u></u>Alignments?align=548&db=core&<u></u>g=ENSG00000171105&r=19%<u></u>3A7112266-7294045</a><br>
<br>
when exported into a phylip format, <br>
<br>
<a href="http://www.ensembl.org/Homo_sapiens/Export/Output/Location/Alignment?align=548;db=core;g=ENSG00000171105;output=alignment;r=19:7112266-7294045" target="_blank">http://www.ensembl.org/Homo_<u></u>sapiens/Export/Output/<u></u>Location/Alignment?align=548;<u></u>db=core;g=ENSG00000171105;<u></u>output=alignment;r=19:7112266-<u></u>7294045</a>;<br>
format=phylip;_format=Text<br>
<br>
I have however seen that the employment of wget is discouraged and the Ensembl API should be used instead.<br>
<br>
Reading the available documentation, I fell that I should use the Compara API Tutorial, <br>
<a href="http://www.ensembl.org/info/docs/api/compara/compara_tutorial.html" target="_blank">http://www.ensembl.org/info/<u></u>docs/api/compara/compara_<u></u>tutorial.html</a><br>
<br>
but I don't really understand it or how to apply it to my objective.<br>
<br>
The question is, can anyone orient me on how to proceed? <br>
<br>
<br>
<br>
Thanks for your attention,<br>
<br>
<br>
<br>
Manuel <br>
<br>
<br>
-- <br>
Dr. Manuel Rodríguez-Pascual<br>
skype: manuel.rodriguez.pascual<br>
phone: (+34) 913466173 // (+34) 679925108<br>
<br>
CIEMAT-Moncloa<br>
Edificio 22, desp. 1.25<br>
Avenida Complutense, 40 <br>
28040- MADRID<br>
SPAIN<br>
<br>
<br>
______________________________<u></u>_________________<br>
Dev mailing list <a href="mailto:Dev@ensembl.org" target="_blank">Dev@ensembl.org</a><br>
Posting guidelines and subscribe/unsubscribe info: <a href="http://lists.ensembl.org/mailman/listinfo/dev" target="_blank">http://lists.ensembl.org/<u></u>mailman/listinfo/dev</a><br>
Ensembl Blog: <a href="http://www.ensembl.info/" target="_blank">http://www.ensembl.info/</a><br>
<br>
<br>
--<br>
Dr Emily Pritchard<br>
Ensembl Outreach Officer<br>
<br>
European Bioinformatics Institute (EMBL-EBI)<br>
European Molecular Biology Laboratory<br>
Wellcome Trust Genome Campus<br>
Hinxton<br>
Cambridge<br>
CB10 1SD<br>
UK<br>
<br>
<br>
______________________________<u></u>_________________<br>
Dev mailing list <a href="mailto:Dev@ensembl.org" target="_blank">Dev@ensembl.org</a><br>
Posting guidelines and subscribe/unsubscribe info: <a href="http://lists.ensembl.org/mailman/listinfo/dev" target="_blank">http://lists.ensembl.org/<u></u>mailman/listinfo/dev</a><br>
Ensembl Blog: <a href="http://www.ensembl.info/" target="_blank">http://www.ensembl.info/</a><br>
<br>
<br>
<br>
<br>
--<br>
Dr. Manuel Rodríguez-Pascual<br>
skype: manuel.rodriguez.pascual<br>
phone: (+34) 913466173 // (+34) 679925108<br>
<br>
CIEMAT-Moncloa<br>
Edificio 22, desp. 1.25<br>
Avenida Complutense, 40 <br>
28040- MADRID<br>
SPAIN<br>
<br>
</blockquote>
</div></div><br>_______________________________________________<br>
Dev mailing list <a href="mailto:Dev@ensembl.org">Dev@ensembl.org</a><br>
Posting guidelines and subscribe/unsubscribe info: <a href="http://lists.ensembl.org/mailman/listinfo/dev" target="_blank">http://lists.ensembl.org/mailman/listinfo/dev</a><br>
Ensembl Blog: <a href="http://www.ensembl.info/" target="_blank">http://www.ensembl.info/</a><br>
<br></blockquote></div><br><br clear="all"><div><br></div>-- <br>Dr. Manuel Rodríguez-Pascual<br>skype: manuel.rodriguez.pascual<br>phone: (+34) 913466173 // (+34) 679925108<br> <br>CIEMAT-Moncloa<br>Edificio 22, desp. 1.25<br>
Avenida Complutense, 40 <br>28040- MADRID<br>SPAIN
</div></div>