蛋白质分类中的数学基础与应用

Author：Shibo Liu
Notice：All rights reserved

1.Abstract

​ 给定蛋白质的氨基酸序列，预测它所属的CATH超家族是一个有趣的问题。我提出了一个神经网络模型，它用来对蛋白质数据进行分类。该模型使用了卷积、残差连接、注意力机制等方法，在seq1024.hdf5数据集上取得了92.01%的准确率。经过调研，我并未发现有其他人在CATH数据集上使用如本文方法的工作，如有雷同，纯属巧合。

2.Introduction

​ 蛋白质是有机大分子，氨基酸是蛋白质的基本组成单位。蛋白质是由氨基酸以“脱水缩合”的方式组成的多肽链、经过盘曲折叠形成的具有一定空间结构的物质。人体内蛋白质的种类很多，性质、功能各异，但都是由20多种氨基酸按不同比例组合而成的。所以，通过对氨基酸编号，我们可以用一个包含氨基酸编号的向量表示不同的蛋白质，向量中前后相邻的元素表示在空间中通过肽键相互连接的氨基酸。

​ 由上述数据结构特性，我们很自然地想到了卷积神经网络[1]。不过，不同于二维的图像，我们应该使用一维卷积。此外，注意力机制[2]对于序列数据也有不错的效果，因此我将其运用到模型中。

3.Related work

​ 一些工作采用传统的机器学习方法来对蛋白质进行检索或分类，如：InterPro[9]数据库使用不同类型的诊断模型（如隐马尔可夫模型，正则表达式等）来搜索蛋白质序列，并将其归入相应的家族或域；Cheng[10]等人提出了一种基于文本文档分类技术的蛋白质分类方法。这种方法类似于文档分类，使用决策树和朴素贝叶斯分类器，结合卡方特征选择，对蛋白质序列中的n-gram（即长度为n的短肽段）进行计数，并根据计数结果对蛋白质序列进行分类。

​ 也有工作采用深度学习的方法对蛋白质进行分类：UDSMProt[11]是一种利用自监督语言模型来对蛋白质序列进行通用的深度表示学习和分类的方法。这种方法使用蛋白质序列的一级结构，即氨基酸序列。

​ 还有的工作结合了机器学习与深度学习的方法：Diplaris[12]等人使用了多种算法，包括支持向量机、随机森林和神经网络等，对蛋白质进行分类，并比较了它们的性能。

​ 综上所述，从方法方面讲：未找到将卷积、残差连接、注意力机制运用到蛋白质分类的工作；从数据方面讲：未找到在CATH数据集上进行蛋白质分类的工作。

4.Method

4.1.Architecture

​ 模型的整体架构如图1(1)所示：嵌入模块对输入进行编码，卷积模块捕获蛋白质的局部特征，注意力模块捕获蛋白质的全局特征，头模块将特征映射到各个蛋白质类别。

​ 需要说明的是，我并非一开始就构造了这样的模型：初始只是简单地增加卷积的层数，达到一定数量后效果不再提升，通过引入残差连接[3]解决了gradient vanishing/gradient explosion的问题；为了捕获氨基酸之间的相关性，加入了注意力层。卷积模块如图1(2)所示，注意力模块如图1(3)所示。此外，我还尝试了Metaformer[4]提出的结构，即在卷积模块和注意力模块中都有MLP，但是它的效果并不好，且全连接层增加了计算量。

​ 在卷积模块中，我们对输入进行第一次卷积，第一次BatchNorm[5]调节分布、并用ReLU[6]激活；然后第二次卷积，卷积后的向量与输入加和，做第二次BatchNorm和ReLU。

​ 在注意力模块中，我们先对输入进行BatchNorm，然后做注意力运算，再与输入加和并用ReLU激活。

​ 嵌入模块对氨基酸序列的两端进行填充（补0），使得所有序列都是相同的长度（400）。

​ 头模块是一个全连接层，输出维度与要预测的蛋白质种类数（6630）相等。

4.2.Batch normalization

​ 这里给出Batch normalization的算法：

​ 输入：$x$的一个小批量： B = { x 1 … m } \Beta=\{x_{1 \dots m}\} B={x1…m​}

​ 需要学习的参数：$\gamma ,\beta$

