BioSemAF-BiLSTM：融合语义与进化信息的蛋白质S-亚磺酰化位点深度预测框架

本综述提出BioSemAF-BiLSTM深度学习框架，整合fastText子词嵌入与位置特异性评分矩阵（PSSM）特征，通过双向长短期记忆网络（BiLSTM）和自适应特征融合模块（AF）提升S-亚磺酰化位点预测性能。该模型在独立测试集上达到89.32%准确率，显著优于传统机器学习方法，为蛋白质翻译后修饰（PTM）研究提供强大计算工具。

　　

1 引言

S-亚磺酰化是一种关键的蛋白质翻译后修饰（PTM），在调节蛋白质功能、氧化还原信号传导和维持细胞稳态中发挥重要作用。准确识别S-亚磺酰化位点对于理解其生物学意义及与疾病的相关性至关重要。然而，通过实验方法专门检测S-亚磺酰化位点仍然具有挑战性，因为这些方法通常耗时且成本高昂。受此问题启发，本研究提出了一个基于深度学习的计算框架，名为BioSemAF-BiLSTM，该框架整合了进化特征和语义特征，以提高S-亚磺酰化位点的预测性能。

传统机器学习方法在S-亚磺酰化位点预测的早期发展中发挥了核心作用。这些方法通常属于蛋白质序列位点预测的更广泛范畴——生物信息学的一个重要领域，旨在识别功能残基。准确的序列表示对于生物信息学中的各种相关预测任务也至关重要，例如确定蛋白质二级结构。然而，这些方法严重依赖复杂的特征工程和先验知识，对大规模数据集的适应性有限，这为深度学习方法的应用铺平了道路。

近年来，几种基于深度学习的模型被开发用于S-亚磺酰化位点的预测，每种模型采用不同的序列表示和特征提取策略。例如，nSSPred结合了n段优化与联合特征编码器，并应用二维卷积神经网络（2D-CNN）进行预测；fastSulf-DNN提取生物子词级特征，并采用自然语言处理技术编码蛋白质序列，用于深度神经网络（DNN）分类；DLF-Sul是一个多模块框架，整合了二进制编码、BLOSUM62和氨基酸指数，然后使用BiLSTM进行序列建模，并进一步应用多头自注意力和CNN进行特征细化，最后通过全连接层进行分类。然而，这些方法存在关键局限性：对专家手工制作特征的强烈依赖限制了生物信息的全面捕获；语义理解能力有限，妨碍了氨基酸功能和相互依赖关系的准确识别；对异构特征的适应能力较弱，限制了动态的、任务感知的特征加权；对局部模式的关注削弱了序列中长程依赖关系的建模。这些问题共同限制了现有方法对S-亚磺酰化位点的预测性能。

2 材料与方法

2.1 基准数据集

本研究使用了iSulf-Cys框架中开发的数据集，该数据集基于Yang等人（2014）的工作。数据集包含7124个非亚磺酰化半胱氨酸位点和1045个亚磺酰化半胱氨酸位点。数据总共源自778个蛋白质序列，包括1105个亚磺酰化半胱氨酸位点。对于每个分析的半胱氨酸位点，Xu等人通过提取上游和下游各10个氨基酸来构建肽段样本，得到一个以目标半胱氨酸为中心的21个氨基酸大小的序列窗口。如果周围序列残基少于10个，则使用占位残基X进行填充。样本的表示如公式（1）所示：Peptide = A_-10A_-9…CA₁A₂A₃…A₁₀，其中A_i ∈ {20种氨基酸} ∪ {X}，i的范围从-10到10，表示相对于中心半胱氨酸的一维坐标。以中心半胱氨酸为亚磺酰化位点的肽段被定义为正样本，所有其他肽段被分类为负样本。为了最小化冗余并减少同源性偏差，去除了序列相似性超过40%的肽段。最终，随机选择145个正样本和268个负样本作为测试集，其余900个正样本和6858个负样本分配给训练集。

2.2 特征提取

2.2.1 PSSM

位置特异性评分矩阵（PSSM）是一种生物信息学表示，用于描述生物序列中的进化保守性。进化上保守的氨基酸残基通常与在蛋白质功能中起关键作用的功能 motif 相关，包括翻译后修饰（PTM）。因此，PSSM作为一种有价值的生物信息学特征，用于表示预测模型中的肽段。PSSM描述和量化了序列中每个位置氨基酸（或核苷酸）的分布，反映了进化保守性和突变倾向。它在各种位点预测任务中得到了广泛应用。例如，Yuan等人（2021）在表征蛋白质-蛋白质相互作用（PPI）节点时引入了PSSM作为关键特征。在PTM位点预测中，PSSM也显示出巨大的实用性。PSSM-Suc使用PSSM编码琥珀酰化位点进行预测建模，而PSSM-Sumo使用称为PsePSSM的变体编码SUMO化位点。这两种方法都展示了强大的预测性能，进一步验证了基于PSSM表示的有效性。

