Ver Tema - Contar y agrupar aminoácidos de diferentes TF

zipetardix · 2012-02-04 16:51 @744

Hola a todos.

Tengo un problema con un script que tengo que hacer en Perl. Llevo muchos días intentando solucionarlo pero no hay manera. Ahí va:

Tengo una lista en .txt formada por una columna de tipo de factor de transcripción (TF) y otra formada por sus aminoácidos. ej:

Syntax: [ Download ] [ Hide ]

Using text Syntax Highlighting

TF1     ATLHLSNIPPSVTMDDLKNLFEAGCSVKAFKFFQKDRKMALIQLGSVEEAIQALIELHNHDLGENHHLRVSFSKS

TF1     ATLHLSNIPPSVTVDDLKNLFEAGCSVKAFKFFQKDRKMALIQLGSVEEAIQALIELHNHDLGENHHLRVSFSKS

TF1     ATLHLSNIPQNVTEDDLRLLFNTGGTVKAFKFFQDRKMALIQMSTVEEAIQALIDLHNYNMGGNQHLRVSFSKS

TF2     GKMFVGGLSWDTSKKDLKDYFTKFGDVSDCTIKFDSDSGRSRGFGFVLFKDPSSVEKVLQQKEHRLDGRLIDPKRAMAMKKEPV

TF2     GKMFVGGLSWDTSKKDLKDYFTKFGEVSDCTIKFDSDSGRSRGFGFVLFKDPSSVKKVLEQKEHRLDGRLIDPQRAMAMKKEPV

TF3     GKMFVGGLSWDTSKKDLKDYFTKFGEVTDCTIKMDPNTGRSRGFGFILFKEPGSVEKVLEQKEHRLDGRLIDPKKAMAMKKDPVKK

[etc]

Lo que necesito es simple: tengo que conseguir que me lea los aminoácidos de cada TF, distinga qué TF es y me los agrupe en ese TF cada a.a. con su frecuencia. Para ello sé que tengo que trabajar con hashes y eso. Mi problema es que de Perl controlo más bien poco, pero bueno, ahí va lo que llevo escrito hasta ahora:

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

#!/usr/bin/perl -w
 
use strict;
 
 
my $arg0 = $ARGV[0];
if ( !$arg0 ){
   print STDERR "Usage: $0 file_name\n";
   exit(0);
}
 
 
open(FILE,"<$arg0") or die("cannot open $arg0");
 
