2010-10-12 19 views
12

Mam dwa data.frames każdy z trzema kolumnami: chrom, start & stop, nazwijmy je zakresami A i zakresamiB. Dla każdego rzędu zakresów A szukam, który (jeśli którykolwiek) wiersz w zakresach B zawiera w całości rząd A rzędów - przez co rozumiem rangesAChrom == rangesBChrom, rangesAStart >= rangesBStart and rangesAStop <= rangesBStop.Znalezienie nakładki w zakresach z R

W tej chwili robię następujące, które po prostu nie lubię bardzo. Zauważ, że zapętlałem się nad rzędami zakresów A z innych powodów, ale żaden z tych powodów nie jest prawdopodobnie wielkim problemem, po prostu sprawia, że ​​rzeczy stają się bardziej czytelne, biorąc pod uwagę to konkretne rozwiązanie.

rangesA:

chrom start stop 
5  100  105 
1  200  250 
9  275  300 

rangesB:

chrom start stop 
    1  200  265 
    5  99  106 
    9  275  290 

dla każdego wiersza w rangesA:

matches <- which((rangesB[,'chrom'] == rangesA[row,'chrom']) && 
       (rangesB[,'start'] <= rangesA[row, 'start']) && 
       (rangesB[,'stop'] >= rangesA[row, 'stop'])) 

I postać nie musi być lepsze (i lepiej, mam na myśli szybszy nad dużymi instancjami zakresów A i zakresów B) sposób, aby to zrobić, niż zapętlenie tego konstruktu. Jakieś pomysły?

Odpowiedz

13

Byłoby to o wiele łatwiejsze/szybsze, jeśli najpierw można scalić oba obiekty.

ranges <- merge(rangesA,rangesB,by="chrom",suffixes=c("A","B")) 
ranges[with(ranges, startB <= startA & stopB >= stopA),] 
# chrom startA stopA startB stopB 
#1  1 200 250 200 265 
#2  5 100 105  99 106 
2

Na swoim przykład dane:

rangesA <- data.frame(
    chrom = c(5, 1, 9), 
    start = c(100, 200, 275), 
    stop = c(105, 250, 300) 
) 
rangesB <- data.frame(
    chrom = c(1, 5, 9), 
    start = c(200, 99, 275), 
    stop = c(265, 106, 290) 
) 

to zrobi to sapply tak, że każda kolumna ma jeden wiersz rangesA a każdy wiersz odpowiada wiersz rangesB:

> sapply(rangesA$stop, '>=', rangesB$start) & sapply(rangesA$start, '<=', rangesB$stop) 
     [,1] [,2] [,3] 
[1,] FALSE TRUE FALSE 
[2,] TRUE FALSE FALSE 
[3,] FALSE FALSE TRUE 
18

Użyj pakietów IRanges/GenomicRanges z Bioconductor, który jest przeznaczony do rozwiązywania tych konkretnych problemów (i masywnie)

source("http://bioconductor.org/biocLite.R") 
biocLite("IRanges") 

Istnieje kilka odpowiednie pojemniki do zakresów na różnych chromosomach, jeden jest RangesList

library(IRanges) 
rangesA <- split(IRanges(rangesA$start, rangesA$stop), rangesA$chrom) 
rangesB <- split(IRanges(rangesB$start, rangesB$stop), rangesB$chrom) 
#which rangesB wholly contain at least one rangesA? 
ov <- countOverlaps(rangesB, rangesA, type="within")>0 
+1

Dobra wskazówka na IRanges, zapomniałem o tym. Nie skończyłem na tym, ponieważ nie pasowało to do mojej sytuacji z różnych powodów, ale nadal dobrze wiedzieć. Mój zabawkowy przykład pomijał kilka kluczowych bitów (moja własna wina), które utrudniały mi pracę z IRLingiem, a rozwiązanie merge() zapewniało ogromne przyspieszenie. Ponadto, podczas gdy skaluje się masowo, widzimy także wiele bardzo małych przypadków i zakładam, że obciążenie S4 spowolniło w tych przypadkach. – geoffjentry

+1

czy istnieje również możliwość policzenia również częściowo zachodzenia na siebie? – Cina

2

RangesA i RangesB są wyraźnie składnia łóżku, można to zrobić na zewnątrz R w linii komend z BEDtools, niezwykle szybki i elastyczne z kilkoma innymi opcjami do pracy z interwałami genomowymi. Następnie umieścić wyniki z powrotem do R.

https://code.google.com/p/bedtools/

9

Pakiet data.table posiada funkcję foverlaps() która jest zdolna do łączenia przez interwał zakresy od v1.9.4:

require(data.table) 
setDT(rangesA) 
setDT(rangesB) 

setkey(rangesB) 
foverlaps(rangesA, rangesB, type="within", nomatch=0L) 
# chrom start stop i.start i.stop 
# 1:  5 99 106  100 105 
# 2:  1 200 265  200 250 
  • setDT() dokonana konwersja data.frame to data.table przez odniesienie

  • setkey() sortuje plik data.table według podanych kolumn (w tym przypadku wszystkich kolumn, ponieważ nie podano żadnych) i oznacza te kolumny jako posortowane, które wykorzystamy później, aby wykonać dołączenie.

  • foverlaps() efektywnie nakłada się na siebie. Zobacz this answer, aby uzyskać szczegółowe wyjaśnienie i porównanie z innymi podejściami.

1

Dodaję rozwiązanie dplyr.

library(dplyr) 
inner_join(rangesA, rangesB, by="chrom") %>% 
    filter(start.y < start.x | stop.y > stop.x) 

wyjściowa:

chrom start.x stop.x start.y stop.y 
1  5  100 105  99 106 
2  1  200 250  200 265 
+0

Wyobraź sobie sytuację, w której zakresy mają 20 tys. Wierszy, a zakresy B - 300 tys. Wierszy -> ogromne łączenie -> niemożliwe do zmieszczenia w pamięci RAM. Rozwiązania z IRanges są lepsze –

+0

IRanges nie działał na PO. Patrz wyżej. – Joe

Powiązane problemy