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:
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
