R - Вложенные для циклов и медленной производительности

Question

Я пытаюсь реализовать функцию для получения значений из одной таблицы на основе другой. Фактические данные имеют> 50 000 наблюдений, поэтому реализация этого вложенного цикла неэффективна. В последние несколько дней я пытаюсь найти SO, чтобы найти что-то, что работает, но не удалось. Мои данные не имеют особого порядка (отдельные лица, сегменты и т. Д.), Поэтому он должен работать, даже если все не в порядке.

Вот игрушечные примеры моих данных для работы с:

region_map <- data.frame(Start = c(721290, 1688193), End= c(1688192, 2926555)) 
individual <- c("Ind1","Ind2","Ind3","Ind4") 
segment <- data.frame(SampleID = c("Ind1","Ind1","Ind2","Ind2","Ind3","Ind3","Ind4","Ind4","Ind4"), 
         Start = c(721290, 1688194, 721290, 1688200, 721290, 2926600, 721290, 1688193, 690), 
         End = c(1688192, 2926555,1688190, 2900000, 2926555, 3000000, 1500000, 2005000, 500000), 
         State = c(1,2,2,5,4,2,2,6,5)) 

И вот упрощенный пример того, что я пытаюсь сделать:

Generate.FullSegmentList <- function(segments, individuals, regionmap){ 
    FullSegments <- data.frame() 
    for(region in 1:nrow(regionmap)){ 

      for(ind in individuals){ 
       # If there is not a segment within that region for that individual 
       if(nrow(
        segments[segments$start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind , ] 
       ) == 0){ 
        Temp <- data.frame(SampleID = ind, 
             Start = regionmap$Start[region], 
             End = regionmap$End[region], 
             State = 3 
        ) 
       } 
       # If there is a segment within that region for that individual 
       if(nrow(
        segments[segments$Start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind , ] 
       ) == 1){ 
        Temp <- data.frame(SampleID = segments$SampleID, 
             Start = regionmap$Start[region], 
             End = regionmap$End[region], 
             State = segments$State[segments$Start >= regionmap$Start[region] & 
                    segments$SampleID == ind ] 
        ) 
       } 
       FullSegments <- list(FullSegments, Temp)    
      } 
    } 
    FullSegments 
} 

В словах, мне нужно посмотреть в каждой области (~ 53 000) и присвойте значение (State, если оно не существует, укажите значение 3) в область для каждого individual, а затем создайте новый data.frame с каждой областью для каждого человека. Чтобы сделать это, я просматриваю регионы, а затем люди, находя segment (их насчитывается ~ 25 000), которые пересекаются с областью, а затем присоединяют ее к таблице.

Вот что вывод из приведенных выше данных игрушечных даст:

SampleID  Start  End  State 
Ind1   721290  1688192  1 
Ind1   1688193  2926555  2 
Ind2   721290  1688192  2 
Ind2   1688193  2926555  5 
Ind3   721290  1688192  4 
Ind3   1688193  2926555  4 
Ind4   721290  1688192  2 
Ind4   1688193  2926555  6 

Эта функция как есть работает точно, как мне это нужно, чтобы, за исключением того, что это займет очень много времени для запуска (с помощью system.time, я понял, что для запуска потребуется более 3 месяцев). Я знаю, что должен быть лучший способ сделать это. Я попытался реализовать функции приложения, и я видел в некоторых других вопросах использование списков вместо data.frame. Я также видел, что для упрощения этого существуют параметры data.table и plyr. Я пробовал их, но не добился успеха в работе с вложенным циклом с операторами if.

Я был бы признателен за объяснение любых полученных ответов, так как это первый раз, когда я написал все это сложное.

Вопросы Я думаю актуальны:

• R Grouping functions
• Post discussing problems with data.frame(...) approach

