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

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;

¡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!

¡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?

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ó.