[Vari] Contest 4: DNA

[gugoXX]

Siano date 2 stringhe molto lunghe, di lunghezza diversa, contenenti solo i caratteri CGAT

CGTGTCGCGATACATAGTACaggagtaaagagatGaggaTacgATAGCGATACGATACAGATGTCGGAGAATC
TTCCTACGCAGaggagtaaagagatCaggaCacgGAGCTGATCGGTCGACCATCATTATA

Domanda 1:
Quanto e' lunga la piu' lunga sottosequenza di DNA comune alle 2 stringhe, e l'indice di partenza della sottosequenza in entrambe le stringhe.
Nell'esempio la risposta e' (dovrebbe) 14, con indici 20 e 11, in quanto la sottostringa aggagtaaagagat, nell'esempio riportata per facilita' in minuscolo, dovrebbe essere la piu' lunga comune.

Domanda 2:
Durante una DNA polimerasi possono generarsi errori di trascrizione.
Detto ERR il numero di errori massimi ammessi, trovare la lunghezza della piu' lunga sottosequenza di DNA comune alle 2 stringhe, e l'indice di partenza in ciascuna di esse.

Gli errori NON potranno essere consecutivi. 2 errori consecutivi invalideranno comunque la prosecuzione della sequenza.
Inoltre una sottosequenza non puo' iniziare o finire con un errore.

Nell'esempio precedente, dato ERR=1 la risposta e' (dovrebbe) 19, sempre con indici 20 e 11, in quanto aggagtaaagagatGaggaTacg differisce da aggagtaaagagatCagga per un solo carattere.
Con ERR=2 la riposta e' (dovrebbe) 23 (analogo), sempre con 20 e 11 come indici.

Le 2 stringhe di TEST arriveranno presto, e si suppone che siano molto, molto lunghe.

Ecco le 2 stringhe di TEST
http://www.usaupload.net/d/s4b7z8fpakd
Sono da circa 100.000 nucleotidi la prima e da 80.000 nucleotidi la seconda.
E i 2 test da provare sono quindi:
1) Trovare la piu' lunga sottosequenza identica (prototipo: Funzione (stringa1, stringa2))
2) Trovare la piu' lunga sottosequenza comprendente un massimo di 2 errori (Prototipo: Funzione (stringa1, stringa2, ERR) )

Proviamo cosi', se poi i risultati non ci soddisfano cercheremo 2 nuove sequenze di DNA.
Per il secondo esercizio consideriamo che il numero di errori sara' comunque sempre ragionevolmente basso.

[cionci]

Ci puoi dare un'idea di quanto siano lunghe ?

[gugoXX]

Originariamente inviato da cionci

Ci puoi dare un'idea di quanto siano lunghe ?

Devo fare prima una stima sulle tempistiche.
Direi che 200.000 caratteri dovrebbe essere una buona partenza, soprattutto per il secondo algoritmo.

[shinya]

Non vorrei rompere le uova nel paniere... ma questo purtroppo è un problema noto e risolto in letteratura.

Non vorrei sfociasse in "la mia versione in C ha un bug ma è velocissima" vs. "L'ho risolto in 3 righe in Blurp!"...

[gugoXX]

Originariamente inviato da shinya

Non vorrei rompere le uova nel paniere... ma questo purtroppo è un problema noto e risolto in letteratura.

Non vorrei sfociasse in "la mia versione in C ha un bug ma è velocissima" vs. "L'ho risolto in 3 righe in Blurp!"...

Quello senza errori lo so... ma anche con gli errori?

[shinya]

Originariamente inviato da gugoXX

Quello senza errori lo so... ma anche con gli errori?

Non sono decisamente un esperto nel campo... ma se faccio una ricerca su citeseer vengono fuori mille-mila paper sull'argomento (con e senza errori)...
Esistono ovviamente tantissime varianti, quindi se non ci si concentra su una sola soluzione facendo a gara di velocità/bug potrebbe venir fuori un thread interessante

[Giullo]

appena rientro a casa gli dò un'occhiata seria.

