• Publicidad

Alineamiento múltiple

Perl aplicado a la bioinformática

Alineamiento múltiple

Notapor Marne » 2014-01-13 10:16 @470

Buenas tardes. Llevo una semana dándole vueltas a la cabeza y como soy biólogo y no informático la verdad es que me encuentro bastante perdido, a ver si me pudierais echar un cable.

El programa que intento hacer abre un fichero con un máximo de 8 cadenas que hay que meter previamente en un array. Luego la segunda función del programa se supone que coge de 2 en 2 esas cadenas, siempre la 1 con la 2, la 1 con la 3, la 1 con la 4, etc...

Por algún motivo no me deja subiros el fichero de las cadenas con lo que he hecho una pequeña modificación al programa para que en vez de abrir un fichero lea unas cadenas que he escrito yo al azar.

El caso es que en vez para alinearme las cadenas siempre me devuelve un return de 0, seguro que he cometido más de un error, pero el caso es que me falta el conocimiento para determinar dónde reside...

El programa debería alinear y meterme los resultados en una matriz, y una vez en ella, debería comprobar todas las columnas y hacerme un print con el número de posiciones de la matriz donde coincide en todas las cadenas el mismo carácter (nucleótido).

Si pudierais ayudarme en cualquiera de mis dudas, que siento que sean tantas... os estaría eternamente agradecido.

El código es este:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #################################################################
  2. ## Programa para alinear por pares 2 secuencias introduciendo  ##
  3. ## gaps para optimizar el número de coincidencias              ##
  4. #################################################################
  5.  
  6. #my @lista_cad;
  7. #my $nombrearchivo= "datos8-500";
  8.  
  9. #sub abrir_archivo {my ($nombrearchivo) = shift;
  10. #       unless(open(Archivo, $nombrearchivo)) {
  11. #           print "Se ha detectado un error al abrir archivo, \n compruebe que el nombre del archivo y el formato son correctos,";
  12. #           exit;
  13. #       }
  14. #       @lista_cad = <Archivo>;
  15. #       close Archivo;
  16. #       return @lista_cad;
  17. #    }
  18.  
  19. my @lista_cad = ( "AGTGTA", "AGTCTA", "AGTCAA" );
  20. my @pareja;
  21. for ( my $i = 1; $i < int(@lista_cad); $i++ ) {
  22.     @pareja = emparejamelas( $array[0], $array[$i] );
  23.     print join( "\n", @pareja );
  24.  
  25. }
  26.  
  27. sub emparejamelas {
  28.     my ( $sec1, $sec2 ) = @ARGV;
  29.  
  30.     # Marcador del esquema
  31.     my $basecoincide   = 1;            # +1 si las bases están emparejadas
  32.     my $basenocoincide = -1;           # -1 si las bases están desparejadas
  33.     my $GAP            = -1;           # -1 para cualquier gap
  34.  
  35.     # Inicializa
  36.     #abrir_archivo($nombrearchivo);
  37.     my @lista_cad;
  38.     $lista_cad[0][0]{score}   = 0;
  39.     $lista_cad[0][0]{pointer} = "none";
  40.     for ( my $j = 1; $j <= length($sec1); $j++ ) {
  41.         $lista_cad[0][$j]{score}   = $GAP * $j;
  42.         $lista_cad[0][$j]{pointer} = "left";
  43.     }
  44.     for ( my $i = 1; $i <= length($sec2); $i++ ) {
  45.         $lista_cad[$i][0]{score}   = $GAP * $i;
  46.         $lista_cad[$i][0]{pointer} = "up";
  47.     }
  48.  
  49.     # Rellenar
  50.     for ( my $i = 1; $i <= length($sec2); $i++ ) {
  51.         for ( my $j = 1; $j <= length($sec1); $j++ ) {
  52.             my ( $marc_diagonal, $marc_izq, $marc_up );
  53.  
  54.             # Calcula el marcador de emparejamiento
  55.             my $letter1 = substr( $sec1, $j - 1, 1 );
  56.             my $letter2 = substr( $sec2, $i - 1, 1 );
  57.             if ( $letter1 eq $letter2 ) {
  58.                 $marc_diagonal = $lista_cad[ $i - 1 ][ $j - 1 ]{score} + $basecoincide;
  59.             }
  60.             else {
  61.                 $marc_diagonal = $lista_cad[ $i - 1 ][ $j - 1 ]{score} + $basenocoincide;
  62.             }
  63.  
  64.             # Calcula el marcador de gaps
  65.             $marc_up  = $lista_cad[ $i - 1 ][$j]{score} + $GAP;
  66.             $marc_izq = $lista_cad[$i][ $j - 1 ]{score} + $GAP;
  67.  
  68.             # Elige el mejor resultado
  69.             if ( $marc_diagonal >= $marc_up ) {
  70.                 if ( $marc_diagonal >= $marc_izq ) {
  71.                     $lista_cad[$i][$j]{score}   = $marc_diagonal;
  72.                     $lista_cad[$i][$j]{pointer} = "diagonal";
  73.                 }
  74.                 else {
  75.                     $lista_cad[$i][$j]{score}   = $marc_izq;
  76.                     $lista_cad[$i][$j]{pointer} = "left";
  77.                 }
  78.             }
  79.             else {
  80.                 if ( $marc_up >= $marc_izq ) {
  81.                     $lista_cad[$i][$j]{score}   = $marc_up;
  82.                     $lista_cad[$i][$j]{pointer} = "up";
  83.                 }
  84.                 else {
  85.                     $lista_cad[$i][$j]{score}   = $marc_izq;
  86.                     $lista_cad[$i][$j]{pointer} = "left";
  87.                 }
  88.             }
  89.         }
  90.     }
  91.  
  92.     #   Rastrear el origen
  93.  
  94.     my $alineam1 = "";
  95.     my $alineam2 = "";
  96.  
  97.     #  Empieza en la última "cell" de la matriz
  98.     my $j = length($sec1);
  99.     my $i = length($sec2);
  100.  
  101.     while (1) {
  102.         last if $lista_cad[$i][$j]{pointer} eq "none";    # Termina en la última "cell" de la matriz
  103.  
  104.         if ( $lista_cad[$i][$j]{pointer} eq "diagonal" ) {
  105.             $alineam1 .= substr( $sec1, $j - 1, 1 );
  106.             $alineam2 .= substr( $sec2, $i - 1, 1 );
  107.             $i--;
  108.             $j--;
  109.         }
  110.         elsif ( $lista_cad[$i][$j]{pointer} eq "left" ) {
  111.             $alineam1 .= substr( $sec1, $j - 1, 1 );
  112.             $alineam2 .= "-";
  113.             $j--;
  114.         }
  115.         elsif ( $lista_cad[$i][$j]{pointer} eq "up" ) {
  116.             $alineam1 .= "-";
  117.             $alineam2 .= substr( $sec2, $i - 1, 1 );
  118.             $i--;
  119.         }
  120.     }
  121.     $alineam1 = reverse $alineam1;
  122.     $alineam2 = reverse $alineam2;
  123.     return $alineam1;
  124.     return $alineam2;
  125. }
