使用ESM-2和持久化景观进行更快的持久同调比对和蛋白质复合物聚类

在本文中，我们将讨论一种更快、计算效率更高的方法，用于计算蛋白质的持久同调比对（PHA），它可以取代多序列比对，适用于黄昏区蛋白质和孤立蛋白质。我们还将讨论如何将其用于蛋白质序列相似性和蛋白质复合物的同源建模。这两种方法都使用了蛋白质语言模型ESM-2。

黄昏区蛋白

"黄昏区蛋白"是生物信息学和分子生物学中使用的术语，用于描述相互之间序列同一性非常低的蛋白质序列，通常低于20-35%。这个术语来源于蛋白质序列比对的语境，指的是序列相似性极低，以至于仅凭序列比较难以确定进化关系或结构相似性的区域。

序列比对挑战：传统的序列比对算法可能难以生成准确的比对结果，因为共同特征稀疏且信噪比低。低序列同一性意味着序列已显著分化，使得识别同源区域（源自共同祖先序列的区域）变得困难。

结构保守性：尽管序列同一性低，黄昏区蛋白仍可能共享相似的三维结构和功能。这是因为蛋白质结构通常比其序列更保守。在序列相似性低的情况下结构保守性可以为关键功能域和进化过程提供见解。

进化研究的意义：理解黄昏区蛋白之间的关系对进化生物学至关重要，因为它有助于追踪物种的分化和蛋白质新功能的开发。

高级分析技术：为了分析黄昏区蛋白，科学家通常依赖于超越简单序列比较的先进方法。这些方法可能包括基于结构的比对技术和机器学习方法，它们可以考虑更复杂的模式和关系，而不仅仅是直接的序列相似性。我们将用持久同调比对（PHA）取代MSA。

"黄昏区"一词突出了分析远缘序列的复杂性和模糊性。它强调了整合多种生物信息来源（包括结构、功能和进化数据）以全面理解这些蛋白质的重要性。

持久化景观

持久化景观：数学概述

持久化景观是拓扑数据分析（TDA）中的一个关键工具，它提供了持久同调的一种稳定且易于解释的表示形式，持久同调是一种从数据中提取多尺度拓扑特征的方法。本节将深入探讨持久化景观的数学框架，阐明其构造、性质以及在TDA中的实用性。

拓扑数据分析（TDA）是数据科学中一个相对较新的领域，专注于数据形状（拓扑）的研究。TDA旨在揭示数据集中固有的几何和拓扑结构，这些结构在传统分析方法下通常是隐藏的。TDA中的一个关键概念是持久同调，它跟踪拓扑特征（如连通分量、孔、空隙）在多个尺度上的演变。

持久同调量化数据集的多尺度拓扑特征。给定点云数据，我们构建一个单纯复形的过滤，通常通过Vietoris-Rips或Čech复形，从而得到一系列嵌套的拓扑空间。随着尺度参数的增加，新的拓扑特征出现，现有特征可能合并或消失。每个特征由其生成时间和消亡时间表征，它们分别是特征出现和消失的尺度参数。形式上，对于过滤参数$tt$，持久对表示为$(b,d)(b,d)$，其中$bb$是拓扑特征的生成时间，$dd$是消亡时间。持久化图是平面中此类点的一个多重集，通常以散点图的形式可视化。

持久化景观是持久化图的功能性表示。它将点的多重集转换为一系列连续的、分段线性的函数，这些函数易于统计分析。给定一个持久化图，其中包含点$\{(b_i,d_i)\}\backslash \{(b\_i,\; d\_i)\backslash \}$，持久化景观由一系列函数$\{{\lambda }_k\}\backslash \{\backslash lambda\_k\backslash \}$组成，其中每个${\lambda }_k:\mathbb{R}\to {\mathbb{R}}_{\ge 0}\backslash lambda\_k:\; \backslash mathbb\{R\}\; \backslash to\; \backslash mathbb\{R\}\_\{\backslash geq\; 0\}$定义如下：

$${\lambda }_k(t)=\sup \{\min (t-b_i,d_i-t):\text{such that }t\in [b_i,d_i]\text{ and }(b_i,d_i)\text{ is the }k\text{-th largest point}\}\backslash lambda\_k(t)\; =\; \backslash sup\backslash \{\backslash min(t-b\_i,\; d\_i-t)\; :\; \backslash text\{such\; that\; \}\; t\; \backslash in\; [b\_i,\; d\_i]\; \backslash text\{\; and\; \}\; (b\_i,\; d\_i)\; \backslash text\{\; is\; the\; \}\; k\backslash text\{-th\; largest\; point\}\backslash \}$$