ah, bella gugo, ormai i tuoi contest sono sono uno dei motivi principali di visita di questo forum (per me chiaramente)

[rеpne scasb]

Originariamente inviato da gugoXX

Siano date ...[CUT]

Vedo un nuovo contest. Spero che cio' non voglia dire che ti sei ritirato dal contest 3.

[gugoXX]

Dai, inizio a puntare io
Per quanto riguarda il secondo esercizio, punto
una stringa lunga 29 caratteri, con 2 errori
trovata nel vergognoso tempo di 174 secondi = 3 minuti circa

[sottovento]

Originariamente inviato da shinya

Non vorrei rompere le uova nel paniere... ma questo purtroppo è un problema noto e risolto in letteratura.

Non vorrei sfociasse in "la mia versione in C ha un bug ma è velocissima" vs. "L'ho risolto in 3 righe in Blurp!"...

Hai ragione, ma... perche' "purtroppo"? Una volta tanto che esiste un'ottima soluzione, e possiamo metterla nel nostro cassetto degli attrezzi...

L'algoritmo e' l'LCS, funziona benissimo ed e' adatto a questo tipo di problemi. Fra l'altro, e' facilmente adattabile anche per risolvere il problema successivo, cioe' quello di trovare la sequenza piu' lunga contenente un numero massimo di errori (basta un intero).
L'LCS e' l'algoritmo usato per la diff di Unix, per intenderci, ed e' ampiamente documentato in rete

[Vincenzo1968]

Per la domanda 1 sto utilizzando l'algoritmo di Knuth-Morris-Pratt (il quale si basa sul concetto di automi a stati finiti che adoro ).

L'idea è quella di partire da una stringa lunga la metà di quella di lunghezza minore. Se la stringa non viene trovata divido ancora per due e così via fino a quando non la trovo(se, invece, la trovo, non divido ma aumento la lunghezza della stringa di ricerca).
La ricerca prosegue sulle stringhe di lunghezza maggiore(minore) di quest'ultima e minore(maggiore) dell'ultima ricerca infruttuosa.

Per ora l'ho provato sulle stringhe di esempio e funziona:



questo è il codice:

#include <stdio.h>
#include <string.h>
#include <malloc.h>

int next[1000];

void InitNext(char *p)
{
	int i, j;
	int M = strlen(p);

	next[0] = -1;

	for ( i = 0, j = -1; i < M; i++, j++, next[i] = j )
	{
		while ( (j >= 0) && (p[i] != p[j]) )
			j = next[j];
	}
}

int KnuthMorrisPratt_search(char *p, char *a)
{
	int i, j;
	int M = strlen(p);
	int N = strlen(a);

	InitNext(p);

	for ( i = 0, j = 0; j < M && i < N; i++, j++ )
	{
		while ( (j >= 0) && (a[i] != p[j]) )
		{
			j = next[j];
		}
	}

	if ( j == M )
		return i - M;
	else
		return i;
}

int main()
{
	int p1_len, p2_len;
	char * strSearch;
	size_t len;

	int pos1, pos2;

	int bTrovato = 0;

	char * p1 = "CGTGTCGCGATACATAGTACaggagtaaagagatGaggaTacgATAGCGATACGATACAGATGTCGGAGAATC";
	char * p2 = "TTCCTACGCAGaggagtaaagagatCaggaCacgGAGCTGATCGGTCGACCATCATTATA";

	p1_len = strlen(p1);
	p2_len = strlen(p2);

	//len = p2_len/2;
	len = p2_len;

	strSearch = (char*)malloc(sizeof(char)*len);
	if ( !strSearch )
	{
		printf("Memoria insufficiente.\n");
		return -1;
	}
	strSearch[0] = '\0';

	memcpy(strSearch, p2, len);
	*(strSearch + len) = '\0';

	bTrovato = 0;
	pos1 = 0;
	pos2 = 0;
	while ( len > 0 )
	{
		while ( p2_len - (pos2 + len) > 0 )
		{
			pos1 = KnuthMorrisPratt_search(strSearch, p1);
			if ( pos1 < p1_len )
			{
				bTrovato = 1;
				printf("\nStringa '%s' di lunghezza %d trovata nelle posizioni %d e %d\n", strSearch, len, pos1, pos2);
				break;
			}
			else
			{
				pos2++;
				memcpy(strSearch, p2 + pos2, len);
				*(strSearch + len) = '\0';
			}
		}
		if ( bTrovato )
			break;
		len--;
		pos1 = 0;
		pos2 = 0;
	}

	free(strSearch);

	return 0;
}

