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?