Эффективная скоординировать матч между data.frames в R с sapply

Question

Я пытаюсь получить вектор рассказывал мне, какие строки в data.frame (transcriptcoords)

   chr start end 
NONHSAT000001 chr1 11868 14409 
NONHSAT000002 chr1 11871 14412 
NONHSAT000003 chr1 11873 14409 
NONHSAT000004 chr1 12009 13670 
NONHSAT000005 chr1 14777 16668 
NONHSAT000006 chr1 15602 29370 

имеют начало/конец координат свободно содержащийся (с 10 +/- о толерантности) в другом data.frame (genecoords)

   chr start end 
NONHSAG000001 chr1 11869 14412 
NONHSAG000002 chr1 14778 29370 
NONHSAG000003 chr1 29554 31109 
NONHSAG000004 chr1 34554 36081 
NONHSAG000005 chr1 36273 50281 
NONHSAG000006 chr1 62948 63887 

Чтобы сделать это, я зацикливание с sapply над рядами индексов первого data.frame, соответствующие координаты с любой строки в второй data.frame. У меня есть решение (как описано ниже), но это, кажется, довольно медленно (около шести секунд с ломтиком 2000 строк):

user system elapsed 
    6.02 0.00 6.04 

Я пытаюсь понять, какие части sapply могут быть оптимизированы. Это блок if/else? Или сравнения строк (==, < =,> =)? Или, проще говоря, является ли он внутренне медленным алгоритмом?

Спасибо! Код я сгенерировал ниже:

load(url("http://www.giorgilab.org/stuff/data.rda")) 

# Pre-vectorize the data frames 
g0<-rownames(genecoords) 
g1<-genecoords[,1] 
g2<-as.integer(genecoords[,2]) 
g3<-as.integer(genecoords[,3]) 

t0<-rownames(transcriptcoords) 
t1<-transcriptcoords[,1] 
t2<-as.integer(transcriptcoords[,2]) 
t3<-as.integer(transcriptcoords[,3]) 

system.time(gs<-sapply(1:2000,function(i){ 
      t<-t0[i] 
      chr<-t1[i] 
      start<-t2[i] 
      end<-t3[i] 

      # Find a match (loose boundaries +/- 10) 
      right1<-which(g1==chr) 
      right2<-which(g2<=start+10) 
      right3<-which(g3>=end-10) 
      right<-intersect(right3,intersect(right1,right2)) 
      right<-g0[right] 

      if(length(right)==1){ 
       g<-right 
      } else if(length(right)>1){ 
       # Get the smallest match 
       heregenecoords<-genecoords[right,] 
       size<-apply(heregenecoords,1,function(x){abs(as.numeric(x[3])-as.numeric(x[2]))}) 
       g<-names(which.min(size)) 
      } else { 
       g<-t 
      } 
      return(g)   
     } 
)) 

Answer 1

с данными

