使用EVOLVEpro对蛋白质热稳定性进行定向进化
对EVOLVEpro还不熟的读者可以查看之前的博客:Rapid in silico directed evolution by a protein language model with EVOLVEpro 文献解读。
本文通过Kaggle竞赛中提供的酶热稳定性数据集对EVOLVEpro的性能进行验证。
数据集为test.csv和test_labels.csv,这是经过实验测量的某种酶的各个突变体的tm值。
野生型的酶序列
VPVNPEPDATSVENVALKTGSGDSQSDPIKADLEVKGQSALPFDVDCWAILCKGAPNVLQRVNEKTKNSNRDRSGANKGPFKDPQKWGIKALPPKNPSWSAQDFKSPEEYAFASSLQGGTNAILAPVNLASQNSQGGVLNGFYSANKVAQFDPSKPQQTKGTWFQITKFTGAAGPYCKALGSNDKSVCDKNKNIAGDWGFDPAKWAYQYDEKNNKFNYVGK
生成突变体文库,生成的突变体文库的大小为4200
from evolvepro.src.process import generate_wt, generate_single_aa_mutants
generate_wt('VPVNPEPDATSVENVALKTGSGDSQSDPIKADLEVKGQSALPFDVDCWAILCKGAPNVLQRVNEKTKNSNRDRSGANKGPFKDPQKWGIKALPPKNPSWSAQDFKSPEEYAFASSLQGGTNAILAPVNLASQNSQGGVLNGFYSANKVAQFDPSKPQQTKGTWFQITKFTGAAGPYCKALGSNDKSVCDKNKNIAGDWGFDPAKWAYQYDEKNNKFNYVGK', output_file='output/kaggle_WT.fasta')
generate_single_aa_mutants('output/kaggle_WT.fasta', output_file='output/kaggle.fasta')
构建突变图文库序列的embedding,PLM使用esm1b_t33_650M_UR50S
python evolvepro/plm/esm/extract.py esm1b_t33_650M_UR50S output/kaggle.fasta output/kaggle_esm1b_t33_650M_UR50S --toks_per_batch 512 --include mean --concatenate_dir output
随机选择12个突变体,作为第一轮进化的训练集
从kaggle test.csv和test_labels.csv找出相应的实验测量数据,其中有3个突变体没有找到相应的实验测量值,所以第一轮训练集数据大小为9.
from evolvepro.src.process import suggest_initial_mutants
suggest_initial_mutants('output/kaggle.fasta', num_mutants=12, random_seed=42)
suggest = [
"L92R",
"L116N",
"A91W",
"Y176N",
"A16Q",
"P97Q",
"K212H",
"T19K",
"P175K",
"A130I",
"A55E",
"I123A",
]
import pandas as pd
def mutation(seq):
for i in range(len(seq)):
if seq[i] != WT[i]:
return WT[i] + str(i+1) + seq[i]
WT = 'VPVNPEPDATSVENVALKTGSGDSQSDPIKADLEVKGQSALPFDVDCWAILCKGAPNVLQRVNEKTKNSNRDRSGANKGPFKDPQKWGIKALPPKNPSWSAQDFKSPEEYAFASSLQGGTNAILAPVNLASQNSQGGVLNGFYSANKVAQFDPSKPQQTKGTWFQITKFTGAAGPYCKALGSNDKSVCDKNKNIAGDWGFDPAKWAYQYDEKNNKFNYVGK'
test = pd.read_csv('test.csv', index_col=0)
test_labels = pd.read_csv('test_labels.csv', index_col=0)
df = pd.concat([test, test_labels], axis=1)[['protein_sequence', 'tm']]
df['length'] = df.apply(lambda x: len(x['protein_sequence']), axis=1)
df = df.loc[df['length'] == 221, ['protein_sequence', 'tm']]
df['mutation'] = df.apply(lambda x: mutation(x['protein_sequence']), axis=1)
df = df.loc[df['mutation'].isin(suggest), :]
df['Variant'] = df.apply(lambda x: x['mutation'][1:], axis=1)
df = df[['Variant', 'tm']]
df.columns = ['Variant', 'activity']
df.to_excel('kaggle/Round1.xlsx', index=False)
进行第一轮定向进化
from evolvepro.src.evolve import evolve_experimental
protein_name = 'kaggle'
embeddings_base_path = 'output'
embeddings_file_name = 'kaggle_esm1b_t33_650M_UR50S.csv'
round_base_path = 'kaggle'
wt_fasta_path = "output/kaggle_WT.fasta"
number_of_variants = 12
output_dir = 'output'
rename_WT = False
round_name = 'Round1'
round_file_names = ['Round1.xlsx']
this_round_variants, df_test, df_sorted_all = evolve_experimental(
protein_name,
round_name,
embeddings_base_path,
embeddings_file_name,
round_base_path,
round_file_names,
wt_fasta_path,
rename_WT,
number_of_variants,
output_dir
)
选择预测值最高的12个突变体,进行第二轮进化
从kaggle test.csv和test_labels.csv找出相应的实验测量数据,其中有5个突变体没有找到相应的实验测量值,所以第二轮训练集数据大小为9 + 7.
P94Q
P93S
P42Y
G88M
N133H
K78H
P80H
V45H
Q158M
D32I
N77H
G118K
选择预测值最高的12个突变体,进行第三轮进化
从kaggle test.csv和test_labels.csv找出相应的实验测量数据,其中有2个突变体没有找到相应的实验测量值,所以第三轮训练集数据大小为9 + 7 + 10.
Q150M
A125M
Q117D
F164W
N140Q
N14Q
P93M
T159F
S115D
F112W
N214Q
K160H
选择预测值最高的12个突变体,进行第四轮进化
从kaggle test.csv和test_labels.csv找出相应的实验测量数据,其中有7个突变体没有找到相应的实验测量值,所以第四轮训练集数据大小为9 + 7 + 10 + 5.
Q158I
A55Q
N14K
P153W
A113T
Q158V
E211W
P107Q
Q150G
F81W
D201Q
D197M
选择预测值最高的12个突变体,进行第五轮进化
从kaggle test.csv和test_labels.csv找出相应的实验测量数据,其中有5个突变体没有找到相应的实验测量值,所以第四轮训练集数据大小为9 + 7 + 10 + 5 + 8.
K192Q
A55Q
D197A
P202Q
S115N
G37H
K178H
N193V
K192M
Q158L
A111H
I123H
最终我们分别选取模型预测和数据集中排名前100,50,10的的突变体,查看它们之间的交集
import pandas as pd
def mutation(seq):
for i in range(len(seq)):
if seq[i] != WT[i]:
return WT[i] + str(i+1) + seq[i]
WT = 'VPVNPEPDATSVENVALKTGSGDSQSDPIKADLEVKGQSALPFDVDCWAILCKGAPNVLQRVNEKTKNSNRDRSGANKGPFKDPQKWGIKALPPKNPSWSAQDFKSPEEYAFASSLQGGTNAILAPVNLASQNSQGGVLNGFYSANKVAQFDPSKPQQTKGTWFQITKFTGAAGPYCKALGSNDKSVCDKNKNIAGDWGFDPAKWAYQYDEKNNKFNYVGK'
pre_df = pd.read_csv('output/kaggle/Round5/df_sorted_all.csv', index_col=0)
test = pd.read_csv('test.csv', index_col=0)
test_labels = pd.read_csv('test_labels.csv', index_col=0)
df = pd.concat([test, test_labels], axis=1)[['protein_sequence', 'tm']]
df['length'] = df.apply(lambda x: len(x['protein_sequence']), axis=1)
df = df.loc[df['length'] == 221, ['protein_sequence', 'tm']]
df['mutation'] = df.apply(lambda x: mutation(x['protein_sequence']), axis=1)
df = df.sort_values(by='tm', ascending=False)
pre_100 = set(pre_df['variant'][:100])
real_100 = set(df['mutation'][:100])
print("top100的交集大小:", len(pre_100 & real_100))
pre_50 = set(pre_df['variant'][:50])
real_50 = set(df['mutation'][:50])
print("top50的交集大小:", len(pre_50 & real_50))
pre_10 = set(pre_df['variant'][:10])
real_10 = set(df['mutation'][:10])
print("top10的交集大小:", len(pre_10 & real_10))
top100的交集大小: 3
top50的交集大小: 2
top10的交集大小: 1
值得一提的是,实验中测定的tm值最高的突变体F112W正是模型预测的tm值最高的突变体,这表明EVOLVEpro对定向进化湿实验中具有一定的指导意义。
不过也可以看到,通过5轮进化,构建的训练集大小为9 + 7 + 10 + 5 + 8=39,在tm值top100,top50和top10预测中并没有太高的命中率。
在未来可以通过结合蛋白质结构的embedding或者更有效的隐变量表示方法,使得EVOLVEpro在定向进化中可以发挥更强大的能力。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 清华大学推出第四讲使用 DeepSeek + DeepResearch 让科研像聊天一样简单!
· 推荐几款开源且免费的 .NET MAUI 组件库
· 实操Deepseek接入个人知识库
· 易语言 —— 开山篇
· Trae初体验