Матрица дизайна для генотипов в R

Question

Я ищу эффективный способ создания «параметризованной» матрицы дизайна для генотипов в R. У меня есть большой файл (около 3 ГБ) с животными и их генотипами. Пример данных выглядит так:

      snp id a1 a2 code
snp1 an1 A A 0
snp1 an2 A B 1
snp1 an3 B B -1
snp2 an1 A B 1
snp2 an2 A A 0
snp2 an3 B B -1

snp - это имя snp (каждое животное имеет каждый snp), id - идентификатор животного (каждое животное имеет уникальный идентификатор), a1 - аллель 1, a2 - аллель 2, code обозначает генотип на основе аллелей, если у животного есть два A, это код 0, если у животного AB, это код -1, а если BB - код 1.

Теперь мне нужно создать на основе этой матрицы дизайна, которая в строке будет иметь животное (столбец id в данных), а в столбцах SNP (столбец snp в данных) и в «ячейке» (на пересечении строки и столбца) Мне нужно значение из столбца кода. В конце концов, это должно выглядеть так:

      an1 0 1
an2 1 0
an3 -1 -1

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

Answer 1

Обычно пакет data.table довольно эффективен в таких случаях. Пример ниже:

      library(data.table)
#> Warning: package 'data.table' was built under R version 4.1.1

df <- fread(text = "snp id a1 a2 code
snp1 an1 A A 0
snp1 an2 A B 1
snp1 an3 B B -1
snp2 an1 A B 1
snp2 an2 A A 0
snp2 an3 B B -1")

dcast(df, id ~ snp, value.var = "code")
#>     id snp1 snp2
#> 1: an1    0    1
#> 2: an2    1    0
#> 3: an3   -1   -1

^{Создано REPEX (v2.0.1)}

Если вам нужен вывод в виде матрицы, вы можете использовать:

      cast <- dcast(df, id ~ snp, value.var = "code")
mat <- as.matrix(cast[, -"id"])
rownames(mat) <- cast$id
mat
#>     snp1 snp2
#> an1    0    1
#> an2    1    0
#> an3   -1   -1

Для файла ~ 3Gb вы можете ожидать, что это будет работать около 10 секунд:

      library(data.table)
#> Warning: package 'data.table' was built under R version 4.1.1

# Setting up larger data
df <- expand.grid(
  snp = paste0("snp", 1:10000),
  id  = paste0("an", 1:10000)
)
df$a1 <- sample(c("A", "B"), nrow(df), replace = TRUE)
df$a2 <- sample(c("A", "B"), nrow(df), replace = TRUE)
df$code <- with(df, dplyr::case_when(
  a1 == "A" & a2 == "A" ~ 0,
  a1 == "B" & a2 == "B" ~ -1,
  TRUE ~ 1
))
setDT(df)

# How big is this data?
format(object.size(df), "Gb")
#> [1] "3 Gb"

# How fast does the function run?
bench::mark(
  dcast(df, id ~ snp, value.var = "code")
)
#> Warning: Some expressions had a GC in every iteration; so filtering is disabled.
#> # A tibble: 1 x 6
#>   expression                                   min   median `itr/sec` mem_alloc
#>   <bch:expr>                              <bch:tm> <bch:tm>     <dbl> <bch:byt>
#> 1 dcast(df, id ~ snp, value.var = "code")    9.32s    9.32s     0.107    6.71GB
#> # ... with 1 more variable: gc/sec <dbl>

^{Создано 2021-10-13 пакетом2021-10-13 пакетом REPEX (v2.0.1)}