[VICIUS]

Carino questo problema sul dna
Per ora sto cercando di implementare l'algoritmo che si trova su wikipedia che purtroppo non è il massimo della velocità. Ho provato con stringhe di 1000 caratteri è piuttosto lentino e occupa una vagonata di ram. Con i file originali non ci ho neanche provato perché verrebbe fuori una array di 8 miliardi di elementi.

def initialize_lcs_matrix(s, t)
  l = Array.new(s.size + 1) { Array.new(t.size + 1, 0)}
  z = 0
  for i in 1..s.size
    for j in 1..t.size
      if s[i-1] == t[j-1]
        l[i][j] = l[i-1][j-1] + 1
        if l[i][j] > z
          z = l[i][j]
        end
      end
    end
  end
  return l
end

def get_index_from_matrix(m)
  x = y = z = 0
  for i in 0..m.size-1
    for j in 0..m[i].size-1
      if m[i][j] > z
        z = m[i][j]
        x = i
        y = j
      end
    end
  end
  return [z, x, y]
end

def get_string(a, s, e)
  string = ""
  for i in s..(s+e)
    string += a[i]
  end
  return string
end

def read_dna(file_name)
  dna = Array.new
  File.open(file_name, 'r') do |f|
    f.each do |s|
      s.scan(/./).each do |c|
        dna.push c
      end
    end
  end
  return dna
end

ary_a = read_dna('DNA1.txt')
ary_b = read_dna('DNA2.txt')

lcsm = initialize_lcs_matrix(ary_a, ary_b)
z, x, y = get_index_from_matrix(lcsm)
seq = get_string(ary_a, x, z)

puts "Trovata la sequenza #{seq}"
puts "Che è lunga #{z} caratteri"
puts "Indice nella stringa a #{x}"
puts "Indice nella stringa b #{y}"

E questo è l'output:

Trovata la sequenza ATAGCTAACACGC
Che è lunga 12 caratteri
Indice nella stringa a 832
Indice nella stringa b 767

Tempo di esecuzione circa 6 secondi.

[shinya]

Originariamente inviato da VICIUS

Carino questo problema sul dna
Per ora sto cercando di implementare l'algoritmo che si trova su wikipedia che purtroppo non è il massimo della velocità. Ho provato con stringhe di 1000 caratteri è piuttosto lentino e occupa una vagonata di ram.

Stai usando l'algoritmo O(N1 * N2). Per andare forte bisogna usare i suffix tree.

[cionci]

Mi ero fatto un algoritmo carino, però appena si sale sopra 8000 elementi esaurisce la memoria. Il modo per terminare la ricerca lo so già, suddividendo il problema, ma credo che mi ci voglia un bel po' a mettere su tutto.
Per 8000 elementi:

time ./Contest4
8000 8000
First string position: 7782 Matched string: ttaggtggtggtaaag
Second string position: 6381 Matched string: ttaggtggtggtaaag

real 0m4.230s
user 0m1.828s
sys 0m1.072s

[DanieleC88]

Ho provato a risolvere il primo pezzo del problema: su questo arcaico computer ci mette la bellezza di 13 minuti...

Ecco il codice del programma:

from time import clock

def leggiFile(percorso):
	f = open(percorso)
	
	stringa = ""
	for linea in f:
		stringa += linea
	f.close()
	return stringa.capitalize()

