• Publicidad

Script para EnsEMBL

¿Apenas comienzas con Perl? En este foro podrás encontrar y hacer preguntas básicas de Perl con respuestas aptas a tu nivel.

Script para EnsEMBL

Notapor bioinfo » 2009-01-31 07:02 @335

Necesito un poco de ayuda con este script para el Perl API de EnsEMBL.

Cuando lo ejecuto en su forma original me devuelve el listado completo de homologías para la lista de genes y funciona aparentemente bien. Cuando sustituyo la línea correspondiente para que solo me liste los ortólogos de mamíferos no lo admite.

A ver si hay por aquí alguien con experiencia en esta API que me pueda ayudar. Gracias de antemano, majos.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
    use strict;
    use warnings;

    use Bio::EnsEMBL::Registry;

    ## Load the registry automatically
    my $reg = "Bio::EnsEMBL::Registry";
    $reg->load_registry_from_url('mysql://[email protected]');

    ## Get the human gene adaptor
    my $human_gene_adaptor = $reg->get_adaptor("Homo sapiens", "core", "Gene");

    ## Get the compara member adaptor
    my $member_adaptor = $reg->get_adaptor("Compara", "compara", "Member");

    ## Get the compara homology adaptor
    my $homology_adaptor = $reg->get_adaptor("Compara", "compara", "Homology");

    ## Lista de genes humanos para buscar homologos
    my @BreastCANgenes = ('ABCA3','ABCB10','ABCB8','ACADM');

    open ("archivo", ">homologos.txt");
    my $gen;

    foreach $gen (@BreastCANgenes) {

        ## Get all existing gene object
        my $ctdp1_genes = $human_gene_adaptor->fetch_all_by_external_name($gen);

        ## For each of these genes...
        foreach my $ctdp1_gene (@$ctdp1_genes) {

            ## Get the compara member
            my $member = $member_adaptor->fetch_by_source_stable_id("ENSEMBLGENE", $ctdp1_gene->stable_id);

            ## Get all the homologues
            my $all_homologies = $homology_adaptor->fetch_all_by_Member($member);

            ## For each homology
            foreach my $this_homology (@$all_homologies) {

                ## print the description (type of homology) and the
                ## subtype (taxonomy level of the event: duplic. or speciation)
                print $this_homology->description, " [", $this_homology->subtype, "]\n";
                print archivo $this_homology->description, " [", $this_homology->subtype, "]\n";

                ## print the members in this homology
                my $members = $this_homology->get_all_Members();

                foreach my $this_member (@$members) {
                    print archivo   $this_member->source_name, " ",
                                    $this_member->stable_id, " (",
                                    $this_member->genome_db->name, ")\n"
                }
                print "\n";
            }
        }
    }
    close ("archivo");
Coloreado en 0.006 segundos, usando GeSHi 1.0.8.4



Esta es la versión modificada para obtener solo ortólogos en mamíferos:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
    use strict;
    use warnings;

    use Bio::EnsEMBL::Registry;

    ## Load the registry automatically
    my $reg = "Bio::EnsEMBL::Registry";

    $reg->load_registry_from_url('mysql://[email protected]');

    ## Get the human gene adaptor
    my $human_gene_adaptor = $reg->get_adaptor("Homo sapiens", "core", "Gene");
    ## Get the compara member adaptor
    my $member_adaptor = $reg->get_adaptor("Compara", "compara", "Member");
    ## Get the compara homology adaptor
    my $homology_adaptor = $reg->get_adaptor("Compara", "compara", "Homology");

    use Bio::EnsEMBL::Compara::MethodLinkSpeciesSet;

    my @BreastCANgenes = ('ABCA3','ABCB10','ABCB8','ACADM');

    open ("archivo", ">homologos.txt");
    my $gen;

    foreach $gen (@BreastCANgenes) {

        ## Get all existing gene object
        my $ctdp1_genes = $human_gene_adaptor->fetch_all_by_external_name($gen);

        ## For each of these genes...
        foreach my $ctdp1_gene (@$ctdp1_genes) {

            ## Get the compara member
            my $member = $member_adaptor->fetch_by_source_stable_id("ENSEMBLGENE", $ctdp1_gene->stable_id);

            ## Get all the homologues
            my $all_homologies
                = $homology_adaptor
                ->fetch_all_by_MethodLinkSpeciesSet_orthology_type_subtype($member,'ortholog_one2one','Mammalia',);

            ## For each homology
            foreach my $this_homology (@$all_homologies) {

                ## print the description (type of homology) and the
                ## subtype (taxonomy level of the event: duplic. or speciation)
                print $this_homology->description, " [", $this_homology->subtype, "]\n";
                print archivo $this_homology->description, " [", $this_homology->subtype, "]\n";

                ## print the members in this homology
                my $members = $this_homology->get_all_Members();
                foreach my $this_member (@$members) {
                    print archivo   $this_member->source_name, " ",
                                    $this_member->stable_id, " (",
                                    $this_member->genome_db->name, ")\n"
                }
                print "\n";
            }
        }
    }
    close ("archivo");
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4
bioinfo
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2009-01-31 06:50 @326

