持久同源对齐（PHA）：使用 ESM-2 和持久同源替代多序列对齐

在论文 利用蛋白质语言模型进行精确多序列比对 和 向量聚类多序列比对：使用蛋白质语言模型比对进入蛋白质序列相似性的“微光区” 中，作者讨论了使用蛋白质语言模型 (pLM) 的嵌入来替代传统多序列比对 (MSA) 的潜力。本文将讨论这种方法的改编，使用 pLM ESM-2 和一种替代方法来确定序列相似性，该方法使用我们称为“持久同源比对”（PHA）的持久同源。

引言

以下方法适用于“微光区”蛋白质。“微光区”在蛋白质序列的语境中指的是一个序列相似性范围，在该范围内，确定两个蛋白质序列是否在进化上相关变得尤为困难。这个概念在生物信息学中至关重要，尤其是在蛋白质序列比对和同源性检测领域。

蛋白质序列通常通过比较来推断其进化关系。这种比较通常使用序列比对技术进行，该技术识别可能由序列之间的功能、结构或进化关系引起的相似区域。相似程度通常用序列同一性来量化，即比对序列中相同氨基酸的百分比。

“微光区”通常指序列同一性在 20% 到 35% 左右的范围。在此范围内：

1. 低序列同一性：序列分歧如此之大，以至于它们的相似性可能是由于遥远的进化关系，也可能仅仅是偶然发生。这使得仅基于序列同一性难以自信地推断同源性（进化相关性）。

2. 结构保守性：尽管序列相似性低，但微光区中的蛋白质通常保留相似的三维结构和功能。这种尽管序列分歧但结构保守性是蛋白质进化中一个显著的方面。

3. 比对挑战：传统的序列比对方法和评分系统（如替换矩阵）在微光区可能难以奏效。它们通常需要额外的信息或复杂的技术，例如那些结合结构或功能数据的技术，才能准确比对序列并推断关系。

4. 进化研究的意义：微光区在进化生物学和生物信息学中具有重要意义，因为研究这些低相似性序列可以深入了解蛋白质结构和功能在漫长进化时间尺度上是如何保守或进化的。

总的来说，比对微光区内的蛋白质提出了生物信息学独特的挑战和机遇，需要先进的方法来辨别其进化和功能关系。

我们将使用的方法基于以下步骤：

1. 嵌入生成：利用蛋白质语言模型（ESM-2），我们为每个蛋白质生成逐位置嵌入，捕获氨基酸的独特上下文和特征。

2. 序列聚类：我们根据蛋白质的序列级嵌入对其进行聚类，将其分组为具有高序列相似性的集合。此步骤有助于识别具有潜在功能或进化关系的蛋白质。此步骤采用持久同源来为每个蛋白质构建条形码图。然后使用 Wasserstein 距离度量和 DBSCAN 对条形码进行聚类，以获得相似蛋白质的簇。

3. 氨基酸相似性和 RBH 识别：在每个簇内，我们测量序列之间氨基酸向量的余弦相似性。我们专注于识别互惠最佳匹配 (RBH)，在不同序列中的氨基酸之间建立对应关系。

4. RBH 网络和引导点聚类：构建一个基于 RBH 关系的网络，然后进行聚类以识别引导点氨基酸。这些引导点作为序列中可靠的比对标记。

5. 将簇排序为列：使用有向无环图 (DAG) 和拓扑排序，我们对簇进行排序以形成 PHA 的列，确保正确的序列排列。

6. 完成比对：算法迭代地将未放置的氨基酸分配到 PHA 列，完善引导点放置。最后一步涉及合并来自各个簇的子比对以创建全面的 PHA。

PHA 基于另一种称为 vcMSA 的方法。有关 vcMSA 方法的视频，请查看此视频。我们将只专注于实现步骤 (1) 和 (2)，因为其余部分已在 vcMSA Github 中实现。用基于持久同源的蛋白质序列聚类来替换 vcMSA 中的步骤 (1) 和 (2) 将更稳健，并且能更好地捕获序列级信息。PHA 方法的其余部分与 vcMSA 技术保持不变。

