• Publicidad

Cálculo de la correlación entre datos de una matriz

Perl aplicado a la bioinformática

Cálculo de la correlación entre datos de una matriz

Notapor pacoparedes » 2012-08-20 18:02 @793

Buenas tardes, estoy trabajando con una matriz de datos que tiene la siguiente estructura:

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
        Muestra 1       Muestra 2       Muestra3
Gen 1   3.92191         4.03648          4.0114
Gen 2   11.93265        12.4678          11.67605
Gen 3   4.22257         3.99756          4.17728
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Los valores numéricos indican la expresión de cada uno de los genes en cada una de las diferentes muestras, por lo que yo necesito sacar la correlación (r) existente entre la expresión del Gen 1 con el Gen 2, Gen 1 con Gen 3, Gen 2 con Gen 3, por lo que he escrito el siguiente código:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl -w
  2. #ejercicio1.pl
  3. my $arch = $ARGV[0];
  4. open( ARCH, $arch );
  5. my @array = <ARCH>;
  6. close ARCH;
  7. foreach $i (@array) {                  #cálculo de los valores de X, media y desviación estándar
  8.     chomp($i);
  9.     @array2 = split( /\t+/, $i );
  10.     $sumatoria_x = 0;
  11.     for ( $j = 1; $j < @array2; $j++ ) {
  12.         $sumatoria_x += $array2[$j];
  13.     }
  14.     $num_elem = @array2;
  15.     $media_x = $sumatoria_x / ( $num_elem - 1 );
  16.  
  17.     #print "$array2[0]\t$media_x\n";
  18.     $sumatoria_x = 0;
  19.     for ( $k = 1; $k < @array2; $k++ ) {
  20.         $sumatoria_x += ( $array2[$k] - $media_x )**2;
  21.     }
  22.     $des_stnd_x = sqrt( $sumatoria_x / ( $num_elem - 2 ) )
  23.         ;                              #le resto dos por la presencia del identificador en la posc. 0 y 1 q indica la fórmula n-1.
  24.  
  25.     #print "$array2[0]\t$media\t$des_stnd_x\n";
  26.  
  27.     foreach $l (@array) {              #cálculo de los valores de Y, media y desviación estándar
  28.         chomp($l);
  29.         @array3 = split( /\t+/, $l );
  30.         $sumatoria_y = 0;
  31.         for ( $m = 1; $m < @array3; $m++ ) {
  32.             $sumatoria_y += $array3[$m];
  33.         }
  34.         $num_elem = @array3;
  35.         $media_y = $sumatoria_y / ( $num_elem - 1 );
  36.  
  37.         #print "$array2[0]\t$media_y\n";
  38.         $sumatoria_y = 0;
  39.         for ( $n = 1; $n < @array3; $n++ ) {
  40.             $sumatoria_y += ( $array3[$n] - $media_y )**2;
  41.         }
  42.         $des_stnd_y = sqrt( $sumatoria_y / ( $num_elem - 2 ) )
  43.             ;                          #le resto dos por la presencia del identificador en la posc. 0 y 1 q indica la fórmula n-1.
  44.  
  45.         #print "$array2[0]\t$media\t$des_stnd_y\n";
  46.         $suma_cov = 0;
  47.         for ( $o = 1; $o < @array2; $o++ ) {    #utilizo una sola variable por que ambos arreglos son del mismo tamaño
  48.             $suma_cov += ( $array2[$o] - $media_x ) * ( $array3[$o] - $media_y );
  49.         }
  50.         $cov         = $suma_cov / ( $num_elem - 1 );
  51.         $correlacion = $cov /      ( $des_stnd_y * $des_stnd_x );
  52.         $co = sprintf( "%.2f", $correlacion );
  53.         print "$array2[0]\t$array3[0]\t$co\n";
  54.     }
  55. }
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4



Este código funciona bien, sin embargo calcula la correlación existente entre
Gen 1-Gen 1, Gen 1-Gen 2, Gen 1-Gen 3,
Gen 2-Gen 1, Gen 2-Gen 2, Gen 2-Gen 3, etc...

Pero como podemos notar el cálculo de la correlación existente entre Gen 1-Gen 2, es lo mismo que el cálculo entre Gen 2-Gen 1, es decir, me está repitiendo cálculos.

Es por eso que pido su ayuda para optimizar mi código y evitar estas situaciones.

Gracias...
Última edición por explorer el 2012-08-20 18:15 @802, editado 1 vez en total
Razón: Formateado de código con Perltidy y poner marcas Perl
pacoparedes
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2012-08-10 17:03 @752

Publicidad

Re: Cálculo de la correlación entre datos de una matriz

Notapor explorer » 2012-08-22 15:34 @690

Según los cálculos que muestras, creo que lo que estás intentando hacer es calcular los coeficientes de correlación de Pearson. ¿Es así?

En ese caso, creo que el error está en la línea 27, ya que si suponemos que los valores x son los leídos en el bucle de la línea 7 (las filas), en el 27 los valores y deberían ser los de las columnas (las muestras), pero en lugar de eso estás volviendo a leer las líneas.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Cálculo de la correlación entre datos de una matriz