标准PSSM的矩阵维度为20行（每行对应20种标准氨基酸之一）和列数等于肽段或序列的长度。PSSM的规范格式定义为公式（2）：PSSM = [a₁₁ a₁₂ … a_1L; a₂₁ a₂₂ … a_2L; …; a_20,1 a_20,2 … a_20,L]，其中a_ij表示第i种氨基酸在第j个位置突变为另一种氨基酸的概率，i的范围从1到20（代表氨基酸类型），j的范围从1到L（代表序列长度，L是序列的总长度）。

在本研究中，PSSM是使用位置特异性迭代BLAST（PSI-BLAST）针对NCBI非冗余（NR）蛋白质数据库生成的。按照常见做法，我们使用了三次迭代，期望值（E-value）为0.001。需要注意的是，PSI-BLAST仅用于执行迭代比对以生成PSSM，而不是用于从头同源序列搜索。其他多序列比对工具如HHblits或Clustal Omega也可用于此目的，但PSI-BLAST在相关工作中仍然被广泛使用。

2.2.2 词嵌入

语言是最直观的表示形式之一，人类长期以来一直试图通过语言表达来描述自然世界。通过将特定字母分配给氨基酸，蛋白质序列可以直接编码为字符串。这些序列包含丰富的语义内容，可以使用适当的自然语言处理（NLP）技术进行挖掘。本节重点介绍局部蛋白质相关任务，特别是从序列数据预测残基级特征。基于窗口的模型在各种PTM位点预测任务中已被证明是有效的。例如，Brandes等人（2015）引入了一个框架，通过将蛋白质序列分割成固定大小的窗口来预测PTM位点，突出了语义特征在蛋白质特征工程中的潜力。

随着NLP技术的不断进步，Word2Vec已成为该领域的重大突破。通过利用上下文信息，它为单词（或字符）生成连续的低维向量表示，捕获语言中嵌入的潜在语义关系。这一进展激发了蛋白质序列分析的新方法。FastText是Word2Vec的增强版本，进一步扩展了这一思想。与Word2Vec相比，fastText引入了子词建模，有效捕获单词的内部形态和结构特征。这特别适用于蛋白质序列建模，因为此类序列经常包含可变和局部化的 motif。FastText不仅在训练期间捕获局部n-gram类信息，还解决了词汇表外（OOV）问题——这在分析复杂生物序列时是一个关键优势。

fastText的核心思想是通过引入子词建模来增强传统的词嵌入。与为每个单词分配一个独特向量的Word2Vec模型不同，fastText将一个单词表示为其组成子词的组合，从而捕获单词内部的细微变化。例如，单词“protein”可以拆分为子词，如“pro”、“ote”和“ein”，所有这些都参与模型的训练。通过对这些子词进行编码和嵌入，fastText学习了更细致的语义表示，使其在处理多样化和不规则的蛋白质序列时特别有效。

从数学角度来看，fastText建立在Skip-gram架构之上，其目标是最大化给定周围上下文预测目标单词的条件概率。给定一个由单词序列w₁, w₂, …, w_T组成的训练语料库，fastText的目标函数是通过最小化以下损失函数来学习单词及其子词的嵌入：L = -∑_t=1^T ∑_{-c≤j≤c,j≠0} log p(w_t+j|w_t)（公式3），其中w_t是当前单词，w_t+j是上下文单词，c是上下文窗口的大小，p(w_t+j|w_t)是预测目标单词w_t+j的条件概率。在fastText中，p(w_{t+jt)不再是简单单词的条件概率，而是子词的条件概率。具体来说，对于每个单词w，它被拆分为几个子词s1, s2, …, sk，目标单词的概率计算为公式（4）：p(wt+j|wt) = ∏s∈S(wt+j) p(s|wt)，其中S(w)表示单词w的子词集合，p(s|wt)表示给定上下文单词wt的子词s的条件概率。}

在实际应用中，fastText通过其自身的嵌入及其组成子词的嵌入来表示单词的语义内容。对于蛋白质序列，由氨基酸残基衍生的二肽或三肽可以被视为子词。使用这种表示，fastText可以将蛋白质序列的短片段（例如氨基酸对或三联体）转换为向量嵌入，捕获子结构之间的潜在上下文关系。在本研究中，每个肽序列被分解为3-g作为子词。对于一个由n个氨基酸组成的序列，这产生n-2个三连词，因此生成n-2个相应的词向量，这些词向量嵌入了序列的隐含顺序和语义内容。这个过程构成了本研究中使用的基于fastText的蛋白质序列嵌入。