Publicidad

Notapor explorer » 2009-01-31 13:33 @606

Bienvenido a los foros de Perl en Español, bioinfo.

Si no me equivoco, para usar el fetch_all_by_MethodLinkSpeciesSet_orthology_type_subtype debes hacerlo sobre un adaptor de tipo MethodLinkSpeciesSet.

Es decir, a la hora de crear el adaptor, estás usando
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
my $homology_adaptor = $reg->get_adaptor("Compara", "compara", "Homology");
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

pero debería ser
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
my $homology_adaptor = $reg->get_adaptor('Multi', 'compara', 'MethodLinkSpeciesSet');
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

(humm... creo que el 'Multi' es para acceder a varias especies a la vez... no estoy seguro).

Aquí tienes un ejemplo de uso de Bio::EnsEMBL::Compara::MethodLinkSpeciesSet (busca por "example script").
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Notapor bioinfo » 2009-01-31 14:41 @653

Muchas gracias, sobre todo por la rapidez. Tenía en mente que era algo así, pero como soy bastante pardillo tanto en Perl como en la API... lo probaré a ver si lo consigo hacer funcionar y ya te cuento.
bioinfo
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2009-01-31 06:50 @326

Notapor bioinfo » 2009-02-01 04:28 @228

He probado y le he dado vueltas a la documentación pero no me funciona. Lo he dejado así:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
    use strict;
    use warnings;
    use Bio::EnsEMBL::Registry;
    ## Load the registry automatically
    my $reg = "Bio::EnsEMBL::Registry";
    $reg->load_registry_from_url('mysql://[email protected]');
    ## Get the human gene adaptor
    my $human_gene_adaptor =
        $reg->get_adaptor('Homo sapiens', 'core', 'Gene');
    ## Get the compara member adaptor
    my $member_adaptor =
        $reg->get_adaptor('Compara', 'compara', 'Member');
    ## Get the compara homology adaptor
    my $homology_adaptor =
        $reg->get_adaptor ('Multi', 'compara', 'MethodLinkSpeciesSet');


    my @BreastCANgenes = ('ABCA3','ABCB10','ABCB8','ACADM');

    open ('archivo', '>homologos.txt');
    my $gen;

    foreach $gen (@BreastCANgenes) {


    ## Get all existing gene object
    my $ctdp1_genes = $human_gene_adaptor->fetch_all_by_external_name($gen);
    ## For each of these genes...
    foreach my $ctdp1_gene (@$ctdp1_genes) {
      ## Get the compara member
      my $member = $member_adaptor->fetch_by_source_stable_id(
          "ENSEMBLGENE", $ctdp1_gene->stable_id);
      ## Get all the orthologues
      my $all_homologies = $homology_adaptor->fetch_all_by_MethodLinkSpeciesSet_orthology_type_subtype($member,'ortholog_one2one','Mammalia');
      ## For each homology
      foreach my $this_homology (@$all_homologies) {
        ## print the description (type of homology) and the
        ## subtype (taxonomy level of the event: duplic. or speciation)
        print $this_homology->description, " [", $this_homology->subtype, "]\n";
        print archivo $this_homology->description, " [", $this_homology->subtype, "]\n";
        ## print the members in this homology
        my $members = $this_homology->get_all_Members();
        foreach my $this_member (@$members) {
          print archivo $this_member->source_name, " ",
        $this_member->stable_id, " (",
        $this_member->genome_db->name, ")\n"


        }
        print "\n";
      }
    }
    }
    close ('archivo');
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4


Me da el siguiente error:
Sintáxis: [ Descargar ] [ Ocultar ]
Using bash Syntax Highlighting
Can't locate object method "fetch_all_by_MethodLinkSpeciesSet_orthology_type_subtype" via package "Bio::EnsEMBL::Compara::DBSQL::MethodLinkSpeciesSetAdaptor" at ortologos.pl line 34.
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4
bioinfo
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2009-01-31 06:50 @326

Notapor explorer » 2009-02-01 10:37 @484

¿Seguro que existe ese método? En la documentación no aparece.

El que sí viene es el fetch_all_by_MethodLinkSpeciesSet_Slice().

