具有 200Gb 的 linux 集群中的 64x R：plyr 包不足以加入表

Posted 2023-04-18

【中文标题】具有 200Gb 的 linux 集群中的 64x R：plyr 包不足以加入表

【英文标题】：64x R in linux cluster with 200Gb: plyr package not enough to join tables

【发布时间】：2017-10-02 11:11:58

【问题描述】：

我在尝试使用 R 中的 plyr 包连接多个 data.frame 时遇到问题。

也许这不是大型数据集连接表的最佳工具/语言？如果有人可以提出替代方案，将不胜感激（例如替代 R 包、unix，或者可能是 MapReduce/Hadoop 或 python [注意：我没有 python 或 mapreduce/hadoop 经验]。

我在一个有高达 200Gb 可用内存的 Linux 集群上运行 64 位版本的 R。

我有 22 个样本/data.frame，每个都有 2 列（"ensembl_gene_id" 和 "FPKM" 值），需要合并约 60,000 行一起。因此，所需的 data.frame 将包含一列 ensembl_gene_id 和 22 列 FPKM 值。

在某些样本中缺少基因，因此我不能简单地对数据框进行排序并逐行连接它们。所以，我的想法是以迭代方式使用 plyr 包的 join() 函数；但是，由于需要内存，我现在质疑这是否是解决此问题的理想方法。

我提供了 3 个 data.frames 和我尝试过的代码。 df1, df2 和 df3（我没有足够的信誉发布两个以上的链接，所以我将把它放在下面的评论中）

library(plyr)

df1 <- read.csv("HAP1_dBet6_1h_1.csv")
df2 <- read.csv("HAP1_dBet6_1h_2.csv")
df3 <- read.csv("HAP1_dBet6_6h_1.csv")

HAP1 <- join(df1, df2, by="ensembl_gene_name")
Error: cannot allocate vector of size 1.5 Gb

关于分配向量失败的要点是强调两点：

1) R 给出了非常令人困惑的错误。如果我打开另一个终端窗口并使用top 命令检查 R 实际使用的是什么，我可以看到我超过了我可以访问的默认 8Gb！因此，我使用salloc 请求 80 Gb 并尝试了下面的代码。 2) 第二点是演示当我以迭代方式执行 plyr join() 命令时内存需求如何增加。

HAP2 <- join(HAP1, df3, by="ensembl_gene_name")
Error: cannot allocate vector of size 34282.1 Gb

理论上，即使我在进行更多连接时可以访问 Tb 的内存，我认为这个向量会增加，从而使这个问题变得难以处理。例如，

HAP3<-join(HAP2, df4, by="ensembl_gene_id")
HAP4<-join(HAP3, df5, by="ensembl_gene_id")
etc.

谁能建议我接下来应该尝试什么？

我注意到每个 data.frame 中有很多没有 ensembl id 的基因（~20,000 NA's）。这些基本上没有意义，所以如果我在加入之前删除它们可能会有所帮助。我想到的另一件事是在plyr 连接之间使用垃圾清理功能gc()。但是，我仍然怀疑这是否会有所帮助。

##########################工作更新

由于友好的用户解决了这个问题，我只在最初的问题中添加了最少的代码，我想我会编写我的工作流程来帮助其他人

library(plyr)
library(dplyr)
library(data.table)

# Load CSV files
df1 <- read.csv("HAP1_dBet6_1h_1.csv")
df2 <- read.csv("HAP1_dBet6_1h_2.csv")
df3 <- read.csv("HAP1_dBet6_6h_1.csv")
df4 <- read.csv("HAP1_dBet6_6h_2.csv")
...
df28 <- read.csv("HAP1_DMSO_6h_2.csv")

# Select only gene column and FPKM columns 
df1 <- select(df1,ensembl_gene_name,FPKM)
...

# Rename FPKM column to something meaningful (i.e. sample name) and omit NA's
df1 <- plyr::rename(df1,c("FPKM"="HAP1_dBet6_1h_1.csv"))
df1 <- na.omit(df1)
...

# Select unique gene name values and summarize FPKM counts
df1 <- df1 %>% group_by(ensembl_gene_name) %>% summarise(sum(K562_dBet6_1h_1))
...

df1 <- data.table(df1) # Turn data.frame into a data.table
setkey(df1, "ensembl_gene_name") # setkey() sorts a data.table and marks it as sorted
...

# Join data.tables
HAP1 <- merge(df1, df2, "ensembl_gene_name")
HAP2 <- merge(HAP1, df3, "ensembl_gene_name")
HAP3 <- merge(HAP2, df4, "ensembl_gene_name")
...
HAP27 <- merge(HAP26, df28, "ensembl_gene_name")

谢谢！

【问题讨论】：

df3 链接gist.github.com/moldach/6d238adcad538269f4fe658f81bf4865

我想知道您是否可以尝试 dplyr 而不是 plyr。它应该更快等。但除此之外，我会选择 data.table （尽管它需要习惯不同/新的语法）。

今晚晚些时候我会试一试，但到目前为止，我的直觉是，连 dplyr 都倒下了。

@MatthewJ.Oldach Matt，有什么对你有用的吗？

我还为所有三个 data.frames 尝试了 plyr packages join_all() 函数，但没有运气。下一个去试试dplyr和data.table

【参考方案1】：

您绝对应该考虑从 data.frame 中删除 NAs。您可能会在右侧 data.frame 中复制每个 NA * NA 的数量。查看nrow(HAP1) 或str(HAP1) 的输出。它是否包含比df1 或df2 更多的行？一个 60,000 x 22 的数据框并没有那么大，但如果您的数据框在每次连接后扩展，它就会呈指数级增长。

【讨论】：