Coloreado en 0.008 segundos, usando GeSHi 1.0.8.4
Última edición por explorer el 2014-01-13 11:04 @502, editado 2 veces en total
Razón: Formateado de código con Perltidy y poner marcas Perl
Marne
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2013-11-12 06:46 @324

Publicidad

Re: Alineamiento múltiple

Notapor explorer » 2014-01-13 13:48 @617

El return de la línea 124 no se ejecuta nunca: cuando Perl llega a la línea 123, sale de la subrutina, devolviendo el valor que contiene la variable $alineam1 a la línea donde se llamó a la subrutina.

Si quieres devolver dos valores, debes escribirlo así:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.     return $alineam1, $alineam2;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Otra cosa, la línea 28 está mal: @ARGV almacena los argumentos pasados al programa, no los argumentos pasados a la subrutina. Debes cambiar @ARGV por @_
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Alineamiento múltiple

Notapor Marne » 2014-01-13 15:27 @685

Muchas gracias, el problema es que aun cambiándolo me sigue devolviendo
returned 0, con lo cual no sé qué debo cambiar para
que me leas las secuencias que le he introducido. Aún
soy un novato en este campo y muchas veces me pierdo en
los algoritmos...
Marne
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2013-11-12 06:46 @324

Re: Alineamiento múltiple

Notapor explorer » 2014-01-13 19:03 @835

En la línea 22 haces referencia a elementos del array @array, pero esa variable no ha sido inicializada.

Consejo: pon

use strict;
use warnings;
use diagnostics;


al principio del programa, para que Perl te ayude con estos casos de variables que salen de la nada.

Otra cuestión, importante... ¿qué algoritmo estás implementando, el Needleman–Wunsch o el Smith–Waterman? Porque... parece una mezcla de los dos :)
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Alineamiento múltiple

Notapor Marne » 2014-01-14 08:09 @381

explorer escribiste:use strict;
use warnings;
use diagnostics;


La verdad es que estuve mirando el algoritmo de Needleman-Wunsch. El otro no lo conozco.
Intenté que funcionara basándome en el primero, pero como siempre me sale todo mal xDD,
supongo que es normal en los comienzos de la programación, como en todo... He cambiado,
lo que me has comentado pero la verdad es que sigo igual... Estoy un poco desesperado
por que llevo sin poder avanzar días enteros por que no entiendo por qué no me devuelve
nada, si sabes de algún algoritmo ya escrito que haga lo mismo te estaría muy agradecido,
por que el de needleman-wunsch como solo está explicado a nivel matemático pero la implantación
en script no la encuentro por ningún lado pues haciéndolo yo solo pasan estas cosas...

Muchas gracias por aguantar lo coñazo que soy...
Marne
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2013-11-12 06:46 @324

