Erhöhung der Effizienz der Tabelle Schnittpunkt

Question

Ich habe 2 Tabellen.Sie sind beide in der Form Chromosom, Anfang und Ende Koordinaten auf diesem Chromosom. Die erste Tabelle enthält Gene und die zweite Tabelle enthält kurze Sequenzen, die in diese Gene fallen können oder nicht. In meinem realen Datensatz sind Gene ungefähr 50.000 Zeilen und Sequenzen sind ungefähr 7.000.000 Zeilen und beide Tabellen haben verschiedene zusätzliche Spalten. Ich würde gerne Überschneidungen zwischen den beiden Tabellen finden.

chromosome=as.character(rep(c(1,2,3,4,5), each=10000)) 
start=floor(runif(50000, min=0, max=50000000)) 
end=start+floor(runif(10000, min=0, max=10000)) 
genes=cbind(chromosome, start, end) 

startseq=floor(runif(7000000, min=0, max=50000000)) 
endseq=startseq+4 
sequences=cbind(chromosome, startseq, endseq) 

Ich versuche, alle intersects zu finden mit:

for (g in 1:nrow(sequences)) { 
    seqrow=as.vector(sequences[g,]) 
    rownr=which(genes[,1]==seqrow[1] & genes[,2] < seqrow[2] & genes[,3] > seqrow[3]) 
    print(rownr) 
} 

Ich beabsichtige, diese Zeilennummern zu verwenden, um Aktionen auf den zusätzlichen Spalten ich in meinem realen Datensatzes zu erfüllen haben. Das Problem im Moment ist, dass der beschriebene Prozess eher langsam ist. Auf welche Weise könnte ich diesen Schnitt beschleunigen?

Answer 1

Sie möchten bioconductor für diese Aufgabe und speziell das GenomicRanges Paket verwenden. Dies gibt ein Objekt der Klasse "Hits" zurück, das die Indizes der Überlappungen enthält. Sie können auch die Funktion intersect verwenden, aber dies gibt das Intervall des Schnittpunkts und nicht die ID des Schnittpunkts zurück. Kurz gesagt, hat Bioconductor und GenomicRanges viele nützliche Set-Funktionen, die ziemlich schnell sind.

## try http:// if https:// URLs are not supported 
source("https://bioconductor.org/biocLite.R") 
biocLite() 
biocLite("GenomicRanges") ## I think genomicranges is part of the standard bioconductor install but if not this will install it. 


library(GenomicRanges) 

set.seed(8675309) 
chromosome <- as.character(rep(c(1,2,3,4,5), each=10000)) 
start <- floor(runif(50000, min=0, max=50000000)) 
end <- start+floor(runif(10000, min=0, max=10000)) 
genes <- cbind(chromosome, start, end) 

startseq <- floor(runif(7000000, min=0, max=50000000)) 
endseq <- startseq+4 
chromosome <- sample(c(1,2,3,4,5), size = 7000000, replace=T) 
sequences=cbind(chromosome, startseq, endseq) 

genes <- GRanges(seqnames = chromosome, ranges = IRanges(start = start, end = end)) 
seqs <- GRanges(seqnames = chromosome, ranges = IRanges(start = startseq, end = endseq)) 

x <- findOverlaps(seqs, genes) 
head(x) 

#Hits object with 6 hits and 0 metadata columns: 
#  queryHits subjectHits 
#  <integer> <integer> 
# [1]   2  41673 
# [2]   2  47476 
# [3]   3  20048 
# [4]   4  9624 
# [5]   4  5662 
# [6]   4  1531 
# ------- 
# queryLength: 7000000 
# subjectLength: 50000