128-R茶话会21-R读取及处理大数据

前言

最近要处理一个100K*1M 左右大小的矩阵，这个矩阵的行为病人记录，列则是每个突变位点的突变信息，记录为0,1,2。

这个矩阵单纯大小就有300多G，我该如何去读取它、处理它呢？

1-如何读取它

首先。毫无疑问的指向data.table 包中的fread。

它有两个优点：

• 效率飞速，自带多线程操作；
• data.table 格式很好地节约内存。

可是，300多G 对我来说还是有些大了。我可不可以分批处理这些数据呢？

1.1-逐行读取数据

使用命令readLines，该函数通过与文件建立某种连接，并设置参数n控制每次读取的行数。通过设置循环，每次固定读取一定行数的文件，并设置循环退出条件为读取结果为零即可：

while( TRUE ){
  # read genotype
  tmp <- readLines(genotype.file, n = N)
  if(length(tmp)==0) break
  tmp <- strsplit(tmp, " ")
  tmp <- do.call("rbind", tmp)
  genotype <- tmp[,-c(1:6)]
  genotype <- apply(genotype, 2, as.numeric)
  genotype <- t(genotype)
  }

但是，这个readLines并不能指定该读哪一行，因为它是逐行读取。也就意味着我必须串行的完成整个文件的处理，排队依次进行。

而如snowfall 等并行处理的包，似乎无法处理readLines 这种文件链接，在我的测试中，每次并行循环都会重建链接，也就是若干个前N 行的文件。

1.2-将数据拆分

那么该如何来并行呢？

不好意思，这里我“作弊”了。除了split命令，我实在想不到其他的办法。也就是非常暴力的将文件拆分：

split -l 1000 -a 2 ../Input/xx.raw ../Input/split/xx_raw_
# -l 设置拆分文件的行数
# -a 用于设置后缀长度，后缀使用字母a-z
# -a 2 则后缀为 aa,ab,ac ...ba,bb ... zz

使用脚本同时处理若干个文件即可。批量处理这些脚本，会在后面的步骤介绍。

2-优化处理过程

首先，我的矩阵是从数据框得到的，而它们读入时被定义为了字符串型，我需要对他们使用转型。

使用apply？来点多线程，mapply? no,no,no。还记得[[125-R编程19-请珍惜R向量化操作的特性]] 吗？

我们将它们直接转型成对应矩阵就好，相当于重新创建了矩阵，接着将矩阵设计成和原矩阵相同的长宽属性。

genotype <- matrix(as.numeric(genotype), ncol = ncol(genotype))

很显然，大部分的记录值都是0，因为纯合野生型占多数，而这样的稀疏矩阵，R 包Matrix 正好满足我的需求。

这个包也非常厉害，可以探索它的说明书：Matrix: Sparse and Dense Matrix Classes and Methods (r-project.org)[1]

关于稀疏矩阵还有一些有趣的探索：(13条消息) R语言的稀疏矩阵学习记录_徐洲更hoptop的博客-CSDN博客[2]

3-写成脚本分别投递

在[[98-R茶话会17-在后台执行R命令]] 我们提过用脚本执行R 命令。

其实脚本非常好写，也就是配置输入与输出：

args <- commandArgs(T)

genotype <- fread(args[1])
genotype.names <- fread(args[2])
# ped.file <- file(args[2], open = "rt")
output.dir <- args[3]

最终的脚本可以写成：

ls ../Input/split | while read id;
do Rscript ../Input/split/${id} ../Input/genotype_name.txt ../Out/${id} ; done

但显然，这样并不能达到我实现并行的目的。

我一共拆分成了100个文件，如何做到同时并行10个脚本呢？也就是1..10,11..20等等，10个为一组。

比如这样的脚本：

for i in `seq 10 10 100`
do 
 cat ./id.text | head -$i | tail -10 | while read id
 do 
  echo "Rscript ./xx.R ../Input/split/${id} ../Input/genotype_name.txt ../Out/${id}" >> script/script_${i}
 done
done

十个这样的脚本，批量执行它们即可：

head -4 script_10
Rscript ./xx.R ../Input/split/id.text ../Input/genotype_name.txt ../Out/id.text 
Rscript ./xx.R ../Input/split/test1 ../Input/genotype_name.txt ../Out/test1
Rscript ./xx.R ../Input/split/test10 ../Input/genotype_name.txt ../Out/test10
Rscript ./xx.R ../Input/split/test100 ../Input/genotype_name.txt ../Out/test100

我先前还写过一个通过取余数来拆分的策略：005. 并行串行的新思路 · 语雀[3]

4-如果我要输出结果

参考：4 Wrangling big data | Exploring, Visualizing, and Modeling Big Data with R[4]

不难发现，data.table::fwrite 又快又省空间。

拓展读物

比如：Exploring, Visualizing, and Modeling Big Data with R (okanbulut.github.io)[5]

就提到了包括data.table 在内的许多处理big data 的方法。

其中The sparklyr package 似乎很有意思，也有一本对应的书：Mastering Spark with R (therinspark.com)[6]

当然，私以为如果是本地几百G 大小的数据处理，R 就可以搞定了，结合脚本实现更加自由的多线程。

如果更大规模的数据量呢？至少我暂时还没有遇到。而且简单的数据处理，linux 中的sed 或awk 也是不错的选择，csvtk 也是一个很好用的软件。

ps：感觉我的这期翻译味好重，奇怪了。

参考资料

[1]

Matrix: Sparse and Dense Matrix Classes and Methods (r-project.org): https://cran.r-project.org/web/packages/Matrix/Matrix.pdf

[2]

(13条消息) R语言的稀疏矩阵学习记录_徐洲更hoptop的博客-CSDN博客: https://blog.csdn.net/u012110870/article/details/115511976

[3]

005. 并行串行的新思路 · 语雀: https://www.yuque.com/mugpeng/sequence/lbpqn5#89Ubp

[4]

4 Wrangling big data | Exploring, Visualizing, and Modeling Big Data with R: https://okanbulut.github.io/bigdata/wrangling-big-data.html#readingwriting-data-with-data.table

[5]

Exploring, Visualizing, and Modeling Big Data with R (okanbulut.github.io): https://okanbulut.github.io/bigdata/

[6]

Mastering Spark with R (therinspark.com): https://therinspark.com/