Re: Alineamiento múltiple

Notapor explorer » 2014-01-14 11:49 @534

El algoritmo Smith–Waterman es una mejora del Needleman-Wunsch. En las páginas de Wikipedia tienes la definición y algún pseudocódigo, junto con la definición formal, que es muy simple, la verdad.

Información sobre estos algoritmos (y otros) y Perl dedicado a la bioinformática, lo tienes en el libro

«Genomic Perl-From Bioinformatics Basics to Working Code», de Rex A. Dwyer, editado por Cambridge University Press, diciembre 2002.

978-0-511-06339-8 (eBook) $ 82.00
978-0-521-80177-5 (libro) Se dejó de imprimir en mayo del año pasado.

Si no eres capaz de encontrarlo por tu cuenta, mándame un mensaje privado.

Algo que sí debe quedar claro es que no es necesaria la matriz 'pointer': la construcción del alineamiento se puede hacer con solo mirar las cantidades puestas en 'score'. De hecho, tampoco es necesario 'score': solo necesitamos una matriz de datos.

He estado ejecutando el código, y hay un problema a la hora de recorrer la matriz para crear los dos alineamientos. Corrigiendo todos los errores que te he comentado antes, me salen como alineamiento las mismas secuencias que le damos de entrada, así que hay un problema a la hora de meter los 'gap'.

Si tengo un rato, lo reviso. Quizás use Text::ASCIITable o Text::TabularDisplay para presentar la matriz de alineamiento.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Alineamiento múltiple

Notapor Marne » 2014-01-14 13:19 @597

Oye: muchísimas gracias por todo. ¿Sabes qué? Voy a empezar de 0 con el libro
que me has dicho, es una pasada, no pierdas más el tiempo, este programa es
un caos pues, a veces es empezar de cero con nuevas ideas y corrigiendo errores
previos, me han ayudado mucho tus consejo, me pongo de nuevo, si lo consigo te cuelgo
el programa funcional por si lo pudiera necesitar alguien en el futuro. Saludos.
Marne
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2013-11-12 06:46 @324

Re: Alineamiento múltiple

Notapor explorer » 2014-01-14 22:57 @997

