Cerca de subseqüències úniques

PROGRAMAS

COMPARACIÓN DE LA COMPLEJIDAD CON EL ALGORITMO SIMPLE DE CORRESPONDECIA EXACTA

Con el objetivo de comparar la eficiencia del algoritmo simple de correspondencia exacta con el algoritmo de Knuth-Morris-Pratt, se hizo un programa en Perl que implentaba el algoritmo simple de correspondencia exacta. La lógica fundamental de este algoritmo está detallada en el apartado de Metodología.

La función primordial del programa es calcular el número de operaciones que realiza el algoritmo de correspondencia exacta, y así poderlo comparar con el algoritmo de Knuth-Morris-Pratt. Por tanto, en el programa se le pasa como argumento el mRNA, a partir del cual genera todos lo pratrones posibles y busca los que son únicos en este mRNA. Además, contiene un contador a partir de la variable $x, en que se guarda el número de operaciones que realiza el algoritmo. Para poder ver el código de este programa, descargar correspondencia.pl.

Para llevar a cabo esta comparación de la complejidad algorítmica, desarrollamos un programa similar, es decir, que también genera patrones a partir de un mRNA y busca si son únicos o no en éste, pero para el algoritmo de Knuth-Morris-Pratt. Los resultados de esta comparación se pueden ver en el apartado de Resultados.

BÚSQUEDA DE SUBSECUENCIAS ÚNICAS DE UN MRNA CONCRETO EN UN GENOMA DADO

DESCARGA DEL CÓDIGO PROGRAMA cercasubsequencies.pl

PARÁMETROS NECESARIOS PARA EJECUTARLO
[Argumento0]	Nombre del fichero que contiene el mRNA en formato FASTA del gen en cuestión
[Argumento1]	Nombre del directorio donde se encuentra el genoma
[Argumento2]	Longitud específica con la se quiere empezar a buscar patrones, si no el programa por defecto iniciará a longitud 14 (parámetro opcional)

La función principal del programa consiste en buscar subsecuencias únicas y específicas de un mRNA concreto en todo un genoma (humano, ratón,...).

La forma como se ha desarrollado el programa para conseguir este objetivo ha sido creando subestructuras óptimas del programa para cada problema a resolver:

Construcción de la next table

La next table es un factor sine qua non del algoritmo de Knuth-Morris-Pratt y, por tanto, uno de los primeros pasos a realizar. La base sobre la que se fundamenta esta tabla de coincidencias dentro del algoritmo KMP está explicada en el apartado de Metodología.

El código de programación, para el cual se calcula esta next table en el programa, se encapsula dentro de una función llamada nexttable:

sub nexttable{

    my $patro = $_[0]; 

    @patro = split(//,$patro);
    $m = scalar(@patro);

    my $i = 2;
    my $c = 0;

    $coin[0] = -1;
    $coin[1] = 0;

    while ($i < $m){

	$j = 1;
    
	while ($j < $i){
	
	    $w = 0;
	    while ($w < $j && $j < $i){
		if ($patro[$j] eq $patro[$w]){
		    $c = $c + 1;
		    $j = $j + 1;
		    $w = $w + 1;
		}else{    
		    if ($c > 0){
			$w = 0;
			$c = 0;
		    }else{
			$j = $j + 1;
			$c = 0;
		    }
		}
	    }
	}
	$coin[$i] = $c;
	$i = $i + 1;
	$c = 0;
    }
    return @coin;
}

Esta función recibe como primero y único parámetro el patrón para el cual se construirá la next table. Ésta se basa en la construcción de una tabla de coincidencias del patrón en un vector llamado @coin (de coindencias) indexado para las mismas posiciones del patrón. Este vector siempre tendrá el valor de -1 en la primera posición (0) y el valor de 0 en la segunda posición (1). A partir de aquí, se recorre el resto de posiciones del patrón comparando los prefijos máximos anteriores a esa posición con el patrón y anotando las coincidencias en cada posición. Si se quiere consultar una explicación detallada de su funcionamiento, ver el código del programa.

Una vez obtenida la next table, aunque es de vital importancia, todavía faltarán resolver muchas otras cuestiones a tener en cuenta para cumplir los objetivos del programa.

Algoritmo de Knuth-Morris-Pratt

La lógica fundamental de este algoritmo y el por qué de su ganancia en eficiencia respecto al algoritmo simple de correspondencia exacta está detallada en el apartado Metodología.

Su implementación en el seno del programa es en dos funciones: una función llamada knuth, para comprobar que los patrones sean únicos dentro del propio mRNA; y otro función, knuth4genome, para hacer la misma comprobación en todo el genoma.

La primera función, knuth, recibe tres parámetros:la next table, ya calculada, la secuencia mRNA y el patrón. Y devuelve 1 (cierto) ó 0 (falso) según si el patrón a comprobar es único o no sobre la secuencia de mRNA.

La segunda función, knuth4genome, recibe cuatro parámetros: la next table, la dirección del directorio donde se localiza el genoma, el patrón a comprobar y la longitud del patrón. Y devuelve 1 (cierto)* ó 0 (falso) según si el patrón a comprobar es único o no sobre el genoma, y en el caso de ser único, también devuelve el desplazamiento y el cromosoma donde se localiza.

* La función devolverá 1 cuando el contador de patrones únicos sea igual o más pequeño que 1. Ésto se hace de esta manera para aquellos patrones únicos en el mRNA pero que no se encuentren en el genoma, debido a que se puede encontrar separado en dos exones diferentes con intrones entremedio.

Asímismo, en esta segunda función, antes de implementar el algoritmo propio KMP, se han de tener presentes ciertas preparaciones:

Para poden comprobar si el patrón es único en todo el genoma, es necesario que este último se encuentre en un directorio. Lo que se hace es abrir el directorio a partir de la función opendir del Perl, y almacenar todos los ficheros de cromosomas que contenga en un vector, excepto los ficheros ocultos "." y "..":

my @fitxerscromosomes = ();

opendir(GENOMA,$genoma);
@fitxerscromosomes = grep { $_ ne '.' and $_ ne '..' } readdir(GENOMA);
closedir(GENOMA);

Después, mediante una iteración con un while, se comprueba que el patrón sea único cromosoma por cromosoma. En cada caso, se abre el cromosoma que toca, y se realiza la comprobación línea por línea. También hay que recordar concatenar al inicio de cada línea los últimos nucleótidos de la línea anterior (según la longitud del patrón).

while ($z < scalar(@fitxerscromosomes && $c<=1){

    open(CROMOSOMA,"< $genoma/$fitxerscromosomes[$z]"); 

    my $trosset = "";
    my $defline = <CROMOSOMA>;
    while (<CROMOSOMA>) {
        chomp;
        my $seq = $trosset . $_;   
        $seq = uc($seq);
        $trosset = substr($_,length($_)-$llargada);
        *******
    }
    close(CROMOSOMA);
    $z = $z + 1;
}

Finalmente, algunas subsecuencias únicas pueden localizarse en el strand negativo de la secuencia de los cromosomas. Para solucionar este hecho, se hace el "reverse complement" de la secuencia del patrón antes de empezar a buscar:

my $patrorc = reverse($patro);     
$patrorc =~ tr/acgtACGT/tgcaTGCA/;

Y cuando se está buscando línea por línea, si en la línea actual no se encuentra, entonces se repite la misma búsqueda pero con el patrón en reverse complement, y si tampoco se encuentra se sigue con el procedimiento habitual del algoritmo.

De todas formas, tanto en una función como en la otra, el algoritmo propio del KMP es el siguiente:


my @mrna = split(//,$mrna);
my $n = scalar(@mrna);

my @patro = split(//,$patro);
my $m = scalar(@patro);

my $s = 0;
my $c = 0;
my $i = 1; 
    
while ($s <= $n-$m){

    $i = $coin[$i];

    $x = $x + 1;
	
    while ($i >= 0 && $i < $m && $mrna[$s+$i] eq $patro[$i]){
	$i = $i + 1;
	$x = $x + 1;
    }
	
    if ($i == $m){
	$c = $c + 1;
	$i = $i - 1;
	$s = $s + ($i - $coin[$i]);
	$i = 1;
    }else{
	if ($s == 0){
            if ($i == 0){
		$s = 1;
		$i = 1;
	    }else{
		$s = $i - $coin[$i];
	    }
	}else{
	    if ($i == 0){
		$s = $s + 1;
		$i = 1;
	    }else{
		$s = $s + ($i - $coin[$i]);
	    }
	}
    }
}

El fundamento de éste es recorrer la secuencia en que queremos comprobar si el patrón es único o no, comparándola con el patrón. Pero a la hora de modificar el desplazamiento según si hemos encontrado el patrón exacto o no, se hace mediante la información proporcionada por la next table. Así directamente se consigue desplazar suficientemente el patrón sobre la secuencia, hasta el punto de ahorrarse posiciones de la secuencia que no hace falta comparar y también, iniciar la siguiente comparación en la posición concreta del patrón (saltándose las posiciones del patrón que tampoco haga falta comparar). Para tener una explicación detallada de su funcionamiento, ver el código del programa.

Generación de patrones

El programa es capaz, por sí mismo, de abrir el fichero FASTA del mRNA y generar todos los posibles patrones que se puedan derivar, por medio de la función substr del Perl, para una longitud determinada. De esta manera, se asegura que el programa analizará todas las subsecuencias susceptibles de ser únicas de un gen determinado.

$inici = 0;
$llargada = 14;
while ($inici <= $n - $llargada){
    $patro = substr($mrna, $inici, $llargada);
    *****
    *****
    $inici = $inici + 1;  
}

Estrategia automática de modificar la longitud

Una vez comprobado si exiten patrones únicos para una longitud concreta en el mRNA, hay dos posibles opciones: que se hayan encontrado o que no.

En el caso que no se hayan encontrado patrones únicos, sumaremos una unidad a la longitud y se volverá a hacerse la búsqueda de subsecuencias únicas a partir del mRNA.

Si se han encontrado, almacenaremos los patrones, su desplazamiento y el cromosoma donde se encuentran, con los tres vectores resultado. Y procederemos a volver a hacer la búsqueda de subsecuencias únicas pero para una longitud de una unidad inferior. Además, en este caso no se efectuará la búsqueda a partir de todo el mRNA entero, sino a partir de los patrones únicos encontrados con la longitud de la cual procedemos, ya que cualquier patrón único de una longitud inferior (por ejemplo de 13) estará contenido dentro de algunos de los patrones únicos encontrados con longitud 14.

Consideraciones:

En el caso de haber ido a una longitud inferior y encontrar patrones, continuamos con longitudes inferiores. Pero en el momento que no encontremos patrones únicos, nos quedamos con los resultados anteriores, previamente almacenados en los tres vectores resultados. Y saldremos de la iteración, mediante la evaluación de la variable $tinc_patro_unic como cierta (1), ya que en este caso se entrará en el "elsif" detallado más abajo.
La otra opción es haber ido a una longitud superior. Si no encontramos patrones, continuaremos con longitudes superiores. En el momento que encontremos patrones únicos, los almacenaremos en los tres vectores resultado, y como la longitud anterior será siempre inferior a la actual (porque se habrá incrementado la longitud), saldremos de la iteración evaluando también como cierta la variable $tinc_patro_unic (1), en el segundo "if" detallado a continuación.

if (scalar(@unics) > 0) {
    @resultat = @unics;
    @desplacaments_resultat = @desplacament;
    @cromosoma_resultat = @cromosoma;

    if ($llargada_anterior < $llargada) {
        $tinc_patro_unic = 1;

    } else {
        $llargada = $llargada - 1;
	@possiblespatrons = @unics;
    }

} elsif (scalar(@resultat) > 0) {
    $tinc_patro_unic = 1;

} else {
    $llargada =$llargada = $llargada + 1;
}