Algoritmo per una ricerca veloce

1

Sto cercando un algoritmo (pseudocodice, codice R) per il mio problema qui. Sembra molto da leggere ma in realtà è molto facile da leggere e sto solo cercando di essere dettagliato.

Ho una lista (diciamo coords ) con 8 milioni di elementi. Sembrano così:

> head(coords)
  cstart  cend
1  35085 35094
2  35089 35098
3  35374 35383
4  35512 35521
5  35513 35522
6  35514 35523
...
98 12600309 12600318 

Ciò che questi elementi rappresentano sono le coordinate su una stringa di DNA. Ad esempio, 35374 si riferisce alla posizione sulla stringa DNA. Non mi interessa se è un A C G T per quello che vale. E TUTTI i coords sono di lunghezza 10. Quindi il campo cend può essere calcolato. (In realtà non uso questo cend nel mio codice)

Ho una lista secondaria (diciamo alignments ) (~ 100 elementi) che ha anche una coordinata start e end . Questo definisce area di interesse . Sembra questo

 chr_name       chr_strand      start_index        end_index
  "chr1"              "+"       "12600311"          "12600325"

Le colonne chr_name, chr_strand non sono importanti. In questo caso, la nostra "area di interesse" potrebbe essere maggiore di 10.

Il mio problema

Ora, quello che voglio fare è esaminare tutti gli 8 milioni di record e vedere se ogni record è contenuto nell'elenco interests (fino all'80% di sovrapposizione, lo spiegherò). Per l'esempio precedente, i primi 6 elementi non sono "vicini" all'area di interesse.

La riga 98 avvia 2 "coppie di basi" prima della nostra area di interesse, MA si sovrappone all'area di interesse fino all'80%. Vale a dire, l'80% del "sito" è contenuto nell'area di interesse.

Risultato

Fondamentalmente, sto cercando una lista (posso preallocarla), dove ogni elemento è un VERO / FALSO. Quindi da sopra avrò una lista come la seguente

FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, ..., TRUE (97'th index)

Il mio codice R attuale Ho documentato il codice R che è molto facile da leggere

  ## mdply is a function that acts essentially like a loop
  ## the .data is the data I am working with.. the 8 million records
  ## for each row, it executes a function where the signature of the function
  ## matches the column names. 
  ## the best part about this .. IS THAT ITS PARALLELIZED
  ## and I can use up to 24 cores to do this.
  b <- mdply(.data = coords, 
             function(cstart, cend){   

               #the function is executed for each row. 
               #loop through our "area of interests"
               for(j in 1:nrow(alignments)){
                 chr_info <- extract_information_from_alignment(alignments[j, ])
                 interesting <- 0
                 uninteresting <- 0
                 hits_index_start <- cstart
                 hits_length <- 10 ## the length of our matrix               
                 aoe_start <- as.numeric(chr_info[["start_index"]])
                 aoe_end <- chr_info[["end_index"]]    

                 ## go through EACH INDEX AND CALCULATE THE OVERLAP
                 for(j in hits_index_start:(hits_index_start + hits_length-1)){
                   if((aoe_start <= j) && (j <= aoe_end)){
                     interesting <- interesting + 1
                   } else {
                     uninteresting <- uninteresting + 1
                   }          
                 } #end for - interesting calc
                 substr_score = interesting/sum(interesting + uninteresting)
                 if(substr_score >= substr_score_threshold){
                   return(TRUE)
                 } else { 
                   return(FALSE)
                 }      
               } #end for
             }, .parallel = TRUE)

Ho eseguito questo codice la scorsa notte per 13 ore e non è stato completato. Stimo che ci vogliono circa 16 ore ..

idee Penso che il collo di bottiglia sia che per ogni riga, ha un ciclo for (100 elementi) e all'interno di quello for loop, c'è un secondo che calcola la sovrapposizione.

Ma sto pensando, non ho bisogno di passare attraverso OGNI allineamento. Dovrei prima verificare se il mio cstart è vicino all'area di interesse.

Forse una sorta di ricerca binaria?

    
posta masfenix 12.02.2015 - 21:36
fonte

1 risposta

0

Come hai notato, i loop annidati possono diventare costosi. 8 milioni di iterazioni del ciclo esterno moltiplicano 100 iterazioni dei cicli centrali (non sono sicuro) le iterazioni del ciclo interno portano a un sacco di elaborazione.

Penso che tu sia sulla strada giusta però - cerca di eliminare il maggior numero possibile di quegli 8 milioni di coord prima che eseguono qualsiasi elaborazione su di essi. Ricerca binaria è sicuramente tuo amico.

L'approccio che prenderei sarebbe approssimativamente simile a

sort the coords list // (done once, use a library function, happens fast)
create your output array, initialize all to 'false'
for each alignment in alignments // this loop should iterate 100x
  // the first alignment value that would fit your 80% coverage
  // would be alignment.start_index - 2 (since all coords are of length 10
  // and you need an 80% overlap)
  let searchStart = alignment.start_index - 2
  call a binary search function for searchStart in the coords list
  // you'll either find the first area of interest, or the location
  // in the list where that first one would go.  

  // do the same for the search end (since the alignment appears to be variable length)
  // and because iterating over the list of coords from the starting point
  // would probably be slower

Mi aspetto che ciò avvenga in modo non parallelo abbastanza rapidamente - di certo non richiede ore.

Tutto ciò che sta facendo è chiamare la ricerca binaria su un elenco di voci da 8 milioni, 100 o 200 volte (a seconda se si cerca end_index o iterato).

    
risposta data 12.02.2015 - 22:04
fonte

Leggi altre domande sui tag