Foro - Perl en Español

#!/usr/bin/perl -w
# Author: Damien Mattei C.N.R.S / U.N.S.A - UMR 6549
# example: ./idfetch.pl AP001266

use Bio::DB::GenBank;

$gb = new Bio::DB::GenBank();

# this returns a Seq object :
$seq1 = $gb->get_Seq_by_acc($ARGV[0]);
print $seq1->display_id() . "\n" ;


foreach $feat ($seq1->all_SeqFeatures()) {

  #print $feat->primary_tag . " " . $feat->source_tag() . "\n" ;

  print "Feature from ", $feat->start, " to ",
   $feat->end, " Primary tag  ", $feat->primary_tag,
   ", produced by ", $feat->source_tag(), "\n";

  if( $feat->strand == 0 ) {
    print "Feature applicable to either strand\n";
  } else {
    print "Feature on strand ", $feat->strand,"\n"; # -1,1
  }

  foreach $tag ( $feat->all_tags() ) {
    print "Feature has tag ", $tag, " with values, ",
    join(' ',$feat->each_tag_value($tag)), "\n";
  }

  print "new feature\n" if $feat->has_tag('new');

}


exit;

__END__

It will display something like that:

[dmattei@pclgmch2 gmap]$ ./idfetch.pl AP001266
AP001266
Feature from 1 to 168978 Primary tag  source, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag chromosome with values, 11
Feature has tag map with values, 11q13
Feature has tag clone with values, RP11-770G2
Feature has tag organism with values, Homo sapiens
Feature has tag db_xref with values, taxon:9606
Feature from 1 to 31550 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 31651 to 48510 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 48611 to 64044 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 64145 to 78208 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 78309 to 89008 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 89109 to 99704 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 99805 to 107965 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 108066 to 116032 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 116133 to 124010 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 124111 to 130494 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 130595 to 136072 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 136173 to 139649 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 139750 to 144590 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 144691 to 148482 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 148583 to 152279 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 152380 to 153632 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment clone_end:T7
vector_side:left
Feature from 153733 to 155746 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 155847 to 156405 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment clone_end:SP6
vector_side:right
Feature from 156506 to 158398 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 158499 to 161333 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 161434 to 163304 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 163405 to 164604 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 164705 to 166693 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment
Feature from 166794 to 168978 Primary tag  misc_feature, produced by
EMBL/GenBank/SwissProt
Feature on strand 1
Feature has tag note with values, assembly_fragment

use Bio::SeqIO;
 
my $seqio_object = Bio::SeqIO->new(-file => $gb_file);
my $seq_object = $seqio_object->next_seq;
 
for my $feat_object ($seq_object->get_SeqFeatures) {
   if ($feat_object->primary_tag eq "CDS") {
      print $feat_object->spliced_seq->seq,"\n";
      # e.g. 'ATTATTTTCGCTCGCTTCTCGCGCTTTTTGAGATAAGGTCGCGT...'
      if ($feat_object->has_tag('gene')) {
         for my $val ($feat_object->get_tag_values('gene')){
            print "gene: ",$val,"\n";
            # e.g. 'NDP', from a line like '/gene="NDP"'
         }
      }
   }
}

#!/usr/bin/perl
#
# Joaquín Ferrero. Diciembre 2007.
#
# Extracción de secuencias de nucleótidos de un fichero GenBank.
#
# El fichero a tratar se pasa como primer argumento
#

use strict;
use warnings;

# Si no hay argumentos, informamos al usuario
@ARGV or die "Uso: $0 <fichero genbank a usar>\n";

my $fichero = shift @ARGV;

# Abrimos y leemos el fichero
open my $gbh, '<', $fichero
     or die "ERROR: No pude abrir el fichero $fichero: $!\n";

my @genbank = <$gbh>;

close $gbh;

# Tratamiento del contenido
imprime_secuencias_del(\@genbank);

