前言

在日常科研和工作中，常常需要处理大规模数据，诸如GWAS、scRNA的数据均可被称为大数据，其大小常常大于500Mb，甚至超过目前大部分家用电脑的内存上限（>16G）。但是，目前高校内对于编程语言的教授主要集中于基础语法与小规模数据处理的方法，而对于生信方法的教授又仅停留于程序包的调用。但是，大数据处理的思路与方法在实施过程中与小数据有巨大区别，使用小数据处理的方法处理大数据将会大大降低运行速度，以图像处理为例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np
import time

arr = np.random.rand(256, 256)

start_time = time.perf_counter()
for i in range(arr.shape[0]):
    for j in range(arr.shape[1]):
        arr[i, j] = arr[i, j] * 2  # 对每个元素进行简单的操作
end_time = time.perf_counter()
print("数组形状:", arr.shape)
print(f"耗时: {(end_time - start_time)*1000:.4f} ms")


start_time = time.perf_counter()
arr = np.random.rand(256, 256)
arr = arr * 2  # 对整个数组进行操作
end_time = time.perf_counter()
print(f"耗时: {(end_time - start_time)*1000:.4f} ms")


输出为

1
2
3
数组形状: (256, 256)
耗时: 16.2331 ms
耗时: 0.4746 ms

可以看到，使用for循环对于一个$256\times256$的图像进行操作耗时远远长于直接采用numpy自带的矩阵点乘。这是因为numpy其内核由C语言编写，其涉及array的操作均经过优化，所以numpy的矩阵点乘效率远远高于使用for循环。
同样的，R语言中也存在类似的情况，在进行两个表格的配对时（这在GWAS meta分析中非常常见），如果用常规的思路对于A表格中的每行在B表格中进行检索，然后再赋值，则时间非常之慢，而如果直接使用 merge 函数，则可大大降低匹配速度。
正因如此，我书写这篇博客以总结一些常用的大数据处理方法。相较于目前市面上已有的介绍，这篇blog在理论与应用间将更加关注应用；在普适性与个性化之间将更加关注个性化（生物学中的使用）。

数据导入与存储

生物信息学中有多种多样的文件类型：vcf, fasta, bed…幸运的是：绝大多数的文件类型都能用记事本打开（可能会费亿点时间），相较于某些不用特定应用就不能打开的文件而言，其可读性相当好。但是与高可读性对应的常是更大的文件大小与更慢的读取速度（电脑读取速度）。
这里，我们随意下载一个gwas文件，文件格式为tsv格式，大小为1.09Gb，这里tsv, csv文件均可以采用data.table包中的fread函数读取。这里使用fread的好处就是 有一个好看的进度条 更快的读取速度与处理速度。 这里并不考虑fst, parquet等保存与读取方法，因为这些方法虽然有更好的读取速度与压缩功能，但是并不兼容目前生信常见的csv, tsv

1
2
3
library("data.table")

df <- fread("D:/data/bioinfo/knee_pain/3773.v1.0.fastGWA.tsv")

上面为读取代码，保存则使用fwrite函数即可。

1
fwrite(df, "D:/data/bioinfo/knee_pain/3773.v1.0.fastGWA.csv", sep = ",")

这就从tsv文件变成了csv文件（tsv采用tab分割不同列，csv则采用comma）。fwrite的速度比write.csv约快30倍（1 minute to 2 second）

基本操作

数据整理的基本操作模式

数据操作的基本范式并不单指任何一个平台或软件，而是一种概念上的模型，其内部包含多种基本操作（创建、删除、检索、插入、排序、过滤、汇总、分组和连接），保证用户在进行数据处理时可以通过对于这些基本操作的排列组合便达到目的。该博客并不会全部涉及，仅讨论与实践密不可分的基本操作。

基本操作的实现

创建

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# Method 1
# data.table()
DT = data.table(
  ID = c("b","b","b","a","a","c"),
  a = 1:6,
  b = 7:12,
  c = 13:18
)
# Method 2
# as.data.table()
as.data.table(iris) -> iris.dt
# Method 3
# fread()
df = fread("abc.csv")

检索

data.table检索使用”行序号, 列名, 分组“进行，先不考虑分组，我们看几个检索的例子。

