基于文本模式列表高效筛选大型VCF文件的最优方法

问题背景 #

我有一个大型制表符分隔的VCF（遗传变异文件），名为file.vcf，包含数百万行数据，格式示例如下：

locus1    1    15    0    0/0,21,2,2,;0
locus1    2    17    0    0/0,21,2,1,;0
locus2    1    10    0    0/1,21,2,2,;0
locus3    1    2     0    0/1,21,2,1,;0
...
locus123929    1    3    0    1/0,22,2,1,;0
locus123929    2    4    0    1/2,1,1,3,;0

需要根据search-file.txt中的基因座列表筛选原文件，保留所有匹配基因座的行。search-file.txt内容示例：

locus1
locus3
locus123929

筛选后预期结果：

locus1    1    15    0    0/0,21,2,2,;0
locus1    2    17    0    0/0,21,2,1,;0
locus3    1    2     0    0/1,21,2,1,;0
locus123929    1    3    0    1/0,22,2,1,;0
locus123929    2    4    0    1/2,1,1,3,;0

解决方案 #

1. Bash 高效处理（推荐，速度最快） #

对于超大型文件，awk是最优选择——它将搜索列表加载为哈希表，逐行扫描VCF文件，内存占用极低且处理速度极快：

awk 'NR==FNR { loci[$1]=1; next } $1 in loci' search-file.txt file.vcf > filtered.vcf

• 逻辑说明：
  • NR==FNR：处理第一个输入文件（search-file.txt）时，把每个基因座存入哈希表loci
  • next：跳过后续逻辑，继续处理下一行
  • $1 in loci：处理第二个文件（file.vcf）时，检查第一列是否在哈希表中，匹配则输出该行

如果search-file.txt包含重复行，可以先去重优化：

sort -u search-file.txt > unique-search.txt
awk 'NR==FNR { loci[$1]=1; next } $1 in loci' unique-search.txt file.vcf > filtered.vcf

2. R 低内存处理 #

避免一次性读取整个大文件，以下两种方法均能控制内存占用：

方法一：逐行读取匹配

# 安装并加载hash包（哈希表查询远快于向量匹配）
install.packages("hash")
library(hash)

# 加载基因座列表到哈希表
loci <- readLines("search-file.txt")
loci_hash <- hash(loci, rep(TRUE, length(loci)))

# 逐行处理VCF文件，匹配则写入输出
con <- file("file.vcf", "r")
out_con <- file("filtered.vcf", "w")
while (length(line <- readLines(con, n = 1)) > 0) {
  locus <- strsplit(line, "\t")[[1]][1]
  if (has.key(locus, loci_hash)) {
    writeLines(line, out_con)
  }
}
close(con)
close(out_con)

方法二：data.table分块读取

如果能接受少量内存占用，data.table的分块读取+键连接效率极高：

library(data.table)

# 加载基因座列表并设置主键
loci_dt <- fread("search-file.txt", col.names = "locus")
setkey(loci_dt, locus)

# 分块读取VCF并匹配输出
chunk_size <- 1e6  # 可根据内存情况调整块大小
con <- file("file.vcf", "r")
out_con <- file("filtered.vcf", "w")
while (TRUE) {
  chunk <- fread(con, nrows = chunk_size, sep = "\t", header = FALSE)
  if (nrow(chunk) == 0) break
  setkey(chunk, V1)
  matched <- loci_dt[chunk, nomatch = 0]
  fwrite(matched, out_con, sep = "\t", col.names = FALSE)
}
close(con)
close(out_con)

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