Notapor pacoparedes » 2012-08-24 17:48 @783

Asi es, explorer. Estoy tratando de calcular el coeficiente de correlación de Pearson. Revisaré tus observaciones. Gracias.
pacoparedes
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2012-08-10 17:03 @752

Re: Cálculo de la correlación entre datos de una matriz

Notapor pacoparedes » 2012-08-25 19:36 @858

Buenas tardes explorer. Revisé mi código y lo que hice fue sustituir los foreach de las líneas 7 y 27 por for, y en la condición de inicio del for para el cálculo de Y utilicé $l=$i+1 y está funcionando bien.

Gracias...
pacoparedes
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2012-08-10 17:03 @752

Re: Cálculo de la correlación entre datos de una matriz

Notapor explorer » 2012-08-28 18:05 @795

Esta es mi versión, aunque no he comprobado si los datos que salen son correctos...
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/env perl
  2. use Acme::Tools;                                # Caja de herramientas marca ACME
  3. use Data::ShowTable;
  4.  
  5. ### Leemos el archivo
  6. my @filas;
  7. readfile('code_32086.txt', \@filas);
  8.  
  9. ### Procesamos las cabeceras de las columnas
  10. my @nombres_columnas = split /\t+/, shift @filas;
  11.  
  12. ### Procesamos las cabeceras de las filas
  13. @filas = map { [ split /\t+/, $_ ] } @filas;    # conversión a array de array
  14. my @nombres_filas = map { shift @$_} @filas;
  15.  
  16. ### Extraemos las columnas
  17. my @columnas;
  18. for my $i (0 .. $#{$filas[0]}) {
  19.     $columnas[$i] = [ map { $_->[$i] } @filas ];
  20. }
  21.  
  22. ### Comprobación
  23. die "ERROR: no coinciden el número de filas con el de columnas.\n"
  24.     if @filas != @columnas;
  25.  
  26. ### Cálculo de las medias
  27. my @medias_filas    = map { avg @$_    } @filas;
  28. my @medias_columnas = map { avg @$_    } @columnas;
  29.  
  30. ### Cálculo de las desviaciones estándar
  31. my @desv_filas      = map { stddev @$_ } @filas;
  32. my @desv_columnas   = map { stddev @$_ } @columnas;
  33.  
  34. ### Recorrer la tabla, calculando la correlación
  35. my @correlaciones;
  36.  
  37. for my $i (0 .. $#filas) {
  38.     for my $j (0 .. $#columnas) {
  39.         my $suma_cov = 0;
  40.         for my $k (0 .. $#filas) {
  41.             $suma_cov += ( $filas   [$i][$k]  -  $medias_filas   [$i] )
  42.                         *
  43.                          ( $columnas[$j][$k]  -  $medias_columnas[$j] )
  44.             ;
  45.         }
  46.  
  47.         my $cov         = $suma_cov / @filas;
  48.         my $correlacion = $cov / ( $desv_filas[$i] * $desv_columnas[$j] );
  49.        
  50.         $correlaciones[$i][$j] = $correlacion;
  51.     }
  52. }
  53.  
  54. ### Presentación de la tabla
  55. my $corr_i = 0;
  56.  
  57. ShowTable({
  58.     titles      => \@nombres_columnas,
  59.     types       => [ 'titulo', ('valor') x @columnas ],
  60.     widths      => [ 7,        ("12.2")  x @columnas ],
  61.     row_sub     => sub {
  62.                     my $rebobinar = shift;
  63.                    
  64.                     if ($rebobinar == 0) {              # nos piden una fila
  65.                         if  ($corr_i < @correlaciones) {
  66.                             my $fila_ref = $correlaciones[$corr_i];
  67.                             my @row = ($nombres_filas[$corr_i], @{$fila_ref});
  68.                             $corr_i++;
  69.                             return @row;
  70.                         }
  71.                         else {
  72.                             return;                     # fin de tabla
  73.                         }
  74.                     }
  75.                     else {                              # rebobinar
  76.                         $corr_i = 0;
  77.                         return 1;
  78.                     }
  79.                 },
  80.     fmtsub      => sub {                                # Formateo de cada celda
  81.                     my($valor, $tipo, $max_width, $width, $precision) = @_;
  82.                     if ($tipo eq 'titulo') {
  83.                         return sprintf "%*s", $width, $valor;
  84.                     }
  85.                     else {
  86.                         return sprintf "%*.*f", $width, $precision, $valor;
  87.                     }
  88.                 },
  89.     max_width   => 80,
  90.     show_mode   => 'Box',
  91. });
  92.  
  93. __END__
  94.  
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
 +-------+-----------+-----------+----------+
 |       | Muestra 1 | Muestra 2 | Muestra3 |
 +-------+-----------+-----------+----------+
 | Gen 1 |      0.46 |      0.44 |     0.46 |
 | Gen 2 |      0.62 |      0.63 |     0.63 |
 | Gen 3 |     -0.66 |     -0.65 |    -0.66 |
 +-------+-----------+-----------+----------+
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España


Volver a Bioinformática

¿Quién está conectado?

Usuarios navegando por este Foro: No hay usuarios registrados visitando el Foro y 0 invitados

cron