高效过滤超大型压缩VCF文件的技术方案咨询

针对8亿+行的压缩VCF文件，直接全量加载到内存不现实，下面是基于pandas的分块流式处理方案，能高效完成你要求的三项过滤：

核心思路 #

通过分块读取压缩文件，逐块执行前两步过滤，将有效数据存入临时文件；最后合并临时文件，完成去重操作，全程避免内存溢出。

具体实现步骤 #

1. 准备工作与参数定义 #

导入依赖库，设置关键参数（可根据硬件配置调整）：

import pandas as pd
import os

# 自定义参数
AC_THRESHOLD = 2
CHUNKSIZE = 1_000_000  # 每个分块的行数，内存充足可调大（如5_000_000）
INPUT_VCF = "input.vcf.gz"  # 你的压缩VCF文件路径
TEMP_DIR = "./vcf_temp_chunks"
OUTPUT_VCF = "filtered_result.vcf"

# 创建临时目录存储分块结果
os.makedirs(TEMP_DIR, exist_ok=True)

2. 提取VCF表头并分块处理 #

VCF文件包含#开头的注释行，先提取表头，再逐块读取数据执行过滤：

# 提取VCF表头（跳过注释行，找到#CHROM开头的表头行）
header = None
with open(INPUT_VCF, "rb") as f:
    for line in f:
        line_str = line.decode().strip()
        if line_str.startswith("#CHROM"):
            header = line_str.split("\t")
            break

# 分块读取并过滤数据
for chunk_idx, chunk in enumerate(pd.read_csv(
    INPUT_VCF,
    sep="\t",
    header=None,
    names=header,
    comment="#",  # 跳过所有#开头的注释行
    compression="gzip",  # 若为其他压缩格式，可改为'zip'或'bz2'
    chunksize=CHUNKSIZE,
    dtype={
        "#CHROM": str,
        "POS": int,
        "REF": str,
        "ALT": str,
        "QUAL": float
    }  # 指定数据类型，避免自动推断浪费资源
)):
    # 第一步过滤：仅保留REF和ALT长度均为1的行
    chunk = chunk[(chunk["REF"].str.len() == 1) & (chunk["ALT"].str.len() == 1)]
    if chunk.empty:
        continue
    
    # 第二步过滤：提取INFO中的AC值，保留AC≥阈值的行
    # 若INFO中AC是第一个键值对，用split比正则更快：
    # chunk["AC"] = chunk["INFO"].str.split(";", n=1).str[0].str.split("=").str[1].astype(int)
    # 通用正则提取方式：
    chunk["AC"] = chunk["INFO"].str.extract(r"AC=(\d+)", expand=False).astype(int)
    chunk = chunk[chunk["AC"] >= AC_THRESHOLD]
    if chunk.empty:
        continue
    
    # 保存过滤后的分块到临时文件
    chunk.to_csv(f"{TEMP_DIR}/chunk_{chunk_idx}.csv", sep="\t", index=False)

3. 合并临时文件并去重 #

合并所有有效分块，按#CHROM和POS去重，保留QUAL最高的行：

# 读取所有临时分块文件
filtered_chunks = []
for file_name in os.listdir(TEMP_DIR):
    if file_name.startswith("chunk_"):
        chunk_df = pd.read_csv(f"{TEMP_DIR}/{file_name}", sep="\t")
        filtered_chunks.append(chunk_df)

combined_df = pd.concat(filtered_chunks, ignore_index=True)

# 按#CHROM和POS分组，保留QUAL最大的行（高效实现）
max_qual_indices = combined_df.groupby(["#CHROM", "POS"])["QUAL"].idxmax()
final_df = combined_df.loc[max_qual_indices]

# 输出结果（恢复VCF格式，可添加注释行如果需要）
final_df.to_csv(OUTPUT_VCF, sep="\t", index=False)

# 清理临时文件
for file_name in os.listdir(TEMP_DIR):
    os.remove(f"{TEMP_DIR}/{file_name}")
os.rmdir(TEMP_DIR)

优化建议 #

• 调整CHUNKSIZE：根据可用内存大小调整，内存越大，分块行数越多，处理速度越快
• AC提取优化：若INFO列中AC始终是第一个键值对，用split方法替代正则表达式，能显著提升处理速度
• 去重方式选择：使用groupby+idxmax比先排序再去重更高效，尤其是数据量极大时
• 压缩格式适配：若你的VCF是其他压缩格式，修改compression参数即可（如zip、bz2）

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