2.3 重采样

如第2.1节所述，本研究中使用的数据集在正样本和负样本之间表现出明显的不平衡，因此需要应用重采样技术。常见的重采样策略包括下采样和上采样。下采样涉及减少多数类实例的数量以实现平衡的类分布，而上采样增加少数类样本的数量以达到同样的效果。为了最小化样本空间的信息丢失并保留蛋白质残基级特征，优先考虑上采样。合成少数过采样技术（SMOTE）是最广泛使用的上采样方法之一。随后研究提出了许多SMOTE的变体，本研究采用支持向量机合成少数过采样技术（SVMSMOTE）进行上采样。重采样过程概述如下。

1. 1.

   训练SVM模型：首先，使用少数类样本训练支持向量机（SVM）分类器以确定决策边界。通过最大化类之间的间隔，SVM有效地学习了训练数据的最佳分离边界。

2. 2.

   识别支持向量：支持向量是那些位于决策边界最近的训练样本。这些点在定义分类边界时至关重要。SVMSMOTE利用这些支持向量来指导新合成实例的生成。

3. 3.

   生成合成样本：对于每个少数类支持向量，SVMSMOTE通过在支持向量与其最近邻之间进行插值来生成合成样本。这个过程类似于SMOTE，插值公式由公式（5）给出：x_new = x_sv + λ ? (x_neigh - x_sv)，其中x_new是生成的合成样本，x_sv是当前支持向量，x_neigh是其最近的少数类邻居，λ是范围[0, 1]内的随机数，用于确定插值比例。

4. 4.

   生成所需数量的合成样本：基于上述方法，SVMSMOTE生成足够的合成样本，以实现训练数据集中平衡的类分布。

2.4 Bi-LSTM

考虑到本研究中用于蛋白质表示的两种特征空间的序列性质，采用了基于Bi-LSTM的方法。LSTM在处理蛋白质序列方面的优势在于其有效捕获长程依赖关系的能力。在蛋白质预测任务中，特别是在翻译后修饰（PTM）位点的识别中，长程上下文依赖关系至关重要，因为蛋白质序列中的相关信息通常超出局部邻域。因此，准确的功能预测依赖于模型提取这种扩展上下文的能力。传统的循环神经网络（RNN）在这些场景中挣扎，因为梯度消失和爆炸等挑战阻碍了跨长时序序列的学习。相比之下，Bi-LSTM结合了序列上前向和后向的上下文信息，使其能够捕获双向依赖关系，这对于准确的PTM位点预测至关重要。

标准的Bi-LSTM由一个前向LSTM和一个后向LSTM组成。具体来说，LSTM架构引入了记忆细胞，决定是否保留或丢弃过去输入的信息。这些记忆细胞由三个门控制：输入门、遗忘门和输出门。输入门控制当前输入有多少存储到记忆细胞中，遗忘门决定丢弃先前记忆的哪些部分，输出门确定传递给网络中下一层的信息。

在本研究中，Bi-LSTM被用作核心模型来处理蛋白质序列数据。通过结合前向和后向传递，Bi-LSTM有效地捕获了蛋白质序列内的双向上下文信息，从而提高了基于序列的预测的准确性。

设f表示前向LSTM过程，b表示后向LSTM过程。对于每个时间步t，前向LSTM计算以下公式（6）-（10）：

i_t = σ(W_i ? [h_t-1^(f), x_t] + b_i) (6)

f_t = σ(W_f ? [h_t-1^(f), x_t] + b_f) (7)

o_t = σ(W_o ? [h_t-1^(f), x_t] + b_o) (8)

c_t^(f) = f_t ? c_t-1^(f) + i_t ? tanh(W_c ? [h_t-1^(f), x_t] + b_c) (9)

h_t^(f) = o_t ? tanh(c_t^(f)) (10)

其中σ是sigmoid激活函数，tanh是双曲正切函数，x_t是当前输入，h_t-1^(f)是前一时间步的隐藏状态，c_t-1^(f)是前一时间步的记忆状态。