​ 输出： { y i = B N γ , β ( x i ) } \{y_i=BN_{\gamma,\beta}(x_i)\} {yi​=BNγ,β​(xi​)}

​ ${\mu }_{\mathrm{B}}\leftarrow \frac{1}{m}{\sum }_{i=1}^m{x}_i$ //小批量均值

​ ${\sigma }_{\mathrm{B}}^2\leftarrow \frac{1}{m}{\sum }_{i=1}^m(x_i-{\mu }_{\mathrm{B}}{)}^2$ //小批量方差

​ ${\hat{x}}_i\leftarrow \frac{x_i-{\mu }_{\mathrm{B}}}{\sqrt{{\sigma }_{\mathrm{B}}^2+ϵ}}$ //正则化

​ $y_i\leftarrow \gamma {\hat{x}}_i+\beta \equiv B{N}_{\gamma ,\beta }(x_i)$ //放缩和平移

​ 参数$\gamma$和$\beta$分别用来学习最优的均值和方差的缩放和偏移参数。

4.3.Activation function

​ 在头模块之前的模块中，我们使用的是ReLU激活函数，如公式(1)所示，$f(x)$表示输出分量，$x$表示输入分量。
$$f(x)=\{\begin{array}{l}0,x\le 0\\ x,x>0\end{array}$$
​ 对于头模块的输出，我们使用的激活函数是softmax[7]激活函数。如公式(2)所示，对输出向量$Z$的每个分量$z_i$，我们通过计算它的softmax值来进行归一化。令归一化后的向量为$\tilde{y}$。
$${\tilde{y}}_i=softmax(z_i)=\frac{e^{z_i}}{{\sum }_{j=1}^n{e}^{z_j}}$$

4.4.Loss function

​ 在完成softmax归一化之后，我们采用交叉熵损失函数来衡量预测值与真实值之间的误差。如公式(3)所示，$y$表示样本的真实标签，$\tilde{y}$表示样本的预测标签。
$$Loss=-ylog(\tilde{y})+(y-1)log(1-\tilde{y})$$

4.5.Back propagation

​ 依据上述内容，我们开始反向传播计算梯度。首先是$Loss$对$\tilde{y}$求导，如公式(4)所示。
$$\frac{dLoss}{d\tilde{y}}=-\frac{y}{\tilde{y}}-\frac{y-1}{1-\tilde{y}}=-\frac{y}{\tilde{y}}+\frac{1-y}{1-\tilde{y}}$$
​ 然后是$\tilde{y}$对求$z_i$导，如公式(5)所示。
$$\frac{d\tilde{y}}{d{z}_i}=\frac{e^{z_i}{\sum }_{j=1}^n{e}^{z_j}-e^{z_i}{e}^{z_i}}{({\sum }_{j=1}^n{e}^{z_j}{)}^2}=\frac{e^{z_i}({\sum }_{j=1}^n{e}^{z_j}-e^{z_i})}{({\sum }_{j=1}^n{e}^{z_j})({\sum }_{j=1}^n{e}^{z_j})}=\tilde{y}(1-\tilde{y})$$
​ 根据链式法则，可得$Loss$对$z_i$的导数，如公式(6)所示。
$$\frac{dLoss}{d{z}_i}=-y(1-\tilde{y})+(1-y)\tilde{y}=-y+\tilde{y}$$
​ 我们惊喜地发现：当最后一个激活函数采用softmax，且损失函数采用交叉熵损失函数时，第一阶段反向传播产生的梯度的形式非常好，避免了精度的损失。

​ 对于头模块的全连接层，其数学形式如公式(7)所示，其中$n$表示上一个模块的输出维度，$m$表示蛋白质种类数。
$$F(X_n)=W_{m\times n}{X}_n$$
​ 对于卷积模块，可以看作是对于输入进行分块（可能是有重叠的分块，取决于卷积核的大小和步长），对每一个分块进行类似全连接的操作，其中$b$为偏置，如公式(8)所示。
$$F(X)=WX+b$$
​ 对于注意力模块，其$Q,K,V$向量也是由输入向量$X$经过线性变换得到的，如公式(9)所示。
$$Q,K,V=W^{Q}X,W^{K}X,W^{V}X$$
​ 所以无论是全连接、卷积还是注意力，在计算权重$W$的梯度时，都是用反向传回的梯度乘以输入$X$。于是权重$W$可以进行更新，如公式(10)所示，$▽$表示回传的梯度，$lr$表示学习率。
$$W-=lr\ast X▽$$

​ 通过反向传播，误差传递到各个参数，参数据此进行更新，以更好地逼近目标。

5.Experiments

5.1.Dataset

​ CATH数据库是 蛋白质数据库中 蛋白质结构的 分层域分类[8]。此层次结构中有四个主要级别：

• Class：根据其二级结构组成进行分类（主α，主β，混合α/β或少量的二级结构）。
• Architecture：根据其总体形状进行分类，由三维空间中二级结构的方向决定，但忽略了它们之间的连接。
• Topology：根据二级结构的整体形状和连接进行分类。
• Homologous：这一层次将被认为有共同祖先的蛋白质结构域分组在一起，因此可以被描述为同源的。

​ 可以将CATH理解为一棵树，第一层是Class，第二层是Architecture，第三层是Topology，第四层是Homologous。每增加一层，都是类别的进一步细化，如下图所示。

​ 我们使用python的h5py库，数据集的路径为：/data/cath/hdf5/seq1024.hdf5。该数据集共有3924783个蛋白质序列，分为6630种蛋白质。蛋白质序列是不等长的向量（最大长度为399），向量的分量是数字1~20，表示20种氨基酸。

​ 举个例子，我们用[19,17,18,15,4,10,11,15,18,1,8,4,7,10,14,14,15,6,18,3,4,16,15,10,16,18,4,10,10,10,1,3,1,10,9,2,13,15,8,15,10,20,17,12,5,3,18,13,18,17,1,1,4,10,6,17,5,15,17,10,11,14,15,15,10,7,6,4,13,18,14,20,8,18,6,16,17,7,5,11,6,10,4,5]来表示某种CATH类别为1.10.8.10的蛋白质。不过，在实际的实验中，我采用长度为6630的one-hot编码表示蛋白质所属的类别。

​ 测试集包含76800个蛋白质序列，剩下的作为训练集。

5.2.Setup

​ 卷积模块的数量为6，所有模块的卷积核大小都是3，其中前三个模块的两次卷积步长都为1，后三个模块的第一次卷积步长为2、第二次卷积的步长为1。在注意力模块中，使用2头自注意。

​ 输入的维度变化如表1所示。

• loss function：CrossEntropyLoss
• Optimizer：Adam
• learning rate：0.001
• iterations：10

6.Results

​ 训练过程中的准确率和损失如图2所示，最高准确率达到了92.01%。我们发现，虽然加入注意力机制并未对模型最终的准确率有较大的提升，但在启动阶段表现良好，刚开始就有很高的准确率。

7.Conclusion

​ 总体来说，我的模型取得了不错的效果，但可能是因为训练集太大，测试集太小，所以效果较好。后续可以尝试增大测试集的比例，并划分出验证集。此外，python的h5py库中还有包含了α碳原子空间坐标的蛋白质数据集，路径为：/data/cath/hdf5/struct256.hdf5，后续工作可以使用图神经网络对其进行分类。

8.Reference

[1]

KRIZHEVSKY A, SUTSKEVER I, HINTON G E. ImageNet Classification with Deep Convolutional Neural Networks[J/OL]. Communications of the ACM, 2017: 84-90. http://dx.doi.org/10.1145/3065386. DOI:10.1145/3065386.

[2]

DEVLIN J, CHANG M W, LEE K, et al. BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding[C/OL]//Proceedings of the 2019 Conference of the North, Minneapolis, Minnesota. 2019. http://dx.doi.org/10.18653/v1/n19-1423. DOI:10.18653/v1/n19-1423.

[3]

HE K, ZHANG X, REN S, et al. Deep Residual Learning for Image Recognition[C/OL]//2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), Las Vegas, NV, USA. 2016. http://dx.doi.org/10.1109/cvpr.2016.90. DOI:10.1109/cvpr.2016.90.

[4]

YU W, SI C, ZHOU P, et al. MetaFormer Baselines for Vision[J]. 2022.

[5]

IOFFE S, SZEGEDY C. Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift[J]. arXiv: Learning,arXiv: Learning, 2015.

[6]

NAIR V, HINTON GeoffreyE. Rectified Linear Units Improve Restricted Boltzmann Machines[J]. International Conference on Machine Learning,International Conference on Machine Learning, 2010.

[7]

Sutton R S, Barto A G. Reinforcement learning: An introduction[M]. MIT press, 2018.

[8]

Orengo C A, Michie A D, Jones S, et al. CATH–a hierarchic classification of protein domain structures[J]. Structure, 1997, 5(8): 1093-1109.

[9]

McDowall J, Hunter S. InterPro protein classification[J]. Bioinformatics for Comparative Proteomics, 2011: 37-47.

[10]

Cheng B Y M, Carbonell J G, Klein‐Seetharaman J. Protein classification based on text document classification techniques[J]. Proteins: Structure, Function, and Bioinformatics, 2005, 58(4): 955-970.

[11]

Strodthoff N, Wagner P, Wenzel M, et al. UDSMProt: universal deep sequence models for protein classification[J]. Bioinformatics, 2020, 36(8): 2401-2409.

[12]

Diplaris S, Tsoumakas G, Mitkas P A, et al. Protein classification with multiple algorithms[C]//Advances in Informatics: 10th Panhellenic Conference on Informatics, PCI 2005, Volas, Greece, November 11-13, 2005. Proceedings 10. Springer Berlin Heidelberg, 2005: 448-456.

9.Appendix

9.1.Definition

gradient vanishing/gradient explosion

​ 在反向传播计算梯度的过程中，需要对激活函数进行求导，运用链式法则将各部分的导数做乘积。如果很多层的导数值都大于1，那么最终的梯度将以指数形式增加，即发生梯度爆炸；如果很多层的导数值都小于1，那么最终的梯度将会以指数形式衰减，即发生了梯度消失。

MLP

​ 多层感知机（Multilayer Perceptron，MLP）是一种基本的人工神经网络模型，它由多个神经元组成的多个层次构成。每个神经元都与前一层的所有神经元连接，并且每条连接都有一个权重，通过调整权重来实现模型的学习。

​ 多层感知机通常包括输入层、隐藏层和输出层。输入层接收外部数据输入，隐藏层用于处理输入数据并提取特征，输出层则产生模型的输出结果。隐藏层可以包含多个层次，每个隐藏层都可以包含多个神经元。

​ 多层感知机通过反向传播算法来训练模型，通过不断调整权重来最小化损失函数，从而使模型能够适应输入数据的特征并进行准确的预测。

​ 多层感知机在图像识别、语音识别、自然语言处理等领域有广泛的应用，它能够处理复杂的非线性关系，并且具有较强的泛化能力。因此，多层感知机是深度学习的基础，也是目前人工智能领域最常用的模型之一。

Batch normalization

​ Batch normalization是一种用于深度神经网络的技术，它可以加快网络的训练速度并增加网络的稳定性。Batch normalization通过在每个训练批次中对每层的输入进行标准化来实现这一点。这意味着它将每个输入的均值和方差调整为零均值和单位方差。

​ 通过标准化输入，Batch normalization可以减少内部协变量偏移，这是指在训练过程中每一层输入分布的变化。这种偏移可能会减慢网络的收敛速度，甚至导致训练失败。Batch normalization可以减少这种偏移，从而加速网络的训练速度。

​ 此外，Batch normalization还可以增加网络的稳定性，减少对初始参数的敏感度，并且可以允许更高的学习率。这使得网络更容易训练，并且可以更快地收敛到一个较好的解。

​ 总之，Batch normalization是一种用于深度神经网络的重要技术，它可以加快网络的训练速度并增加网络的稳定性，从而使得训练更加高效和可靠。

Activation function

​ 激活函数是神经网络中非常重要的一个组成部分，它在神经元中负责对输入信号进行非线性变换，从而使神经网络具有非线性映射能力。激活函数的作用是引入非线性因素，使得神经网络可以学习和表示更加复杂的函数关系。

​ 常见的激活函数包括Sigmoid函数、Tanh函数、ReLU函数等。Sigmoid函数将输入值映射到0到1之间，Tanh函数将输入值映射到-1到1之间，而ReLU函数则是对输入值进行阈值处理，小于0的部分置为0，大于0的部分保持不变。

​ 激活函数在神经网络中起到了至关重要的作用，它可以帮助神经网络学习并表示复杂的非线性关系，提高神经网络的表达能力和学习能力。因此，选择合适的激活函数对于神经网络的性能和效果至关重要。

Loss function

​ 损失函数是在机器学习和深度学习中用来衡量模型预测结果与真实标签之间差异的函数。它的作用是衡量模型在训练过程中的表现，并且帮助优化算法调整模型参数以最小化预测结果与真实标签之间的差异。

​ 常见的损失函数包括均方误差（Mean Squared Error，MSE）、交叉熵损失函数（Cross-Entropy Loss）、对数损失函数（Log Loss）等。不同的损失函数适用于不同的问题和模型结构。

​ 在训练过程中，优化算法会尝试最小化损失函数的值，以使模型的预测结果更加接近真实标签。因此，选择合适的损失函数对于模型的训练和性能至关重要。

Back propagation

​ 反向传播是一种用于训练神经网络的算法。它通过计算损失函数对每个参数的梯度，并将这些梯度传播回网络的每一层，从而调整网络中的权重和偏置，以最小化损失函数。反向传播算法的主要步骤包括前向传播，计算损失函数，计算梯度，然后使用梯度下降或其他优化算法来更新网络参数。这个过程一直持续到网络的性能达到满意的水平。反向传播算法的发明被认为是深度学习领域的重要里程碑，它使得训练深度神经网络成为可能。

9.2.Partial key codes

9.2.1.One-dimensional convolution

class Conv1d(nn.Module):
    def __init__(self, in_channels:int, out_channels:int, kernel_size:int=3, 
                 stride:int=1, padding:int=0, bias:bool=True) -> None:
        super(Conv1d, self).__init__()
        self.stride = stride
        self.pad = padding
        self.ker_size = kernel_size        
        self.out_ch = out_channels
        self.kernel = nn.Parameter(nn.init.xavier_normal_(pt.randn(out_channels, in_channels, kernel_size)))
        self.bias = None
        if bias:
            self.bias = nn.Parameter(pt.zeros(out_channels))

    def forward(self, x):    
        # x：(batch_size, in_channel, in_len)
        out_len = (x.shape[2] + 2*self.pad - self.ker_size)//self.stride + 1
        x = nn.functional.pad(x, (self.pad, self.pad), "constant", 0)
        
        # (batch_size, in_ch, 1, in_len) -> (batch_size, in_ch*ker_size, out_len)
        x_unfold = nn.functional.unfold(x.unsqueeze(2), kernel_size=(1,self.ker_size), stride=(1,self.stride))
        out = self.kernel.view(self.out_ch,-1).matmul(x_unfold).view(x.shape[0], self.out_ch, out_len)
                 
        if self.bias is not None:
            b = pt.unsqueeze(self.bias, 0).unsqueeze(2)
            out += b
 
        return out

9.2.2.Multi-head self-attention

class SelfAttention(nn.Module):
    def __init__(self, d_m, d_k) -> None:
        super(SelfAttention, self).__init__()
        self.d_m = d_m
        self.d_k = d_k
        self.Wq = nn.Linear(in_features=self.d_m, out_features=self.d_k)
        self.Wk = nn.Linear(in_features=self.d_m, out_features=self.d_k)
        self.Wv = nn.Linear(in_features=self.d_m, out_features=self.d_k)

    def forward(self, x):
        q,k,v = self.Wq(x), self.Wk(x),self.Wv(x)
        score = pt.einsum('nci, ncj -> nc',q,k)
        score /= pt.sqrt(pt.tensor(self.d_k))
        score = nn.functional.softmax(score, dim=-1)
        out = v * score[:,:,None]
        return out
    

class MultiHeadSelfAttention(nn.Module):
    def __init__(self, embed_dim:int, num_heads=1) -> None:
        super(MultiHeadSelfAttention, self).__init__()
        assert embed_dim % num_heads == 0, 'embed_dim must be divided by num_heads'
        self.heads = nn.ModuleList([SelfAttention(d_m=embed_dim, d_k=embed_dim//num_heads) for _ in range(num_heads)])
        self.Wo = nn.Linear(in_features=embed_dim, out_features=embed_dim)
    
    def forward(self, x):
        Z = pt.cat([head(x) for head in self.heads], dim=-1)
        if len(self.heads) > 1:
            return self.Wo(Z)
        else:
            return Z