Применение пользовательской функции к data.table by row возвращает неверное количество значений

Я новичок в data.tables и у меня есть таблица, содержащая геномные координаты ДНК, например:

       chrom   pause strand coverage
    1:     1 3025794      +        1
    2:     1 3102057      +        2
    3:     1 3102058      +        2
    4:     1 3102078      +        1
    5:     1 3108840      -        1
    6:     1 3133041      +        1

Я написал пользовательскую функцию, которую хочу применить к каждой строке в моей таблице из примерно 2 миллионов строк, она использует mapToTranscripts GenomicFeatures для извлечения двух связанных значений в виде строки и новой координаты. Я хочу добавить их в мою таблицу в двух новых столбцах, например:

       chrom   pause strand coverage       transcriptID CDS
    1:     1 3025794      +        1 ENSMUST00000116652 196
    2:     1 3102057      +        2 ENSMUST00000116652  35
    3:     1 3102058      +        2 ENSMUST00000156816 888
    4:     1 3102078      +        1 ENSMUST00000156816 883
    5:     1 3108840      -        1 ENSMUST00000156816 882
    6:     1 3133041      +        1 ENSMUST00000156816 880

Функция следующая:

    get_feature <- function(dt){

      coordinate <- GRanges(dt$chrom, IRanges(dt$pause, width = 1), dt$strand) 
      hit <- mapToTranscripts(coordinate, cds_canonical, ignore.strand = FALSE) 
      tx_id <- tx_names[as.character(seqnames(hit))] 
      cds_coordinate <- sapply(ranges(hit), '[[', 1)

      if(length(tx_id) == 0 || length(cds_coordinate) == 0) {  
        out <- list('NaN', 0)
      } else {
        out <- list(tx_id, cds_coordinate)
      }

      return(out)
    } 

Затем я делаю:

    counts[, c("transcriptID", "CDS"):=get_feature(.SD), by = .I] 

И я получаю эту ошибку, указывающую, что функция возвращает два списка более короткой длины, чем исходная таблица, вместо одного нового элемента в строке:

Warning messages:
    1: In `[.data.table`(counts, , `:=`(c("transcriptID", "CDS"),  ... :
      Supplied 1112452 items to be assigned to 1886614 items of column 'transcriptID' (recycled leaving remainder of 774162 items).
    2: In `[.data.table`(counts, , `:=`(c("transcriptID", "CDS"),  ... :
      Supplied 1112452 items to be assigned to 1886614 items of column 'CDS' (recycled leaving remainder of 774162 items).

Я предполагал, что использование оператора . I применило бы функцию на основе строки и вернуло бы одно значение на строку. Я также убедился, что функция не возвращает пустые значения, используя оператор if.

Затем я попробовал эту ложную версию функции:

    get_feature <- function(dt) {

      return('I should be returned once for each row')

    }

И назвал это так:

    new.table <- counts[, get_feature(.SD), by = .I] 

Создает таблицу данных из 1 строки вместо одной исходной длины. Поэтому я пришел к выводу, что моя функция, или, как я ее называю, каким-то образом сворачивает элементы результирующего вектора. Что я делаю неправильно?

Обновление (с решением): как указано @StatLearner, в этом ответе объясняется, что, как объяснено в ?data.table, .I предназначен только для использования в j (как в DT[i,j,by=]). Следовательно, by=.I эквивалентно by=NULL и правильный синтаксис by=1:nrow(dt) для группировки по номеру строки и применения функции по строкам.

К сожалению, для моего конкретного случая это совершенно неэффективно, и я рассчитал время выполнения 20 секунд для 100 строк. Для моего набора данных из 36 миллионов строк, который занимает 3 месяца.

В моем случае мне пришлось отказаться и использовать mapToTranscripts Функция на всей таблице, как это, которая занимает пару секунд и, очевидно, была предназначена для использования.

    get_features <- function(dt){
      coordinate <- GRanges(dt$chrom, IRanges(dt$pause, width = 1), dt$strand) # define coordinate
      hits <- mapToTranscripts(coordinate, cds_canonical, ignore.strand = FALSE) # map it to a transcript
      tx_hit <- as.character(seqnames(hits)) # get transcript number
      tx_id <- tx_names[tx_hit] # get transcript name from translation table

      return(data.table('transcriptID'= tx_id, 
                       'CDS_coordinate' =  start(hits))
    }

     density <- counts[, get_features(.SD)]

Затем вернитесь к геному, используя mapFromTranscripts от GenomicFeatures пакет, чтобы я мог использовать data.tables присоединиться, чтобы получить информацию из исходной таблицы, что и было целью того, что я пытался сделать.

1 ответ

Решение

То, как я делаю это, когда мне нужно применить функцию для каждой строки в data.table, группирует ее по номеру строки:

counts[, get_feature(.SD), by = 1:nrow(counts)]

Как объясняется в этом ответе,.I не предназначен для использования в by так как он должен возвращать последовательность индексов строк, которая создается путем группировки. Причина по которой by = .I не выдает ошибку в том, что data.table создает объект .I равняется NULL в пространстве имен data.table, следовательно by = .I эквивалентно by = NULL,

Обратите внимание, что с помощью by=1:nrow(dt) группирует по номеру строки и позволяет вашей функции обращаться только к одной строке из data.table:

require(data.table)
counts <- data.table(chrom = sample.int(10, size = 100, replace = TRUE),
                     pause = sample((3 * 10^6):(3.2 * 10^6), size = 100), 
                     strand = sample(c('-','+'), size = 100, replace = TRUE),
                     coverage = sample.int(3, size = 100, replace = TRUE))

get_feature <- function(dt){
    coordinate <- data.frame(dt$chrom, dt$pause, dt$strand)
    rowNum <- nrow(coordinate)
    return(list(text = 'Number of rows in dt', rowNum = rowNum))  
}

counts[, get_feature(.SD), by = 1:nrow(counts)]

создаст таблицу данных с тем же количеством строк, что и в counts, но coordinate будет содержать только одну строку из counts

   nrow                 text rowNum
1:    1 Number of rows in dt      1
2:    2 Number of rows in dt      1
3:    3 Number of rows in dt      1
4:    4 Number of rows in dt      1
5:    5 Number of rows in dt      1

в то время как by = NULL предоставит всю data.table для функции:

counts[, get_feature(.SD), by = NULL]

                   text rowNum
1: Number of rows in dt    100

который предназначен для by работать.

Другие вопросы по тегам