Ошибка в моем коде и улучшить функцию - разделить по запросу - PullRequest
0 голосов
/ 14 февраля 2019

У меня есть этот пример фрейма данных

nucleotide  start  end    strand  block_id  query  pid       
AE002161.1  5537   6724   1       1         0      AAF73616.1
AE002161.1  6714   7727   1       1         0      AAF37902.1
AE002161.1  7687   10839  -1      1         1      AAF37903.1
AE002161.1  10826  13900  -1      1         0      AAF37904.1
AE002161.1  13887  15596  1       1         0      AAF37905.1
AE002161.1  18606  19487  -1      2         0      AAF37910.1
AE002161.1  19822  19998  -1      2         0      AAF37911.1
AE002161.1  19982  21625  1       2         1      AAF37912.1
AE002161.1  21728  22996  1       2         0      AAF37913.1
AE002161.1  23108  25063  1       2         0      AAF37914.1
AE002161.1  36276  36575  -1      3         0      AAF37924.1
AE002161.1  36680  38116  -1      3         0      AAF37925.1
AE002161.1  38120  39928  -1      3         1      AAF37926.1
AE002161.1  40478  41497  1       3         0      AAF37927.1
AE002161.1  41864  42256  1       3         0      AAF37928.1
AE002161.1  45880  46554  1       4         0      AAF37933.1
AE002161.1  46556  47884  1       4         0      AAF37934.1
AE002161.1  47902  48408  1       4         1      AAF37935.1
AE002161.1  48412  49254  1       4         1      AAF37936.1
AE002161.1  49264  50379  1       4         0      AAF73618.1
AE002161.1  50395  51903  1       4         0      AAF73619.1

и эта функция

library(tidyverse)

splitq <- function(data){
  a <- data %>%
    mutate(., block_id = group_indices(., nucleotide, block_id) ) %>%
    group_by(nucleotide, block_id) %>%
    mutate(old=cumsum(query)) %>%
    mutate( query = ifelse( old > 1 , 0,  query ) ) %>%
    ungroup()

  a_max <- max(a$block_id)

  b <- data %>%
    arrange( desc(row_number() ) ) %>%
    mutate(., block_id = group_indices(., nucleotide, block_id) + a_max ) %>%
    group_by(nucleotide, block_id) %>%
    mutate(old=cumsum(query)) %>%
    mutate( query = ifelse( old > 1 , 0,  query ) ) %>%
    ungroup() %>%
    bind_rows(a) %>%
    select(-old)
}

, когда я запускаю эту функцию, я получаю этот результат

nucleotide  start  end    strand  block_id  query  pid         type 
AE002161.1  50395  51903  1       8         0      AAF73619.1  CDS  
AE002161.1  49264  50379  1       8         0      AAF73618.1  CDS  
AE002161.1  48412  49254  1       8         1      AAF37936.1  CDS  
AE002161.1  47902  48408  1       8         0      AAF37935.1  CDS  
AE002161.1  46556  47884  1       8         0      AAF37934.1  CDS  
AE002161.1  45880  46554  1       8         0      AAF37933.1  CDS  
AE002161.1  41864  42256  1       7         0      AAF37928.1  CDS  
AE002161.1  40478  41497  1       7         0      AAF37927.1  CDS  
AE002161.1  38120  39928  -1      7         1      AAF37926.1  CDS  
AE002161.1  36680  38116  -1      7         0      AAF37925.1  CDS  
AE002161.1  36276  36575  -1      7         0      AAF37924.1  CDS  
AE002161.1  23108  25063  1       6         0      AAF37914.1  CDS  
AE002161.1  21728  22996  1       6         0      AAF37913.1  CDS  
AE002161.1  19982  21625  1       6         1      AAF37912.1  CDS  
AE002161.1  19822  19998  -1      6         0      AAF37911.1  CDS  
AE002161.1  18606  19487  -1      6         0      AAF37910.1  CDS  
AE002161.1  13887  15596  1       5         0      AAF37905.1  CDS  
AE002161.1  10826  13900  -1      5         0      AAF37904.1  CDS  
AE002161.1  7687   10839  -1      5         1      AAF37903.1  CDS  
AE002161.1  6714   7727   1       5         0      AAF37902.1  CDS  
AE002161.1  5537   6724   1       5         0      AAF73616.1  CDS  
AE002161.1  5537   6724   1       1         0      AAF73616.1  CDS  
AE002161.1  6714   7727   1       1         0      AAF37902.1  CDS  
AE002161.1  7687   10839  -1      1         1      AAF37903.1  CDS  
AE002161.1  10826  13900  -1      1         0      AAF37904.1  CDS  
AE002161.1  13887  15596  1       1         0      AAF37905.1  CDS  
AE002161.1  18606  19487  -1      2         0      AAF37910.1  CDS  
AE002161.1  19822  19998  -1      2         0      AAF37911.1  CDS  
AE002161.1  19982  21625  1       2         1      AAF37912.1  CDS  
AE002161.1  21728  22996  1       2         0      AAF37913.1  CDS  
AE002161.1  23108  25063  1       2         0      AAF37914.1  CDS  
AE002161.1  36276  36575  -1      3         0      AAF37924.1  CDS  
AE002161.1  36680  38116  -1      3         0      AAF37925.1  CDS  
AE002161.1  38120  39928  -1      3         1      AAF37926.1  CDS  
AE002161.1  40478  41497  1       3         0      AAF37927.1  CDS  
AE002161.1  41864  42256  1       3         0      AAF37928.1  CDS  
AE002161.1  45880  46554  1       4         0      AAF37933.1  CDS  
AE002161.1  46556  47884  1       4         0      AAF37934.1  CDS  
AE002161.1  47902  48408  1       4         1      AAF37935.1  CDS  
AE002161.1  48412  49254  1       4         0      AAF37936.1  CDS  
AE002161.1  49264  50379  1       4         0      AAF73618.1  CDS  
AE002161.1  50395  51903  1       4         0      AAF73619.1  CDS  

