最近在kaggle参加了诺维信举办的酶稳定性预测比赛 ，最终很幸运获得了第三名，这篇文章主要是简单介绍一下解决方案，具体的数据和题目要求可访问上面的链接。
文章目录
- 模型概述
- 特征工程
- 模型
- XGB1
- XGB2
- XGB3
- Ensemble
- NESP 3D geometry
- FoldX
- ThermoNet v2
- Relaxed rosetta scores
- RMSD
- Different pLDDT
- 总结
模型概述
我分别使用了如下三个数据集训练了三个11-fold xgb模型。
- 数据集1
- 数据集2
- 数据集3
和以下公共笔记本ensemble后得到最终结果。
- 3D geometry
- ThermoNet v2
- rosetta
- rmsd
- different pLDDT
特征工程
受这篇论文[1]的启发，我选择了离突变残基最近的21个残基来构建结构特征。
- 计算21个相邻残基和突变残基之间的相对距离。 r C r_{C} rC表示突变残基的位置， r i r_{i} ri表示第 i i i个相邻残基。 D c i D_{ci} Dci表示突变残基和最近邻残基之间的相对距离。 x i ， y i ， z i x_{i}，y_{i}，z_{i} xi，yi，zi表示残基的三维坐标，由PDB文件计算。我们可以按原子进行分组，在对每组求平均得到残基坐标。
D c i = r c − r i m a x i ∈ K ( r c − r i ) , i = 1 , 2 , . . . , 21 D_{ci}=\frac{r_{c}-r_{i}}{max_{i∈K}{(r_{c}-r_{i})}},i=1,2,...,21 Dci=maxi∈K(rc−ri)rc−ri,i=1,2,...,21
r c i = ( x c − x i ) 2 + ( y c − y i ) 2 + ( z c − z i ) 2 r_{ci}=\sqrt{{(x_{c}-x_{i})}^{2}+{(y_{c}-y_{i})}^{2}+{(z_{c}-z_{i})}^{2}} rci=(xc−xi)2+(yc−yi)2+(zc−zi)2
- 计算平均相对距离(aver dist)。
- 为了更好地代表蛋白质的几何构型，还需计算了突变残基和21个相邻残基之间的夹角。
α
1
=
a
r
c
c
o
s
(
z
c
−
z
i
)
r
c
i
\alpha_{1}=arccos\frac{(z_{c}-z_{i})}{r_{ci}}
α1=arccosrci(zc−zi)
α
2
=
a
r
c
c
o
s
(
x
c
−
x
i
)
r
c
i
\alpha_{2}=arccos\frac{(x_{c}-x_{i})}{r_{ci}}
α2=arccosrci(xc−xi)
α
3
=
a
r
c
c
o
s
(
y
c
−
y
i
)
r
c
i
\alpha_{3}=arccos\frac{(y_{c}-y_{i})}{r_{ci}}
α3=arccosrci(yc−yi)
β
4
=
a
r
c
t
a
n
(
y
c
−
y
i
)
(
x
c
−
x
i
)
\beta_{4}=arctan\frac{(y_{c}-y_{i})}{(x_{c}-x_{i})}
β4=arctan(xc−xi)(yc−yi)
β
5
=
a
r
c
t
a
n
(
z
c
−
z
i
)
(
y
c
−
y
i
)
\beta_{5}=arctan\frac{(z_{c}-z_{i})}{(y_{c}-y_{i})}
β5=arctan(yc−yi)(zc−zi)
β
6
=
a
r
c
t
a
n
(
x
c
−
x
i
)
(
z
c
−
z
i
)
\beta_{6}=arctan\frac{(x_{c}-x_{i})}{(z_{c}-z_{i})}
β6=arctan(zc−zi)(xc−xi)
-
计算突变残基到蛋白质中心的距离(location3d)。
-
B_factor, blosum40,45,60,62,80,100, demask, 溶剂可及表面积, Apolar, Backbone, Sidechain, Ratio, In/Out等等。
-
B_factor：由AlphaFold2(https://alphafold.com/) 生成的PDB文件提供。PDB文件包含代表AlphaFold2模型预测置信度(每个残基)的列b_factor，已发现该置信度与蛋白质稳定性相关。
-
Blosum(Substitution Matrix)：替换矩阵是生物学中的一个概念，用于确定两个序列之间的变化的影响。使用不同的比对数据库存在若干组Blosum矩阵，以数字命名。具有高数目的Blosum矩阵被设计用于比较紧密相关的序列，而那些具有低数目的Blosum矩阵被设计用于比较遥远的相关序列。
-
Demask：DemASK(https://demask.princeton.edu/) 使用一个简单的线性模型来预测单一氨基酸替换对功能的影响。
-
溶剂可及表面积：该特征表示暴露在溶剂中的残留物的面积。一般来说，埋藏的蛋白质残基更有可能改变热稳定性。我们可以使用Bioython来获得这个特征，或者我们可以使用getarea(http://curie.utmb.edu/area_man.html) 网站来获得这个功能。如果我们使用getArea网站，我们将获得以下形式的.txt文件。名为TATAL的栏目代表SASA。Backbone和Sidechain代表主链和侧链原子的贡献。Ratio列出了侧链表面积与每个残基的“无规卷曲”值的比率。图2中列出的20种氨基酸的“随机卷曲”值，如果比率值超过50%，则被认为是溶剂暴露，如果比率低于20%，则被掩埋，在In/Out栏中分别标记为“o”和“i”。在处理此功能时，我将‘o’更改为0，将‘i’更改为1，并将None值设置为0.5。
-
-
使用blosum62 乘sasa。
当我们使用SASA特征时，相同位置的所有突变都被赋予相同的SASA值，这意味着所有突变是等价的。为了区分同一位置的不同突变，我们可以使用 1 − 1 1 + e − x s f 1-\frac{1}{1+e^{\frac{-x}{s_{f}}}} 1−1+esf−x1将它们缩放到[0，1]范围，然后乘以SASA值。
-
Transformer ESM。
为了将氨基酸序列转换为有意义的特征，我们将使用Facebook预先训练的蛋白质transform ESM。我们将每个训练和测试野型输入到我们的转换器中，并提取最后的隐藏层激活。对于每个蛋白质，其形状为
(1，len_Protein_seq，1280)。我们将保存完整d的embeddings和pooled embeddings以供使用。transform embedding的维度为1280。由于我们只有几千行列车数据，这些特征太多了，无法全部输入我们的XGB模型，否则会造成过拟合。此外，我们希望使用局部、21个相邻残基、池化和增量embedding。最后，利用RAPIDS PCA对嵌入数据进行降维处理。 -
其他特征。
我们可以使用biopython得到 stability, aromaticity, isoelectric 等特征，但这些特征帮助不大，public score 和private score只有小幅提升。
下面来介绍一下三个XGB模型。
模型
XGB1
该模型的public score为0.52132，private score为0.46642。
我使用上面提到的第一个数据集来训练我的第一个XGBoost。这个数据集是从Kaggle的训练集中选择的野生型和突变体，有4000多行训练数据。首先，我对数据集做了一些简单的处理，去除了重复行和删除突变，删除了没有匹配PDB的行，保留了突变样本>20的组。最后，我们得到了一个1800+行的训练集。
利用上述特征，以rankdata™为目标对模型进行训练。
然后，我们使用WEIGHT和GAIN特征重要性来获得前30个重要特征。我们可以看到，ESM提取特征，B_factor、相对距离等更为重要，但我最终保留了所有的特征，因为我没有一个可靠的交叉验证集来选择特征。
XGB2
该模型的public score为0.40512，private score为0.30495。
这个模型产生了严重的过拟合，实际上在ensemble后降低了我的最终得分。
我们使用上面提到的第二个数据集 来训练我的第二个XgBoost模型。这个数据集是从Kaggle的训练集中选择的野生型和突变体。但这个数据集被清理得更详细，清理后的数据集只包含来自30个组的986个突变，所有组都是同源、相同pH的。
我们使用的功能和目标与XGB1相同。以下是功能的重要性：
XGB3
该模型的公开得分为0.40435，私人得分为0.43294。
我们使用第三个数据集 来训练我的第三个XgBoost模型。由于dTm和ddg具有很强的相关性，我们也可以使用ddg作为目标，以增加训练数据。这个外部数据集包含14656个独特的突变。由于野生型比大多数突变体更稳定，我们只保留了不稳定的突变，并确保每组的pH值和来源都是相同的。经过简单的清理，训练集包含6000多行样本。
与前两个模型不同，我只使用了一些结构特征来训练这个模型，因为我发现单独使用这些特征可以获得很好的public score。为了防止严重的过拟合和增加多样性，我没有添加其他特征。
以下是特征重要性，可以发现平均相对距离是很重要的。
Ensemble
大部分ensemble来自于公共笔记本，下面我简单介绍一下：
NESP 3D geometry
这个笔记本是由Elias创建的。
Public score: 0.32046, Private score: 0.33897。
此笔记本没有使用训练集，直接利用创建的3D几何特征应用于测试集。这种方法在这次比赛中是很常见的。
我们将以下两个位置之间的距离进行rankdata：
- 发生突变的残基的中心坐标。
- 蛋白质中心。
pdb 文件(https://www.kaggle.com/datasets/roberthatch/nesp-kvigly-test-mutation-pdbs)包含所有的原子坐标，我们只是取平均值来获得蛋白质中心坐标。为了得到突变残基的坐标，我们可以按残基对原子进行分组，使用平均值。如果突变是缺失的，我们使用下一个残基位置。
D i s t = ( x m u t − x c e n ) 2 + ( y m u t − y c e n ) 2 + ( z m u t − z c e n ) 2 Dist = \sqrt{(x_{mut}-x_{cen})^{2}+(y_{mut}-y_{cen})^{2}+(z_{mut}-z_{cen})^{2}} Dist=(xmut−xcen)2+(ymut−ycen)2+(zmut−zcen)2
FoldX
这个笔记本是我创建的。
Public score: 0.41332, Private score: 0.37538。
我从FoldX 网站(https://foldxsuite.crg.eu/documentation) 下载了Foldx5Win32.Zip。按照文档说明，我首先运行命令来修复野生型PDB。
然后，我创建了.cfg文件并运行了positionscan命令。根据该文件的指导，.cfg的格式如下：
命令=位置扫描。
pdb=wildtype_structure_prediction_af2_Repair.pdb。
位置=LA17E、LA17K、KA18C、KA18F…。(野生型+链+位置+突变)。
最后，我们将在完成时获得repair_scanning_output.txt。文件的第二列是ddg。我将其与测试集合。
ThermoNet v2
这个笔记(https://www.kaggle.com/code/vslaykovsky/nesp-thermonet-v2/notebook) 由Vladimir Slaykovskiy创建。
Public score: 0.49424, Private score: 0.46078。
这个模型只使用不稳定的突变进行训练。因为训练集同时包含ddg和dTm目标，并且ddg目标多于dTm目标，所以我们可以使用dTm作为辅助目标。loss使用的是
L = ( y d d G − y ^ d d G ) 2 + C ∗ ( y Δ d T m − y ^ d T m ) 2 \mathcal{L}=(y_{ddG}−\widehat{y}_{ddG})^2+C∗(y_{ΔdTm}−\widehat{y}_{dTm})^2 L=(yddG−y ddG)2+C∗(yΔdTm−y dTm)2， C = 0.01 C=0.01 C=0.01是一个决定dTm贡献的超参。
这本笔记本使用pdb(https://www.kaggle.com/code/vslaykovsky/14656-unique-mutations-voxel-features-pdbs) 与Acellera HTMD库生成体素特征。体素特征包括一组7个属性打包到一个大小(16，16，16)的网格中。七个特征分别是：'hydrophobic', 'aromatic', 'hbond_acceptor', 'hbond_donor', 'positive_ionizable', 'negative_ionizable', 'occupancies'对于每个训练样本，特征来自野生型和突变型PDB，因此每个样本的最终形状是：(14，16，16，16)。
最后用一个简单的VGG型CNN对回归模型进行训练。一个由10个模型组成的集合被用来预测分数。只有ddg输出用于预测，dTm仅作为损失函数的一部分用于训练。
模型架构如下图所示：
Relaxed rosetta scores
这个笔记本 由greySnow创建。
Public score: 0.471, Private score: 0.438。
这本笔记本使用relax结构 ，使用计算得到的的能量分数。Rosetta是一个软件套件，包括用于计算建模和分析蛋白质结构的算法(参见此处)。注册后即可免费从这里下载。此外，它还有一个很好的Python接口，即pyrosetta。我们可以用它来计算能量分数。
RMSD
这个笔记 由scar Villarreal Escamilla创建。
Public score: 0.39329, Private score: 0.37927。
本笔记本根据通过VMD和NAMD软件进行的分子动力学(MD)模拟的轨迹计算出的每个残基的平均均方根偏差(RMSD)。
VMD软件允许进行许多有趣的分析形式。在这场比赛中，均方根偏差(RMSD)似乎是最相关的分析，因为它直接衡量残基结构的稳定性。基本上，如果残基的RMSD非常小(即它移动很少)，这意味着该残基在其当前位置非常稳定，因此如果它发生突变，新突变的残基可能不像原始残基那样稳定。因此，这种突变残基更有可能具有较低的热稳定性等级。
Different pLDDT
这个笔记 由Chris Deotte创建。
Public score: 0.29709, Private score: 0.28793。
这本笔记本我们利用这里的pdb创建了不同的特征，对于测试集的每个突变，我们使用突变型的pLDDT减去野生型pLDDT。
我的最终解决方案是三个XGB模型和上面这些公共笔记本进行ensemble。下表是ensemble的权重：
| Model | Weight |
|---|---|
| NESP 3D geometry | 1 |
| Foldx | 1.5 |
| Thermonet v2 | 1.5 |
| Relaxed rosetta scores | 2 |
| RMSD | 1.5 |
| XGB3 | 1 |
| XGB1 | 1.5 |
| Different pLDDT | 1.5 |
| XGB2 | 1.5 |
总结
第一点是利用三维坐标来构建突变残基的环境，例如相邻残基和突变残基之间的相对距离。第二点是，现有的服务器或模型，如Rosetta、FoldX、ESM等，在这场比赛中取得了不错的成绩。第三点是找到方法来区分野生型和突变型。在这次比赛中，很多公共笔记本都使用了野生型，但注意野生型和突变型的区别可能也会有所帮助。这个方案有很多不成熟的地方，感谢各位批评指正。
参考资料：
[1]蛋白质三级结构嵌入编码及其在蛋白质工程中的应用研究
[2] How to use the information in PDB file
[3] Important Clarification Regarding b_factor
[4] NOVO ESP – ELI5 - Performant Approaches LB=0.451
[5] Surface area of the amino acids in the model structure
[6] XGBoost - 5000 Mutations 200 PDB Files
[7] Accessible surface area scaled by BLOSUM62 substitutions
[8] NOVO-ESP – Residue Depth and More w/ BioPython
[9] NESP 3D geometry
[10] ThermoNet v2
[11] nesp-relaxed-rosetta-scores
[12] rmsd-from-molecular-dynamics
[13] deletion-specific-ensemble
[14] difference-features-lb-0-600