Esta es mi versión, con el código del libro.
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. ##################################################################
  3. ## Programa para alinear por pares 2 secuencias introduciendo   ##
  4. ## gaps para optimizar el número de coincidencias               ##
  5. ##                                                              ##
  6. ## Algoritmo Needleman–Wunsch                                   ##
  7. ## Joaquin Ferrero, 2014-01-15                                  ##
  8. ##################################################################
  9.  
  10. use v5.16;
  11. use Text::ASCIITable;
  12.  
  13. ## Secuencias de prueba
  14. my @array = qw(
  15.     ACACACTA
  16.     AGCACACA
  17. );
  18.  
  19. ## Constantes
  20. my $GAP     = -2;                       # Valor de la no coincidencia
  21. my $similar =  1;                       # Usamos una constante en lugar de una matriz de similaridad
  22.  
  23. ## Variables
  24. my @M;                                  # Matriz bidimensional de similaridad del alineamiento
  25.                                         # Es un array de referencias a arrays
  26. my $tabla;                              # Tabla en ASCII
  27.  
  28.  
  29. ## Proceso: comparamos la primera secuencia con todas las demás
  30. my $seq1 = shift @array;
  31.  
  32. for my $seq2 (@array) {
  33.  
  34.     $tabla = Text::ASCIITable->new();
  35.  
  36.     my @alineamientos = alinea($seq1, $seq2);
  37.  
  38.     printf "%10s : %10s\n", $seq1, $alineamientos[0];
  39.     printf "%10s : %10s\n", $seq2, $alineamientos[1];
  40. }
  41.  
  42.  
  43. sub alinea {
  44.     my($seq1, $seq2) = @_;
  45.  
  46.     similaridad($seq1, $seq2);                  # calcula la matriz de alineamiento de las dos secuencias
  47.  
  48.     return obtenAlineamiento($seq1, $seq2);     # devuelve los dos alineamientos
  49. }
  50.  
  51. sub similaridad {                               # Devuelve el valor máximo de similaridad, y actualiza @M
  52.     my($s, $t) = @_;                                    # s: source (fuente). t: target (objetivo)
  53.     my $ls = length $s;
  54.     my $lt = length $t;
  55.  
  56.     ## Inicialización de la matriz
  57.     undef @M;
  58.  
  59.     for my $i (0 .. $ls) { $M[$i][0] = $GAP * $i }      # la matriz tiene una columna y fila más que el largo
  60.     for my $j (0 .. $lt) { $M[0][$j] = $GAP * $j }      # de las secuencias
  61.  
  62.     ## Relleno de la matriz
  63.     for my $i (1 .. $ls) {
  64.     for my $j (1 .. $lt) {
  65.  
  66.         ## Peso de la coincidencia entre dos letras de las dos secuencias
  67.         my $p = p(substr($s, $i-1, 1), substr($t, $j-1, 1));
  68.  
  69.         ## Actualizamos matriz
  70.         $M[$i][$j] = max(
  71.             $M[$i-1][$j]   + $GAP,      # Caso de insertar letra   (Insert)
  72.             $M[$i][$j-1]   + $GAP,      # Caso de faltar una letra (Delete)
  73.             $M[$i-1][$j-1] + $p         # Caso de coincidir        (Match)
  74.         );
  75.     }}
  76.  
  77.     ## Relleno de la tabla ASCII
  78.     my @s = ('', split //, $s);
  79.     my @t = ('', split //, $t);
  80.     $tabla->setCols('', @t);
  81.  
  82.     for my $i (0 .. $ls) {
  83.         $tabla->addRow($s[$i], @{$M[$i]});
  84.     }
  85.  
  86.     print $tabla;
  87.  
  88.     return $M[$ls][$lt];
  89. }
  90.  
  91. sub obtenAlineamiento {                                                 ## Versión recursiva
  92.     my($seq1, $seq2) = @_;
  93.     my($i, $j) = (length($seq1), length($seq2));
  94.  
  95.     return ("-" x $j, $seq2) if $i == 0;
  96.     return ($seq1, "-" x $i) if $j == 0;
  97.  
  98.     my($ultimo1, $ultimo2) = (substr($seq1, -1), substr($seq2, -1));
  99.  
  100.     if ($M[$i][$j] == $M[$i-1][$j-1] + p($ultimo1,$ultimo2)) {          ## Case 1
  101.         ## last letters are paired in the best alignment
  102.         my ($sa,$ta) = obtenAlineamiento(substr($seq1,0,-1),substr($seq2,0,-1));
  103.  
  104.         return ($sa . $ultimo1, $ta . $ultimo2);
  105.     }
  106.     elsif ($M[$i][$j]==$M[$i-1][$j]+$GAP) {                             ## Case 2
  107.         ## last letter of the first string is paired with a gap
  108.         my ($sa,$ta) = obtenAlineamiento(substr($seq1,0,-1),$seq2);
  109.  
  110.         return ($sa . $ultimo1, $ta . "-");
  111.     }
  112.     else {
  113.         ## Case 3: last letter of the second string is paired with a gap
  114.         my ($sa,$ta) = obtenAlineamiento($seq1,substr($seq2,0,-1));
  115.  
  116.         return ($sa . "-", $ta . $ultimo2);
  117.     }
  118. }
  119.  
  120. sub p {                                 # Devuelve un valor según la igualdad de las letras
  121.     my($a, $b) = @_;
  122.     return ($a eq $b) ? $similar : -$similar;
  123. }
  124.  
  125. sub max {                               # Devuelve el máximo de una serie de argumentos
  126.     my ($m, @l) = @_ ;
  127.     for my $x (@l) { $m = $x  if  $x > $m }
  128.     return $m;
  129. }
  130.  
  131. __END__
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4
Esta es la salida:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
.-----------------------------------------------------------.
|     |     | A   | G   | C   | A   | C   | A   | C   | A   |
+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
|     |   0 |  -2 |  -4 |  -6 |  -8 | -10 | -12 | -14 | -16 |
| A   |  -2 |   1 |  -1 |  -3 |  -5 |  -7 |  -9 | -11 | -13 |
| C   |  -4 |  -1 |   0 |   0 |  -2 |  -4 |  -6 |  -8 | -10 |
| A   |  -6 |  -3 |  -2 |  -1 |   1 |  -1 |  -3 |  -5 |  -7 |
| C   |  -8 |  -5 |  -4 |  -1 |  -1 |   2 |   0 |  -2 |  -4 |
| A   | -10 |  -7 |  -6 |  -3 |   0 |   0 |   3 |   1 |  -1 |
| C   | -12 |  -9 |  -8 |  -5 |  -2 |   1 |   1 |   4 |   2 |
| T   | -14 | -11 | -10 |  -7 |  -4 |  -1 |   0 |   2 |   3 |
| A   | -16 | -13 | -12 |  -9 |  -6 |  -3 |   0 |   0 |   3 |
'-----+-----+-----+-----+-----+-----+-----+-----+-----+-----'
  ACACACTA :  A-CACACTA
  AGCACACA :  AGCACAC-A
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: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Alineamiento múltiple

Notapor Marne » 2014-01-17 06:10 @298

Muchas gracias de nuevo, he conseguido implementar tu código y el que he sacado
del libro que me recomendastes y hacer una comparativa de tiempos, lo cual me
consigue siempre un valor más fiable al hacerlo de varias formas. Estoy ya con
la fase final del programa. De verdad que me han servido mucho tus consejos. :)
Un saludo,
Marne
Marne
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2013-11-12 06:46 @324

Re: Alineamiento múltiple

Notapor Marne » 2014-01-20 11:33 @523

Os dejo el programa funcional hasta la línea 380,
solo me queda un problemilla más, a partir de dicha línea.
Me gustaría que repitiera lo que es el programa en sí pero
haciendo un reverse de cada string para que me dé otro
resultado de alineamiento pero con las cadenas al revés.
Lo he intentado pero no sé cómo hacerlo... Si pudierais
echarme un último cable, yo feliz, y si no, pues haré los
estudios solo con las cadenas del derecho. Muchas gracias
por toda la ayuda prestada hasta conseguir que me funcionara esto.
Saludos,
Marne

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. ##############################################
  2. ## PROGRAMA DE OPTIMIZACIÓN DE ALINEAMIENTO ##
  3. ##############################################
  4.  
  5. use strict;
  6. use Benchmark;
  7.  
  8. $| = 1;                                # De este modo no hace todos los prints al final del programa y los va sacando uno a uno
  9.  
  10. print "1) Introduzca el valor de penalización por GAP: \n";
  11.  
  12. my $castigoGAP = <STDIN>;
  13. my @secuencias;
  14. my $tiempo_inicial;
  15. my $tiempo_final;
  16. my $tiempo;
  17. my @ss;
  18. my $ss;
  19. my $archivo;
  20. my @filedata;
  21. my $archivo_txt;
  22. my @matriz;
  23. my @ssRev;
  24. my $ssRev;
  25. my $linea;
  26.  
  27. ########################################################################################################
  28.  
  29. ##########################################################
  30. ## SUBRUTINAS UTILIZADAS EN EL PROGRAMA DE ALINEAMIENTO ##
  31. ##########################################################
  32.  
  33. sub get_file_data {
  34.     my ($filename) = @_;
  35.  
  36.     unless ( open( GET_FILE_DATA, $filename ) ) {
  37.         print STDERR "No se puede abrir el fichero \"$filename\"\n\n";
  38.         exit;
  39.     }
  40.     @filedata = <GET_FILE_DATA>;
  41.  
  42.     return @filedata;
  43. }
  44.  
  45. sub extract_seq {
  46.     my @archivo = @_;
  47.     foreach my $line (@archivo) {
  48.         if ( $line =~ /^\s*$/ ) {
  49.             next;
  50.         }
  51.         elsif ( $line =~ /^\s*[0-9]/ ) {
  52.             next;
  53.         }
  54.         else {
  55.             $line =~ s/\n//g;          #Elimina el salto de linea
  56.             push( @ss, $line );
  57.         }
  58.  
  59.     }
  60. }
  61.  
  62.        
  63. sub Alineamiento_Fusion {
  64.  
  65.     #####################################################################
  66.     ##  Subrutina que rellena 2 tablas con programación dinámica @M    ##
  67.     ##  y @como, las cuales nos devuelven el mejor resultado posible   ##
  68.     ##  para la fusión de alineamientos de $S y $T, y como añadir a    ##
  69.     ##  estas dos $ los GAPs de la manera más cercana al óptimo.       ##
  70.     ##  Devuelve el mejor marcador, haciendo referencia a @como.       ##
  71.     #####################################################################
  72.  
  73.     my ( $S, $T ) = @_;                # Alineamiento que hace referencia a una lista de secuencias
  74.                                        # con los GAPs
  75.  
  76.     my ( $sLong, $tLong ) = ( length( $$S[0] ), length( $$T[0] ) );
  77.     my @Sgaps = ( ("-") x @$S );
  78.     my @Tgaps = ( ("-") x @$T );
  79.     my @M;
  80.     my @como;
  81.     $M[0][0] = 0;
  82.     foreach my $i ( 1 .. $sLong ) {
  83.         $como[$i][0] = "|";
  84.         my @Scol = map { substr( $_, $i - 1, 1 ) } @$S;
  85.         $M[$i][0] = $M[ $i - 1 ][0] + marcadorCol( @Scol, @Tgaps );
  86.     }
  87.     foreach my $j ( 1 .. $tLong ) {
  88.         $como[0][$j] = "-";
  89.         my @Tcol = map { substr( $_, $j - 1, 1 ) } @$T;
  90.         $M[0][$j] = $M[0][ $j - 1 ] + marcadorCol( @Sgaps, @Tcol );
  91.     }
  92.     foreach my $i ( 1 .. $sLong ) {
  93.         foreach my $j ( 1 .. $tLong ) {
  94.             my @Scol = map { substr( $_, $i - 1, 1 ) } @$S;
  95.             my @Tcol = map { substr( $_, $j - 1, 1 ) } @$T;
  96.             $M[$i][$j] = $M[ $i - 1 ][ $j - 1 ] + marcadorCol( @Scol, @Tcol );
  97.             $como[$i][$j] = "\\";
  98.             my $left = $M[$i][ $j - 1 ] + marcadorCol( @Sgaps, @Tcol );
  99.             if ( $left > $M[$i][$j] ) {
  100.                 $M[$i][$j]    = $left;
  101.                 $como[$i][$j] = "-";
  102.             }
  103.             my $up = $M[ $i - 1 ][$j] + marcadorCol( @Scol, @Tgaps );
  104.             if ( $up > $M[$i][$j] ) {
  105.                 $M[$i][$j]    = $up;
  106.                 $como[$i][$j] = "|";
  107.             }
  108.         }
  109.     }
  110.     return ( $M[$sLong][$tLong], \@como );
  111. }
  112.  
  113. ########################################################################################################
  114.  
  115. sub AñadeColumPrevia {
  116.  
  117.         ###################################################################################
  118.         ##  Añade una nueva columna en el extremo izquierdo de un alineamiento existente ##   
  119.         ###################################################################################
  120.        
  121.     my ($A, @col)= @_;          ## Este alineamiento hace referencia a una lista de cadenas
  122.                                 ## Lista de bases a añadir a cada cadena de alineamiento
  123.        
  124.     foreach (@$A) {$_ = (shift @col).$_};
  125. }
  126.  
  127. ########################################################################################################
  128.  
  129. sub ComparaCadenas {
  130.  
  131.     ################################################################################
  132.     ## Dada una serie de cadenas, las ordena en función de su longitud y compara  ##
  133.     ## las bases T, C, G y A de cada una de las cadenas devolviendo el número     ##
  134.     ## de bases que coinciden en todas las cadenas introducidas.                  ##
  135.     ################################################################################
  136.  
  137.     my @matriz   = @_;
  138.     my $longitud = length( $matriz[1] );
  139.     my @bases    = ( 'T', 'C', 'G', 'A' );
  140.     my $coincid  = 0;
  141.  
  142.     my @similitud;
  143.     for ( my $i = 0; $i < $longitud; $i++ ) {
  144.         for ( my $z = 0; $z < int(@matriz); $z++ ) {
  145.             $similitud[$z] = substr( $matriz[$z], $i, 1 );
  146.         }
  147.         for ( my $j = 0; $j < 4; $j++ ) {
  148.             my $c      = $bases[$j];
  149.             my $contar = 0;
  150.             foreach (@similitud) {
  151.                 if (/$c/) {
  152.                     ++$contar;
  153.                 }
  154.             }
  155.             if ( $contar == int(@similitud) ) {
  156.                 $coincid++;
  157.             }
  158.         }
  159.     }
  160.     return $coincid;
  161. }
  162.  
  163.  
  164. #######################################################################################################
  165.  
  166. sub EstrategiaOptima {
  167.  
  168.     #################################################################################
  169.     ##  Dada una lista de secuencias, aplica una estrategia de fusión que busca    ##
  170.     ##  maximizar el alineamiento, pese a no tener por que ser el óptimo.          ##
  171.     ##  Dicha estrategia consiste en la maximización de manera repetitiva de la    ##
  172.     ##  ganancia inmediata por la fusión de cualquier secuencia individual con     ##
  173.     ##  las alineaciones existentes.                                               ##
  174.     ##  Devuelve el resultado final del alineamiento múltiple.                     ##
  175.     #################################################################################
  176.  
  177.     my @ss = @_;                       ## Argument list is a list of sequences.
  178.  
  179.     my ( $bestI, $bestJ, $bestScore ) = ( -1, -1, -999999999 );    ## Busca el mejor molde para la fusión.
  180.     foreach my $i ( 0 .. $#ss ) {
  181.         foreach my $j ( $i + 1 .. $#ss ) {
  182.             my ( $score, $como ) = Alineamiento_Fusion( [ $ss[$i] ], [ $ss[$j] ] );
  183.             ( $bestI, $bestJ, $bestScore ) = ( $i, $j, $score )
  184.                 if $score > $bestScore;
  185.         }
  186.     }
  187.  
  188.     my $alineamiento = FusionAlineamiento( [ $ss[$bestI] ], [ $ss[$bestJ] ] );    ## Fusión.
  189.     splice( @ss, $bestJ, 1 );
  190.     splice( @ss, $bestI, 1 );
  191.  
  192.     while (@ss) {
  193.  
  194.         my ( $bestI, $bestScore ) = ( -1, -999999999 );                           ## Busca la mejor secuencia para
  195.         foreach my $i ( 0 .. $#ss ) {  ## añadirla al alineamiento.
  196.             my ( $score, $como ) = Alineamiento_Fusion( [ $ss[$i] ], $alineamiento );
  197.             ( $bestI, $bestScore ) = ( $i, $score ) if $score > $bestScore;
  198.         }
  199.  
  200.         $alineamiento = FusionAlineamiento( [ $ss[$bestI] ], $alineamiento );     ## Añade la mejor secuencia al
  201.         splice( @ss, $bestI, 1 );      ## alinemiento, la elimina de @ss
  202.     }
  203.     return $alineamiento;
  204. }
  205.  
  206.  
  207. ########################################################################################################
  208.  
  209. sub FusionAlineamiento {
  210.  
  211.     ########################################################################################
  212.     ##  Fusiona dos alineamientos para hallar el máximo resultado de fusión posible        ##
  213.     ##  Devuelve el resultado del alineaminto, hace referencia a una lista de cadenas.    ##
  214.     ########################################################################################
  215.  
  216.     my ( $S, $T ) = @_;
  217.     my ( $score, $como ) = Alineamiento_Fusion( $S, $T );
  218.     my @result = ( ("") x ( @$S + @$T ) );
  219.     my ( $i, $j ) = ( length( $$S[0] ), length( $$T[0] ) );
  220.     my @Sgaps = ( ("-") x @$S );
  221.     my @Tgaps = ( ("-") x @$T );
  222.     while ( $i > 0 || $j > 0 ) {
  223.         if ( $$como[$i][$j] eq "\\" ) {
  224.             my @Scol = map { substr( $_, $i - 1, 1 ) } @$S;
  225.             my @Tcol = map { substr( $_, $j - 1, 1 ) } @$T;
  226.             AñadeColumPrevia(\@result, @Scol, @Tcol);
  227.             $i--; $j--;
  228.         } elsif ($$como[$i][$j] eq "|") {
  229.             my @Scol = map { substr($_,$i-1,1) } @$S;
  230.             AñadeColumPrevia(\@result, @Scol, @Tgaps);
  231.             $i--;
  232.         } elsif ($$como[$i][$j] eq "-") {
  233.             my @Tcol = map { substr($_,$j-1,1) } @$T;
  234.             AñadeColumPrevia(\@result, @Sgaps, @Tcol);
  235.             $j--;
  236.         }
  237.     }
  238.     return \@result;
  239. }
  240.  
  241.  
  242. ########################################################################################################
  243.  
  244. sub MarcadorAlinMulti {
  245.  
  246.     #####################################################################
  247.     ## Resume los resultados de la suma de pares de todas las columnas ##
  248.     ## de un alineamiento, devolviendo el marcador de la suma de los   ##
  249.     ## pares del alineamiento completo.                                ##
  250.     #####################################################################
  251.  
  252.     my ($alineamiento) = @_;
  253.     my $score;
  254.     foreach my $i ( 0 .. length( $$alineamiento[0] ) - 1 ) {
  255.         $score += marcadorCol( map { substr( $_, $i, 1 ) } @$alineamiento );
  256.     }
  257.     return $score;
  258. }
  259.  
  260. #######################################################################################################
  261.  
  262. sub marcadorCol {
  263.  
  264.     ############################################################
  265.     ##  Dada una lista de bases en una columna, devuelve      ##
  266.     ##  el marcador de la suma de los pares para una columna  ##
  267.     ############################################################
  268.  
  269.     my @col = @_;                      # Las bases de la columna
  270.     my ( $gaps, $aas, $score ) = ( 0, 0, 0 );
  271.     while (@col) {
  272.         my $aa = shift @col;
  273.         ( $gaps++, next ) if $aa eq "-";
  274.         $aas++;
  275.         foreach my $aa1 (@col) {
  276.             next if $aa1 eq "-";
  277.             $score += ( $aa eq $aa1 ) ? +1 : -1;
  278.         }
  279.     }
  280.     return $score + ( $castigoGAP * $gaps * $aas );
  281. }
  282.  
  283.  
  284. #######################################################################################################
  285.  
  286. sub ImprimeAlin {
  287.  
  288.     ########################################
  289.     ##  Imprime el alineamiento múltiple  ##
  290.     ########################################
  291.  
  292.     my ($alineamiento)                 ## Hace referencia a una lista de cadenas
  293.         = @_;
  294.     foreach (@$alineamiento) { push( @secuencias, $_ ); }
  295.     return @secuencias;
  296. }
  297.  
  298. ########################################################################################################
  299.  
  300. sub Exportar_archivo {
  301.  
  302.     my ( $archivo_txt, @exportar ) = @_;
  303.     open( Archivo, ">" . $archivo_txt . '.txt' );
  304.     my $N = int(@exportar);
  305.     for ( my $i = 0; $i < $N; $i++ ) {
  306.         print Archivo $exportar[$i] . "\n";
  307.     }
  308.     close Archivo;
  309. }
  310.  
  311. ###################################################################################################
  312.  
  313. sub Reverse_ss {
  314.     @ssRev = reverse @ss;
  315.     foreach $linea (@ss) {
  316.         push( @ssRev, reverse $linea );
  317.     }
  318.     return "@ssRev";
  319. }
  320.  
  321.  
  322. #######################################################################################################
  323. #######################################################################################################
  324.  
  325.  
  326.         ###########################
  327.         ##  PROGRAMA: EJECUCIÓN  ##
  328.         ###########################
  329.  
  330.  
  331. #######################################################################################################
  332. #######################################################################################################
  333.  
  334.  
  335. print "2) Introduzca el nombre del archivo a analizar\n";
  336. $archivo = <STDIN>;
  337. print("\n");
  338. get_file_data($archivo);
  339. extract_seq(@filedata);
  340.  
  341. {
  342.     print "El resultado de comparar las secuencias previo al alineamiento devuelve:\n\n";
  343.     print join( "\n", @ss ), "\n\n     Cuyo valor de coincidencia en todas las cadenas es:      ", ComparaCadenas(@ss),
  344.         "\n\n\n";
  345. }
  346.  
  347. $tiempo_inicial = new Benchmark;
  348. my $count = 0;
  349. @matriz = ();
  350. print "4.1) Introduzca el tamaño de las particiones a realizar previas al alineamiento; \n";
  351. my $val = <STDIN>;
  352. chomp($val);
  353.  
  354. for ( my $i = 0; $i < length( $ss[1] ); $i += $val ) {    ### Si quieres cambios modifica el valor
  355.     my @secu;
  356.     for ( my $j = 0; $j < int(@ss); $j++ ) {
  357.         $secu[$j] = substr( $ss[$j], $i, $val );          ###  Si quieres cambios modifica el valor
  358.     }
  359.     ImprimeAlin( EstrategiaOptima(@secu) );
  360.  
  361.     $count += ComparaCadenas(@secuencias);
  362.     for ( my $jj = 0; $jj < int(@secuencias); $jj++ ) {
  363.         $matriz[$jj] = $matriz[$jj] . $secuencias[$jj];
  364.     }
  365.     @secuencias = ();                  #Borra todo el array.
  366. }
  367.  
  368. print "\n5.1) El resultado de alinear las secuencias de $val en $val y juntarlas posteriormente es:\n\n";
  369. print join( "\n", @matriz ), "\n\n El valor máximo de alineamiento es: ", $count, "\n\n";
  370.  
  371. $tiempo_final = new Benchmark;
  372. $tiempo = timediff( $tiempo_final, $tiempo_inicial );
  373. print "El tiempo invertido es:", timestr($tiempo), "\n \n \n";
  374.  
  375. print "6.1) Introduzca el nombre del archivo de salida, por defecto se le añadirá la extensión .txt \n";
  376. $archivo_txt = <STDIN>;
  377. chomp($archivo_txt);
  378. Exportar_archivo( $archivo_txt, @matriz );
  379.  
  380.  
  381.  
  382.  
  383. ##############################################################################################################
  384. ##############################################################################################################
  385. #      A partir de aquí me lío y no consigo que funcione
  386.  
  387.  
  388. Reverse_ss(@ss);
  389. {
  390.     print "El resultado de comparar las secuencias previo al alineamiento devuelve:\n\n";
  391.     print join( "\n", @ssRev ), "\n\n     Cuyo valor de coincidencia en todas las cadenas es:   ",
  392.         ComparaCadenas(@ssRev), "\n\n\n";
  393. }
  394.  
  395. my $tiempo_inicial2 = new Benchmark;
  396. my $count2          = 0;
  397. my @matriz2         = ();
  398. print "4.2) Introduzca el tamaño de las particiones a realizar previas al alineamiento; \n";
  399. $val = <STDIN>;
  400. chomp($val);
  401.  
  402. for ( my $i = 0; $i < length( $ssRev[1] ); $i += $val ) {    ### Si quieres cambios modifica el valor
  403.     my @secu;
  404.     for ( my $j = 0; $j < int($ssRev); $j++ ) {
  405.         $secu[$j] = substr( $ssRev[$j], $i, $val );          ###  Si quieres cambios modifica el valor
  406.     }
  407.     ImprimeAlin( EstrategiaOptima(@secu) );
  408.  
  409.     $count += ComparaCadenas(@secuencias);
  410.     for ( my $jj = 0; $jj < int(@secuencias); $jj++ ) {
  411.         $matriz2[$jj] = $matriz[$jj] . $secuencias[$jj];
  412.     }
  413.     @secuencias = ();                  #Borra todo el array.
  414. }
  415.  
  416. print "\n5.2) El resultado de alinear las secuencias de $val en $val y juntarlas posteriormente es:\n\n";
  417. print join( "\n", @matriz ), "\n\n El valor máximo de alineamiento es: ", $count2, "\n\n";
  418.  
  419. $tiempo_final = new Benchmark;
  420. $tiempo = timediff( $tiempo_final, $tiempo_inicial );
  421. print "El tiempo invertido es:", timestr($tiempo), "\n \n \n";
  422.  
  423. print "6.2) Introduzca el nombre del archivo de salida, por defecto se le añadirá la extensión .txt \n";
  424. $archivo_txt = <STDIN>;
  425. chomp($archivo_txt);
  426. Exportar_archivo( $archivo_txt, @matriz );
Coloreado en 0.009 segundos, usando GeSHi 1.0.8.4
Última edición por explorer el 2014-01-20 12:28 @561, editado 1 vez en total
Razón: Formateado de código con Perltidy y poner marcas Perl
Marne
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2013-11-12 06:46 @324

Siguiente

Volver a Bioinformática

¿Quién está conectado?

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