РЕДАКТИРОВАТЬ: КакиеКажется, не приятно, так как он генерирует некоторую избыточность, он должен создать 5 блоков, а не 8.

Я просто хочу разделить на query == 1.Таким образом, для каждого запроса у меня должно быть n строк выше и n ниже (те же строки в том же порядке).Эта операция должна выполняться с помощью block_id.

Когда два соседа query == 1 будут рядом, как

AE002161.1  45880  46554  1       4         0      AAF37933.1
AE002161.1  46556  47884  1       4         0      AAF37934.1
AE002161.1  47902  48408  1       4         1      AAF37935.1
AE002161.1  48412  49254  1       4         1      AAF37936.1
AE002161.1  49264  50379  1       4         0      AAF73618.1
AE002161.1  50395  51903  1       4         0      AAF73619.1

Должно возвращаться

AE002161.1  45880  46554  1       4         0      AAF37933.1
AE002161.1  46556  47884  1       4         0      AAF37934.1
AE002161.1  47902  48408  1       4         1      AAF37935.1
AE002161.1  48412  49254  1       4         0      AAF37936.1
AE002161.1  49264  50379  1       4         0      AAF73618.1
AE002161.1  50395  51903  1       4         0      AAF73619.1
AE002161.1  45880  46554  1       5         0      AAF37933.1
AE002161.1  46556  47884  1       5         0      AAF37934.1
AE002161.1  47902  48408  1       5         0      AAF37935.1
AE002161.1  48412  49254  1       5         1      AAF37936.1
AE002161.1  49264  50379  1       5         0      AAF73618.1
AE002161.1  50395  51903  1       5         0      AAF73619.1

Это означает, что яне волнует, если все block_id изменяются, так как он уникален для каждого блока (не повторять нигде в выходных данных).

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

Но когда я запускаю это для файла 591MB с 2070926 строками, где 332236 - query == 1, где 330409 из них различны, я получаю некоторые ошибки.

Сообщение об ошибке не генерируется, но я пропускаю некоторыевопросы.

Кто-нибудь знает, что происходит?

Заранее спасибо

1 Ответ

0 голосов
/ 14 февраля 2019

Фиксированная функция

splitq <- function(data){
  a <- data %>% filter(query == 1) %>% mutate( old = block_id, new = row_number()) %>% select(pid, new, old)
  b <- data %>% 
    left_join(a, by = c("block_id" = "old")) %>%
    group_by(new) %>%
    mutate( query = ifelse( pid.x == pid.y, 1, 0), block_id = new ) %>%
    arrange(nucleotide, block_id, start, end) %>%
    select(-pid.y, -new) %>%
    rename(pid=pid.x) %>%
    ungroup()
}
...