导入

# Let's start by importing necessary libraries
from transformers import EsmModel, AutoTokenizer
import gudhi as gd
import numpy as np
import torch
from scipy.spatial.distance import pdist, euclidean, squareform
from sklearn.cluster import KMeans, AgglomerativeClustering, DBSCAN
from sklearn.metrics import silhouette_score
from gudhi.hera import wasserstein_distance

辅助函数：

• get_hidden_states：从 ESM-2 模型的特定层提取蛋白质序列的隐藏状态。
• compute_euclidean_distance_matrix_scipy：计算隐藏状态向量之间的成对欧几里得距离。
• compute_persistent_homology：从距离矩阵计算持久同源，生成持久图。
• compute_wasserstein_distances：确定持久图之间的 Wasserstein 距离，这对于比较不同蛋白质的拓扑特征至关重要。

以下 Python 代码封装了一个复杂而全面的工作流，它集成了深度学习模型、蛋白质建模和高级数学技术，用于分析蛋白质序列。该代码涉及多个关键步骤，包括从蛋白质语言模型中提取隐藏状态、计算欧几里得距离矩阵、使用持久同源分析拓扑数据，以及最后对结果进行聚类。让我们详细分解这个工作流的每个部分。

隐藏状态提取

get_hidden_states 函数接收蛋白质序列，对其进行标记化，并将其输入神经网络模型（具体而言，是一个基于 Transformer 的蛋白质语言模型，如 ESM-2），以从指定的层提取隐藏状态。数学上，如果我们有一个输入序列 $XX$ 和一个由权重 $\theta \backslash theta$ 参数化的神经网络模型 $ff$，则在层 $ll$ 的隐藏状态 $HH$ 可以表示为

$$H^{(l)}=f_{\theta }^{(l)}(X)H^\{(l)\}\; =\; f^\{(l)\}\_\{\backslash theta\}(X)$$

其中 $f_{\theta }^{(l)}f^\{(l)\}\_\{\backslash theta\}$ 表示与模型第 $ll$ 层对应的函数。

# Helper function to get the hidden states of a specific layer for a given input sequence
def get_hidden_states(tokenizer, model, layer, input_sequence):
    model.config.output_hidden_states = True
    encoded_input = tokenizer([input_sequence], return_tensors='pt', padding=True, truncation=True)
    with torch.no_grad():
        model_output = model(**encoded_input)
    hidden_states = model_output.hidden_states
    specific_hidden_states = hidden_states[layer][0]
    return specific_hidden_states

欧几里得距离矩阵计算

compute_euclidean_distance_matrix_scipy 函数计算所有向量对（隐藏状态）之间的成对欧几里得距离并返回距离矩阵。如果 $HH$ 是一个矩阵，其中每一行代表一个向量（隐藏状态），则两个向量 $H_{i}H\_i$ 和 $H_{j}H\_j$ 之间的欧几里得距离由下式给出

$$D_{ij}=\mathrm{\parallel }{H}_i-H_j{\mathrm{\parallel }}_{2}D\_\{ij\}\; =\; \backslash |\; H\_i\; -\; H\_j\; \backslash |\_2$$

其中 $\mathrm{\parallel }\cdot {\mathrm{\parallel }}_2\backslash |\; \backslash cdot\; \backslash |\_2$ 表示欧几里得范数。结果矩阵 $DD$ 是一个对称矩阵，其中每个元素 $D_{ij}D\_\{ij\}$ 表示第 $ii$ 个向量和第 $jj$ 个向量之间的距离。

# Helper function to compute the Euclidean distance matrix
def compute_euclidean_distance_matrix_scipy(hidden_states):
    euclidean_distances = pdist(hidden_states.numpy(), metric=euclidean)
    euclidean_distance_matrix = squareform(euclidean_distances)
    return euclidean_distance_matrix

持久同源计算