while (<FILE>){
    chomp;
    
    my ($id, $seq) = split("\t", $_);
    my $i=0;
    my %TF1;
    if ( $id eq 'TF1') {
         
        
        my @TF1_aa = split(//,$seq);
 
#       for(my $i=0; $i < scalar(@TF1_aa); $i++){
#        $TF1{$TF1_aa[$i]}++;     
# esto es lo que nos han dado como ejemplo para trabajar pero con esto el hash se llena con TODOS los aminoácidos de la lista sin distinguir el tipo de TF. Así pues, cambié a if/elsif pero dentro de estas condiciones veo que no puedo declarar un for y por ello pasé a while
 
        while ($i < scalar(@TF1_aa)){ $i++;                            
        $TF1{$TF1_aa[$i]}++};
 
                        my @aa  = keys( %TF1 );
                        my @aa_count  = values( %TF1 );
                        print "@aa\n";
                        print "@aa_count\n";
# estos print los pongo para poder ver cómo está trabajando el script y si está haciendo lo que quiero o no
 
 
    } elsif ($id eq 'TF2') {
 # aquí repito lo anterior pero cambiándolo por $TF2 y así elsif con todos los tipos que tengo en la lista
 

Vale, pues ejecutando este script me da el siguiente output:

Syntax: [ Download ] [ Hide ]

Using text Syntax Highlighting

Use of uninitialized value within @TF1_aa in hash element at test.pl line 31, <FILE> line 1.

 S F T N K E Y V Q M C L A P D R I G

1 6 6 4 2 5 5 4 6 2 1 1 6 5 2 7 7 3 6

Use of uninitialized value within @TF1_aa in hash element at test.pl line 31, <FILE> line 2.

 F S T N K E Y V Q L A W P D R I G

1 4 3 2 4 2 7 3 7 1 4 8 1 3 5 8 2 9

Use of uninitialized value within @TF1_aa in hash element at test.pl line 31, <FILE> line 3.

 F S T K E Y V Q C L A W P D R I G

1 4 4 2 2 6 2 9 1 1 4 7 1 3 7 10 3 7

Use of uninitialized value within @TF1_aa in hash element at test.pl line 31, <FILE> line 4.

 F S T N K E Y V Q L A P H D I R G

1 5 3 1 2 5 6 3 5 1 6 2 1 2 8 5 7 8

[...]

y así hasta que termina con todas las líneas en las que $id eq 'TF1'.

Pues eso, quiero que me las agrupe todas dentro del hash %TFn y no me vaya dando el resultado una a una sin almacenarlas.

Muchas gracias por adelantado.

explorer · 2012-02-04 20:59 @916

Bienvenido a los foros de Perl en español, zipetardix.

Pues... la pista que te han dado es muy buena... es más: es justo lo que puedes usar para resolver el problema:

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

#!/usr/bin/perl
use strict;
use warnings;
 
@ARGV == 1  or  die "Usage: $0 <archivo TF>\n";
 
my $archivo = $ARGV[0];
 
 
# Estructura donde guardaremos la información
# Cada clave de %TF corresponde a un TF
# y su valor correspondiente es una referencia a un hash,
# cuyas claves serán los aminoácidos, y sus valores, las veces que aparecen.
my %TF;
 
open(my $FILE, '<', $archivo) or die "cannot open $archivo: $!\n";
 
while (<$FILE>) {
    chomp;
    
    my ($id, $seq) = split;                             # (¡magia!)
    
    my @TF_aa = split(//, $seq);
 
    for my $aa (@TF_aa) {
        $TF{$id}{$aa}++;                                # ¡Ta chán! 
    }
}
 
close $FILE;
 
for my $id (sort keys %TF) {
    for my $aa (sort keys %{$TF{$id}}) {
        print "[$id][$aa][$TF{$id}{$aa}]\n";
    }
}

sale:

Syntax: [ Download ] [ Hide ]

Using text Syntax Highlighting

[TF1][A][17]

[TF1][C][2]

[TF1][D][12]

[TF1][E][13]

[TF1][F][15]

[TF1][G][10]

[TF1][H][13]

[TF1][I][12]

[TF1][K][16]

[TF1][L][29]

[TF1][M][6]

[TF1][N][14]

[TF1][P][5]

[TF1][Q][11]

[TF1][R][7]

[TF1][S][19]

[TF1][T][9]

[TF1][V][13]

[TF1][Y][1]

[TF2][A][4]

[TF2][C][2]

[TF2][D][19]

[TF2][E][7]

[TF2][F][14]

[TF2][G][16]

[TF2][H][2]

[TF2][I][4]

[TF2][K][24]

[TF2][L][12]

[TF2][M][6]

[TF2][P][6]

[TF2][Q][4]

[TF2][R][10]

[TF2][S][16]

[TF2][T][6]

[TF2][V][12]

[TF2][W][2]

[TF2][Y][2]

[TF3][A][2]

[TF3][C][1]

[TF3][D][8]

[TF3][E][5]

[TF3][F][6]

[TF3][G][9]

[TF3][H][1]

[TF3][I][3]

[TF3][K][15]

[TF3][L][6]

[TF3][M][4]

[TF3][N][1]

[TF3][P][4]

[TF3][Q][1]

[TF3][R][4]

[TF3][S][4]

[TF3][T][5]

[TF3][V][5]

[TF3][W][1]

[TF3][Y][1]

Es decir: usamos una estructura bidimensional para guardar toda la información. En las "filas" de esa estructura están los TF, mientras que las "columnas" son los aminoácidos.

zipetardix · 2012-02-05 05:03 @252

Pues muchísimas gracias. A veces uno empieza a complicarse él solo y acaba con un cacao importante.

Ahora para el siguiente paso mi pregunta es muy sencilla pero no me sale la respuesta por mi mismo: ¿Cómo puedo hacer para sacar la frecuencia relativa de cada aminoácido dentro de cada TF? Algo tal que así:

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

        print STDOUT $TF{$id}{$aa}/ $totaldeaadelTF;
 

EDITO:

¿Sería posible asignar como valores dichas frecuencias relativas? Esto es debido a que en realidad esos son los datos con los que he de trabajar después y entonces si pudiera acceder a cada uno de ellos como una "coordenada" del hash para posteriores operaciones me sería totalmente salvador...

Entiendo que tal vez sea ya pedir demasiado, pero lo pregunto porque quizá sea más sencillo de lo que me parece....

explorer · 2012-02-05 06:52 @327

El número total de aminoácidos por cada TF, es, justamente, scalar(@TF_aa), así que te vale con tener un a variable sumatorio antes:

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

my %totales_aa;

y luego ir sumando:

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

    $totales_aa{$id} += @TF_aa;

Y para guardar las frecuencias, puedes hacerlo en el mismo lugar en que guardas el número de aminoácidos. Pero esto hay que hacerlo una vez leídos todos los TF, así que, por ejemplo, puedes hacerlo todo en el último bucle:

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

for my $id (sort keys %TF) {
    my $n = $totales_aa{$id};
    for my $aa (sort keys %{$TF{$id}}) {
        my $a = $TF{$id}{$aa};
        my $m = $a/$n;
        print "[$id][$aa][$a][$n][$m]\n";
        $TF{$id}{$aa} = $m;
    }
}

De forma más abreviada,

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

for my $id (sort keys %TF) {
    my $n = $totales_aa{$id};
    for my $aa (sort keys %{$TF{$id}}) {
        $TF{$id}{$aa} /= $n;
    }
}

zipetardix · 2012-02-05 08:13 @384

Pues muchísimas gracias por todo, me funciona perfectamente. Voy a ver si puedo seguir yo solo desde aquí... Ahora tengo que, de diferentes secuencias prueba de aminoácidos, predecir a qué grupo de TF pertenecería multiplicando cada aminoácido que aparezca por la frecuencia obtenida en cada uno de los distintos tipos de TF, obteniendo un valor diferente para uno y quedándome con el más alto. ¡A ver si hay suerte!

De todas formas insisto en agradecer toda la ayuda obtenida hasta ahora, me ha salvado.

Publicidad

Socializa