Повышение эффективности таблицы пересекается

Question

У меня есть 2 таблицы. Они оба находятся в форме хромосомы, начальные и конечные координаты на этой хромосоме. Первая таблица содержит гены, а вторая таблица содержит короткие последовательности, которые могут или не могут попадать в эти гены. В моем реальном наборе данных гены составляют около 50 000 строк, а последовательности - около 7.000.000 строк, а обе таблицы имеют дополнительные столбцы. Я хотел бы найти перекрытия между двумя таблицами.

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) 

Я пытаюсь найти все пересекается с помощью:

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

я намерен использовать эти номера строк для выполнения действий дополнительных столбцов есть в моих реальных данных. На данный момент проблема заключается в том, что описанный процесс довольно медленный. Каким образом я мог бы ускорить это пересечение?

Answer 1

Вы хотите использовать bioconductor для этой задачи, и в частности пакет GenomicRanges. Это вернет объект класса «Хиты», который будет содержать индексы перекрытий. Вы также можете использовать функцию intersect, но это возвращает интервал пересечения, а не идентификатор пересекающегося seq. Короче говоря, bioconductor и GenomicRanges имеет много полезных функций набора, которые довольно быстры.

## 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