R.基因型设计矩阵

Question

我正在寻找一种有效的方法来创建一个"参数化";基因型设计矩阵。我有一个大文件(大约3gb)，里面有动物和它们的基因型。示例数据如下所示:

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是动物的id(每个动物都有唯一的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

通常是数据。表包很性能在这些类型的情况下。在下面的例子:

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

^{由reprex包(v2.0.1)在2018-10-13上创建}

如果你需要一个矩阵输出，你可以使用:

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>

^{由reprex包(v2.0.1)在2018-10-13上创建}