Foro - Perl en Español

¡Hola!

Tengo un problema: necesito extraer la secuencia intergénica de solo una cadena de genbank. En un foro anterior lo vi con ambas cadenas, ahora necesito extraerlas pero solo en una cadena, bien sea la cadena molde, o la complementaria.

Espero contar con su colaboración.

Gracias

Este tema ya está comentado en otros hilos de este foro:

* ADN y partes codificantes

En uno de mis primeros mensajes hay enlaces a otros hilos donde hay ejemplos de cómo leer ficheros GenBank, extraer la secuencia genómica, interpretar los CDS, etc.

Mira, tengo este código que haya las secuencias intergénicas pero resulta que mi archivo genbank tiene 3988 genes y él solo me haya 3058 zonas intergénicas. ¿Sabes por qué? Muchas gracias.


use Bio::SeqIO;

# prog1.pl
# recibe hasta 4 argumentos:
# 1) archivo entrada en formato GenBank
# 2) nombre del archivo de salida en formato FASTA
# 3) longitud mínima de las regiones intergénicas (opcional, por defecto toma todas)
# 4) longitud máxima de las regiones intergénicas (opcional)

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

# calcular ORFs vecinos y corta las secuencias intergenicas
for($gen=1;$gen<scalar(@genes);$gen++)
{
next if($genes[$gen-1][1] >= $genes[$gen][0]); #no solapa asume genes en orden
next if( defined($min_intergenic_size) && $min_intergenic_size > 0 &&
$genes[$gen][0]-$genes[$gen-1][1]+1 < $min_intergenic_size); # evita cortas
next if( defined($max_intergenic_size) && $max_intergenic_size > 0 &&
$genes[$gen][0]-$genes[$gen-1][1]+1 > $max_intergenic_size); # evita largas

$n_of_intergenic++;

# calcula bordes de intergenes
$start = $genes[$gen-1][1] + 1;
$end = $genes[$gen][0] - 1;
$length = $end - $start + 1;

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-1,$length)."\n";
}
}

close(FNA);

return $n_of_intergenic;
}

&extract_intergenic_from_genbank;

Pues no, no lo sé.

Prueba a enviarme el fichero, adjuntándolo en un mensaje privado (menú de arriba, debajo de la palabra "Foro").

¡Hola!
¡Qué pena no poder enviar el fichero en un mensaje privado! Lo que paso es que es muy grande y no subió el archivo, al ver eso lo que hice fue subirlo a skydrive y de allí lo puedes descargar si no es mucha molestia. Qué pena y aquí esta el link

https://skydrive.live.com/P.mvc#!/?cid= ... 864506!211

Ahí aparece mi fichero genbank ¡gracias!

El problema está en que la subrutina que has copiado del blog #!/perl/bioinfo, en la línea 68, descarta procesar cualquier secuencia genética que se solape al gen anterior (el número de nucleótico final del segmento anterior es mayor o igual al número inicial del nucleótido actual).

Del fichero que me has pasado, la subrutina descarta 929 genes que se solapan entre sí.

Los puedes ver si cambias la línea 68 y añades estas líneas:

Y aún así falta un gen más (3058+929=3987=3988-1): el primero, pues el fichero de salida comienza pory en el fichero vemos que hay uno de 1524 bases (del 1 al 1524).

¡Ok! Muchísimas gracias. Me funcionó a las mil maravillas. ¡Gracias!

¿Cómo podría saber cuántos son los genes que solapan e imprimir por pantalla el valor o la cantidad de genes que solaparon?

Puedes probar a grabarlo en otro fichero de salida...

Sí, esa salida ya la tengo en otro fichero. Lo que necesito es el número total de genes que se solaparon. Intenté con un contador pero no me funcionó.