是的，我注意到了 NA 和一些重复的基因名称。所以我删除它们：df1 <- read.csv("HAP1_dBet6_1h_1.csv")df1<-unique(df1)df1<-na.omit(df1) df2 也是如此。当我nrow(HAP1) 时，行数仍然比 df1 或 df2 多得多。如果我执行以下操作：HAP1 <- na.omit(HAP1) 并再次执行nrow(HAP1)，则比以前少，但仍比 df1 或 df2 多。不知道这个结果如何？

df1 <- unique(df1) 可能不会返回您所期望的（只有唯一 ID 的数据框）。查看子集后的nrow(df1) 是否与unique(df1$ensembl_gene_id) 的值相同。要获得具有唯一基因 ID 的数据框，您需要执行类似library(dplyr) uniquedf1 <- df1 %>% group_by(ensembl_gene_id) %>% summarise(sum(FPKM)) 之类的操作。 df2 也是如此。然后加入。

请注意，dplyr 也有 dplyr::distinct(ensemble_gene_id)，它根据基因 id 给出不同的行。我在这里发表评论，因为我最近才知道 dplyr 的 distinct 可以实现这一点。一定要添加类似 keep.others = T 的东西（检查 ?dplyr::distinct 的确切名称）以确保其他东西仍然存在。

@PaulLemmens 感谢您的注意。在这种情况下，我并没有真正考虑过distinct，但它是一个有用的动词。 OP 数据的问题是多个基因 ID 可能具有不同的有意义的 FPKM 值。在这种情况下，FPKM 值需要某种聚合。

非常感谢@CPak uniquedf1 <- df1 %>% group_by(ensembl_gene_id) %>% summarise(sum(FPKM)) 工作。这完全是由于unique() 没有按照我认为的那样做的错误（如@PaulLemmens 在下面提到的那样，在三个data.frame 中的每一个中留下重复的键值）。这完全有道理，因为人们期望行号不超过 df1 或 df2。随着您加入更多 data.frame，行数可能会减少（或保持不变），永远不会增加。

【参考方案2】：

你能测试一下这个data.table 解决方案吗：

library(data.table)
df1 <- fread("HAP1_dBet6_1h_1.csv")
df2 <- fread("HAP1_dBet6_1h_2.csv")
setkey(df1, "ensembl_gene_name")
setkey(df2, "ensembl_gene_name")
HAP1 <- merge(df1, df2, "ensembl_gene_name")

在处理较大的数据集时，我更喜欢带有setkey 选项的data.table 对象。

编辑：

如果fread 正在将您的基因名称转换为其他名称，那么您可以使用它：

library(data.table)
df1 <- read.csv("HAP1_dBet6_1h_1.csv")
df2 <- read.csv("HAP1_dBet6_1h_2.csv")
setDT(df1)
setDT(df2)
setkey(df1, "ensembl_gene_name")
setkey(df2, "ensembl_gene_name")
HAP1 <- merge(df1, df2, "ensembl_gene_name")

【讨论】：

当我输入上面的代码并通过write.csv(HAP1, "HAP1.csv") 保存 HAP1 并打开查看它时，一些基因名称已转换为日期。我可以通过head(HAP1) 确认这一点，然后查看 df1 和 df2 以确保它们之前不存在。

@MatthewJ.Oldach 真的吗？并且您的 csv 文件中的基因没有转换为日期？因为这是非常常见的excel问题：genomebiology.biomedcentral.com/articles/10.1186/…

@MatthewJ.Oldach Still.... 你应该测试加入，转换（如果它的 data.table 问题是可以解决的）

哎呀，我检查了原始文件，结果发现这些日期是基因名称的产物 cuffdiff 所以它不是上面的代码。对不起

@MatthewJ.Oldach 发生了，这是一个众所周知的 excel 问题，我还没有听说 R 这样做

【参考方案3】：

我无法对您的数据进行尝试，但我确实有一个建议。转换为矩阵并使用 rbind.fill，可能在前后转置矩阵。

祝你好运！

【参考方案4】：

我认为数据的主要问题/困难在于两个/所有数据框都包含重复的键值（该 data.table 然后还抱怨说它需要笛卡尔积）。所以假设 df1 类似于

A 1
B 2
C 3
B 3
A 4

df2 是

A 10
B 20
C 30
B 30
A 50

那么合并就会变成

A 1 ?
A 4 ?
A ? 10
A ? 50
# etc.

? 是我对 NA 的代表，R 将放在那里。这就是 data.table 所抱怨的（请注意，它建议的 by = .EACHI 不起作用，因为 nrow(df1) != nrow(df2)）

因此，您可以算出 66k 行（我检查了三个示例）大约会以多快的速度爆炸。这些重复的一半。所以你必须决定如何处理重复的键。如果您的最终数据框中的问号正常，那么您可以按照我对 PoGibas 答案的评论中的建议将 allow.cartesian = TRUE 添加到 data.table::merge()。

如果您确实担心?'s，那么我认为您可以尝试找到存在于您拥有的所有 22 个文件中的最大基因名称集，然后手动循环它们并以这种方式构建您的最终数据框。但这会丢失您的数据。

也许您可以尝试首先将所有文件连接成一个长 df（22 x 66k 行），并为 FPKM 标签添加一列，然后使用tidyr::spread() 传播它（可能按基因分组并执行spread()出于内存原因，在do() 循环中）？？？

# pseudocode
myfiles <- as.data.frame(f = list.files(...))
hap <- myfiles %>% group_by(f) %>%
  do(
    d <- read.csv(.$f, ...)
    d$fpkm <- names(d)[2]
  ) %>%
  group_by(gene) %>% 
  do(
    d <- spread(., ...) #something with .$fpkm, and the other column
  )

可能有一种 data.table 方法可以更快地执行此操作。