compute_persistent_homology 使用距离矩阵计算数据的持久同源。持久同源是一种来自拓扑数据分析的方法，通过检查拓扑特征（如连通分量、环等）如何随着阈值参数的变化而出现和消失来研究数据的“形状”。对于给定的距离矩阵 $DD$，构建一个单纯复形，并分析其在不同尺度上的拓扑结构。拓扑特征的持久性通常由持久图或条形码表示。

# Helper function to compute the persistent homology of a given distance matrix
def compute_persistent_homology(distance_matrix, max_dimension=3):
    max_edge_length = np.max(distance_matrix)
    rips_complex = gd.RipsComplex(distance_matrix=distance_matrix, max_edge_length=max_edge_length)
    st = rips_complex.create_simplex_tree(max_dimension=max_dimension)
    persistence = st.persistence()
    return st, persistence

Wasserstein 距离计算

函数 compute_wasserstein_distances 计算持久图之间的成对 Wasserstein 距离。Wasserstein 距离是衡量两个概率分布之间差异的度量，在此上下文中用于量化持久图之间的不相似性。数学上，对于两个持久图 $PP$ 和 $QQ$，Wasserstein 距离 $W(P,Q)W(P,\; Q)$ 定义为

$$W_q^p(P,Q)={(\underset{\gamma \in \mathrm{\Gamma }(P,Q)}{\inf }\sum _{(x,y)\in \gamma }\mathrm{\parallel }x-y{\mathrm{\parallel }}_q^p)}^{1\mathrm{/}p}W\_q^p(P,\; Q)\; =\; \backslash left(\backslash inf\_\{\backslash gamma\; \backslash in\; \backslash Gamma(P,\; Q)\}\; \backslash sum\_\{(x,\; y)\; \backslash in\; \backslash gamma\}\; \backslash |\; x\; -\; y\; \backslash |\_q^p\; \backslash right)^\{1/p\}$$

其中 $\mathrm{\Gamma }(P,Q)\backslash Gamma(P,\; Q)$ 是 $PP$ 和 $QQ$ 中点之间所有可能匹配的集合，$\mathrm{\parallel }\cdot {\mathrm{\parallel }}_q\backslash |\; \backslash cdot\; \backslash |\_q$ 是 $qq$-范数。

# Helper function to compute the Wasserstein distances between all pairs of persistence diagrams
def compute_wasserstein_distances(persistence_diagrams, dimension):
    n_diagrams = len(persistence_diagrams)
    distances = np.zeros((n_diagrams, n_diagrams))
    filtered_diagrams = [[point for point in diagram if point[0] == dimension] for diagram in persistence_diagrams]
    for i in range(n_diagrams):
        for j in range(i+1, n_diagrams):
            X = np.array([p[1] for p in filtered_diagrams[i]])
            Y = np.array([p[1] for p in filtered_diagrams[j]])
            distance = wasserstein_distance(X, Y)
            distances[i][j] = distance
            distances[j][i] = distance
    return distances

加载 ESM-2 模型和分词器

接下来，我们加载要使用的 ESM-2 模型。这里我们使用 Hugging Face 的 facebook/esm2_t6_8M_UR50D。

# Load the tokenizer and model
tokenizer = AutoTokenizer.from_pretrained("facebook/esm2_t6_8M_UR50D")
model = EsmModel.from_pretrained("facebook/esm2_t6_8M_UR50D")

定义用于嵌入的层

现在，我们选择模型的一个层来获取隐藏状态（嵌入向量）。这可以根据 蛋白质语言模型隐藏表示的几何结构 和 大型 Transformer 模型隐藏表示的几何结构 中的结果进行选择。

# Define layer to be used
num_layers = model.config.num_hidden_layers
layer = num_layers - 1  # Index of the last layer

计算单个蛋白质的持久图

我们可以按如下方式计算单个蛋白质序列的持久图：

