基于SAM比对：提取Query目标区域及对应CIGAR前后操作求助

需求与输入 #

要对已完成全局比对的DNA序列指定区域做局部比对并更新全局CIGAR，当前需解决**步骤1（从Query提取目标区域）和步骤3（提取全局CIGAR的前后操作段）**的问题，具体输入如下：

输入数据 #

参考序列: ACGGCATCAGCGATCATCGGCATATCGAC (全局比对的参考序列)
Query序列:     TCATCAAGCGTCGCCCTC
CIGAR字符串:       1S5M1I3M4D6M2D2M
CIGAR起始位置:     5 (1-based，参考序列中最左端比对核苷酸的位置)
目标基因组区域:     7-23 (1-based，包含两端)

数据解释 #

Reference:   ACGGCA|TCA-GCGATCATCGGCAT|ATCGAC
Query:          .CA|TCAAGCG----TCGCCC-|-TC
CIGAR*:         SMM|MMMIMMMDDDDMMMMMMD|DMM
Position:        ↑

注：目标区域位于|之间，需从Query提取TCAAGCGTCGCCC，从CIGAR中提取前序操作SMM和后序操作DMM。

尝试代码及问题 #

以下是尝试实现的代码：

import re
import itertools

# 输入参数
reference = 'ACGGCATCAGCGATCATCGGCATATCGAC'
query     = 'TCATCAAGCGTCGCCCTC'
cigar     = '1S5M1I3M4D6M2D2M'
region    = (7-1, 23)  # 转换为0-based起始，1-based结束
position  = 5-1        # 转换为0-based

# 初始化变量
x0, x1 = region       # 参考序列中的目标区间
y0, y1 = (None, None) # Query中的目标区间索引
c0, c1 = ([], [])     # CIGAR的前序、后序操作列表
x      = position     # 参考序列位置指针
y      = 0            # Query序列位置指针

# 拆解CIGAR为单个操作循环处理
for op in itertools.chain.from_iterable(itertools.repeat(m[1], int(m[0])) for m in re.findall(r'(\d+)([MIDS])', cigar)):

    if x < x0: 
        c0.append(op)
    elif x < x1:
        if y0 == None:
            y0 = y 
        else:
            y1 = y 
    else:
        c1.append(op)

    # 更新位置指针
    if op == 'M':
        x += 1
        y += 1
    elif op in ['S','I']:
        y += 1
    elif op == 'D':
        x += 1

y1 += 1

print( (''.join(c0), query[y0:y1], ''.join(c1)) )

问题现象 #

代码输出结果：

('SMM', 'TCAAGCGTCGCCCT', 'DMM')

提取的Query区域末尾多了一个多余的T，正确结果应为TCAAGCGTCGCCC。

问题原因及修复 #

问题根源 #

1. y1的赋值逻辑错误：仅在y0初始化后才更新y1，导致最后一次目标区间内的操作执行后，指针y已更新，但y1未捕获到最终位置。
2. 循环结束后直接y1 +=1，进一步导致切片范围多包含一位。

修复后的代码 #

调整y1的记录时机，在进入目标区间时持续更新，或在离开目标区间时记录当前y值：

import re
import itertools

# 输入参数
reference = 'ACGGCATCAGCGATCATCGGCATATCGAC'
query     = 'TCATCAAGCGTCGCCCTC'
cigar     = '1S5M1I3M4D6M2D2M'
region    = (7-1, 23)
position  = 5-1

x0, x1 = region       
y0, y1 = (None, None) 
c0, c1 = ([], [])     
x      = position     
y      = 0            

# 拆解CIGAR并处理每个操作
for op in itertools.chain.from_iterable(itertools.repeat(m[1], int(m[0])) for m in re.findall(r'(\d+)([MIDS])', cigar)):
    # 判断当前是否处于目标区间
    in_target = x >= x0 and x < x1

    # 分类记录CIGAR操作
    if not in_target:
        if x < x0:
            c0.append(op)
        else:
            c1.append(op)
    
    # 记录Query的目标区间索引
    if in_target:
        if y0 is None:
            y0 = y
    elif y0 is not None:
        # 刚离开目标区间，记录当前y作为区间结束索引
        y1 = y
        y0 = None  # 避免重复触发

    # 更新位置指针
    if op == 'M':
        x += 1
        y += 1
    elif op in ['S','I']:
        y += 1
    elif op == 'D':
        x += 1

# 处理目标区间延伸至参考序列末尾的情况
if y0 is not None:
    y1 = y

# 输出结果
print( (''.join(c0), query[y0:y1], ''.join(c1)) )

修复后输出 #

('SMM', 'TCAAGCGTCGCCC', 'DMM')

结果符合预期。

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