后向LSTM的操作与前向LSTM类似，但以相反顺序处理序列。在标准的Bi-LSTM模型中，前向和后向LSTM的输出被组合起来以计算每个时间步的最终输出。然而，在本研究中，后向LSTM的初始隐藏状态被修改。具体来说，后向LSTM的初始隐藏状态被设置为前向LSTM初始隐藏状态的加权和，权重是可训练的。这种修改表示为公式（11）：h₀^(b) = σ(W_ih_i^(f))，其中w_i是前向LSTM第i步的权重，h_i是第i步的隐藏状态。这一改进使后向传递能够整合从前向传递学到的信息，从而增强后向传播。

此外，该模型还包含了残差连接，可以表示为公式（12）：h_t+1^(f) = α ? h_t+1^(f) + (1 - α) ? h_t^(f)，其中α是一个超参数，用于调整当前隐藏状态和先前隐藏状态之间的权重。残差连接保留了前一时间步隐藏状态的信息，增强了模型捕获长程依赖关系的能力。

2.5 自适应特征融合模块

本研究的特征提取基于生物信息学特征和来自蛋白质序列片段的语义特征。这两种类型的特征强调信息表示的不同方面：生物信息学特征突出序列的结构和功能特征，而语义特征在自然语言水平上捕获上下文关系。然而，由于这两种模态之间可能存在冗余或不一致，提出了一个自适应特征融合模块来更有效地整合和利用它们。该模块设计有两个主要目标：首先，准确识别和消除冗余信息，以防止不相关或嘈杂的特征损害模型性能；其次，增强对目标任务至关重要的判别性特征，从而提高模型的学习效率和预测准确性。自适应特征融合模块不仅管理不同特征类型的重要性分配，还动态调整融合策略以适应不同任务或数据分布的变化。

本研究提出了一个利用自注意力和交叉注意力机制的融合模块。该模块的目标是使用这些注意力机制动态分配特征权重，从而提取显著信息并实现高效融合。下面详细介绍该模块。

设要融合的特征为H^P（大小N × d₁）和H^W（大小N × d₂）。由于H^P和H^W的维度可能不同，最初建立了投影层，将两个特征映射到共同的维度d。这是通过线性投影后接非线性激活函数σ完成的（公式13，14）：

H₁^P = σ(W₁H^P) (13)

H₁^W = σ(W₂H^W) (14)

其中W₁ ∈ R^d×d₁和W₂ ∈ R^d×d₂。

接下来，自注意力应用于H₁^P和H₁^W。自注意力是一种用于计算序列中每个元素与所有其他元素相关性的机制，使其能够捕获全局上下文。它动态地为序列中的重要元素分配更高的权重，从而提高了模型理解复杂依赖关系的能力。自注意力的计算如公式（15）所示：Self-Attention(Q, K, V) = softmax(QK^T/√d)V，其中Q、K和V是H₁^P和H₁^W的线性变换。具体来说，自注意力应用于提取它们各自的关键特征模式（公式16，17）：

H₂^P = Self-Attention(H₁^P), (16)

H₂^W = Self-Attention(H₁^W), (17)

接下来，使用交叉注意力机制将H₂^P和H₂^W组合起来，以模拟两种类型特征之间的相互作用。交叉注意力通过关注一个特征集相对于另一个特征集的最相关部分来提取重要信息。通过这样做，它有效地融合了来自不同来源的特征，并突出了它们之间的关键关联。交叉注意力的计算类似于自注意力，但Q、K和V的设计不同（公式18）：Cross-Attention(Q, K, V) = softmax(QK^T/√d)V，这里Q来自一个特征集的线性变换，而K和V来自另一个特征集的线性变换。因此，可以从生物信息学特征中提取相关的语义特征信息，如公式（19），（20）所示：

H_PW = Cross-Attention(H₂^P, H₂^W), (19)

H_WP = Cross-Attention(H₂^W, H₂^P), (20)

最后，我们对自注意力和交叉注意力的输出进行加权融合（公式21）：

H^F = α(H_PW + H₂^W) + (1 - α)(H_WP + H₂^P)

其中α是权重参数。融合后的特征表示H^F随后传递给后续的密集层。

由于修饰半胱氨酸周围的位置偏差得到氧化还原化学和化学蛋白质组学的支持，自适应特征融合模块中的注意力权重作为±10窗口内残基级重要性的代理，鼓励模型对与这些生物学上合理上下文一致的位置进行加权。

2.6 交叉验证与性能评估

为了评估我们提出方法的鲁棒性和泛化能力，我们采用了重复10折交叉验证策略。具体来说，整个数据集被随机分成10折，其中9折用于训练，一折用于测试。与其他研究类似，我们采用了常用的评估指标，包括受试者工作特征曲线下面积（auROC）、精确召回曲线下面积（auPRC）、灵敏度（Sn）、特异性（Sp）、准确度（ACC）和马修斯相关系数（MCC），以全面评估模型的预测性能。后四个指标的公式如下（公式22-25）：