def trovaSequenza(stringa1, stringa2):
	x = 0
	coordinate = ()
	lunghezza = 0
	sequenza = ""
	while x + lunghezza < len(stringa1):
		sequenza += stringa1[x + lunghezza]
		try:
			y = stringa2.index(sequenza)
			coordinate = (x, y)
			lunghezza += 1
		except ValueError:
			x += 1
			sequenza = sequenza[1:]
	return lunghezza, coordinate

def msec(t1, t2):
	return (t2-t1)*1000.0

origine = "DNA"
print "Lettura in corso...",
t1 = clock()
stringa1 = leggiFile(origine + "1.txt")
stringa2 = leggiFile(origine + "2.txt")
t2 = clock()
print "completata in", msec(t1, t2), "millisecondi."
print

print "Cerco la stringa di lunghezza massima..."
t1 = clock()
lunghezza, coordinate = trovaSequenza(stringa1, stringa2)
t2 = clock()

x, y = coordinate
print " ->", lunghezza, "caratteri alla posizione", coordinate
print " -> sottosequenza della stringa 1:", stringa1[x:(x+lunghezza)]
print " -> sottosequenza della stringa 2:", stringa2[y:(y+lunghezza)]
print "completato in", msec(t1, t2), "millisecondi."

Ed ecco l'output che ottengo (senza la stampa delle due sottosequenze perché l'ho aggiunta solo ora, di aspettare altri 13 minuti non mi va ):

Lettura in corso... completata in 252.805109036 millisecondi.

Cerco la stringa di lunghezza massima...
 -> 21 caratteri alla posizione (10002, 40002)
completato in 777964.344022 millisecondi.

[sottovento]

Provato anch'io, limitando le stringhe ad 8000 caratteri.

Usato l'algoritmo LCS sotto Java, sul mio laptop DELL Latitude D830. Risultato: 800 millisecondi

Non mi aspetto grosse variazioni di tempo per quanto riguarda la seconda parte del contest

[grigor91]

Originariamente inviato da DanieleC88

Ho provato a risolvere il primo pezzo del problema: su questo arcaico computer ci mette la bellezza di 13 minuti...

Ecco il codice del programma:

code ...

ho provato a eseguirlo sul mio pc e i tempi sono questi:

Lettura in corso... completata in 10.0537917529 millisecondi.

Cerco la stringa di lunghezza massima...
 -> 21 caratteri alla posizione (10002, 40002)
 -> sottosequenza della stringa 1: actgtcctgtcaacaaggagt
 -> sottosequenza della stringa 2: actgtcctgtcaacaaggagt
completato in 38560.2209219 millisecondi.

[gugoXX]

Originariamente inviato da sottovento

Provato anch'io, limitando le stringhe ad 8000 caratteri.

Usato l'algoritmo LCS sotto Java, sul mio laptop DELL Latitude D830. Risultato: 800 millisecondi

Non mi aspetto grosse variazioni di tempo per quanto riguarda la seconda parte del contest

Il problema dell'LCS e' che richeide una matrice rettangolare MxN, e M e N per le sequenze di DMA sono taaaaaanto lunghe
Scherzi a parte, come pensi di conciliare LCS con la trattazione degli errori ammessi nelle stringhe?

[cionci]

Solo per la prima parte con 200000 caratteri nella stringa uno e 200000 nella stringa due:

~$ time Contest4
200000 200000
First string position: 96827 Matched string: aaaaggtagtttggtagtggatt
Second string position: 134120 Matched string: aaaaggtagtttggtagtggatt

real 0m21.040s
user 0m21.001s
sys 0m0.016s

Algoritmo fatto in casa basato su un albero che enumera totalmente tutte le possibili stringhe.
Metto un po' in ordine il codice e lo posto.

[gugoXX]

Originariamente inviato da Vincenzo1968

Per la domanda 1 sto utilizzando l'algoritmo di Knuth-Morris-Pratt (il quale si basa sul concetto di automi a stati finiti che adoro ).

Secondo me e' una strada molto promettente.
Non la sto percorrendo (e infatti pago), ma ritengo che sia una delle migliori. Facci sapere...