短蛋白序列相似性分析技术咨询及CD-Hit替代方案求助

一、现代序列聚类与相似性工具（2020年后仍活跃维护） #

这些工具兼容性更强，无需手动编译（可通过包管理器安装），且针对短序列优化：

• MMseqs2：速度远超CD-Hit，支持短蛋白序列聚类/相似性分析，用CMake构建避免编译冲突。可通过conda直接安装：conda install -c bioconda mmseqs2，核心聚类命令：

  # 先将CSV序列转为FASTA格式（R/Python均可快速实现）
  mmseqs easy-cluster input.fasta output_dir tmp_dir --min-seq-id 0.9
  

• Linclust：MMseqs2内置模块，专门针对短序列优化，聚类速度比MMseqs2常规模式更快，适合大规模短蛋白数据集。
• VSEARCH：原本用于扩增子分析，现已支持蛋白序列聚类，纯C编写兼容性极佳，brew或conda均可安装：brew install vsearch，聚类命令类似CD-Hit：

  vsearch --cluster_fast input.fasta --id 0.9 --centroids centroids.fasta --uc clusters.uc
  

二、自定义两两比对实现（R/Python，适合小数据集） #

如果需要精细的两两比对结果（找重复序列和相似区域），可基于现有语言实现：

R实现（利用Biostrings） #

假设CSV文件包含seq_id（序列ID）和sequence（蛋白序列）两列：

library(Biostrings)
library(tidyverse)

# 读取序列数据
seq_df <- read_csv("protein_seqs.csv")
seqs <- AAStringSet(seq_df$sequence)
names(seqs) <- seq_df$seq_id

# 全局比对计算序列一致性矩阵
similarity_matrix <- pairwiseAlignmentMatrix(
  seqs, seqs,
  substitutionMatrix = "BLOSUM62",
  gapOpening = 10, gapExtension = 1,
  type = "global", scoreOnly = FALSE
)

# 转换为一致性百分比
identity_matrix <- apply(similarity_matrix, c(1,2), function(x) {
  if (inherits(x, "PairwiseAlignmentsSingle")) pid(x) else 0
})

# 筛选高相似配对（阈值设为90%，可调整）
high_similar_pairs <- which(identity_matrix >= 90 & upper.tri(identity_matrix), arr.ind = TRUE)
high_similar_df <- data.frame(
  seq1 = rownames(identity_matrix)[high_similar_pairs[,1]],
  seq2 = colnames(identity_matrix)[high_similar_pairs[,2]],
  similarity = identity_matrix[high_similar_pairs]
)

# 提取完全重复的序列对
duplicate_pairs <- which(identity_matrix == 100 & upper.tri(identity_matrix), arr.ind = TRUE)
duplicates_df <- data.frame(
  seq1 = rownames(identity_matrix)[duplicate_pairs[,1]],
  seq2 = colnames(identity_matrix)[duplicate_pairs[,2]]
)

Python实现（利用Biopython，适合学习阶段） #

from Bio import pairwise2
from Bio.Seq import Seq
from Bio.Align import substitution_matrices
import pandas as pd
import numpy as np

# 读取CSV序列
df = pd.read_csv("protein_seqs.csv")
seq_ids = df["seq_id"].tolist()
sequences = [Seq(seq) for seq in df["sequence"].tolist()]

# 加载BLOSUM62替换矩阵
blosum62 = substitution_matrices.load("BLOSUM62")

# 初始化一致性矩阵
n = len(sequences)
identity_matrix = np.zeros((n, n))

# 两两全局比对计算一致性
for i in range(n):
    for j in range(i+1, n):
        alns = pairwise2.align.globalds(sequences[i], sequences[j], blosum62, -10, -1)
        best_aln = alns[0]
        # 计算一致性百分比
        match_count = sum(a == b for a, b in zip(best_aln.seqA, best_aln.seqB))
        identity = (match_count / max(len(sequences[i]), len(sequences[j]))) * 100
        identity_matrix[i][j] = identity
        identity_matrix[j][i] = identity

# 转换为DataFrame方便筛选
identity_df = pd.DataFrame(identity_matrix, index=seq_ids, columns=seq_ids)

# 筛选高相似配对
high_similar_pairs = identity_df.where(np.triu(np.ones(identity_df.shape), k=1) > 0)\
                       .stack().reset_index()
high_similar_pairs.columns = ["seq1", "seq2", "similarity"]
high_similar_pairs = high_similar_pairs[high_similar_pairs["similarity"] >= 90]

# 提取完全重复序列
duplicate_pairs = high_similar_pairs[high_similar_pairs["similarity"] == 100]

三、关键注意事项 #

• CSV转FASTA：所有命令行工具均以FASTA为标准输入，可通过R的writeXStringSet(seqs, "input.fasta")或Python的Bio.SeqIO.write快速转换。
• 大数据集优先用工具：自定义两两比对时间复杂度为O(n²)，序列数量超过1000时会显著变慢，MMseqs2/Linclust基于k-mer预聚类，速度提升几个数量级。
• CD-Hit编译问题替代方案：如果仍想使用CD-Hit，直接安装conda预编译版本：conda install -c bioconda cd-hit，跳过手动编译步骤。

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