如何用Pandas或替代工具高效查询大型VCF衍生文件

我从16GB的dbSNP VCF文件生成了00-All_relevant.vcf.gz，使用的shell命令如下：

gzcat 00-All.vcf.gz | grep -v ## | awk -v FS='\t' -v OFS='\t' '{print $3, $1, $2, $4, $4}' | gzip > 00-All_relevant.vcf.gz

我需要从这个文件中查询多个rsID（对应文件中的ID列），获取对应的基因组位置和变异信息（#CHROM、POS、REF、ALT列）。项目中使用Pandas，也接受其他Python兼容工具。

直接用Pandas读取全量文件会导致内核崩溃：

import pandas as pd

rsid_df = pd.read_csv('00-All_relevant.vcf.gz',
                      sep='\t')

于是我采用分块读取的方式，实现了如下代码：

rsid_df = pd.read_csv('00-All_relevant.vcf.gz',
                      sep='\t',
                      chunksize=1000000)

rsids = ['rs537152180','rs376204250','rs181326522']
variants = pd.DataFrame()

for data in rsid_df:
    rsid_found = data['ID'].isin(rsids)
    if rsid_found.astype(int).sum() > 0:
        variant = data.loc[rsid_found]
        variants = pd.concat([variants,variant])
        for id in variant['ID'].tolist():
            rsids.remove(id)
        if not rsids:
            break

print(variants)

这个方案可行但速度极慢，尤其是当目标rsID位于最后一个分块时。请问有什么方法可以提速？

1. Shell层提前过滤（最快解决方案） #

既然目标是精准匹配特定rsID，直接用系统级文本工具提前过滤，避免读取全量数据：

1. 先把目标rsID写入文本文件rs_list.txt，每行一个ID
2. 用zgrep快速过滤：

zgrep -w -F -f rs_list.txt 00-All_relevant.vcf.gz > filtered_results.vcf

参数说明：

• -w：匹配完整单词，避免rsID被部分匹配（比如rs123不会匹配rs1234）
• -F：按固定字符串匹配而非正则，大幅提升速度
• -f：从指定文件读取匹配模式

之后直接用Pandas读取小体积的过滤结果：

import pandas as pd
variants = pd.read_csv('filtered_results.vcf', sep='\t', names=['ID', '#CHROM', 'POS', 'REF', 'ALT'])

2. 优化Pandas分块逻辑 #

如果必须用Python处理，针对现有代码做核心优化：

import pandas as pd

# 提前指定列名和数据类型，减少Pandas的推断开销
col_names = ['ID', '#CHROM', 'POS', 'REF', 'ALT']
dtypes = {'ID': str, '#CHROM': str, 'POS': int, 'REF': str, 'ALT': str}

# 用集合存储目标rsID，in操作比列表快100倍以上
target_rsids = {'rs537152180','rs376204250','rs181326522'}
result_chunks = []

for chunk in pd.read_csv('00-All_relevant.vcf.gz',
                         sep='\t',
                         chunksize=1_000_000,
                         names=col_names,
                         dtype=dtypes):
    matched = chunk[chunk['ID'].isin(target_rsids)]
    if not matched.empty:
        result_chunks.append(matched)
        # 批量移除已找到的ID，避免重复匹配
        target_rsids -= set(matched['ID'])
        if not target_rsids:
            break

# 最后一次性合并结果，避免多次concat的性能损耗
variants = pd.concat(result_chunks, ignore_index=True)
print(variants)

核心优化点：

• 集合替代列表存储目标ID，大幅提升匹配效率
• 用列表收集结果块，最后一次合并，减少内存碎片和合并开销
• 显式指定列名和数据类型，跳过Pandas的自动推断步骤

3. 用Dask并行处理 #

如果需要处理更大规模数据，或者利用多CPU核心，用Dask DataFrame（语法与Pandas兼容）：

import dask.dataframe as dd

col_names = ['ID', '#CHROM', 'POS', 'REF', 'ALT']
dtypes = {'ID': str, '#CHROM': str, 'POS': int, 'REF': str, 'ALT': str}

target_rsids = ['rs537152180','rs376204250','rs181326522']

# Dask自动分块并并行读取
ddf = dd.read_csv('00-All_relevant.vcf.gz',
                  sep='\t',
                  names=col_names,
                  dtype=dtypes)

# 筛选后计算出结果
variants = ddf[ddf['ID'].isin(target_rsids)].compute()
print(variants)

Dask会自动利用多CPU核心并行处理分块，比单线程Pandas分块速度提升2-8倍（取决于核心数）。

4. 构建索引文件（适合多次查询） #

如果需要反复查询不同rsID，提前构建索引：

# 生成rsID到行号的映射文件
gzcat 00-All_relevant.vcf.gz | awk '{print $1, NR}' > rsid_index.txt

之后查询时，先从索引文件找到目标ID的行号，再用sed直接读取对应行：

# 比如查rs537152180
row_num=$(grep -w 'rs537152180' rsid_index.txt | awk '{print $2}')
gzcat 00-All_relevant.vcf.gz | sed -n "${row_num}p"

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