Python实现基因区间内测序读段片段匹配标记的代码问题排查

需求说明 #

• 筛选所有读段区间完全落在基因区间（第5列genestart、第6列genestop）内的测序读段
• 首个符合条件且未被其他标记读段覆盖的读段，添加后缀基因名+a；后续符合条件且未被覆盖的读段依次添加基因名+b、基因名+c等
• 不符合条件（读段不在基因内，或已被其他标记读段覆盖）的读段标记为covered

输入数据（制表符分隔） #

contig  230  3888 read1  1 8397 Gene1
contig  3343 3885 read2  1 8397 Gene1
contig  6030 8257 read4  1 8397 Gene1
contig  6030 8257 read10 1 8397 Gene1
contig  6260 8257 read7  1 8397 Gene1
contig  6030 8257 read8  1 8397 Gene1
contig  6190 6345 read10 1 8397 Gene1
contig   0   3124 read27 26 406 Gene12
contig  100  200  read30 26 406 Gene12 

期望输出 #

contig  230  3888 read1  1 8397 Gene1 Gene1a
contig  3343 3885 read2  1 8397 Gene1 covered
contig  6030 8257 read4  1 8397 Gene1 Gene1b
contig  6030 8257 read10 1 8397 Gene1 covered
contig  6260 8257 read7  1 8397 Gene1 covered
contig  6030 8257 read8  1 8397 Gene1 covered
contig  6190 6345 read10 1 8397 Gene1 covered
contig   0   3124 read27 26 406 Gene12 Gene12a
contig  100  200  read30 26 406 Gene12 covered

现有代码 #

with open(input_file, "r") as f_in, open(output_file, "w") as f_out:
    gene_dict = {}
    for line in f_in:
        fields = line.strip().split()
        gene_name = fields[6]
        read_start = int(fields[1])
        read_end = int(fields[2])
        gene_start = int(fields[4])
        gene_end = int(fields[5])
        if gene_name not in gene_dict:
            gene_dict[gene_name] = []
        gene_found = False
        for gene_frag in gene_dict[gene_name]:
            frag_start = gene_frag[0]
            frag_end = gene_frag[1]
            if read_start <= frag_start <= read_end or read_start <= frag_end <= read_end:
                gene_found = True
                break
        if gene_found:
            gene_dict[gene_name].append((gene_start, gene_end))
            gene_count = chr(ord('a') + len(gene_dict[gene_name]) - 1)
            fields.append(f"{gene_name}{gene_count}")
        else:
            gene_dict[gene_name] = [(gene_start, gene_end)]
            fields.append(f"{gene_name}")
        f_out.write("\t".join(fields) + "\n")

实际输出 #

contig  230  3888 read1  1 8397 Gene1 Gene1
contig  3343 3885 read2  1 8397 Gene1 Gene1
contig  6030 8257 read4  1 8397 Gene1 Gene1
contig  6030 8257 read10 1 8397 Gene1 Gene1
contig  6260 8257 read7  1 8397 Gene1 Gene1
contig  6030 8257 read8  1 8397 Gene1 Gene1
contig  6190 6345 read10 1 8397 Gene1 Gene1
contig   0   3124 read27 26 406 Gene12 Gene12
contig  100  200  read30 26 406 Gene12 Gene12

错误排查 #

1. 核心逻辑颠倒：

   • 需求是判断读段是否在基因区间内，但现有代码反向判断基因区间是否在读段内
   • gene_found的逻辑完全错误：现在是检测基因区间是否和已存片段重叠，实际应该检测当前读段是否和已标记的读段重叠

2. 存储对象错误：

   • 现有字典存储的是基因的区间，而实际需要存储已标记的符合条件的读段区间，用来判断后续读段是否被覆盖

3. 标记逻辑错误：

   • 现有代码在gene_found为True时添加新标记，False时重置字典，完全和需求相反；正确逻辑应该是：
     • 先判断读段是否在基因内
     • 若不在，标记covered
     • 若在，再检查是否和已标记读段重叠：无重叠则添加新标记并存储读段区间，有重叠则标记covered

修正后的代码 #

with open(input_file, "r") as f_in, open(output_file, "w") as f_out:
    # 字典结构：{基因名: [已标记的读段区间列表, 当前标记序号]}
    gene_dict = {}
    for line in f_in:
        fields = line.strip().split()
        gene_name = fields[6]
        read_start = int(fields[1])
        read_end = int(fields[2])
        gene_start = int(fields[4])
        gene_end = int(fields[5])
        
        # 初始化基因对应的存储项：[已标记读段区间, 当前标记索引]
        if gene_name not in gene_dict:
            gene_dict[gene_name] = [[], 0]
        
        # 第一步：判断读段是否完全落在基因区间内
        if not (read_start >= gene_start and read_end <= gene_end):
            fields.append("covered")
            f_out.write("\t".join(fields) + "\n")
            continue
        
        # 第二步：判断当前读段是否和已标记的读段有重叠
        overlapped = False
        for seg_start, seg_end in gene_dict[gene_name][0]:
            # 区间重叠判断：两个区间有交集的情况
            if not (read_end < seg_start or read_start > seg_end):
                overlapped = True
                break
        
        if overlapped:
            fields.append("covered")
        else:
            # 生成标记后缀：a, b, c...
            gene_dict[gene_name][1] += 1
            suffix = chr(ord('a') + gene_dict[gene_name][1] - 1)
            fields.append(f"{gene_name}{suffix}")
            # 存储当前读段区间
            gene_dict[gene_name][0].append((read_start, read_end))
        
        f_out.write("\t".join(fields) + "\n")

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