# Define a sample protein sequence
input_sequence = "GLSDGEWQQVLNVWGKVEADIPGHGQEVLIRLFKGHPETLEKFDKFKHLKSEDEMKASEDLKKHGATVLTALGGILKKKGHHEAELKPLAQSHATKHKIPIKYLEFISEAIIHVLHSRHPGDFGADAQGAMNKALELFRKDIAAKYKELGYQG"

# Compute the hidden states matrix
hidden_states_matrix = get_hidden_states(tokenizer, model, layer, input_sequence)

# Compute the distance matrix
distance_matrix = compute_euclidean_distance_matrix_scipy(hidden_states_matrix)

# Compute the persistent homology
# max_dimension = 4 may take several minutes
# max_dimension = bigger number means the computation will take longer
st, persistence = compute_persistent_homology(distance_matrix, max_dimension=3)

gd.plot_persistence_diagram(persistence)

这将返回一个持久图的绘制。

使用持久同源聚类蛋白质

接下来，我们定义一个要使用持久同源聚类的蛋白质列表。我们计算每个蛋白质的持久图，然后计算每对持久图之间的 Wasserstein 距离。最后，我们运行一个持久同源引导的 DBSCAN，使用我们刚刚计算的 Wasserstein 距离矩阵对持久图进行聚类。这为我们提供了基于蛋白质持久同源（即基于它们的持久图）的语义丰富且有意义的蛋白质聚类。请注意，我们也可以使用条形码代替持久图，因为两者之间存在双射，表明它们是等效的。因此，如果您更喜欢条形码，可以将其替换为持久图。

sequences = [

    "MAHMTQIPLSSVKRPLQVRGIVFLGICTQKGTVYGNASWVRDQARH",
     "MKHVTQIPKSSVRRPLQFRGICFLGTCQKGTVYGKASWVHDQARHA",
     "MNHITQVPLSSAKRPLQVRGICFLGITCKNGTVYGKACWVRDQARH",

    "MKLITILGLLALATLVQSTGCVTVNAAHCGVTTGQTVCAGVAKCRAE",
     "MKLITILGALALATLVQSTGCVNVNAAHCVTTGQTVCAGVAKCRAET",
     "MKLITILGALALATLVQSTGCVNVNAAHCVTAGQTVCAGVAKCRAETS",

    "MGSSHHHHHHSSGLVPRGSHMENITVVKFNGTQTFEVHPNVSVGQAGV",
     "MGSSHHHHHHSSGLVPRGSHMENITVVKFNGTQTFEVHPNVSVGQAGVR",
     "MGSSHHHHHHSSGLVPRGSHMENITVVKFNGTQTFEVHPNVSVGQAGVRR",

    "MGGHNGWQILVKGKWTTMDFLRNAVIDQKLRRARRELKLMKAFESLK",
     "MGGHNGWQILVKGKWTTMDFLRNAVIDQKLRRARRELKLMKAFESLKN",
     "MGGHNGWQILVKGKWTTMDFLRNAVIDQKLRRARRELKLMKAFESLKNN",

    "MAQSNISDAMVQLTPAGRSLMLLVQHGSQVAAGVTFQDNQRFPGGRD",
     "MAQSNISDAMVQLTPAGRSLMLLVQHGSQVAAGVTFQDNQRFPGGRDF",
     "MAQSNISDAMVQLTPAGRSLMLLVQHGSQVAAGVTFQDNQRFPGGRDFF"
]


# Initialize list to store persistent diagrams
persistent_diagrams = []

# Compute persistent homology for each sequence
for sequence in sequences:
    hidden_states_matrix = get_hidden_states(tokenizer, model, layer, sequence)
    distance_matrix = compute_euclidean_distance_matrix_scipy(hidden_states_matrix)
    _, persistence_diagram = compute_persistent_homology(distance_matrix)
    
    # Store the persistent diagram
    persistent_diagrams.append(persistence_diagram)