Otra opción para recuperar información a partir de los nombres de los genes podría ser usando Bio::Tools::Run::Ensembl, con la función get_gene_by_name(), indicando además, -use_orthologues para indicar qué ortologías queremos usar.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Notapor bioinfo » 2009-02-01 10:50 @493

Lo intentaré de esta forma que comentas. Por otro lado yo entiendo que lo otro tiene que funcionar de alguna forma ya que está en la documentación de la api. O a lo mejor es que yo no me entero...

¿que te parece?:

http://www.ensembl.org/info/docs/api/Pdoc/ensembl-compara/modules/Bio/EnsEMBL/Compara/DBSQL/HomologyAdaptor.html#POD7
bioinfo
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2009-01-31 06:50 @326

Notapor explorer » 2009-02-01 11:46 @532

Sí, está en la documentación, y además su código. Pero no veo ningún ejemplo en que se llame específicamente a este método. En el ejemplo que te enlacé antes, no se llama ni se ejecuta ningún método de Bio::EnsEMBL::Compara::MethodLinkSpeciesSet, sino que se hace de forma indirecta, a través del Bio::EnsEMBL::Registry.

Lamento no ayudarte más, pero en estos temas estoy muy pez.

Pongo aquí el código, como referencia para los demás.
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
use strict;
use Bio::EnsEMBL::Registry;
use Bio::EnsEMBL::Utils::Exception qw(throw);
use Bio::SimpleAlign;
use Bio::AlignIO;
use Bio::LocatableSeq;
use Getopt::Long;

my $usage = qq{
perl DumpMultiAlign.pl
  Getting help:
    [--help]

  For the query slice:
    [--species species]
        Query species. Default is "human"
    [--coord_system coordinates_name]
        Query coordinate system. Default is "chromosome"
    --seq_region region_name
        Query region name, i.e. the chromosome name
    --seq_region_start start
    --seq_region_end end

  For the alignments:
    [--alignment_type method_link_name]
        The type of alignment. Default is "BLASTZ_NET"
    [--set_of_species species1:species2:species3:...]
        The list of species used to get those alignments. Default is
        "human:mouse". The names should correspond to the name of the
        core database in the registry_configuration_file or any of its
        aliases

  Ouput:
    [--output_format clustalw|fasta|...]
        The type of output you want. "clustalw" is the default.
    [--output_file filename]
        The name of the output file. By default the output is the
        standard output
};

my $species = "human";
my $coord_system = "chromosome";
my $seq_region = "14";
my $seq_region_start = 75000000;
my $seq_region_end = 75010000;
my $alignment_type = "BLASTZ_NET";
my $set_of_species = "human:mouse";
my $output_file = undef;
my $output_format = "clustalw";
my $help;

GetOptions(
    "help" => \$help,
    "species=s" => \$species,
    "coord_system=s" => \$coord_system,
    "seq_region=s" => \$seq_region,
    "seq_region_start=i" => \$seq_region_start,
    "seq_region_end=i" => \$seq_region_end,
    "alignment_type=s" => \$alignment_type,
    "set_of_species=s" => \$set_of_species,
    "output_format=s" => \$output_format,
    "output_file=s" => \$output_file);

# Print Help and exit
if ($help) {
    print $usage;
    exit(0);
}

if ($output_file) {
    open(STDOUT, ">$output_file") or die("Cannot open $output_file");
}

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

# Getting all the Bio::EnsEMBL::Compara::GenomeDB objects
my $genome_dbs;
my $genome_db_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'GenomeDB');

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

foreach my $this_species (split(":", $set_of_species)) {
    my $this_meta_container_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
        $this_species, 'core', 'MetaContainer');

    throw("Registry configuration file has no data for connecting to <$this_species&")
        if (!$this_meta_container_adaptor);

    my $this_binomial_id = $this_meta_container_adaptor->get_Species->binomial;

    # Fetch Bio::EnsEMBL::Compara::GenomeDB object
    my $genome_db = $genome_db_adaptor->fetch_by_name_assembly($this_binomial_id);

    # Add Bio::EnsEMBL::Compara::GenomeDB object to the list
    push(@$genome_dbs, $genome_db);
}

# Getting Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object
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 $set_of_species!")
    if (!$method_link_species_set);

# Fetching the query Slice:
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;
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4
Última edición por explorer el 2009-02-08 07:02 @335, editado 1 vez en total
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Notapor bioinfo » 2009-02-01 12:21 @556

Muchas gracias por todo, me has dado ideas con las que pelear. Cuando lo tenga lo posteo.
bioinfo
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2009-01-31 06:50 @326


Volver a Básico

¿Quién está conectado?

Usuarios navegando por este Foro: No hay usuarios registrados visitando el Foro y 1 invitado

cron