1
2
3
4
5
6
>>> df[1,variant_id,]
[1] "rs2531267"
>>> df[1:3,variant_id,]
[1] "rs2531267"   "rs12238997"  "rs144155419"
>>> df[1,c(variant_id, beta),]
[1] "rs2531267" "-0.155689"

即方括号内第一个值为检索的行号，第二个值为检索的列名，可以使用c()或者1:3检索多行或多列。而如果需要删掉某行，则只需要在前面加一个-即可，即

1
>>> df[-1,,]

如果需要返回一个data.table，则需要输入.()格式，即

1
2
3
4
>>> df[1, .(variant_id),]
   variant_id
       <char>
1:  rs2531267

插入

使用:=来加入一个新列，熟悉python的人可能知道，这被称为”海象运算符“。该操作可能是最常见的操作之一，例如，我们在harmoised等位基因时，需要将 beta 设置为负值。

1
2
3
>>> df[,reverse_beta := -beta,]
>>> df[1, c(reverse_beta, beta),]
[1]  0.155689 -0.155689

又或者说，我们需要创建一个全为0的列

1
df[, all_zero := 0,]

删去某一列也可使用该方法

1
df[, reverse_beta := NULL]

值得注意的是，海象运算符:=操作均在其本身进行操作，所以在对某些重要变量进行操作时，记得拷贝备份。
对于行的合并或者插入，则直接采用rbind（列合并也可以采用cbind）

1
2
dt <- data.table(id = 1:3, value = c(10,20,30))
dt <- rbind(dt, data.table(id=4, value=40))

过滤

过滤操作通常对于行进行，例如，我们需要找到GWAS研究中pvalue<0.05的SNP

1
df[p_value < 0.05,,]

当然，也可以采用一些逻辑运算符，如果我希望找到第1个染色体上的显著SNP：

1
df[p_value < 0.5 & chromosome == 1,,]

汇总

虽然汇总在GWAS中并不常用，但是在其他很多统计过程中非常重要，所以我简单举个例子

1
2
3
4
5
6
7
>>> df[,mean(beta),,]
[1] 0.0003188427
# For several columns
>>> df[,lapply(.SD, mean, na.rm=TRUE), .SDcols = c("beta", "N")]
           beta        N
          <num>    <num>
1: 0.0003188427 95979.12

分组

分组操作在GWAS分析中较少见，但是在其他分析中却很常见，例如我们想查看单细胞测序并标注的结果中，不同类细胞之间的基因表达差异，便可以使用分组。分组采用检索中提到的最后一个参数。例如，我想看看不同染色体中SNP的beta均值（无意义，仅作方法演示）。

1
2
3
4
5
6
>>> df[,lapply(.SD, mean, na.rm=TRUE), .SDcols = c("beta"),by = chromosome]
    chromosome          beta
         <int>         <num>
 1:          1  5.526048e-04
 2:          2  1.851875e-04
 ......

这时候就可以做一些相对复杂的操作了，比如我希望1. 筛选出所有pvalue < 0.05 的SNP；2. 对于 beta 求绝对值；3. 按照染色体求 beta 均值：

1
2
3
4
5
6
df[p_value < 0.05,lapply(.SD, function(x) mean(abs(x), na.rm=TRUE)), by = chromosome, .SDcols = c("beta")]
    chromosome      beta
         <int>     <num>
 1:          1 0.1257330
 2:          2 0.1234975
 ......

这里使用了lapply：对于列表中每列进行操作；匿名函数：对于每列进行取绝对值abs()后取均值mean().

连接

这是一个相当重要的操作，因为在 harmonise 操作中需要将两个数据按照 rsid 连接后再进行后续操作。连接操作使用merge()函数。

1
df <- merge(x = df, y = df2, by = "variant_id")

其中分为左连接、右连接、内连接、全连接，分别代表保留左侧（df）表格的所有行all.x=T，保留右侧（df2）所有行all.y=T，保留均存在(默认)的和保留所有行all=T。如果我只想看两个数据中均存在的部分，我们可以直接使用上面的代码进行，如果希望保留所有，则：

1
df <- merge(x = df, y = df2, by = "variant_id", all=T)

merge操作后的表格可能存在重名的列，这些列被加入后坠以示区分，我们也可以修改列名：

1
2
3
4
cols <- c("T", "SE_T", "N", "P_noSPA", "CONVERGE", "reverse_beta")
setnames(df, cols, paste0(cols, ".x"))
cols <- c("INFO")
setnames(df, cols, paste0(cols, ".y"))

