机器学习是数理统计和计算机科学的延伸，包含许多统计模型和计算机程序算法。1992年，MUGGLETON等［1］运用机器学习算法预测蛋白质二级结构，将蛋白质的二级结构与机器学习的编码相关联。随着机器学习算法的逐渐成熟，出现了用于预测结构、折叠、结合甚至催化活性的新版本机器学习算法，其主要目的是处理有关突变体及其功能特性的积累信息。易华伟等［2］基于蛋白质的氨基酸序列通过机器学习算法预测其稳定性。程淑萍等［3］运用3种不同的机器学习算法预测非编码RNA和蛋白质之间的相互作用，提高预测的准确率。徐开琨等［4］运用特征选择与机器学习算法研究生物蛋白质中的标志物，这对疾病的早期诊断和临床治疗等有重要意义。胡如云等［5］介绍了机器学习在合成生物学领域的典型应用，如启动子预测、酶催化设计、代谢途径构建、基因线路设计等。
以上机器学习算法都是基于大数据集进行的蛋白质工程研究，在只有少量数据集下进行的机器学习算法研究较少。同时，大多研究基于蛋白质的结构，对未知结构的蛋白质分子及其功能特性的研究也很少。基于序列-活性关系（innovative sequence-activity relationship，ISAR）［11］的算法，为蛋白质工程提供了新的研究方向。
ISAR算法解决了数据集少和未知蛋白质结构的问题，有助于探究蛋白质序列和功能特性的相互关系。ISAR算法由FONTAINE等［12］于2018年开发，以氨基酸的物理化学性质为基础，用光谱表示蛋白质。CADET等［13］解释了ISAR算法在改善来自黑曲霉的环氧化物水解酶对映选择性中的应用。该实验基于9个单点突变（29）的组合，预测了这512个变异体的对映选择性，发现了具有更高对映选择性的突变体。并运用ISAR算法对4个不同类型的数据集（热稳定性的细胞色素P450、结合亲和力的TNF-alpha、结合亲和力的GLP-2和热稳定性的肠毒素）进行拟合建模，验证了ISAR算法能有效捕捉突变对蛋白质功能的影响。NICOLAS等［14］提出了用不同组合编码扩展序列的方法建模，测试了4种不同长度的蛋白质（GLP-2、TNFalpha、细胞色素P450和环氧水解酶）及其活性（cAMP激活、结合亲和力、热稳定性和对映选择性），并预测了一个多肽链的适应度值。OSTAFE等［15］通过不同pH的不同介质用ISAR算法提高了葡萄糖氧化酶的活性，使葡萄糖氧化酶突变体（P1）对介质二茂铁-甲醇（12倍）和硝基苯胺（8倍）具有更大的特异性，并在3个pH调整缓冲区中表现出更好的性能，在存在二茂铁甲醇的情况下，与pH为5.5的野生型酶相比，P1的kcat/KM比率增加了121倍。
质子泵型视紫红质（proton-pumping rhodopsin，PPR）［16］ 蛋白是一种结构比较简单且广泛存在的光能转换系统。结合视黄醛分子在吸收光子后令其构象改变，引起整体蛋白的变构等，在一定程度上能促进三磷酸腺苷酶（adenosine triphosphate，ATP）的合成，将光能转化为化学能［17-18］。PPR通过将光能转化为化学能，从而提高细胞的积累生物量。紫色球杆菌视紫红质（gloeobacter violaceus rhodopsin，GR）是一种PPR，通过取代GR的不同氨基酸获得不同的GR突变体，影响其吸收光谱的峰值。吸收光谱的峰值越大，光能转化率越高，微生物的生产和细胞的耐受性越强。因此，需对GR光谱特性进行研究。
通过对GR视网膜结合口袋进行视紫红质的氨基酸替换，组合成新的GR突变体，以调节其光谱最大吸收波长λmax［19］。从数据集中提取81条少量GR突变体，用ISAR算法对数据进行数字信号处理（digital signal processing，DSP），再与实验获得的最大光谱吸收波长进行偏最小二乘回归（partial least squares regression，PLSR）建模，采用交叉验证方法获得符合GR实验数据要求的最佳模型和蛋白质氨基酸序列与光谱特性的关系。
1 ISAR算法
图1
1.1 编码阶段
在氨基酸指数索引（AAindex）数据库中找到符合蛋白质特性的最佳索引。AAindex是含566组索引数字指数的数据库，每组索引包括20种标准氨基酸，代表各自的物理化学和生物化学性质，同时给出了这些索引之间的相关性及其他信息［21］。依据AAindex数据库，将蛋白质的氨基酸序列数字化，再用FFT将数字信息转换为能量对频率的蛋白质光谱形式。其可通过
实现，其中，N为输入信号的长度，S为输出光谱（复数），k为光谱的频率，s（n）为有限长连续信号，n为输入信号的位置，i为虚数单位。
DSP是一种数据预处理技术，对信号进行分解和处理，以不同方式进行数据表达。FFT是DSP的一种方法，BENSON［22］很早就将傅里叶变换方法应用于生物序列的比较。FFT中一个点的变换将影响整条数据的光谱，因此蛋白质位点中任何一个突变点的变化，都会令蛋白质光谱整体变换，这也是ISAR算法的亮点之一。已有很多方法研究单点突变的变换，将每个特性值累加，忽略了氨基酸之间的相互作用。
1.2 建模阶段
对蛋白质光谱与蛋白质突变体的生物化学特性值进行PLSR建模，并用留一交叉验证法（leave-one-out cross validation，LOOCV）验证模型和提高模型的鲁棒性。蛋白质的生物化学特性值除通过实验获得外，还可将生物、数学和计算机相结合，通过机器学习拟合蛋白质的数据预测得到。对无法获得晶体结构的蛋白质，通过已知蛋白质的氨基酸序列特征属性，用ISAR算法便可获得蛋白质的生物化学特性值。
由机器学习算法进行建模，通常需要大量数据集进行训练和验证。在ISAR算法中，可用PLSR解决数据集小的问题。为防止模型过拟合，通过LOOCV优化模型参数和验证模型，同时用确定系数R2和均方误差E（
其中，
根据AAindex数据库中566组索引与实验数据进行建模，获得多个模型。在同一组索引下，由不同参数得到若干个不同模型。在参数不变的情况下，566组索引通过LOOCV得到n×566个模型。不同参数和不同组索引得到不同的确定系数R2和均方误差E。根据最小E和较优R2找到最佳索引，从而获得最佳模型。
1.3 预测阶段
将建模阶段选出的最佳索引经实验数据预处理后放入模型，预测与实验数据相对应的特征属性值，同时验证氨基酸序列与蛋白质功能属性的相关性。因此，使用有限的实验数据通过ISAR算法拟合可获得适应其特征属性的最佳模型，进一步预测其他可能存在的更好的突变体。
2 结果与讨论
2.1 实验数据
GR通过定向进化获得81条突变体吸收波长数据集，包含单点突变和组合突变（位点：W122，V126，M158，G162，E166，G178，S181，F185和A256）［23］。GR的野生型光谱最大吸收波长λmax为539 nm，通过替换不同的氨基酸，获得的突变体的λmax为50~650 nm。因此，组合突变体和氨基酸之间的相互作用会影响GR的λmax。
数据预处理是建模前的关键步骤。首先将81条数据进行氨基酸编码，根据AAindex数据库，每组索引有20种不同的氨基酸，每种氨基酸用不同的数值代替。GR野生型（wide type，WT）是由20个标准氨基酸组成的长为298的氨基酸序列。首先根据AAindex数据库进行数值化处理，然后由FFT转换为蛋白质光谱。图2为GR的WT、单点突变体W122A和组合突变体G162L/E166W/F185A/A256S的编码过程。其他突变体的光谱转换类似。W122A为GR的第122个位点，是色氨酸（tryptophan，W）突变成丙氨酸（alanine，A）的单点突变体。组合突变体G162L/E166W/F185A/A256S由4个位点突变成其他氨基酸得到。3种蛋白质通过FFT得到蛋白质光谱，图3分别为野生型、突变体W122A和组合突变体G162L/E166W/F185A/A256S在不同频率下的振幅。由图3可知，虽然W122A为单点突变、G162L/E166W/F185A/A256S只有4个位点突变，但影响的是整个蛋白的光谱振幅，与WT得到的蛋白质光谱完全不同，这也是ISAR算法的优势。因此，81条数据通过FFT得到完全不同的81条蛋白质光谱，为后续的吸收光谱建模打下良好基础。
图2
图2
用ISAR算法由GR数据得到的不同蛋白质光谱
Fig. 2
Different protein spectra are obtained from GR data by ISAR method
图3
图3
3个蛋白质通过FFT转变为蛋白质光谱
Fig.3
Three proteins are transformed into protein spectra by FFT
2.2 选择最佳索引
GR是由20种标准氨基酸组成的蛋白质，氨基酸的数值化和相互作用均会影响GR的吸收光谱特性。AAindex数据库有566组氨基酸索引，其中16组因部分氨基酸数值缺失被剔除。利用剩下的550组索引在固定参数下与PLSR建立模型，并测试其E值，如图4所示。不同索引的E值不同，550组索引的E值差异较大，对模型的建立和预测有很大影响。由
图4
图4
根据不同参数获得550个λmax模型的E值
Fig.4
The E of λmax 550 models according to different parameters
表1 不同索引下的R2和E
Table 1
2.3 建模和预测
由最小E选择最佳索引LEVM760106，选出最佳模型参数，进行建模与预测。运用ISAR算法预测GR的81条实验数据，得到λmax的实验值和预测值、E=11.64、确定系数R2=0.944，模型的拟合度很高。E和R2是衡量模型优劣的标准，E用于选择最佳索引和最佳模型，R2体现模型的预测能力。由图5知，λmax的预测值和实验值在直线附近聚集，偏离不大，模型鲁棒性较好。
图5
图5
R2= 0.944时GR及突变体的LOOCV预测
Fig. 5
Prediction of GR and mutants by LOOCV when R2=0.944
2.4 方法比较
ISAR算法中的一个重要环节是用FFT进行数据处理。FFT将单个位点或多个位点的突变体转换为影响整个蛋白质的变换，改进了以往生物特性功能单纯累加的缺陷。同时本研究也运用小波变换预处理数据，得到的结果如表2所示。小波变换也是DSP的一种方法，是对短傅里叶变换的继承和发展，通过时间或空间的局部变换，突出某地方的特征。
表2 不同方法对GR的验证结果
Table 2
| 索引 | 数据处理方法 | 交叉验证 方法 | 样本数 | E | R2 |
|---|---|---|---|---|---|
| LEVM760106 | FFT | LOOCV | 81 | 11.64 | 0.944 |
| LEVM760106 | FFT | 十折交叉 验证法 | 81 | 14.76 | 0.908 |
| LEVM760106 | 小波变换 | LOOCV | 81 | 12.12 | 0.940 |
| LEVM760107 | 小波变换 | 十折交叉 验证法 | 81 | 17.35 | 0.871 |
除用LOOCV外，还用十折交叉验证法验证模型的优劣。十折交叉验证法是将样本数n分成10份，选择9份作为训练集，1份作为验证集。通过不同的数据预处理和不同的交叉验证方法，得到不同的E和R2。E越小、R2越接近于1，模型的拟合度越高、鲁棒性越好。由表2可知，用FFT选择的最佳索引均为LEVM760106，且十折交叉验证法增加了验证集数量，令E和R2均有降低。小波变换采用的小波基为db1，得到的E值均较FFT方法大。虽然用小波变换和LOOCV得到的R2与ISAR算法的相差不多，但是E值高了约0.5，模型拟合度不及ISAR算法。实验证明ISAR算法中的FFT非常适合GR的λmax拟合，可获得较优模型，且方便后期的拟合和预测。由表2可知，用4种方法寻找最佳索引，3种方法得到的最佳索引均为LEVM760106。
表3 索引LEVM760106中20种氨基酸的数值表示
Table 3
| 氨基酸名称 | 数值 | 氨基酸名称 | 数值 |
|---|---|---|---|
| 丙氨酸 | 5.2 | 亮氨酸 | 7.0 |
| 精氨酸 | 6.0 | 赖氨酸 | 6.0 |
| 天冬酰胺 | 5.0 | 甲硫氨酸 | 6.8 |
| 天冬氨酸 | 5.0 | 苯丙氨酸 | 7.1 |
| 半胱氨酸 | 6.1 | 脯氨酸 | 6.2 |
| 谷氨酰胺 | 6.0 | 丝氨酸 | 4.9 |
| 谷氨酸 | 6.0 | 苏氨酸 | 5.0 |
| 甘氨酸 | 4.2 | 色氨酸 | 7.6 |
| 组氨酸 | 6.0 | 酪氨酸 | 7.1 |
| 异亮氨酸 | 7.0 | 缬氨酸 | 6.4 |
3 结 论
运用AAindex数据库中最佳索引的数字化蛋白质氨基酸序列，用FFT预处理得到影响整体蛋白质光谱的部分变化，与GR的光谱最大吸收波长λmax进行PLSR建模，得到拟合度较高的模型。最佳模型选取的索引为LEVM760106，其中确定系数R2为0.944，均方误差E为11.64。由该索引的物理化学性质，可得到在突变前后的紫色球杆菌视紫红质蛋白序列中，其范德华参数R0变化所引起的物理化学性质与吸收光谱波长的相关性最大。这证明了经傅里叶变换后的频率（蛋白质谱）可有效描述和预测氨基酸序列的蛋白质活性（序列-活性关系用蛋白质谱建模）。此外，此频率考虑了突变对光谱的整体影响，而不是局部影响。ISAR算法的优点之一是不依赖于蛋白质的空间结构，只要获得一定长度的蛋白质氨基酸序列就可与蛋白质的生物化学特性进行建模。优点之二是考虑了氨基酸之间的相互作用力，而不是单纯功能属性值的累加。另外，ISAR算法中的LOOCV和PLSR有效解决了机器学习中数据集少的问题，在数据集量较少时也可进行机器学习。
研究发现，统计预测模型基于氨基酸的物理化学特性建立，并研究了蛋白质氨基酸序列与功能特性的相互关系。用基于已有数据的ISAR算法做预测能否得到比实验数据更好的突变体有待进一步研究。如能通过ISAR算法预测更好的突变体，这将大大节约生物化工实验时间，降低实验成本。统计模型的建立依赖于机器学习，可以预见，机器学习算法将适应更多的生物工程需求，应用于更广泛的生物化工领域。
http://dx.doi.org/10.3785/j.issn.1008-9497.2022.03.003
参考文献
Protein secondary structure prediction using logic-based machine learning
［J］.
基于氨基酸序列和模拟结构预测蛋白质稳定性的研究进展
［J］.
Research progress on the prediction of protein stability based on amino acid sequence and simulated structure
［J］.
基于机器学习方法的非编码RNA-蛋白质相互作用的预测
［J］.
Prediction of ncRNA-protein interactions based on machine learning methods
［J］.
基于质谱的蛋白质生物标志物发现中的特征选择与机器学习方法研究进展
［J］.
Research progress of feature selection and machine learning methods for mass spectrometry-based protein biomarker discovery
［J］.
面向合成生物学的机器学习方法及应用
［J］.
Machine learning for synthetic biology： Methods and applications
［J］.
Design and evolution of enzymes for non-natural chemistry
［J］.
Solubilization and iterative saturation mutagenesis of α1，3-fucosyltransferase from helicobacter pylori to enhance its catalytic efficiency
［J］.
蛋白质工程： 从定向进化到计算设计
［J］.
Protein engineering： From directed evolution to computational design
［J］.
机器学习助力酶定向 进化
［J］.
Machine learning assisted enzyme directed evolution
［J］.
Introduction to machine learning
［J］.
A machine learning approach for reliable prediction of amino acid interactions and its application in the directed evolution of enantioselective enzymes
［J］.
Method and electronic system for predicting at least one fitness value of a protein， related computer program product
Application of fourier transform and proteochemometrics principles to protein engineering
［J］.
Novel descriptors and digital signal processing-Based method for protein sequence activity relationship study
［J］.
One-shot optimization of multiple enzyme parameters： Tailoring glucose oxidase for pH and electron mediators
［J］.
Bacterial rhodopsin： Evidence for a new type of phototrophy in the sea
［J］.
Bacteriorhodopsin-like proteins of eubacteria and fungi： The extent of conservation of the haloarchaeal proton-pumping mechanism
［J］.
Potential of proton-pumping rhodopsins： Engineering photosystems into microorganisms
［J］.
Directed evolution of Gloeobacter violaceus rhodopsin spectral properties
［J］.
An algorithm for the machine calculation of complex Fourier series
［J］.
AAindex： Amino acid index database， progress report 2008
［J］.
Digital signal processing methods for biosequence comparison
［J］.
Learned protein embeddings for machine learning
［J］.
Digital signal processing techniques： Calculating biological functionalities
［J］.
A simplified representation of protein conformations for rapid simulation of protein folding
［J］.
Relation between amino acid composition and cellular location of proteins
［J］.
Physical reasons for secondary structure stability： Alpha-helices in short peptides
［J］.
