Cerca de subseqüències úniques

PROGRAMES

COMPARACIÓ DE LA COMPLEXITAT AMB L'ALGORISME SIMPLE DE CORRESPONDÈNCIA EXACTA

Amb l'objectiu de comparar l'eficiència de l'algorisme simple de correspondència exacta amb l'algorisme de Knuth-Morris-Pratt, es va fer un programa en Perl que implementava l'algorisme simple de correspondència exacta. La lògica fonamental d'aquest algorisme està detallada a l'apartat de Metodologia.

La funció primordial del programa és calcular-nos el nombre d'operacions que realitza l'algorisme de correspondència exacta, i així poder-ho comparar amb l'algorisme de Knuth-Morris-Pratt. Per tant, en el programa se li passa com a argument el mRNA, del qual en genera tots els patrons possibles i cerca els que són únics en aquest mRNA. A més, conté un comptador a partir de la variable $x, en què, s'hi guarda el nombre d'operacions que realitza l'algorisme. Per poder veure el codi d'aquest programa, descarregar correspondencia.pl.

Per dur a terme aquesta comparació de la complexitat algorísmica, vam fer un programa similar, és a dir, que també genera patrons a partir d'un mRNA i busca si són únics o no en aquest, però per l'algorisme de Knuth-Morris-Pratt. Els resultats d'aquesta comparació es poden veure a l'apartat de Resultats.

CERCA DE SUBSEQÜÈNCIES ÚNIQUES D'UN MRNA CONCRET EN UN GENOMA DONAT

DESCÀRREGA DEL CODI PROGRAMA cercasubsequencies.pl

PARÀMETRES NECESSARIS PER A EXECUTAR-LO
[Argument0]	Nom del fitxer que conté el mRNA en format FASTA del gen en qüestió
[Argument1]	Nom del directori on es troba el genoma
[Argument2]	Llargada específica amb la que es vol que comenci a buscar patrons, si no el programa per defecte iniciarà a llargada 14 (paràmetre opcional)

La funció principal del programa consisteix en cercar subseqüències úniques i específiques d'un mRNA concret en tot un genoma (humà, ratolí,...).

La forma com s'ha desenvolupat el programa per aconseguir aquest objectiu ha estat creant subestructures òptimes del programa per a cada problema a resoldre:

Construcció de la next table

La next table és un factor sine qua non de l'algorisme de Knuth-Morris-Pratt i, per tant, un dels primers passos a realitzar. La base sobre la que es fonamenta aquesta taula de coincidències dins l'algorisme KMP està explicada a l'apartat de Metodologia.

El codi de programació, pel qual es calcula aquesta next table en el programa, s'encapsula dins una funció anomenada 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;
}

Aquesta funció rep com a primer i únic paràmetre el patró pel que es construirà la next table. Aquesta es basa en la construcció d'una taula de coincidències del patró en un vector anomenat @coin (de coincidències) indexat per les mateixes posicions del patró. Aquest vector sempre tindrà el valor de -1 en la primera posició (0) i el valor de 0 en la segona posició (1). A partir d'aquí, es recorre la resta de posicions del patró tot comparant els prefixs màxims anteriors a aquella posició amb el patró i anotant les coincidèndies a cada posició. Si es vol tenir una explicació detallada del seu funcionament, veure el codi del programa.

Un cop obtingutda la next table, tot i ser de vital importància, encara faltaran resoldre moltes altres qüestions a tenir en compte per a complir els objectius del programa.

Algorisme de Knuth-Morris-Pratt

La lògica fonamental d'aquest algorisme i el perquè del seu guany en eficiència respecte a l'algorisme simple de correspondència exacta està detallat en l'apartat de Metodologia.

La seva implementació en el si del programa és en dues funcions: una funció anomenada knuth, per comprovar que els patrons siguin únics dins del propi mRNA; i un altra funció, knuth4genome, per tal de fer la mateixa comprovació en el genoma sencer.

La primera funció, knuth, rep tres paràmetres: la next table ja calculada, la seqüència mRNA i el patró. I retorna 1 (cert) o 0 (fals) segons si el patró a comprovar és únic o no sobre la seqüència de mRNA.

La segona funció, knuth4genome, rep quatre paràmetres: la next table, l'adreça del directori on es localitza el genoma, el patró a comprovar i la llargada del patró. I retorna 1 (cert)* o 0 (fals) segons si el patró a comprovar és únic o no sobre el genoma, i en el cas de ser únic, també retorna el desplaçament i el cromosoma en què es localitza.