Sn = TP / (TP + FN) (22)

Sp = TN / (TN + FP) (23)

ACC = (TP + TN) / (TP + TN + FP + FN) (24)

MCC = (TP?TN - FP?FN) / √((TP+FP)(TP+FN)(TN+FP)(TN+FN)) (25)

其中TP、TN、FP和FN分别表示真阳性、真阴性、假阳性和假阴性。由于我们使用不平衡数据集，将auROC和auPRC作为评估指标是非常合适的。auPRC更好地反映了模型在正类上的性能，尤其是在类别不平衡的情况下。

除了模型性能评估之外，我们还使用基于序列压缩的Kolmogorov复杂度（KC）近似来定量评估特征充分性。简单来说，Kolmogorov复杂度问的是“能够再现给定序列的最短描述是什么？”。如果特征表示保留了原始序列中的大部分规律和模式，那么压缩特征（或特征与原始序列的连接）将几乎与压缩原始序列本身一样有效。因此，小的估计信息丢失率表明特征保留了大多数与预测相关的序列级生化信号，而大的丢失表明重要的模式可能被丢弃。信息丢失率η_loss定义为公式（26）：η_loss = K(X|f(X)) / K(X) = (K(X, f(X)) - K(f(X))) / K(X)，其中K(X)表示字符串X的（前缀）Kolmogorov复杂度，正式地（公式27）：K(X) = min{|p| : 程序p在固定的通用图灵机U上输出X}。换句话说，K(X)是生成X的最短程序的长度（以位为单位）。这个定义由于不变性定理而不依赖于加性常数（独立于X）。同样，K(X, f(X))是X和f(X)连接描述的联合复杂度，K(X|f(X))是条件复杂度（给定f(X)重建X所需的额外信息）。

在本研究中，X表示原始氨基酸序列片段，f(X)表示特征提取映射（fastText子词嵌入和PSSM）。由于精确的Kolmogorov复杂度是不可计算的（由于停机问题的不可判定性），我们使用现成的无损压缩器来近似它，这提供了K的实用上界。具体地，让C(?)是在固定ZIP/DEFLATE压缩器下的压缩长度（字节）；我们估计（公式28）：K?(X) = C(X), K?(f(X)) = C(f(X)), K?(X|f(X)) ≈ C(f(X)‖X) - C(f(X))，其中‖表示与分隔符的连接以避免边界伪影。这个估计器与归一化压缩距离框架密切相关。对于生物符号序列，基于压缩的复杂度估计——通常基于Lempel-Ziv变体——有着悠久的历史，并已成功应用于遗传文本。

这个度量对于我们的亚磺酰化位点预测任务特别关键，我们在没有结构信息的情况下处理21个残基序列片段。鉴于我们的特征将符号序列转换为用于BiLSTM分类的数值表示，量化信息保存确保提取的特征保持生化相关模式同时减少噪声。

3 结果与讨论

3.1 样本序列内容分析

为了更好地说明亚磺酰化位点样本和非亚磺酰化位点样本在残基水平上的差异，本研究采用了Two-Sample Logo技术。该方法利用统计分析来识别正样本和负样本中中心半胱氨酸残基及其周围氨基酸的显著区别。Two-Sample Logo是序列内容分析中各种PTM常用的可视化工具，提供了围绕修饰位点的氨基酸分布模式的清晰见解，并作为位点预测的解释支持。

图2描绘了本研究中使用的正样本集和负样本集之间的Two-Sample Logo比较，显著性阈值（P值）设置为0.5。该图揭示了多个正样本和负样本之间的显著差异。例如，非中心半胱氨酸残基在负样本中更普遍，特别是在位置+1、+4、+5、+6、+7、+8、+9、+10、+14、+15、+17和+21。这一观察结果与双向门控循环单元网络与自注意力（BiGRUD-SA）模型的发现一致，尽管数据集不同，但突出了本研究中提出的BioSemAF-BiLSTM模型的强大泛化能力。此外，其他氨基酸在两个样本组之间也表现出统计学上的显著差异，例如：

? 赖氨酸（K）：经常出现在正样本中的位置+1、+3、+4、+5、+6、+7、+9、+15、+18、+19、+20；

? 谷氨酸（E）：经常出现在正样本中的位置+1、+4、+6、+7、+8、+12、+13、+14、+15、+16、+18；

? 组氨酸（H）：经常出现在负样本中的位置+3、+6、+7、+8、+9、+10、