# Subrutina que imprime las secuencias del genbank
sub imprime_secuencias_del {
    my $genbank_ref = shift;

    # Tratamiento del fichero: por líneas
    # Vamos leyendo línea a línea hasta que encontramos
    # un final de sección ORIGIN o fin de fichero, y
    # procesamos esa parte
    my $parte;
    foreach ( @{$genbank_ref} ) {
        if ( m{^//} ) {         # Si es final de ORIGIN
            procesa($parte);
            $parte = '';
        }
        else {                  # Si no lo es, la guardamos
            $parte .= $_;
        }
    }

    # Fin de fichero. Vemos si nos queda algo más que procesar
    if ( $parte ) { procesa($parte) }
}

sub procesa {
    my $parte = shift;
    my @parte = split "\n", $parte;

    # Procesamiento de un fichero GenBank.
    # Tenemos que leer las partes CDS para saber los rangos,
    # y el ORIGIN, con las secuencias de nucleótidos

    my @features_cds;           # Aquí guardaremos las features
    my $origin = '';            # y aquí, el ORIGIN

    for ( my $i=0; $i < @parte; $i++ ) {        # para todas las líneas

        my $linea = $parte[$i];                 # leemos la línea

        # Encontrada una sección de Coding Sequence
        if ( $linea =~ /^     CDS/ ) {
            my %secuencia_codificante;                  # almacena una secuencia

            # Debemos interpretar las siguientes líneas
            # hasta la siguiente sección
            do {
                if    ( $linea =~ m/(\d+)\.\.(\d+)/ ) {         # rango de la secuencia
                    $secuencia_codificante{ parte  } = [ $1, $2 ];
                }
                elsif ( $linea =~ m{/gene="(.+)"} ) {           # nombre de la secuencia
                    $secuencia_codificante{ nombre } = $1        ;
                }

                $linea = $parte[++$i];          # leemos siguiente línea

            } while ( $linea =~ /^ {21}/ );     # mientras estemos dentro del CDS

            # Almacenamos la secuencia encontrada (en un array de hash)
            push @features_cds, \%secuencia_codificante;

            redo;                               # Volvemos a chequear la misma
                                                # línea (reiniciamos el bucle
                                                # sin mirar la condición del
                                                # for, no sea que venga otra
                                                # sección CDS a continuación)
        }

        # Encontrada sección ORIGIN
        elsif ( $linea =~ /^ORIGIN/ ) {

            # leemos la primera línea de nucleótidos
            $i++;
            $linea = $parte[$i];

            # para todas las líneas hasta el final
            for ( ; $i < @parte; $linea = $parte[++$i] ) {

                # División de la línea en dos partes
                my ( $offset, $nucleotidos ) = unpack "A9 A*", $linea;

                # Tratamiento de la línea
                $offset      = 0+$offset;       # Lo pasamos a número
                $nucleotidos =~ s/ //g;         # fuera espacios

                # Ver si tiene un offset numérico
                if ( $offset ) {
                    # colocamos los nucleótidos en el lugar indicado por el offset
                    substr( $origin, $offset-1 ) = $nucleotidos;
                }
                # sino, vamos añadiendo los nucleotidos al final
                else {
                    $origin .= $nucleotidos;
                }
            }
        }
    }

    # Extracción de las secuencias codificantes
    foreach my $secuencia_codificante ( @features_cds ) {

        my $nombre_de_la_secuencia = ( defined $secuencia_codificante->{nombre} )
                                   ? $secuencia_codificante->{nombre}
                                   : ''
                                   ;

        my $rango_inicio = $secuencia_codificante->{parte}->[0];
        my $rango_fin    = $secuencia_codificante->{parte}->[1];

        print qq("$nombre_de_la_secuencia"($rango_inicio..$rango_fin): );

        print substr $origin, $rango_inicio -1, $rango_fin - $rango_inicio +1;

        print "\n";
    }
}

__END__