1 de junio de 2010

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

9 comentarios:

  1. como hago para imprimir el nombre de cada gen o el lucus tag de cada gen antes de la zona intergenica de este gracias

    ResponderEliminar
  2. El código anterior ya anota el locus tag del gen anterior ($genes[$gen-1][2]) y posterior ($genes[$gen][2]) a la región intergénica.

    Si quisieras el nombre del gen, deberías anotarlo al leer el archivo GenBank de esta forma:
    $gename = join(',',sort $f->each_tag_value('gene'));

    posteriormente guardarlo en el array @genes:
    push(@genes,[$start,$end,$gi,$strand,$gename]);

    e imprimirlo junto a las regiones intergénicas:
    "length:$length|neighbours:$genes[$gen-1][4]($genes[$gen-1][3]),$genes[$gen][2]($genes[$gen][4])|$taxon\n";

    ResponderEliminar
  3. mira tengo otro incoveniente mi archivo genbank tiene 3988 genes y el solo me haya 3058 zonas intergenicas sabes por que? muchas gracias

    ResponderEliminar
  4. Si el DNA es vírico o bacteriano puede ser por genes solapados, pero por el número de genes parece que es eucariota por lo cual imagino que será debido a que en los 3988 están anotadas las variantes de splicing de los genes: http://en.wikipedia.org/wiki/Alternative_splicing

    ResponderEliminar
  5. no mi DNA es de un procariota que podria estar pasando
    aqui esta mi archivo gb

    https://skydrive.live.com/P.mvc#!/?cid=4ba9133381864506&permissionsChanged=1&id=4BA9133381864506!211

    la verdad no se por que esta pasando te agradeceria muchisimo si le encuentras la solucion a mi problema gracias

    ResponderEliminar
  6. no ya se soliciono eso no te preocupes un compañero del foro perl en español me colaboro

    http://perlenespanol.com/foro/analizar-solo-una-cadena-de-genbank-t6278.html

    pero ahora tengo un problema

    como hago para hayar la secuencia intergenica que ahi entre cadena complementaria y cadena complementaria y asu vez la que ahi entre cadena molde y cadena molde

    ResponderEliminar
  7. Pues imagino que es tan fácil como coger la secuencia completa de nucleótidos y extraer los trozos indicados por las coordenadas. Para la secuencia complementaria sería darle la vuelta una vez extraída (las coordenadas se refieren a la directa)

    ResponderEliminar
  8. dentro del genbank ahi genes expresados de la siguiente forma

    gene complement(11874..12311)
    /locus_tag="MT0010"
    /db_xref="GeneID:922447"

    como hago para que el resultado de mis zonas intergenicas omita o no cuente estos genes si no que los pase por alto solo nesesito el resultado de los que se expresan sin la palabra complement


    gene 10887..10963
    /locus_tag="MTt01"
    /db_xref="GeneID:922440"

    ResponderEliminar
  9. mira tu me dijiste
    Pues imagino que es tan fácil como coger la secuencia completa de nucleótidos y extraer los trozos indicados por las coordenadas. Para la secuencia complementaria sería darle la vuelta una vez extraída (las coordenadas se refieren a la directa)

    como podria implementar eso detro del codigo e intentado pero no he podido es que soy muy nuevo en esto y me urge saber gracias

    ResponderEliminar