最近在kaggle参加了诺维信举办的酶稳定性预测比赛 ,最终很幸运获得了第三名,这篇文章主要是简单介绍一下解决方案,具体的数据和题目要求可访问上面的链接。
模型概述
我分别使用了如下三个数据集训练了三个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文章来源:https://www.toymoban.com/news/detail-419428.html
[14] difference-features-lb-0-600文章来源地址https://www.toymoban.com/news/detail-419428.html
到了这里,关于kaggle酶稳定性预测第三名解决方案分享的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!