如何用Hamming距离对FASTQ提取的DNA序列进行聚类（Python）

嘿，我来帮你搞定这个DNA序列聚类的问题！你提到想用Hamming距离来分组相似序列，确实Hamming距离是衡量等长序列差异的好工具，但你纠结的「参考序列」问题其实可以通过聚类算法来绕开——不需要固定某一条当参考，而是让算法自动把相似的序列归到一起。

首先敲个核心前提：Hamming距离只适用于长度完全相同的序列，所以第一步先确保你的所有DNA序列长度一致，不然得先做截断/补全处理，或者换用编辑距离（比如Levenshtein）替代。

1. 先实现Hamming距离计算函数 #

这是整个聚类的基础，先写一个能计算两条等长序列差异碱基数量的函数：

def hamming_distance(seq1, seq2):
    # 先检查序列长度是否一致，不一致直接抛出错误
    if len(seq1) != len(seq2):
        raise ValueError("Hamming distance requires sequences of equal length")
    # 统计对应位置碱基不同的数量
    return sum(c1 != c2 for c1, c2 in zip(seq1, seq2))

2. 小数据集方案：单链聚类（简单直接） #

如果你的序列数量不多（比如几千条以内），单链聚类是个易理解、易实现的方法：只要两个序列的Hamming距离小于你设定的阈值（比如2，代表最多2个碱基差异），就把它们归为同一类。

具体实现代码：

def cluster_sequences(sequences, threshold=2):
    clusters = []
    for seq in sequences:
        # 遍历现有聚类，看当前序列能不能加入
        added = False
        for cluster in clusters:
            # 只要和聚类里任意一条序列的距离小于阈值，就加入这个聚类
            if any(hamming_distance(seq, s) <= threshold for s in cluster):
                cluster.append(seq)
                added = True
                break
        # 如果没找到匹配的聚类，新建一个
        if not added:
            clusters.append([seq])
    # 把聚类列表转换成你需要的元组格式
    return [tuple(cluster) for cluster in clusters]

代码说明： #

• threshold：你可以根据实验需求调整，比如允许1个碱基差异就设为1，这个值取决于测序误差容忍度或者你的研究目标。
• 单链聚类的逻辑是「只要沾边就归为一类」，能把所有互相相似（或间接相似）的序列都拉进同一个组里。

3. 大数据集优化：DBSCAN结合自定义距离 #

如果你的序列数量很大（几万条以上），上面的单链聚类效率会很低，这时候可以用sklearn的DBSCAN算法，自定义Hamming距离作为度量：

from sklearn.cluster import DBSCAN
from sklearn.metrics.pairwise import pairwise_distances

def cluster_with_dbscan(sequences, threshold=2):
    # 先检查所有序列长度是否一致
    seq_lengths = set(len(s) for s in sequences)
    if len(seq_lengths) > 1:
        raise ValueError("All sequences must have the same length for Hamming distance")
    
    # 计算所有序列两两之间的Hamming距离矩阵
    distance_matrix = pairwise_distances(
        [[s] for s in sequences],  # 包装成二维数组适配sklearn格式
        metric=lambda x, y: hamming_distance(x[0], y[0])
    )
    
    # 用DBSCAN聚类：eps是距离阈值，min_samples设为1允许单序列聚类
    dbscan = DBSCAN(eps=threshold, min_samples=1, metric="precomputed")
    labels = dbscan.fit_predict(distance_matrix)
    
    # 根据聚类标签分组，转成元组列表
    clusters = {}
    for idx, label in enumerate(labels):
        if label not in clusters:
            clusters[label] = []
        clusters[label].append(sequences[idx])
    return [tuple(cluster) for cluster in clusters.values()]

为什么选DBSCAN？ #

DBSCAN是基于密度的聚类算法，能自动识别噪声（和其他序列都不相似的单独序列），而且处理大数据集的效率比单链聚类高得多。

4. 额外注意事项 #

• 序列去重：如果你的序列里有完全重复的，可以先去重再聚类，能大幅减少计算量：unique_sequences = list(set(sequences))
• 阈值选择：这个要结合测序数据质量定，比如如果测序错误率是0.1%，100bp的序列就可以把阈值设为1。
• 长度不一致处理：如果序列长度不同，要么统一截断到最短序列的长度，要么改用Levenshtein距离（编辑距离）来衡量差异。

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