查找同染色体区间重叠位点及Python代码索引报错修复

技术需求 #

查找与指定区间重叠且位于同一染色体上的位点。

位点文件示例（实际含10万行） #

with open(file_path, 'r') as f:
    lines = [l for l in f if not l.startswith('#')]
print(lines)
# 输出结果：
['chr1\t36931696\t.\tT\t.\t100\tPASS\tDP=839\tGT:GQ:AD:DP:VF:NL:SB:NC\t0/.:100:830:839:0.0107:24:-100.0000:0.0071\n', 
 'chr2\t25457280\t.\tA\t.\t100\tPASS\tDP=1410\tGT:GQ:AD:DP:VF:NL:SB:NC\t0/0:19:1403:1410:0.0050:24:-100.0000:0.0014\n', 
 '\n', '\n']

# 这里仅保留了2行有效数据，实际文件包含10万行

区间文件示例（实际含500个元组） #

print(bedregions)
# 输出结果：
[('chr1', 36931694, 36931909, 'CSF3R.exon.17.line.1.chr1.36931697.36932509--tile--1.probe--coordinates(36931694-36931909)'), 
 ('chr2', 25466989, 25467211, 'DNMT3A.CDS.17.line.57.merged--with.DNMT3A.CDS.16.li.probe--coordinates(25466989-25467211)')]

# 这里仅保留了2个元组，实际包含500个元组

当前问题与修复方案 #

问题描述 #

编写roi2函数实现上述需求时，触发IndexError: list index out of range错误。

错误原因分析 #

常见触发该错误的场景：

• 处理位点文件时未过滤空行（如示例中的\n），分割空行时会得到空列表，访问索引会报错
• 位点行格式不符合预期，分割后字段数量不足，导致索引越界

修正后的代码实现 #

def roi2(vcf_lines, bed_regions):
    matched_sites = []
    for line in vcf_lines:
        # 跳过空行
        line = line.strip()
        if not line:
            continue
        # 分割位点字段
        fields = line.split('\t')
        # 过滤格式异常的行
        if len(fields) < 2:
            continue
        try:
            chrom = fields[0]
            pos = int(fields[1])
        except (ValueError, IndexError):
            continue
        # 匹配区间
        for region in bed_regions:
            reg_chrom, reg_start, reg_end, _ = region
            if chrom == reg_chrom and reg_start <= pos <= reg_end:
                matched_sites.append(line)
                break
    return matched_sites

# 使用示例
result = roi2(lines, bedregions)
for site in result:
    print(site)

代码说明 #

1. 空行过滤：对每行做strip()处理，空行直接跳过
2. 格式校验：检查分割后的字段数量，转换位置为整数时捕获异常，避免格式错误的行导致崩溃
3. 重叠判断：染色体匹配后，判断位点位置是否在区间起始和结束之间（包含边界）
4. 效率优化：找到匹配区间后立即跳出循环，减少不必要的遍历

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