многие другие вопросы на вложенные в петлях вовлекают делать расчеты, которые хорошо работают для этого функцию применения (например apply(df, 1, function(x){ mean(x) }) , но я не смог принять это значение для отображения значений из data.frame в data.frame.

Answer 1

Пакет Bioconductor IRanges работает на «целое число в диапазоне», как область и сегмент начала и конца координат. Установите пакет с

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

загрузить его и создать представление диапазонов интереса

library(IRanges) 
r <- with(region_map, IRanges(Start, End)) 
s <- with(segments, IRanges(Start, End)) 

В результате до сих пор

> r 
IRanges object with 2 ranges and 0 metadata columns: 
      start  end  width 
     <integer> <integer> <integer> 
    [1] 721290 1688192 966903 
    [2] 1688193 2926555 1238363 
> s 
IRanges object with 9 ranges and 0 metadata columns: 
      start  end  width 
     <integer> <integer> <integer> 
    [1] 721290 1688193 966904 
    [2] 1688194 2926555 1238362 
    [3] 721290 1688190 966901 
    [4] 1688200 2900000 1211801 
    [5] 721290 2926555 2205266 
    [6] 2926600 3000000  73401 
    [7] 721290 1500000 778711 
    [8] 1688193 2005000 316808 
    [9]  690 500000 499311 

Вы заинтересованы в поиске совпадений между сегменты «запроса» и «субъект» region_map

olaps <- findOverlaps(s, r) 

давая

> olaps 
Hits object with 9 hits and 0 metadata columns: 
     queryHits subjectHits 
     <integer> <integer> 
    [1]   1   1 
    [2]   1   2 
    [3]   2   2 
    [4]   3   1 
    [5]   4   2 
    [6]   5   1 
    [7]   5   2 
    [8]   7   1 
    [9]   8   2 
    ------- 
    queryLength: 9/subjectLength: 2 

Это хорошо масштабироваться до миллионов наложений.

Вы сказали, что вас интересует состояние всех людей во всех регионах, и из вашего кода это похоже на то, что человек не в регионе имеет состояние 3.Я создал матрицу со всего государства 3

state <- matrix(3, nrow(region_map), length(individual), 
       dimnames=list(NULL, individual)) 

затем созданный индекс в два столбца в матрице на основе наложений мы нашли

idx <- matrix(c(subjectHits(olaps), 
       match(segments$SampleID[queryHits(olaps)], individual)), 
       ncol=2) 

и использовал матрицу индекса для обновления состояния

state[idx] <- segments$State[queryHits(olaps)] 

Это фактически суммирует ваш желаемый результат - состояние в каждом регионе x индивидуальная комбинация. Одна из возможных проблем заключается в том, что два сегмента одного и того же лица перекрывают один регион, а сегменты имеют разное состояние; будет назначено только одно состояние.

> state 
    Ind1 Ind2 Ind3 Ind4 
[1,] 1 2 4 2 
[2,] 2 5 4 6 

Cast это как data.frame с, например,

data.frame(SampleID=colnames(state)[col(state)], 
      Start=region_map[row(state), "Start"], 
      End=region_map[row(state), "End"], 
      State=as.vector(state)) 

Answer 2

Я не думаю, что вам нужно что-либо «этот комплекс». Вы можете сделать все, что вам нужно, с помощью нескольких подключений. В этом случае я буду использовать data.table.

Вы попросили объяснить какой-либо ответ, однако для этого я не могу сделать лучше, чем указать вам в направлении data.table homepage. Будет важно понять, что делают команды set* и := и как работает «обновление по ссылке».

Задайте свои данные data.table s.

library(data.table) 

dt_individual <- data.table(SampleID = individual) 
dt_region <- data.table(region_map) 
dt_segment <- data.table(segment) 

Просто присоединиться все это вместе

## Change some column names of `dt_segment` so we can identify them after the joins 
setnames(dt_segment, c("Start", "End"), c("seg_Start", "seg_End")) 

## create a 'key_col' to join all the individuals to the regions 
dt_join <- dt_individual[, key_col := 1][ dt_region[, key_col := 1], on="key_col", allow.cartesian=T][, key_col := NULL] 
# SampleID Start  End 
# 1:  Ind1 721290 1688192 
# 2:  Ind2 721290 1688192 
# 3:  Ind3 721290 1688192 
# 4:  Ind4 721290 1688192 
# 5:  Ind1 1688193 2926555 
# 6:  Ind2 1688193 2926555 
# 7:  Ind3 1688193 2926555 
# 8:  Ind4 1688193 2926555 

Теперь используйте функцию foverlaps, чтобы найти пересекающиеся регионы

setkey(dt_join, SampleID, Start, End) 
setkey(dt_segment, SampleID, seg_Start, seg_End) 

foverlaps(dt_join, 
      dt_segment, 
      type="any") 

# SampleID seg_Start seg_End State Start  End 
# 1:  Ind1 721290 1688192  1 721290 1688192 
# 2:  Ind1 1688194 2926555  2 1688193 2926555 
# 3:  Ind2 721290 1688190  2 721290 1688192 
# 4:  Ind2 1688200 2900000  5 1688193 2926555 
# 5:  Ind3 721290 2926555  4 721290 1688192 
# 6:  Ind3 721290 2926555  4 1688193 2926555 
# 7:  Ind4 721290 1500000  2 721290 1688192 
# 8:  Ind4 1688193 2005000  6 1688193 2926555 

---

Чтобы посмотреть все данные (то есть и те, которые подпадают в регионах и тех, которые этого не делают), вы можете сделать соединение cartesian, а затем assi дп значения для тех, кто в этом регионе, и те, за ее пределами, как вы хотели

dt_join[dt_segment, on="SampleID", nomatch=0, allow.cartesian=T] 

Answer 3

У вас есть много строк в коде, которые читают nrow(some-subset-of-your-data). Вы бы быстро увеличили производительность, если бы переключили их на sum(the-conditions). Например:

Turn:

nrow(segments[segments$start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind , ]) == 0 

В

sum(segments$start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind) == 0 

Таким образом, R не хранит subsetted кадр данных в памяти каждый раз.

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

isEmpty <- sum(segments$start >= regionmap$Start[region] & 
            segments$End <= regionmap$End[region] & 
            segments$SampleID == ind) == 0 

if(isEmpty){ 
### do something 
} else if(!isEmpty) { 
### do something else 
}