Extraer secuencias intergénicas de archivos GenBank

La tarea de hoy consiste en leer un archivo GenBank, por ejemplo de un genoma circular bacteriano, con el fin de extraer las secuencias intergénicas. Para esta tarea vamos a utilizar código de BioPerl y antes necesitaremos descargar al menos un fichero GenBank, como se explica por ejemplo en el taller de (bio)perl. Una vez tengamos los archivos de entrada necesarios, podemos resolver el problema de las secuencias intergénicas con una subrutina como ésta:

1:  use Bio::SeqIO;
2:
3:  # prog1.pl
4:  # recibe hasta 4 argumentos:
5:  # 1) archivo entrada en formato GenBank
6:  # 2) nombre del archivo de salida en formato FASTA
7:  # 3) longitud mínima de las regiones intergénicas (opcional, por defecto toma todas)
8:  # 4) longitud máxima de las regiones intergénicas (opcional)
9:
10:  sub extract_intergenic_from_genbank
11:  {
12:    my ($infile,$out_intergenic_file,$min_intergenic_size,$max_intergenic_size) = @_;
13:
14:    my ($n_of_intergenic,$gi,$start,$end,$length,$strand,$taxon) = (0);
15:    my $in = new Bio::SeqIO(-file => $infile, -format => 'genbank' );
16:
17:    open(FNA,">$out_intergenic_file") ||
18:      die "# extract_intergenic_from_genbank : cannot create $out_intergenic_file\n";
19:
20:    while( my $seq = $in->next_seq) # scan all sequences found in $input
21:    {
22:      my ($gbaccession,$sequence,$gen,@genes) = ( $seq->accession() );
23:      $sequence = $seq->primary_seq()->seq() ||
24:         'no encuentro la secuencia primaria, necesito un fichero GenBank completo!';
25:      $taxon = '';
26:
27:      for my $f ($seq->get_SeqFeatures)
28:      {
29:        if($f->primary_tag =~ /CDS|rRNA|tRNA/) # campos de 'genes'
30:        {
31:          $gi = '';
32:          if($f->has_tag('db_xref'))
33:          {
34:            my $crossrefs = join(',',sort $f->each_tag_value('db_xref'));
35:            if($crossrefs =~ /(GI\:\d+)/){ $gi = $1 }
36:          }
37:          elsif($f->has_tag('locus_tag'))
38:          {
39:            if($gi eq '' && $f->has_tag('locus_tag'))
40:            {
41:              $gi = "ID:".join(',',sort $f->each_tag_value('locus_tag'));
42:            }
43:          }
44:
45:          next if($gi eq '');
46:
47:          $start = $f->start();
48:          $end  = $f->end();
49:          $strand = $f->location()->strand();
50:          push(@genes,[$start,$end,$gi,$strand]);
51:        }
52:        elsif($f->primary_tag() =~ /source/)
53:        {
54:          if($f->has_tag('organism'))
55:          {
56:            foreach my $element ($f->each_tag_value('organism'))
57:            {
58:              $taxon .= "[$element],";
59:            }
60:            chop $taxon;
61:          }
62:        }
63:      }
64:
65:      # calcular ORFs vecinos y corta las secuencias intergenicas
66:      for($gen=1;$gen<scalar(@genes);$gen++)
67:      {
68:        next if($genes[$gen-1][1] >= $genes[$gen][0]); #no solapa asume genes en orden
69:        next if( defined($min_intergenic_size) && $min_intergenic_size > 0 &&
70:          $genes[$gen][0]-$genes[$gen-1][1]+1 < $min_intergenic_size); # evita cortas
71:        next if( defined($max_intergenic_size) && $max_intergenic_size > 0 &&
72:          $genes[$gen][0]-$genes[$gen-1][1]+1 > $max_intergenic_size); # evita largas
73:
74:        $n_of_intergenic++;
75:
76:        # calcula bordes de intergenes
77:        $start = $genes[$gen-1][1] + 1;
78:        $end  = $genes[$gen][0] - 1;
79:        $length = $end - $start + 1;
80:
81:        print FNA ">intergenic$n_of_intergenic|$gbaccession|coords:$start..$end|".
82:          "length:$length|neighbours:$genes[$gen-1][2]($genes[$gen-1][3]),".
83:          "$genes[$gen][2]($genes[$gen][3])|$taxon\n";
84:        print FNA substr($sequence,$start-1,$length)."\n";
85:      }
86:    }
87:
88:    close(FNA);
89:
90:    return $n_of_intergenic;
91:  }