Foro - Perl en Español

Publicado: **2011-08-23 11:19 @513**

Hola.
Necesito hacer un script que extraiga regiones al azar de un archivo genbank. La región debe tener 500 nt de longitud.

Mi código es el siguiente:

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

use Bio::SeqIO;
 
sub extract_intergenic_from_genbank {
    my ( $infile, $out_intergenic_file, $min_intergenic_size, $max_intergenic_size ) = @_;
 
    my ( $n_of_intergenic, $gi, $start, $end, $length, $strand, $taxon ) = (0);
    my $in = new Bio::SeqIO( -file => "NC_002755.gb", -format => 'genbank' );
 
    open( FNA, ">$salida.FASTA" )
        || die "# extract_intergenic_from_genbank : cannot create $out_intergenic_file\n";
 
    while ( my $seq = $in->next_seq )  # scan all sequences found in $input
    {
        my ( $gbaccession, $sequence, $gen, @genes ) = ( $seq->accession() );
        $sequence = $seq->primary_seq()->seq()
            || 'no encuentro la secuencia primaria, necesito un fichero GenBank completo!';
        $taxon = '';
 
        for my $f ( $seq->get_SeqFeatures ) {
            if ( $f->primary_tag =~ /CDS|rRNA|tRNA/ )    # campos de 'genes'
            {
                $gi = '';
                if ( $f->has_tag('db_xref') ) {
                    my $crossrefs = join( ',', sort $f->each_tag_value('db_xref') );
                    if ( $crossrefs =~ /(GI\:\d+)/ ) { $gi = $1 }
                }
                elsif ( $f->has_tag('locus_tag') ) {
                    if ( $gi eq '' && $f->has_tag('locus_tag') ) {
                        $gi = "ID:" . join( ',', sort $f->each_tag_value('locus_tag') );
                    }
                }
 
                next if ( $gi eq '' );
 
                $start  = $f->start();
                $end    = $f->end();
                $strand = $f->location()->strand();
                push( @genes, [ $start, $end, $gi, $strand ] );
            }
            elsif ( $f->primary_tag() =~ /source/ ) {
                if ( $f->has_tag('organism') ) {
                    foreach my $element ( $f->each_tag_value('organism') ) {
                        $taxon .= "[$element],";
                    }
                    chop $taxon;
                }
            }
        }
 
        for ( $gen = 1; $gen < scalar(@genes); $gen++ ) {
            $n_of_intergenic++;
 
            $start  = int( rand(4403249) );
            $end    = $start + 500;
            $length = $end - $start;
 
            print FNA ">intergenic$n_of_intergenic|$gbaccession|coords:$start..$end|"
                . "length:$length|neighbours:$genes[$gen-1][2]($genes[$gen-1][3]),"
                . "$genes[$gen][2]($genes[$gen][3])|$taxon\n";
 
            print FNA substr( $sequence, $start, $length ) . "\n";
 
        }
    }
 
    close(FNA);
 
    return $n_of_intergenic;
}
 
&extract_intergenic_from_genbank;
 
Coloreado en 0.004 segundos,  usando GeSHi 1.0.8.4

En esta parte del código $start = int(rand(4403249));, 4403249 es igual a la última coordenada del genbank pero eso yo lo elaboré solo para este archivo genbank. Necesito que en vez de ese número esté en una variable que contenga la última coordenada del genbank. ¿Cómo podría hacer eso?

Gracias

Publicado: **2011-08-23 11:31 @521**

Se lo puedes pasar como un argumento más, a la subrutina. Solo tienes que recibir el valor en una nueva variable que pondrás en la línea 4.

Publicado: **2011-08-23 11:42 @529**

Sí, pero lo que necesito en sí es poder hacer este script para cualquier archivo genbank donde la última coordenada nunca va a ser igual.

En está posición del arreglo, $genes[$gen][1], ahí hay una lista de coordenadas. ¿Cómo extraigo la última de ellas en una variable?

Publicado: **2011-08-23 11:53 @536**

En $genes[$gen][1] solo hay una coordenada: la posición final del gen número $gen.

Si quieres extraer la posición final del último gen, sería:

my $ultimo = $genes[-1][1];

Publicado: **2011-08-23 11:59 @541**

Muchísimas gracias: ¡esa es la última coordenada!

Foro - Perl en Español

Extraer regiones al azar de un genoma (genbank)

Extraer regiones al azar de un genoma (genbank)

Re: Extraer regiones al azar de un genoma (genbank)

Re: Extraer regiones al azar de un genoma (genbank)

Re: Extraer regiones al azar de un genoma (genbank)

Re: Extraer regiones al azar de un genoma (genbank)