在这里，我们需要选择持久同源要使用的维度。我们可以选择维度 0 来根据零维持久同源特征进行聚类。我们也可以选择更高的维度来根据更高维的拓扑特征进行聚类。该方法尚未经过广泛测试以确定哪种方法最好，因此我们暂时将此选择留给您。这里我们选择维度零以保持可解释性和简单性。

# Compute the Wasserstein distances between all pairs of persistence diagrams
wasserstein_distances = compute_wasserstein_distances(persistent_diagrams, 0)
# Compute the persistent homology of the Wasserstein distance matrix
st_2, persistence_2 = compute_persistent_homology(wasserstein_distances)
# Plot the persistence diagram
gd.plot_persistence_diagram(persistence_2)

这将输出另一个持久图，我们从中为 DBSCAN 选择 epsilon。特别是，我们在红点之间找到大的间隙，并选择一个落入这些间隙中的 epsilon，直到我们为 DBSCAN 获得高轮廓系数。

DBSCAN 聚类

最后，应用 DBSCAN 聚类，根据 Wasserstein 距离对持久图进行聚类。DBSCAN（基于密度的带噪声空间聚类应用）是一种算法，它将紧密聚集的点分组，同时将低密度区域中的点标记为离群值。使用轮廓系数来评估聚类的质量。在这里，我们根据上面的持久图选择 eps=40.0。

# Perform DBSCAN clustering of persistence diagrams
dbscan = DBSCAN(metric="precomputed", eps=40.0, min_samples=1).fit(wasserstein_distances)
labels_dbscan = dbscan.labels_
if len(set(labels_dbscan)) > 1:  # More than 1 cluster
    print("Silhouette Coefficient for DBSCAN: %0.3f" % silhouette_score(wasserstein_distances, labels_dbscan))
else:
    print("Cannot compute Silhouette Coefficient for DBSCAN as there is only one cluster.")

# Print the clusters for DBSCAN
print("\nClusters for DBSCAN:")
dbscan_clusters = {i: [] for i in set(labels_dbscan)}
for sequence, label in zip(sequences, labels_dbscan):
    dbscan_clusters[label].append(sequence)
for label, cluster in dbscan_clusters.items():
    print(f"Cluster {label}: {cluster}")

Clusters for DBSCAN:
Cluster 0: ['MAHMTQIPLSSVKRPLQVRGIVFLGICTQKGTVYGNASWVRDQARH', 'MNHITQVPLSSAKRPLQVRGICFLGITCKNGTVYGKACWVRDQARH']
Cluster 1: ['MKHVTQIPKSSVRRPLQFRGICFLGTCQKGTVYGKASWVHDQARHA']
Cluster 2: ['MKLITILGLLALATLVQSTGCVTVNAAHCGVTTGQTVCAGVAKCRAE']
Cluster 3: ['MKLITILGALALATLVQSTGCVNVNAAHCVTTGQTVCAGVAKCRAET', 'MKLITILGALALATLVQSTGCVNVNAAHCVTAGQTVCAGVAKCRAETS']
Cluster 4: ['MGSSHHHHHHSSGLVPRGSHMENITVVKFNGTQTFEVHPNVSVGQAGV', 'MGSSHHHHHHSSGLVPRGSHMENITVVKFNGTQTFEVHPNVSVGQAGVR', 'MGSSHHHHHHSSGLVPRGSHMENITVVKFNGTQTFEVHPNVSVGQAGVRR']
Cluster 5: ['MGGHNGWQILVKGKWTTMDFLRNAVIDQKLRRARRELKLMKAFESLK', 'MGGHNGWQILVKGKWTTMDFLRNAVIDQKLRRARRELKLMKAFESLKN', 'MGGHNGWQILVKGKWTTMDFLRNAVIDQKLRRARRELKLMKAFESLKNN', 'MAQSNISDAMVQLTPAGRSLMLLVQHGSQVAAGVTFQDNQRFPGGRD', 'MAQSNISDAMVQLTPAGRSLMLLVQHGSQVAAGVTFQDNQRFPGGRDF', 'MAQSNISDAMVQLTPAGRSLMLLVQHGSQVAAGVTFQDNQRFPGGRDFF']

