Bajar de ensembl todas las secuencias de proteína de los ortólogos de un gen (o de varios)

Dada una secuencia de proteína, puede interesarnos bajar de una sola vez todas las secuencias de proteínas ortólogas conocidas. Si aceptamos el criterio de ensembl para determinar qué es y qué no es ortólogo (una buena opción en una decisión a veces complicada), las siguientes instrucciones pueden servir de guía para esta tarea.

i) pasos previos

Tendréis que tener instalado el ensemble-api (ver instrucciones de instalación) y el módulo DBD::Mysql. En Mac OS X he conseguido instalar este módulo usando macports. En 5 líneas, sería:

~% sudo port selfupdate

~% sudo port install mysql5-server

~% sudo port install p5-dbd-mysql

~% sudo port install p5-dbi

~% sudo port -f activate p5-dbi

ii) bajar los ortólogos

En un fichero de texto, ponéis una lista con los identificadores de genes para los que queréis hacer la búsqueda.

(por ejemplo)

~% cat IDs.txt

ENSMUSG00000031628

ENSMUSG00000027997

ENSG00000003400

Ejecutáis get_orthologues.pl dándole este fichero como $ARGV[0]

~% perl get_orthologues.pl IDs.txt

Si veis algo así en pantalla, todo OK.

.

.

--

Obtaining 1-to-1 and 1-to-many orthologues for ENSMUSG00000031628

main IDs for each orthologue written to ENSMUSG00000031628.mainIDs.txt

all IDs for each orthologue written to ENSMUSG00000031628.allIDs.txt

sequeences for the main IDs dumped to ENSMUSG00000031628.ort.fa

--

Obtaining 1-to-1 and 1-to-many orthologues for ENSMUSG00000027997

main IDs for each orthologue written to ENSMUSG00000027997.mainIDs.txt

all IDs for each orthologue written to ENSMUSG00000027997.allIDs.txt

sequeences for the main IDs dumped to ENSMUSG00000027997.ort.fa

--

Obtaining 1-to-1 and 1-to-many orthologues for ENSG00000003400

main IDs for each orthologue written to ENSG00000003400.mainIDs.txt

all IDs for each orthologue written to ENSG00000003400.allIDs.txt

sequeences for the main IDs dumped to ENSG00000003400.ort.fa

.

.

(en mi ordenador, para estos 3 IDs tarda unos 5 minutos)

El resultado, para cada identificador ENSG.., son como véis 3 ficheros:

ENSG*.mainIDs.txt --> listado con los protein-IDs del transcript principal de cada ortólogo

ENSG*.allIDs.txt --> listado con los protein-IDs de todas las secuencias que ensembl guarda para cada ortólogo

ENSG*.ort.fa --> secuencias de proteina (formato fasta) correspondientes al transcript principal de cada ortólogo

Aclaración por si alguien se ha liado con esto de transcript_principal / todas_las_secuencias. Como sabéis, para cada gen (y para cada uno de sus ortólogos) hay descritos varios transcripts normalmente. ¿Cuál utiliza ensembl como transcript principal? El más largo. La secuencia de proteína de ese se graba en OUT3, su ID en OUT1 y las IDs de todos en OUT2.

El código de get_orthologues.pl es el siguiente:

#!/usr/bin/perl

# Script to retrieve protein sequences for 1-to-1(and 1-to-many)-orthologues for a given gene.
# (adaptado del que me pasó Bert Overduin de helpdesk@enseml.org)

# Definicion de librerias y modulos

use strict;
use warnings;

use lib '/Users/jramon/Soft/ensembl_API/ensembl/modules/';
use lib '/Users/jramon/Soft/ensembl_API/ensembl-compara/modules/';
use lib '/Users/jramon/Soft/ensembl_API/ensembl-variation/modules/';
use lib '/Users/jramon/Soft/ensembl_API/ensembl-functgenomics/modules/';
use lib '/Users/jramon/Soft/BioPerl-1.6.1/';
use Bio::EnsEMBL::Registry;

my $reg = "Bio::EnsEMBL::Registry";
$reg->load_registry_from_url('mysql://anonymous@ensembldb.ensembl.org');

my $member_adaptor = $reg->get_adaptor("Multi", "compara", "Member");
my $homology_adaptor = $reg->get_adaptor("Multi", "compara", "Homology");


# Leer identificadores de genes y buscar ortologos para cada uno
# guardando IDs en OUT1,2 y secuencias en OUT3

open (INPUT_IDs,"$ARGV[0]");
while () {
my $gene_stable_id = $_;
chop $gene_stable_id;

my $member = $member_adaptor->fetch_by_source_stable_id("ENSEMBLGENE", $gene_stable_id);

open (OUT1,">$gene_stable_id.mainIDs.txt");
open (OUT2,">$gene_stable_id.allIDs.txt");
open (OUT3,">$gene_stable_id.ort.fa");
print "--\nObtaining 1-to-1 and 1-to-many orthologues for $gene_stable_id\n";
print "main IDs for each orthologue written to $gene_stable_id.mainIDs.txt\n";
print "all IDs for each orthologue written to $gene_stable_id.allIDs.txt\n";
print "sequences for the main IDs dumped to $gene_stable_id.ort.fa\n";

my $all_homologies = $homology_adaptor->fetch_all_by_Member($member);

foreach my $homology (@$all_homologies) {
if($homology->description eq (('ortholog_one2one')||('ortholog_one2many'))){
my $members = $homology->get_all_Members();
foreach my $member (@$members) {
 if($member->stable_id ne $gene_stable_id){
  my $canonical_peptide = $member->get_canonical_peptide_Member;
  print OUT1
   $canonical_peptide-> stable_id, "\n";

  my $gene_adaptor = $reg->get_adaptor($member->genome_db->name, "core", "Gene");
  my $gene = $gene_adaptor->fetch_by_stable_id($member->stable_id);
  my @transcripts = @{$gene->get_all_Transcripts};
  foreach my $transcript(@transcripts){
   if($transcript->translation){
     print OUT2 $transcript->translation->stable_id, "\n";
   }
  }

  print OUT3
   ">", $canonical_peptide->stable_id, "\n",
   $canonical_peptide->sequence, "\n";
  }
}
}
}
close OUT1;
close OUT2;
close OUT3;
}
close INPUT_IDs;

De todas formas, para el volumen de datos que ha de bajar, se me hace un poco lento el código. Se agradecen mejoras en este sentido.