对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 pddef 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_experimentalprotein_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 = Falseround_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 pddef 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在定向进化中可以发挥更强大的能力。