* La funció retornarà 1 quan el comptador de patrons únics sigui igual o més petit que 1. Això es fa d'aquesta manera per aquells patrons únics en el mRNA però que no es troben en el genoma, degut a que es pot trobar separat en dos exons diferents amb introns entremig.

Tanmateix, en aquesta segona funció, abans d'implementar l'algorisme propi KMP, cal tenir present certes preparacions:

Per tal de comprovar si el patró és únic en tot el genoma, cal que aquest últim es trobi contingut dins d'un directori. El què es fa és obrir el directori a partir de la funció opendir del Perl, i emmagatzemar tots els fitxers de cromosomes que contingui en un vector, excepte els fitxers ocults "." i "..":

my @fitxerscromosomes = ();

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

Després, mitjançant una iteració amb un while, es comprova que el patró sigui únic cromosoma per cromosoma. En cada cas, s'obre el cromosoma que toca, i es realitza la comprovació línia per línia. També cal recordar concatenar a l'inici de cada línia els últims nucleòtids de l'anterior línia (segons la llargada del patró).

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;
}

Finalment, algunes subseqüències úniques poden localitzar-se a l'strand negatiu de la seqüència dels cromosomes. Per tal de solucionar aquest fet, es fa el "reverse complement" de la seqüència del patró abans de començar a buscar:

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

I quan s'està fent la cerca línia per línia, si a la línia actual no es troba, aleshores es prova de repetir la mateixa cerca però amb el patró en reverse complement, i si tampoc el troba es segueix amb el procediment habitual de l'algorisme.

De totes maneres, tant en una funció com en l'altre, l'algorisme propi del KMP és el següent:


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 fonament d'aquest és recórrer la seqüència en què volem comprovar si el patró és únic o no, tot comparant-la amb el patró. Però a l'hora de modificar el desplaçament segons si hem trobat el patró exacte o no, es fa mitjançant la informació proporcionada per la next table. Així directament s'aconsegueix desplaçar suficientment el patró sobre la seqüència, fins al punt d'estalviar-se posicions de la seqüència que no cal comparar i també, iniciar la següent comparació en la posició concreta del patró (saltant-se les posicions del patró que tampoc cal comparar). Si es vol tenir una explicació detallada del seu funcionament, veure el codi del programa.

Generació de patrons

El programa és capaç per ell mateix d'obrir el fitxer FASTA del mRNA i generar tots les patrons possibles que se'n puguin derivar, mitjançant la funció substr del Perl, per a una llargada determinada. D'aquesta manera, s'assegura que el programa analitzarà totes les subseqüències susceptibles de ser úniques d'un gen determinat.

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

Estratègia automàtica de moure la llargada

Un cop comprovat si hi ha patrons únics per una llargada concreta en el mRNA, hi ha dues opcions possibles: que se n'hagin trobat o que no.

En el cas que no s'hagin trobat patrons únics, sumarem una unitat a la llargada i es tornarà a fer la cerca de subseqüències úniques a partir del mRNA.

Si se n'han trobat, emmagatzemarem els patrons, el seu desplaçament i el cromosoma, amb tres vectors resultat. I procedirem a tornar fer la cerca de subseqüències úniques però per una llargada d'una unitat inferior. A més a més, en aquest cas no s'efectuarà la cerca a partir de tot el mRNA sencer, sinó a partir dels patrons únics trobats amb la llargada de la qual provenim, ja que qualsevol patró únic d'una llargada inferior (per exemple de 13) estarà contingut dins d'alguns dels patrons únics trobats amb llargada 14.

Consideracions:

En el cas que haguem anat a una llargada inferior i trobem patrons, continuarem anant a llargades inferiors. Pero en el moment que no trobem patrons únics, ens quedarem amb els resultats anteriors, prèviament emmagatzemats als tres vectors resultat. I sortirem de la iteració, mitjançant l'avaluació de la variable $tinc_patro_unic com a certa (1), ja que en aquest cas s'entrarà a l'"elsif" detallat més avall.
L'altre opció és que haguem anat a una llargada superior. Si no trobem patrons, continuarem anant a llargades superiors. En el moment que trobem patrons únics, els emmagatzemarem als tres vectors de resultats, i com que la llargada anterior serà sempre inferior a l'actual (perquè haurem incrementat la llargada), sortirem de la iteració també avaluant com a certa la variable $tinc_patro_unic (1), en el segon "if" detallat a continuació.

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;
}