实战操作

刚刚我们在介绍基础功能的过程中便已经实现数据的读取、配对的操作，接下来就是需要进行harmonise与分割，最终形成一个可以直接进行Meta分析的结果：

1
2
3
4
5
6
7
8
9
10
11
12
library("data.table")

df <- fread("D:/data/bioinfo/knee_pain/3773.v1.0.fastGWA.tsv")
# fwrite(df, "D:/data/bioinfo/knee_pain/3773.v1.0.fastGWA.csv", sep = ",")

df2 <- fread("D:/data/bioinfo/knee_pain/kneP.tsv2")
df <- merge(x = df, y = df2, by = "variant_id")

cols <- c("T", "SE_T", "N", "P_noSPA", "CONVERGE")
setnames(df, cols, paste0(cols, ".x"))
cols <- c("INFO")
setnames(df, cols, paste0(cols, ".y"))

接下来进行harmonise操作，我们希望比较两个数据库的每个SNP，其效应与非效应基因是否相同：1. 相同则不需要进行后续操作；2. 相反则交换并将beta设置为-beta；3. 既不相同也不相反则删去

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Change the name
cols <- c("T", "SE_T", "N", "P_noSPA", "CONVERGE")
setnames(df, cols, paste0(cols, ".x"))
cols <- c("INFO")
setnames(df, cols, paste0(cols, ".y"))

# Find the same and the flip SNP
df[,same := FALSE,]
df[effect_allele.x == effect_allele.y & other_allele.x == other_allele.y,
   same := TRUE,]
df[,flip := FALSE,]
df[effect_allele.x == other_allele.y & other_allele.x == effect_allele.y,
   flip := TRUE,]

# Change the alleles and betas of the flip
df <- df[same | flip,,]
df[flip == TRUE, `:=`(
  effect_allele.y = effect_allele.x,
  other_allele.y  = other_allele.x,
  beta.y = -beta.y
),]
# or
# df[flip == TRUE, c("effect_allele.y", "other_allele.y", "beta.y") :=
# .(effect_allele.x, other_allele.x, -beta.y)]

这里其实我们已经得到了用于Meta分析的data.table，当然，我们也可以将其分别导出到其他软件中使用：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
select_cols <- c("variant_id", "chromosome.x", "base_pair_location.x",
                 "effect_allele.x", "other_allele.x", "N.x",               
                 "effect_allele_frequency.x", "T.x", "SE_T.x",
                 "P_noSPA.x", "beta.x", "standard_error.x",
                 "p_value.x", "CONVERGE.x")
df1 <- df[, .SD, .SDcols = select_cols]
select_cols <- c(
  "chromosome.y", "base_pair_location.y", "other_allele.y", "effect_allele.y",
  "effect_allele_frequency.y", "beta.y",  "standard_error.y",
  "p_value.y", "INFO.y"
)
df2 <- df[, .SD, .SDcols = select_cols]
fwrite(df1, "D:/data/bioinfo/knee_pain/3773.v1.0.fastGWA_harmonised.csv")
fwrite(df2, "D:/data/bioinfo/knee_pain/kneP.tsv2_harmonised.csv")

这里的.SD代表当前分组中被选择的列的子集，.SDcols则用于指定.SD中的列，这样就可以接受在df的索引外定义的select_cols，如果不使用该方法，直接使用列索引，则列名不能带引号。

参考资料

  1. 实战大数据：基于R语言
  2. Matter, Ulrich. Big data analytics: a guide to data science practitioners making the transition to big data. CRC Press, 2023.
  3. Wickham, Hadley, Mine Çetinkaya-Rundel, and Garrett Grolemund. R for data science. O’Reilly Media, 2023.
  4. Luraschi, Javier, Kevin Kuo, and Edgar Ruiz. Mastering Spark with R: the complete guide to large-scale analysis and modeling. O’Reilly Media, 2019.
  5. 刘艺非. R语言高效能实战:更多数据和更快速度. 人民邮电出版社, 2022.
  6. Larger-Than-Memory Data Workflows with Apache Arrow
  7. Apache Arrow R Cookbook
  8. R interface to Apache Spark
  9. An Introduction to Polars from R
  10. tidypolars
  11. fastverse

Cover image icon: Big data icons created by Eucalyp - Flaticon