如何优化SNPs频率统计循环？高效处理百万级基因组区域计数

我需要统计每100,000个碱基位置区间内的SNPs频率，目前使用已预处理的VCF文件，教授提供的参考代码如下：

inputfile=open("bcftools_snps.txt", 'r')
X=0
for line in inputfile:
  A=line.split()
  if float(A[1]) < 100000:
    X=X+1
print("0-100000=", X)

目前我通过定义新变量逐个统计200,000-300,000等区间，但因基因组长度达7,000,000个碱基，需定义700个变量，效率极低。请问是否有更高效的实现方法？

补充说明：VCF文件即Variant Call Format（变异调用格式），用于记录组装基因组与参考基因组（测序研究时使用的同物种参考基因组）的差异；SNPs即Single Nucleotide Polymorphisms（单核苷酸多态性），指染色体上与参考基因组存在碱基差异的位点，这类突变是基因组进化的重要指标，因此统计其整体及特定区域的频率十分关键。

不用逐个定义变量，用字典或者列表就能批量统计所有区间的SNP数量，以下是两种可行的实现方式：

方法一：使用字典统计 #

字典的键设为区间的起始值（比如0、100000、200000...），值对应该区间的SNP数量。遍历文件时，通过碱基位置计算所属区间的键，直接累加计数：

interval_size = 100000
max_position = 7000000
# 初始化所有区间的计数为0
snp_counts = {i: 0 for i in range(0, max_position, interval_size)}

with open("bcftools_snps.txt", 'r') as inputfile:
    for line in inputfile:
        # 跳过VCF文件的注释行（如果有的话）
        if line.startswith('#'):
            continue
        parts = line.split()
        pos = int(float(parts[1]))  # 碱基位置转整数避免精度问题
        # 计算所属区间的起始值
        interval_start = (pos // interval_size) * interval_size
        # 处理超过最大基因组长度的位点
        if interval_start >= max_position:
            interval_start = max_position - interval_size
        snp_counts[interval_start] += 1

# 输出每个区间的结果
for start in snp_counts:
    end = start + interval_size
    print(f"{start}-{end}={snp_counts[start]}")

方法二：使用列表统计 #

如果区间是连续且固定步长的，用列表更高效，列表索引对应区间的序号（比如索引0对应0-100000，索引1对应100000-200000...）：

interval_size = 100000
max_position = 7000000
# 计算总区间数
total_intervals = max_position // interval_size
# 初始化列表，每个元素对应一个区间的计数
snp_counts = [0] * total_intervals

with open("bcftools_snps.txt", 'r') as inputfile:
    for line in inputfile:
        if line.startswith('#'):
            continue
        parts = line.split()
        pos = int(float(parts[1]))
        # 计算区间索引
        interval_idx = pos // interval_size
        # 处理超过最大位置的位点
        if interval_idx >= total_intervals:
            interval_idx = total_intervals - 1
        snp_counts[interval_idx] += 1

# 输出结果
for idx in range(total_intervals):
    start = idx * interval_size
    end = start + interval_size
    print(f"{start}-{end}={snp_counts[idx]}")

注意事项 #

• 代码中加入了跳过VCF注释行的判断，如果你的预处理文件已经去掉注释，可以删除这部分逻辑。
• 将碱基位置转换为整数，能避免浮点运算的精度误差，统计更准确。
• 两种方法只需遍历一次文件即可完成所有区间的统计，无需手动定义大量变量，效率大幅提升。

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