如何在Python中基于汉明距离聚类相似DNA序列？

阿华AIGC实验室

2026-5-20

优化DNA序列汉明距离聚类的思路与实现

看来你当前用的“取首个序列为中心、阈值内归为一类”的思路是聚类的雏形，但确实存在两个核心问题：一是初始中心的随机性（直接取第一个）会影响聚类结果的合理性，二是只做一轮划分，没有迭代优化的过程，很容易把相似序列分到不同类，或者把不那么相似的硬凑到一起。

下面分享几个实用的优化方向，附Python代码示例（生物信息学场景常用）：

一、先优化汉明距离的计算效率

如果你的序列数量多、长度长，循环计算两两距离会非常慢。可以把DNA序列转成数值向量，用numpy的向量化操作大幅提速：

import numpy as np

def dna_to_vector(seq):
    # 把DNA字符映射为数值：A=0, T=1, C=2, G=3
    mapping = {'A':0, 'T':1, 'C':2, 'G':3}
    return np.array([mapping[c] for c in seq])

def hamming_distance(vec1, vec2):
    # 向量化计算汉明距离，比普通循环快10-100倍
    return np.sum(vec1 != vec2)

二、改进聚类策略：迭代式K-medoids聚类

因为汉明距离是离散的非欧氏距离，K-means并不适用，K-medoids（基于真实序列的中心点聚类）更贴合你的场景——它用列表中真实存在的序列作为中心，避免了K-means的“虚拟中心”问题。完整流程如下：

从序列列表中智能选择初始中心（比如选两两距离最远的K个，避免初始中心过于相似）
计算每个序列到所有中心的汉明距离，将序列分到距离最近的中心对应的类
对每个类，重新选择新中心（类内所有序列两两距离总和最小的那个序列）
重复步骤2-3，直到中心不再变化，或达到最大迭代次数

完整实现示例

def k_medoids_dna_clustering(seq_list, k, max_iter=100):
    # 把所有序列转成数值向量
    vec_list = [dna_to_vector(seq) for seq in seq_list]
    n = len(vec_list)
    
    # 初始化中心：随机选k个不重复的序列索引
    centers_idx = np.random.choice(n, k, replace=False)
    centers = [vec_list[i] for i in centers_idx]
    
    for _ in range(max_iter):
        # 步骤1：将每个序列分配到最近的中心
        clusters = [[] for _ in range(k)]
        for idx, vec in enumerate(vec_list):
            distances = [hamming_distance(vec, c) for c in centers]
            closest_center = np.argmin(distances)
            clusters[closest_center].append(idx)
        
        # 步骤2：更新每个类的中心（选类内两两距离总和最小的序列）
        new_centers_idx = []
        for cluster in clusters:
            if not cluster:
                # 若类为空，重新随机选一个未被选为中心的序列
                available_idx = [i for i in range(n) if i not in new_centers_idx]
                new_idx = np.random.choice(available_idx)
                new_centers_idx.append(new_idx)
                continue
            
            min_total_dist = float('inf')
            best_idx = cluster[0]
            for idx_candidate in cluster:
                total_dist = 0
                for idx_member in cluster:
                    total_dist += hamming_distance(vec_list[idx_candidate], vec_list[idx_member])
                if total_dist < min_total_dist:
                    min_total_dist = total_dist
                    best_idx = idx_candidate
            new_centers_idx.append(best_idx)
        
        # 若中心无变化，停止迭代
        if set(new_centers_idx) == set(centers_idx):
            break
        centers_idx = new_centers_idx
        centers = [vec_list[i] for i in centers_idx]
    
    # 将聚类结果转回原DNA序列
    result_clusters = []
    for cluster in clusters:
        result_clusters.append([seq_list[i] for i in cluster])
    
    return result_clusters

三、不确定聚类数量K？试试基于阈值的层次聚类

你提到每个列表的序列数量不同，可能没法提前确定K值。这种情况下，基于距离阈值的层次聚类更适合——它不需要指定聚类数，完全按照你设定的相似度阈值划分：

计算所有序列的两两汉明距离矩阵
从每个序列单独成类开始，不断合并距离最小的两个类
直到所有类之间的距离都超过设定阈值，停止合并

简单实现思路

from scipy.cluster.hierarchy import linkage, fcluster
from scipy.spatial.distance import pdist

def threshold_based_clustering(seq_list, threshold):
    # 转换为数值向量并计算两两汉明距离矩阵
    vec_list = [dna_to_vector(seq) for seq in seq_list]
    # pdist计算两两绝对汉明距离（scipy默认hamming是比例，这里自定义计算逻辑）
    distance_matrix = pdist(vec_list, lambda x,y: np.sum(x!=y))
    # 用单链接聚类（也可替换为平均链接、完全链接等）
    Z = linkage(distance_matrix, method='single')
    # 根据阈值划分聚类
    cluster_labels = fcluster(Z, threshold, criterion='distance')
    
    # 整理聚类结果
    clusters = {}
    for label, seq in zip(cluster_labels, seq_list):
        if label not in clusters:
            clusters[label] = []
        clusters[label].append(seq)
    return list(clusters.values())