基于Python基础包识别容错DNA序列的唯一序列及计数

解决思路 #

要搞定这个问题，核心是用**编辑距离（Levenshtein距离）**来衡量序列间的差异——这个距离正好对应测序时的插入、删除、替换三种误差，把距离≤2的序列归为同一簇，再为每个簇选一个代表作为“原始序列”，最后统计簇的大小就是对应的reads数。

下面是基于Python基础包（pandas、纯Python实现的工具函数）的完整方案：

完整代码实现 #

先定义必要的工具函数和数据结构：

import pandas as pd
from collections import defaultdict

# 实现Levenshtein编辑距离，计算两个序列的插入/删除/替换操作次数
def levenshtein_distance(s1, s2):
    if len(s1) == 0:
        return len(s2)
    if len(s2) == 0:
        return len(s1)
    
    # 动态规划表：(len(s1)+1)行 x (len(s2)+1)列
    dp = [[0]*(len(s2)+1) for _ in range(len(s1)+1)]
    
    # 初始化边界：空序列到目标序列的距离就是目标长度
    for i in range(len(s1)+1):
        dp[i][0] = i
    for j in range(len(s2)+1):
        dp[0][j] = j
    
    # 填充DP表
    for i in range(1, len(s1)+1):
        for j in range(1, len(s2)+1):
            # 字符相同则替换成本为0，否则为1
            cost = 0 if s1[i-1] == s2[j-1] else 1
            dp[i][j] = min(
                dp[i-1][j] + 1,    # 删除操作
                dp[i][j-1] + 1,    # 插入操作
                dp[i-1][j-1] + cost # 替换操作
            )
    return dp[-1][-1]

# 并查集（Union-Find）数据结构，高效管理序列分组
class UnionFind:
    def __init__(self, size):
        self.parent = list(range(size))
    
    def find(self, x):
        # 路径压缩，提升查询效率
        if self.parent[x] != x:
            self.parent[x] = self.find(self.parent[x])
        return self.parent[x]
    
    def union(self, x, y):
        # 合并两个集合
        fx = self.find(x)
        fy = self.find(y)
        if fx != fy:
            self.parent[fy] = fx

接下来处理你的DataFrame：

# 示例数据（替换成你的实际DataFrame）
df = pd.DataFrame({
    'sequence': ['12344', '12344', '12334', '1234', '123444', 
                 'ATGCTAG', 'ATGCTAGG', 'ATGCTAGCG', 'ATGCTAGC']
})

# 1. 提取序列并统一转为字符串（如果你的序列是数字列表/整数，可先转成字符串）
sequences = df['sequence'].astype(str).tolist()
total_sequences = len(sequences)

# 2. 优化版：先按序列长度分组，只比较长度差≤2的序列（减少不必要的计算）
groups_by_length = defaultdict(list)
for idx, seq in enumerate(sequences):
    groups_by_length[len(seq)].append((idx, seq))

# 3. 初始化并查集，开始分组
uf = UnionFind(total_sequences)

for length in groups_by_length:
    current_group = groups_by_length[length]
    # 组内序列两两比较
    for i in range(len(current_group)):
        idx1, seq1 = current_group[i]
        for j in range(i+1, len(current_group)):
            idx2, seq2 = current_group[j]
            if levenshtein_distance(seq1, seq2) <= 2:
                uf.union(idx1, idx2)
    # 和长度+1的组比较
    if length + 1 in groups_by_length:
        for idx1, seq1 in current_group:
            for idx2, seq2 in groups_by_length[length + 1]:
                if levenshtein_distance(seq1, seq2) <= 2:
                    uf.union(idx1, idx2)
    # 和长度+2的组比较
    if length + 2 in groups_by_length:
        for idx1, seq1 in current_group:
            for idx2, seq2 in groups_by_length[length + 2]:
                if levenshtein_distance(seq1, seq2) <= 2:
                    uf.union(idx1, idx2)

# 4. 整理每个簇的序列，并选定代表原始序列
clusters = defaultdict(list)
for idx, seq in enumerate(sequences):
    root_idx = uf.find(idx)
    clusters[root_idx].append(seq)

# 生成结果：这里选择「最短序列→出现次数最多」的序列作为原始序列（符合你的示例逻辑）
result_list = []
for cluster_seqs in clusters.values():
    # 统计簇内各序列的出现次数
    seq_counts = pd.Series(cluster_seqs).value_counts()
    # 排序规则：先按长度升序，再按出现次数降序
    sorted_seqs = sorted(seq_counts.index, key=lambda x: (len(x), -seq_counts[x]))
    original_seq = sorted_seqs[0]
    total_reads = len(cluster_seqs)
    result_list.append({
        'original_sequence': original_seq,
        'reads': total_reads
    })

# 转成最终的结果DataFrame
result_df = pd.DataFrame(result_list)
print(result_df)

代码说明 #

1. Levenshtein距离：精准计算两个序列间的插入/删除/替换操作次数，完美匹配你提到的测序误差类型，距离≤2就视为同一原始序列。
2. 并查集：高效管理序列分组，避免重复标记同一簇的序列，比暴力分组效率高很多。
3. 长度分组优化：两个序列长度差超过2时，编辑距离必然大于2，所以只比较长度差≤2的序列，能大幅减少计算量。
4. 代表序列选择：示例中选最短的序列作为原始序列，若有多个长度相同的序列，选出现次数最多的。你可以根据需求修改排序规则（比如直接选出现次数最多的序列）。

运行结果 #

对于示例数据，输出结果如下：

original_sequence  reads
0              1234      5
1           ATGCTAG      4

完全符合你给出的示例预期。

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