本质上，对于持久化图中的每个点，都会构建一个“帐篷”函数，其中心位于其持久区间的中间点，高度等于持久度（消亡时间与生成时间之差）的一半。第$kk$个景观函数${\lambda }_k\backslash lambda\_k$是这些帐篷函数在每个时间$tt$处的第$kk$个最大值。

• 稳定性：持久化景观在输入数据的微小扰动下是稳定的，这使得它们成为稳健的分析特征。
• 向量化：每个${\lambda }_k\backslash lambda\_k$可以离散化并视为欧几里德空间中的向量，这有助于应用标准统计和机器学习技术。
• 组合结构：分段线性的性质简化了计算和解释。

持久化景观已成功应用于各种领域，包括形状分析、信号处理和生物数据分析。它们提供了拓扑特征的紧凑和稳定摘要，允许根据其拓扑特征有效比较、聚类和分类数据集。持久化景观提供了TDA中强大而通用的工具。通过将复杂的拓扑信息转换为更易于理解和分析的形式，它们弥合了抽象拓扑概念与实际数据分析之间的鸿沟，为深入了解复杂数据集的内在结构铺平了道路。

下面我们有一个脚本，用于生成三个包含25个点的随机持久化图。然后我们为每个持久化图计算持久化景观，并计算它们的成对$L_{2}L\_2$-距离以获得距离矩阵。从这个距离矩阵，我们计算一个二级持久化图。

import numpy as np
import matplotlib.pyplot as plt
from gudhi.representations import Landscape

# Number of diagrams and points
num_diagrams = 3
num_points = 25

# Generate three random persistence diagrams each with 25 points
# Assuming birth and death times are uniformly distributed between 0 and 10
diagrams = [np.random.uniform(0, 10, (num_points, 2)) for _ in range(num_diagrams)]

# Ensure that death time is greater than birth time
for d in diagrams:
    d[:, 1] = np.maximum(d[:, 1], d[:, 0] + np.random.uniform(0, 1, num_points))

# Compute landscapes for each diagram
landscape = Landscape(num_landscapes=3, resolution=100)
landscapes = landscape.fit_transform(diagrams)

# Plot the landscapes
plt.figure(figsize=(15, 5))
for i, l in enumerate(landscapes, start=1):
    plt.subplot(1, num_diagrams, i)
    plt.plot(l[:100], label="Landscape 1")
    plt.plot(l[100:200], label="Landscape 2")
    plt.plot(l[200:], label="Landscape 3")
    plt.title(f"Diagram {i} Landscapes")
    plt.xlabel("Sample Points")
    plt.ylabel("Amplitude")
    plt.legend()

plt.tight_layout()
plt.show()

景观应该类似于以下景观

成对的$L_{2}L\_2$-距离矩阵将看起来像这样

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.spatial.distance import cdist
import gudhi as gd
from gudhi import RipsComplex
from gudhi.representations import DiagramSelector

# Assuming 'landscapes' variable contains the landscapes computed previously
# If 'landscapes' is not available, we can use a placeholder or regenerate them

# Compute L2 distances between each pair of landscapes
dist_matrix = cdist(landscapes, landscapes, metric='euclidean')

# Print the distance matrix as a heatmap
plt.figure(figsize=(6, 6))
sns.heatmap(dist_matrix, annot=True, cmap='viridis')
plt.title("L2 Distance Matrix Heatmap")
plt.xlabel("Landscape Index")
plt.ylabel("Landscape Index")
plt.show()

# Compute the persistence diagram of the distance matrix with Gudhi's Rips complex
rips_complex = RipsComplex(distance_matrix=dist_matrix)
simplex_tree = rips_complex.create_simplex_tree(max_dimension=2)
diag = simplex_tree.persistence()

# Plot the persistence diagram
plt.figure()
gd.plot_persistence_diagram(diag)
plt.title("Persistence Diagram of the Rips Complex")
plt.xlabel("Birth")
plt.ylabel("Death")
plt.show()

一旦我们得到了$L_{2}L\_2$-距离矩阵，我们就会得到一个二级持久化图，它看起来像这样

