使用EVOLVEpro對蛋白質熱穩定性進行定向進化

crazypigf發表於2024-11-25

對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在定向進化中可以發揮更強大的能力。

相關文章