如何在Python中比较含空位的不等长序列碱基相似度？

首先，我们需要修改原函数，让它能接收两个序列作为输入，同时处理长度差异和空位符号-的规则：只要某位置任一序列是-，该位置视为不相似；仅当两个位置都是A/T/G/C且碱基相同时，才算匹配。

修改后的代码 #

def calculate_sequence_similarity(seq1, seq2):
    matches = 0
    valid_positions = 0
    
    # 遍历两个序列的对应位置，自动取较短序列的长度作为遍历上限
    for base1, base2 in zip(seq1, seq2):
        # 跳过有空位的位置，不计入有效比较
        if base1 == '-' or base2 == '-':
            continue
        # 统计有效比较位置
        valid_positions += 1
        # 统计匹配的碱基
        if base1 == base2:
            matches += 1
    
    # 处理无有效比较的边界情况
    if valid_positions == 0:
        return 0.0
    
    # 计算相似度百分比
    similarity_percent = (matches / valid_positions) * 100
    return {
        "similarity_percent": round(similarity_percent, 2),
        "total_matches": matches,
        "valid_comparison_positions": valid_positions
    }

# 测试示例序列
seq1 = "AAAATCCCTAGAAAGGTCAT"
seq2 = "AAGATC---TTTCTACT"
result = calculate_sequence_similarity(seq1, seq2)
print(result)

代码说明 #

1. 参数调整：函数改为接收seq1和seq2两个参数，不再只处理单个序列。
2. 遍历逻辑：用zip()同时遍历两个序列，自动以较短序列的长度为上限，解决长度不同的问题。
3. 空位处理：只要当前位置任一碱基是-，直接跳过该位置，不纳入有效统计。
4. 计数逻辑：
   • valid_positions统计两个碱基都不是-的位置数（有效比较的总次数）
   • matches统计有效位置中碱基相同的次数
5. 边界处理：如果没有有效比较位置（比如两个序列全是-），返回0避免除以0错误。
6. 返回结果：返回包含相似度百分比、匹配数、有效位置数的字典，结果更直观，还对百分比做了四舍五入保留两位小数。

示例输出 #

运行测试代码后，会得到：

{'similarity_percent': 42.86, 'total_matches': 6, 'valid_comparison_positions': 14}

对应示例序列的计算：有效比较位置共14个，其中6个位置碱基匹配，相似度约为42.86%。

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