现在我们可以使用这个持久化图，通过选择一个距离阈值来捕获二级持久化图中的部分红点，对持久化景观执行DBSCAN聚类。在下一节中，我们将对蛋白质-蛋白质复合物进行此操作，并选择一个阈值，该阈值捕获二级持久化图中80%的点。换句话说，我们为DBSCAN选择一个距离阈值$ϵ\backslash epsilon$，使得图中80%的点落在此阈值以下。

蛋白质序列和蛋白质复合物聚类

如果您想运行以下脚本，请首先下载此数据集，其中包含相互作用的人类蛋白质对（最初从UniProt获取），或者前往UniProt自行管理相互作用蛋白质对的数据集。请务必在运行脚本之前调整文件路径。此脚本将使用上述方法，利用蛋白质语言模型ESM-2对相互作用的蛋白质对进行聚类。与此帖子中描述的方法相比，这在速度方面有了显著改进。特别是，计算1000个蛋白质-蛋白质复合物对之间的景观$L_{2}L\_2$-距离大约需要30分钟，而计算持久化图对之间的Wasserstein距离需要一整天，或者瓶颈距离需要数小时。

import pandas as pd
import numpy as np
from transformers import EsmModel, AutoTokenizer
import torch
from scipy.spatial.distance import pdist, squareform
from gudhi import RipsComplex
from gudhi.representations.vector_methods import Landscape
from sklearn.cluster import DBSCAN
# import matplotlib.pyplot as plt
from tqdm import tqdm

# Define a helper function for hidden states
def get_hidden_states(sequence, tokenizer, model, layer):
    model.config.output_hidden_states = True
    encoded_input = tokenizer([sequence], return_tensors='pt', padding=True, truncation=True, max_length=1024)
    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.numpy()

# Define a helper function for Euclidean distance matrix
def compute_euclidean_distance_matrix(hidden_states):
    euclidean_distances = pdist(hidden_states, metric='euclidean')
    euclidean_distance_matrix = squareform(euclidean_distances)
    return euclidean_distance_matrix

# Define a helper function for persistent homology
def compute_persistent_homology(distance_matrix, max_dimension=0):
    max_edge_length = np.max(distance_matrix)
    rips_complex = RipsComplex(distance_matrix=distance_matrix, max_edge_length=max_edge_length)
    st = rips_complex.create_simplex_tree(max_dimension=max_dimension)
    st.persistence()
    persistence_pairs = np.array([[p[1][0], p[1][1]] for p in st.persistence() if p[0] == 0 and p[1][1] < np.inf])  # Filter out infinite death times
    return st, persistence_pairs

# Define a helper function for persistent homology
def compute_persistent_homology2(distance_matrix, max_dimension=0):
    max_edge_length = np.max(distance_matrix)
    rips_complex = RipsComplex(distance_matrix=distance_matrix, max_edge_length=max_edge_length)
    st = rips_complex.create_simplex_tree(max_dimension=max_dimension)
    st.persistence()
    return st, st.persistence()

# Define a helper function for Landscape transformations with tqdm
def compute_landscapes(persistence_diagrams, num_landscapes=5, resolution=10000):
    landscape_transformer = Landscape(num_landscapes=num_landscapes, resolution=resolution)
    landscapes = landscape_transformer.fit_transform([d for d in persistence_diagrams if len(d) > 0])  # Filter out empty diagrams
    return landscapes

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

# Define layer to be used
layer = model.config.num_hidden_layers - 1

# Load the TSV file
file_path = 'filtered_protein_interaction_pairs.tsv'
protein_pairs_df = pd.read_csv(file_path, sep='\t')

# Only process the first 1000 proteins
protein_pairs_df = protein_pairs_df.head(1000)

# Extract concatenated sequences
concatenated_sequences = protein_pairs_df['Protein1'] + protein_pairs_df['Protein2']

# Initialize list to store persistent diagrams
persistent_diagrams = []

# Loop over concatenated sequences to compute their persistent diagrams
for sequence in tqdm(concatenated_sequences, desc="Computing Persistence Diagrams"):
    hidden_states_matrix = get_hidden_states(sequence, tokenizer, model, layer)
    distance_matrix = compute_euclidean_distance_matrix(hidden_states_matrix)
    _, persistence_diagram = compute_persistent_homology(distance_matrix)
    persistent_diagrams.append(persistence_diagram)

# Compute landscapes
landscapes = compute_landscapes(persistent_diagrams)