tx0 <- read.table(textConnection("chr start end 
NONHSAT000001 chr1 11868 14409 
NONHSAT000002 chr1 11871 14412 
NONHSAT000003 chr1 11873 14409 
NONHSAT000004 chr1 12009 13670 
NONHSAT000005 chr1 14777 16668 
NONHSAT000006 chr1 15602 29370")) 

gene0 <- read.table(textConnection("chr start end 
NONHSAG000001 chr1 11869 14412 
NONHSAG000002 chr1 14778 29370 
NONHSAG000003 chr1 29554 31109 
NONHSAG000004 chr1 34554 36081 
NONHSAG000005 chr1 36273 50281 
NONHSAG000006 chr1 62948 63887")) 

GenomicRanges пакет в Bioconductor делает это легко и эффективно (для миллионов наложений).

library(GenomicRanges) 
tx <- with(tx0, GRanges(chr, IRanges(start, end))) 
gene <- with(gene0, GRanges(chr, IRanges(start, end))) 

## increase width by 10 on both sides of the center of the gene range 
gene <- resize(gene, width=width(gene) + 20, fix="center") 
## find overlaps of 'query' tx and 'subject' gene, where query is within subject 
olaps <- findOverlaps(tx, gene, type="within") 

Показано, например, что 'запрос' (ТХ) 1, 2, 3 и 4 находятся в пределах 'субъект' (генов) 1.

> findOverlaps(tx, gene, type="within") 
Hits of length 6 
queryLength: 6 
subjectLength: 6 
    queryHits subjectHits 
    <integer> <integer> 
1   1   1 
2   2   1 
3   3   1 
4   4   1 
5   5   2 
6   6   2 

и что ген 1 перекрывается 4 транскрипта, транскрипты гена 2 на 2.

> table(subjectHits(olaps)) 

1 2 
4 2 

См. Также publication. Используя больший набор данных:

tx <- with(transcriptcoords, GRanges(V1, IRanges(V2, V3, names=rownames(tx0)))) 
gene <- with(genecoords, GRanges(V1, IRanges(V2, V3, names=rownames(gene0)))) 

с некоторыми задержками:

system.time(gene <- resize(gene, width=width(gene) + 20, fix="center")) 
## user system elapsed 
## 0.056 0.000 0.057 
system.time(findOverlaps(tx, gene, type="within")) 
## user system elapsed 
## 2.248 0.000 2.250 

Я думаю, что это примерно время для data.table решения от @ danas.zuokos с только 1000 транскриптов

system.time({ 
    dt <- genecoords[transcriptcoords, allow.cartesian = TRUE]; 
    res <- dt[start <= start.1 + tol & end >= end.1 - tol, 
     list(gene = gene[which.min(size)]), by = transcript] 
}) 
## user system elapsed 
## 2.148 0.244 2.400 

Answer 2

Ha! Мартин победил меня с лучшим ответом. Практически всегда лучше использовать проверенный на работоспособность код в хорошо установленной библиотеке, а не сворачивать свои собственные. Определенно используйте решение Мартина, а не это.

Но, просто для смеха, вот еще один способ сделать это.

Во-первыхи, сделать некоторые гены и распечатки:

gs = 1:10*500 
genes = data.frame(start=gs, end=gs+400) 
rownames(genes) = sprintf('g%05d', 1:nrow(genes)) 

ts = sample(1:max(genes$end), size=10) 
transcripts = data.frame(start=ts, end=ts+60) 
rownames(transcripts) = sprintf('t%05d', 1:nrow(transcripts)) 

Мы можем векторизацию сравнения с использованием outer который применяет функцию к каждой комбинации два векторных аргументов.

overlaps = function(genes, transcripts, min_overlap=1) { 
    b1 = outer(genes$end, transcripts$start, min_overlap=min_overlap, 
      function(e,s,min_overlap) e-s+1>min_overlap) 
    b2 = outer(genes$start, transcripts$end, min_overlap=min_overlap, 
      function(s,e,min_overlap) e-s+1>min_overlap) 
    result = b1 & b2 
    rownames(result) = rownames(genes) 
    colnames(result) = rownames(transcripts) 
    return(result) 
} 

Для наших генов и транскриптов, мы могли бы получить что-то вроде:

> genes 
     start end 
g00001 500 900 
g00002 1000 1400 
g00003 1500 1900 
g00004 2000 2400 
g00005 2500 2900 
g00006 3000 3400 
g00007 3500 3900 
g00008 4000 4400 
g00009 4500 4900 
g00010 5000 5400 

> transcripts 
     start end 
t00001 535 595 
t00002 2502 2562 
t00003 4757 4817 
t00004 3570 3630 
t00005 3094 3154 
t00006 1645 1705 
t00007 5202 5262 
t00008 13 73 
t00009 788 848 
t00010 4047 4107 

o1 = overlaps(genes, transcripts, 1) 

Результата является булевой матрицей, которая говорит вам, перекрывает ли каждая запись с каждым геном.

> o1 
     t00001 t00002 t00003 t00004 t00005 t00006 t00007 t00008 t00009 t00010 
g00001 TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE 
g00002 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
g00003 FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE 
g00004 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
g00005 FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
g00006 FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE 
g00007 FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE 
g00008 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE 
g00009 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
g00010 FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE 

Answer 3

Я пользуюсь data.table библиотекой.

rm(list = ls()) 
load(url("http://www.giorgilab.org/stuff/data.rda")) 
library(data.table) 
tol <- 10 # tolerance 
id <- 1:2000 # you can comment this out, but make sure you have big RAM 

Преобразование в data.table формат. Дополнительный вычислительный размер (я не уверен, почему вы принимаете abs, а не конец всегда больше, чем старт?).

genecoords <- data.table(genecoords, keep.rownames = TRUE) 
setnames(genecoords, c("gene", "chr", "start", "end")) 
genecoords[, size := end - start] 
transcriptcoords <- data.table(transcriptcoords, keep.rownames = TRUE)[id] # comment out [id] when running on all trascripts 
setnames(transcriptcoords, c("transcript", "chr", "start", "end")) 

И это дает результат.

setkeyv(genecoords, "chr") 
setkeyv(transcriptcoords, "chr") 
dt <- genecoords[transcriptcoords, allow.cartesian = TRUE] 
res <- dt[start <= start.1 + tol & end >= end.1 - tol, list(gene = gene[which.min(size)]), by = transcript] 

Понимая, что это не будет включать транскрипты, которые не удовлетворяют условию (g<-t в коде).