Submit gene with unknown intron to GenBank

Hola de nuevo,
el 31 de diciembre conseguí finalmente enviar a GenBank unas secuencias parciales de genes de cebada utilizadas por mi colega Ana Casas en un estudio. Éste es un paso necesario para publicar en casi cualquier revista seria, pero además es la manera de asegurar que tus secuencias van a ser útiles para otras personas en el futuro.

Antes de que se me olvidé comparto aquí como lo hice, teniendo en cuenta que uno de los genes tiene un intrón de longitud desconocida y se secuenció en dos amplicones que abarcan respectivamete los exones 1-9 y 10-14:
 

Para poder archivar una secuencia así en GenBank utilice tbl2asn en varios pasos:

1.  Obtención de una plantilla .sbt con https://submit.ncbi.nlm.nih.gov/genbank/template/submission . El fichero resultante especifica los autores y otros metadatos, y lo puedes usar para distintas secuencias. 
2. Confección de un fichero FASTA con extensión .fsa (unk_intron_gene.fsa) que contenga las secuencias de ambos amplicones separadas por un tramo de 100 Ns. Para facilitar el siguiente paso intrones y exones deben ir en mayúsculas y minúsculas, respectivamente, como en este fragmento que comprende los exones 8, 9 y 10:
   

   ...TGTCAGatactatgcaattgccacaccaagtgctacacaaagattgctttttggtct
   TGTCAGatactatgcaattgccacaccaagtgctacacaaagattgctttttggtcttct
   tgaagcaccaccatcatgggctccagatgcacttgatgcagcagttcagcttgttgaact
   ccttcgggcagctgaagattatgctactggcatgcggGTATGACATACTGCATGCTGGCT
   GTTGTTTCAGTCCTGTTAGTTGTGATGCCTCACGATACAAAATTTCCATATTCGTATGTT
   TTGGGTGTGCATGTTTATTAATCTTGGTAACTTTAAATTCCTGTTCAGcttccaaaaaat
   tggttgcatcttcatttcttgcgtgcgattggaactgcaatgtctatgagggctggtatt
   gctgccgatacagctgctgcgttgctttttcgcatactatcccaaccaacgttgcttttt
   cctccactaaggcatgctgaaggagttgaagtgcaacatgaaccactgggtggctatgta
   tcatcatacaaaagacagGTATGCAGTAGTTTCTGCATCTAGTTAATTTTTCATTATCTG
   TTCTTCTTTAGTAAAGACTCAANNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
   NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
   NNGGATCCATGTTTTAGTCTTCTTGGTTTTACTGATTGTTGCCTTATGTCTGCATGACTA
   ATTTACCTGCTTGCACTTTGAACTATTCACAGctggaagttcctgcatctgaaaccacaa
   ttgatgccactgcacaaggcattgcttccttgctgtgtgctcatggtcctgatgttgagt
   ggagaatatgtaccatctgggaagctgcctatggtttgttacctctgaattcatcagcag
   ttgatttgcccgaaatcgttgtagctgctccgcttcagccacctactttgtcatggagcc
   tatacttgccactgttgaaagtattcgagtatctacctcgtggaagtccatctgaagcat
   gccttatgagaatatttgtggcaacagttgaagctatactcagaagaactttcccttcgg
   aaacctctgaatcatctaaaagaccaagaagtcaatccaagaaccttgctgttgctgaac
   tccgtacaatgatacattcactctttgttgaatcatgtgcttcaatgaaccttgcttccc
   ggttgttgtttgttgtattaactgtttgcgtcagtcatcaagctttgccagggggcagca
   aaagaccaacgggtagtgaaaaccattcttctgaggaggccactgaggacccaagattaa
   ccaatggaagaaataaggtcaagaagaaacaagggcctgttggtacatttgactcgtatg
   tgctggctgctgtttgtgccttatcttgtgagcttcagctgttccctatcctttgcaaga
   gtgcaacaaactcaaaagtaaaagactctataaagatcctgaagcctggaaaaaacaatg
   ggatcagtaatgagctacagaatagcattagctcagcaattctccatactcgtagaattc
   ttggcatcctggaagctcttttctccttgaagccatcatcagttggtacctcctggaact
   atagttcaaatgagatagttgcagcggctatggttgccgctcatgtttctgagttatttc
   gccggtcgaggccatgcctaaatgcactatcttcactgaagcgatgtaagtgggatgctg
   agatttctaccagggcatcatccctttaccatttgatcgatttgcatggtaaaactgtgt
   cctccatcgtgaacaaagctgagcctctagaagctcacctgacttttacatcagtaaaga
   gagatggtcaacaacacattgaggaaaacagcaccagctcatcgggtaatggcaacttgg
   aaaagaagaatgcttcagcctcacacatgaaaaatggtttttcaagaccactcttgaaat
   gctcagaagaggctaggcgaaatggtaatgttgcaagtacatccgggaaagttcctgcaa
   ctttacaggctgaagcatctgatttggctaacttccttaccatggatagaaatgggggtt
   atcgaggctctcagactctcctaagttctgttatctcagaaaaacaggaattatgcttct
   ctgttgtctcattgctctggcataagcttattgcatctcctgaaacgcagatgtctgcag
   aaagtacatcagctcatcaaggttggagaaagGTA...

3. Obtención de las coordenadas de los exones con ayuda del siguiente script Perl que usa variables especiales de expresiones regulares:

   perl -lne 'if(/^>/){} else{ while(/[a-z]+/g){ printf("%d\t%d\n",$-[0]+1,$+[0])} }' unk_intron_gene.fsa

4. Prepara y completa un fichero .tbl (feature table) por gen. Verás que los campos van separados por tabuladores. Si es un gen parcial deberás indicarlo usando ">" en las coordenadas de inicio y fin, tal como se explica en https://www.ncbi.nlm.nih.gov/Sequin/table.html . Para el gen ejemplo yo obtuve el siguiente .tbl:

   >Feature HvGI Table1
   <1    >7751 gene
       gene HvGI
   1383 1469 mRNA
   1564 1641
   1738 1792
   1881 2023
   2120 2376
   2811 2869
   3079 3319
   3390 3540
   3672 3881
   4821 6320
   6420 6488
   6636 6854
   7014 7106
   7195 7500
   1383 1469 CDS
   1564 1641
   1738 1792
   1881 2023
   2120 2376
   2811 2869
   3079 3319
   3390 3540
   3672 3881
   4821 6320
   6420 6488
   6636 6854
   7014 7106
   7195 7500
       product HvGI
   

5. Guarda en una carpeta propia los ficheros .fsa y .tbl, un par por gen, por ejemplo genes/
6. Descarga del binario tbl2asn adecuado para tu sistema operativo desde ftp://ftp.ncbi.nih.gov/toolbox/ncbi_tools/converters/by_program/tbl2asn
7. Haz el binario ejecutable si fuera necesario 
8. Haz la conversión, por ejemplo en Linux:

   linux.tbl2asn -t template.sbt -p genes/ -V vb -a r10u

9. Comprueba errores (errorsummary.val) y corrige los respectivos ficheros .tbl y vuelve al paso 8.
10. Envía los ficheros .sqn resultantes por medio de http://www.ncbi.nlm.nih.gov/LargeDirSubs/dir_submit.cgi
     

Hasta pronto,
Bruno

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