如何快速且内存高效地将27GB的unitigs.rtab导入R？

问题背景 #

需导入约27GB的unitigs.rtab文件，该文件为二进制矩阵：2806行（对应基因组名称）、5682556列（对应unitig名称），值为0/1表示unitig在基因组中的存在/缺失状态。使用data.table::fread时，即使在925GB内存、8核环境下仍触发段错误。

文件示例（前10行前5列） #

用以下命令查看示例：

head -n 10 "unitigs.rtab" | cut -f 1-5

输出结果：

尝试的导入代码及错误 #

导入代码（scriptA.R）：

library(data.table)
unitig_file <- fread("unitigs.rtab", verbose = TRUE)

触发的段错误信息（翻译后）：

OpenMP版本 (_OPENMP)       201511
omp_get_num_procs()            8
R_DATATABLE_NUM_PROCS_PERCENT  未设置（默认50）
R_DATATABLE_NUM_THREADS        未设置
R_DATATABLE_THROTTLE           未设置（默认1024）
omp_get_thread_limit()         2147483647
omp_get_max_threads()          8
OMP_THREAD_LIMIT               未设置
OMP_NUM_THREADS                未设置
RestoreAfterFork               true
data.table使用4线程，throttle==1024。详见?setDTthreads。
输入不含\n，将其视为文件名打开
[01] 检查参数
  使用4线程（omp_get_max_threads()=8，nth=4）
  NAstrings = [<<NA>>]
  无NA字符串看起来像数字。
  显示进度=0
  0/1列将被读取为整数
[02] 打开文件
  打开文件unitigs.rtab
  文件已打开，大小=27.27GB（29279214177字节）。
  内存映射成功
[03] 检测并跳过BOM
[04] 调整mmap以\0结尾
  在输入中发现\n，不同行可能以不同换行符结尾（如同一文件混合\n和\r\n）。这很常见且符合预期。
[05] 若需要则跳过初始行
  定位到第1行，起始内容：<<Unitig_sequence  AAAAGTTCGATTTA>>
[06] 检测分隔符、引用规则和列数
  自动检测分隔符...
  sep=0x9，100行均有5682556个字段，引用规则0
  在第1行检测到5682556列。该行可能是列名或首行数据。行起始内容：<<Unitig_sequence AAAAGTTCGATTTA>>
  选择引用规则0
  fill=false，找到的最大列数为5682556
[07] 检测列类型、估算行数及判断首行是否为列名
  采样跳转点数量=10，因为（从第1行到文件末尾的29279214176字节）/(2 * 1423294172 jump0size) ==10
  类型代码（跳转000）: C5555555555555555555555555555555555555555555555555555555555555555555555555555555...5555555555  引用规则0
  类型代码（跳转010）: C5555555555555555555555555555555555555555555555555555555555555555555555555555555...5555555555  引用规则0
  判定'header'为true，因为第2列在第1行是字符串，而在1062个采样行的其余部分是更低类型（int32）
  =====
  在11个跳转点采样1062行（处理引号内的\n）
  从第2行首数据行到最后一行末尾的字节数：28981067180
  行长度：均值=11365124.05 标准差=-nan 最小值=11365116 最大值=11365134
  估算行数：28981067180 / 11365124.05 =2551
  初始分配=2806行（2551 +9%），使用bytes/max(mean-2*sd,min)，限制在[1.1*estn, 2.0*estn]之间
  =====
[08] 分配列名
[09] 应用用户对列类型的覆盖
  0个类型覆盖和0个删除覆盖后：C5555555555555555555555555555555555555555555555555555555555555555555555555555555...5555555555
[10] 为数据表分配内存
  分配5682556个列槽（5682556 -0个被删除），共2806行
[11] 读取数据
  jumps=[0..2), chunk_size=14490533590, total_size=28981067180

 *** 捕获段错误 ***
地址0x7f6115d70ebe，原因'memory not mapped'

回溯：
 1: fread("unitigs.rtab",     verbose = TRUE)
发生不可恢复的异常。R正在终止...
job3362275/slurm_script:第12行:14265 段错误 (核心已转储) Rscript scriptA.R

内存高效的导入方案 #

1. 使用稀疏矩阵存储（最推荐） #

由于数据是二进制0/1，且从示例看大部分值为0，稀疏矩阵可将内存占用降低几个数量级，仅存储非零值的位置：

library(Matrix)

# 读取列名（unitig名称）
col_names <- scan("unitigs.rtab", nlines = 1, what = character(), sep = "\t")[-1]
# 读取行名（基因组名称）
row_names <- read.table("unitigs.rtab", skip = 1, colClasses = c("character", rep(NULL, length(col_names))), sep = "\t")[[1]]

# 初始化存储非零值位置的容器
i <- integer()
j <- integer()
x <- integer()

# 逐行读取，仅记录值为1的位置
con <- file("unitigs.rtab", "r")
readLines(con, n = 1)  # 跳过表头
line_idx <- 1
while (length(line_content <- readLines(con, n = 1)) > 0) {
  values <- as.integer(strsplit(line_content, "\t")[[1]][-1])
  one_positions <- which(values == 1)
  if (length(one_positions) > 0) {
    i <- c(i, rep(line_idx, length(one_positions)))
    j <- c(j, one_positions)
    x <- c(x, rep(1, length(one_positions)))
  }
  line_idx <- line_idx + 1
}
close(con)

# 构建稀疏矩阵
unitig_sparse <- sparseMatrix(i = i, j = j, x = x, dimnames = list(row_names, col_names))

2. 指定最小数据类型 #

默认fread会将数值列识别为int32（占4字节），而0/1可用**逻辑型（logical，占1字节）**存储，直接减少75%内存占用，同时关闭多线程避免内存分配冲突：

library(data.table)

# 关闭多线程，避免内存竞争
setDTthreads(1)

# 获取总列数，指定列类型
total_cols <- length(scan("unitigs.rtab", nlines = 1, what = character(), sep = "\t"))
col_types <- c("character", rep("logical", total_cols - 1))

# 导入文件
unitig_dt <- fread("unitigs.rtab", colClasses = col_types)

3. 转置后导入（适配R的列优先存储） #

R采用列优先存储，行少列多的矩阵存储效率低。先用命令行转置文件，再导入R后转回原结构：

# 使用awk转置大型Tab分隔矩阵
awk 'BEGIN{FS=OFS="\t"} {for(i=1;i<=NF;i++) a[i]=a[i] $i OFS; n=NF} END{for(i=1;i<=n;i++) print substr(a[i],1,length(a[i])-length(OFS))}' unitigs.rtab > unitigs_transposed.rtab

导入转置后的文件：

library(data.table)
transposed_dt <- fread("unitigs_transposed.rtab", colClasses = c("character", rep("logical", 2805)))
# 转回原矩阵（可选）
unitig_dt <- t(transposed_dt[, -1, with=FALSE])
colnames(unitig_dt) <- transposed_dt[[1]]
rownames(unitig_dt) <- colnames(transposed_dt)[-1]

4. 使用vroom包（延迟加载） #

vroom基于延迟加载机制，不会一次性将所有数据加载到内存，适合超大型文件：

library(vroom)
# 指定列类型：第一列为字符，其余为逻辑型
unitig_vroom <- vroom("unitigs.rtab", col_types = c(col_character(), rep(col_logical(), 5682555)))

内容的提问来源于stack exchange，提问作者Daisy238