现在我们已经有了蛋白质的簇，我们可以使用 利用蛋白质语言模型进行精确多序列比对 和 向量聚类多序列比对：使用蛋白质语言模型比对进入蛋白质序列相似性的“微光区” 中描述的其余方法来比对蛋白质序列。特别是，到目前为止，我们已经用基于持久同源的蛋白质聚类替换了下图中步骤 (A) 和 (B)，这有望表现得更好。

我们注意到，此时我们可以从簇中删除离群值，就像在 vcMSA 中所做的那样，或者我们可以使用 Fréchet 均值持久图来确定离群值。有关持久图的 Fréchet 均值（重心）的更多信息，请参阅此处，以及提供的示例笔记本。其余的实现可以按照论文的 Github 仓库进行。例如，我们现在可以使用以下函数实现步骤 (3) 以在每个簇中找到互惠最佳匹配：

from scipy.spatial.distance import cosine

def identify_rbh(hidden_states_matrices, threshold=0.8):
    """
    Identifies Reciprocal Best Hits (RBHs) based on cosine similarity.
    :param hidden_states_matrices: List of hidden states matrices for each protein.
    :param threshold: Cosine similarity threshold for considering a hit.
    :return: Dictionary of RBH pairs for each protein pair.
    """
    num_proteins = len(hidden_states_matrices)
    rbh_pairs = {}

    for i in range(num_proteins):
        for j in range(i + 1, num_proteins):
            max_similarity = {}
            protein_i, protein_j = hidden_states_matrices[i], hidden_states_matrices[j]

            # Compute cosine similarity between all pairs of amino acids
            for idx_i, amino_acid_i in enumerate(protein_i):
                for idx_j, amino_acid_j in enumerate(protein_j):
                    similarity = 1 - cosine(amino_acid_i, amino_acid_j)
                    if similarity > threshold:
                        max_similarity[(idx_i, i)] = max(similarity, max_similarity.get((idx_i, i), 0))
                        max_similarity[(idx_j, j)] = max(similarity, max_similarity.get((idx_j, j), 0))

            # Identify RBHs
            rbh_pairs[(i, j)] = [(idx_i, idx_j) for (idx_i, _), sim_i in max_similarity.items() 
                                 for (idx_j, _), sim_j in max_similarity.items() 
                                 if sim_i == sim_j and sim_i > threshold]

    return rbh_pairs

这将提供每个簇内氨基酸嵌入向量之间的余弦相似度，并允许我们聚类具有高相似度的氨基酸。

结论

这种新方法（PHA）将取代传统的多序列比对（MSA），并且将对序列相似性低或长插入/缺失的蛋白质序列更具鲁棒性和适用性。我们注意到，这避免了对 MSA 算法的需求，从而避免了初始引导树构建、中间成对比对、间隙罚分和替换矩阵。此外，上下文嵌入提供的额外信息使得对结构相似但氨基酸相似性低的蛋白质比对具有更高的准确性。这标志着在用更鲁棒的方法实现相同目标方面，取代 MSA 迈出了重要一步。

现在出现一个问题，如果像 ESM-2 和 ESMFold 这样的蛋白质语言模型（它们在没有 MSA 的单个序列上进行训练），结合像 vcMSA 或 PHA 这样的几何方法可以取代通常的 MSA，为什么我们不尝试将其构建到我们的深度学习模型中，以完全消除对 MSA 的需求呢？这是 AlphaFold2 比 ESMFold 慢得多的主要原因，并且在计算和时间方面是一个巨大的瓶颈。虽然 MSA 通常在某些任务上提供更高的准确性，但对于许多蛋白质来说，它们很难或不可能构建。此外，对于涉及蛋白质复合物与小分子、DNA 或 RNA 相互作用的更复杂结构，它们甚至更难构建。因此，拥有一个使用替代方法（如 PHA 或 vcMSA）的模型，或者完全规避对它们的需求的模型，将非常有益。