# Compute the L2 distances between landscapes
with tqdm(total=len(landscapes)*(len(landscapes)-1)//2, desc="Computing Pairwise L2 Distances") as pbar:
    l2_distances = np.zeros((len(landscapes), len(landscapes)))
    for i in range(len(landscapes)):
        for j in range(i+1, len(landscapes)):
            l2_distances[i, j] = l2_distances[j, i] = np.linalg.norm(landscapes[i] - landscapes[j])
            pbar.update(1)

# Compute the second-level persistent homology using the L2 distance matrix
with tqdm(total=1, desc="Computing Second-Level Persistent Homology") as pbar:
    st_2, persistence_2 = compute_persistent_homology2(l2_distances)
    pbar.update(1)

# Function to calculate the epsilon for DBSCAN
def calculate_epsilon(persistence_diagrams, threshold_percentage, max_eps=np.inf):
    lifetimes = [p[1][1] - p[1][0] for p in persistence_diagrams if p[0] == 0]
    lifetimes.sort()
    threshold_index = int(threshold_percentage * len(lifetimes))
    epsilon = lifetimes[threshold_index]
    # Ensure epsilon is within a reasonable range
    epsilon = min(epsilon, max_eps)
    return epsilon

# Calculate epsilon with a maximum threshold
threshold_percentage = 0.8  # 80%
max_epsilon = 5000.0  # Example maximum threshold
epsilon = calculate_epsilon(persistence_2, threshold_percentage, max_eps=max_epsilon)

# Perform DBSCAN clustering
with tqdm(total=1, desc="Performing DBSCAN Clustering") as pbar:
    dbscan = DBSCAN(metric="precomputed", eps=epsilon, min_samples=1)
    dbscan.fit(l2_distances)  # Use L2 distances here
    labels = dbscan.labels_
    pbar.update(1)

# Add the cluster labels to the DataFrame
protein_pairs_df['Cluster'] = labels

# Save the DataFrame with cluster information
output_file_path = 'clustered_protein_interaction_pairs_l2_distances.tsv'
protein_pairs_df.to_csv(output_file_path, sep='\t', index=False)

print(f"Clustered data saved to: {output_file_path}")

我们注意到，调整`num_landscapes`参数和`resolution`为适当值仍在实验阶段，但这将帮助您开始使用ESM-2和持久同调对蛋白质复合物进行聚类的旅程。

用PHA取代MSA

接下来，我们将讨论用持久同调比对（PHA）取代多序列比对（MSA）的可能性，这参照了这篇博客文章。在我们的设置中，我们现在有了一种方法可以为任何蛋白质构建持久化景观，这与上述过程非常相似，但现在是针对单个蛋白质而不是连接的蛋白质对。我们还有一种方法可以根据持久化景观之间的$L_{2}L\_2$-距离来构建二级持久化图。从此，我们像以前一样执行DBSCAN，然后继续使用这篇博客文章中描述的vcMSA方法。

下一步，您可以尝试使用其他方法对持久化图进行向量化，而不是持久化景观。例如，您可以尝试持久化图像，或者Gudhi 文档中解释的其他方法。您还可以尝试使用$L_{2}L\_2$-范数之外的其他相似性度量。

结论

总而言之，本文介绍了一种创新方法，利用持久同调比对（PHA）和先进的蛋白质语言模型ESM-2，分析蛋白质序列，特别是那些处于“黄昏区”的蛋白质。这里讨论的方法在计算效率和处理低序列同源性蛋白质的能力方面，比传统的多序列比对有了显著改进。

持久化景观在拓扑数据分析中的应用已被证明是提取蛋白质数据中有意义拓扑特征的强大工具。这个数学框架使得能够处理复杂的数据结构，并提供了一种稳定的方法来根据蛋白质的内在拓扑特性进行比较和聚类。此外，将ESM-2集成到蛋白质-蛋白质复合物的聚类中，标志着生物信息学领域的重大进展。这种方法不仅加快了计算过程，而且提高了蛋白质相互作用预测的准确性。PHA与vcMSA方法的结合，进一步展示了这些技术在取代传统序列比对方法方面的潜力，特别是对于序列模糊或独特的蛋白质。

总的来说，计算方法的进步以及ESM-2和PHA等复杂模型在蛋白质序列分析中的整合，代表了生物信息学领域的重大飞跃。这些方法为理解蛋白质、其相互作用及其进化关系的复杂世界开辟了新的途径。它们还为未来的研究铺平了道路，其中进一步的改进和创新可能导致更高效、更准确的蛋白质分析方法。