2024 Datawhale AI夏令营机器学习方向
一、赛事详情
1.赛事链接
https://challenge.xfyun.cn/h5/detail?type=molecular-properties&ch=dw24_3Gp2Ih
2.赛事背景
在当今科技日新月异的时代，人工智能（AI）技术正以前所未有的深度和广度渗透到科研领域，特别是在化学及药物研发中展现出了巨大潜力。精准预测分子性质有助于高效筛选出具有优异性能的候选药物。以PROTACs为例，它是一种三元复合物由目标蛋白配体、linker、E3连接酶配体组成，靶向降解目标蛋白质。本次大赛聚焦于运用先进的人工智能算法预测其降解效能，旨在激发参赛者创新思维，推动AI技术与化学生物学的深度融合，进一步提升药物研发效率与成功率，为人类健康事业贡献智慧力量。通过此次大赛，我们期待见证并孵化出更多精准、高效的分子性质预测模型，共同开启药物发现的新纪元。
3.赛事任务
选手根据提供的demo数据集，可以基于demo数据集进行数据增强、自行搜集数据等方式扩充数据集，并自行划分数据。运用深度学习、强化学习或更加优秀人工智能的方法预测PROTACs的降解能力，若DC50>100nM且Dmax<80% ，则视为降解能力较差（demo数据集中Label=0）；若DC50<=100nM或Dmax>=80%，则视为降解能力好（demo数据集中Label=1）。
4.评审规则
本模型依据提交的结果文件，采用F1-score进行评价。
本赛题均提供下载数据，选手在本地进行算法调试，在比赛页面提交结果。
每支团队每天最多提交3次。
得分从高到低排序，排行榜将选择团队的历史最优成绩进行排名。
二、零基础baseline
# 导入需要用到的相关库
import pandas as pd
import numpy as np
from lightgbm import LGBMClassifier
# 读取训练集和测试集
train = pd.read_excel('./dataset-new/traindata-new.xlsx')
test = pd.read_excel('./dataset-new/testdata-new.xlsx')
# 特征工程
train = train.drop(['DC50 (nM)', 'Dmax (%)'], axis=1)
# 将object类型的数据进行目标编码处理
for col in train.columns[2:]:
if train[col].dtype == object or test[col].dtype == object:
train[col] = train[col].isnull()
test[col] = test[col].isnull()
# 加载决策树模型进行训练
model = LGBMClassifier(verbosity=-1)
model.fit(train.iloc[:, 2:].values, train['Label'])
pred = model.predict(test.iloc[:, 1:].values, )
# 保存结果文件到本地
pd.DataFrame(
{
'uuid': test['uuid'],
'Label': pred
}
).to_csv('submit.csv', index=None)
baseline采用决策树模型进行训练，使用的是lightgbm的LGBMClassifier。
其主要参数：
LGBMClassifier(
boosting_type='gbdt', # 弱学习器类型，默认为 'gbdt'
num_leaves=31, # 叶子节点数目，默认为31
max_depth=-1, # 树的最大深度，默认为-1，表示不限制
learning_rate=0.1, # 学习率，默认为0.1
n_estimators=100, # 弱学习器（决策树）数量，默认为100
objective='binary', # 目标函数，默认为二分类
min_split_gain=0.0, # 分裂节点时的最小增益，默认为0
min_child_samples=20, # 叶子节点最小样本数，默认为20，用于防止过拟合
subsample=1.0, # 每次迭代时用于训练的子样本比例，默认为1.0
subsample_freq=0, # 子样本的频率，默认为0，表示不使用子样本
colsample_bytree=1.0, # 每棵树的特征采样比例，默认为1.0
reg_alpha=0.0, # L1正则化参数，默认为0
reg_lambda=0.0, # L2正则化参数，默认为0
random_state=None, # 随机数种子，默认为None
silent=True # 是否静默模式，默认为True，不打印训练过程信息
)
三、赛题深入解析
赛题字段理解
UUID: 一个唯一的标识符，用于识别数据记录。
Label: 可能指的是化合物或数据集的标签或名称。
Uniprot: 一个蛋白质数据库，提供关于蛋白质序列、结构和功能的信息。
Target: 目标蛋白，即PROTACs设计来降解的特定蛋白质。
E3 ligase: E3连接酶，一种在泛素化过程中扮演角色的酶，帮助标记蛋白质以供降解。
Smiles: 一种表示分子结构的字符串格式。
Assay: 实验方法，这里可能指的是用于测定DC50、Dmax等的实验。
Article DOI: 文章的数字对象标识符，用于唯一标识科学文献。
Molecular Weight, Exact Mass: 分子量和精确质量。
XLogP3: 预测的辛醇/水分配系数，用于预测脂溶性。
Heavy Atom Count, Ring Count: 重原子数和环的数量。
Hydrogen Bond Acceptor/Donor Count: 氢键受体/供体计数。
Rotatable Bond Count: 可旋转键的数量。
Topological Polar Surface Area: 拓扑极性表面积。
Molecular Formula, InChI, InChI Key: 分子式、国际化学标识符及其密钥。
四、进阶baseline
# 导入库
import numpy as np
import pandas as pd
from catboost import CatBoostClassifier
from sklearn.model_selection import StratifiedKFold, KFold, GroupKFold
from sklearn.metrics import f1_score
from rdkit import Chem
from rdkit.Chem import Descriptors
from sklearn.feature_extraction.text import TfidfVectorizer
import tqdm, sys, os, gc, re, argparse, warnings
warnings.filterwarnings('ignore')
train = pd.read_excel('./dataset-new/traindata-new.xlsx')
test = pd.read_excel('./dataset-new/testdata-new.xlsx')
# test数据不包含 DC50 (nM) 和 Dmax (%)
train = train.drop(['DC50 (nM)', 'Dmax (%)'], axis=1)
# 定义了一个空列表drop_cols，用于存储在测试数据集中非空值小于10个的列名。
drop_cols = []
for f in test.columns:
if test[f].notnull().sum() < 10:
drop_cols.append(f)
# 使用drop方法从训练集和测试集中删除了这些列，以避免在后续的分析或建模中使用这些包含大量缺失值的列
train = train.drop(drop_cols, axis=1)
test = test.drop(drop_cols, axis=1)
# 使用pd.concat将清洗后的训练集和测试集合并成一个名为data的DataFrame，便于进行统一的特征工程处理
data = pd.concat([train, test], axis=0, ignore_index=True)
cols = data.columns[2:]
# 将SMILES转换为分子对象列表,并转换为SMILES字符串列表
data['smiles_list'] = data['Smiles'].apply(lambda x:[Chem.MolToSmiles(mol, isomericSmiles=True) for mol in [Chem.MolFromSmiles(x)]])
data['smiles_list'] = data['smiles_list'].map(lambda x: ' '.join(x))
# 使用TfidfVectorizer计算TF-IDF
tfidf = TfidfVectorizer(max_df = 0.9, min_df = 1, sublinear_tf = True)
res = tfidf.fit_transform(data['smiles_list'])
# 将结果转为dataframe格式
tfidf_df = pd.DataFrame(res.toarray())
tfidf_df.columns = [f'smiles_tfidf_{i}' for i in range(tfidf_df.shape[1])]
# 按列合并到data数据
data = pd.concat([data, tfidf_df], axis=1)
# 自然数编码
def label_encode(series):
unique = list(series.unique())
return series.map(dict(zip(
unique, range(series.nunique())
)))
# 对每个类转换为其编码
for col in cols:
if data[col].dtype == 'object':
data[col] = label_encode(data[col])
# 提取data中label行不为空的，将其作为train的数据并更新索引
train = data[data.Label.notnull()].reset_index(drop=True)
# 提取data中label行为空的，将其作为teat的数据并更新索引
test = data[data.Label.isnull()].reset_index(drop=True)
# 特征筛选
features = [f for f in train.columns if f not in ['uuid','Label','smiles_list']]
# 构建训练集和测试集
x_train = train[features]
x_test = test[features]
# 训练集标签
y_train = train['Label'].astype(int)
def cv_model(clf, train_x, train_y, test_x, clf_name, seed=2022):
# 进行5折交叉验证
kf = KFold(n_splits=5, shuffle=True, random_state=seed)
train = np.zeros(train_x.shape[0])
test = np.zeros(test_x.shape[0])
cv_scores = []
# 每一折数据采用训练索引和验证索引来分割训练集和验证集
for i, (train_index, valid_index) in enumerate(kf.split(train_x, train_y)):
print('************************************ {} {}************************************'.format(str(i+1), str(seed)))
trn_x, trn_y, val_x, val_y = train_x.iloc[train_index], train_y[train_index], train_x.iloc[valid_index], train_y[valid_index]
# 配置CatBoost分类器的参数
params = {'learning_rate': 0.05, 'depth': 8, 'l2_leaf_reg': 10, 'bootstrap_type':'Bernoulli','random_seed':seed,
'od_type': 'Iter', 'od_wait': 100, 'random_seed': 11, 'allow_writing_files': False, 'task_type':'CPU'}
# 使用CatBoost分类器训练模型
model = clf(iterations=20000, **params, eval_metric='AUC')
model.fit(trn_x, trn_y, eval_set=(val_x, val_y),
metric_period=100,
cat_features=[],
use_best_model=True,
verbose=1)
val_pred = model.predict_proba(val_x)[:,1]
test_pred = model.predict_proba(test_x)[:,1]
train[valid_index] = val_pred
test += test_pred / kf.n_splits
cv_scores.append(f1_score(val_y, np.where(val_pred>0.5, 1, 0)))
print(cv_scores)
print("%s_score_list:" % clf_name, cv_scores)
print("%s_score_mean:" % clf_name, np.mean(cv_scores))
print("%s_score_std:" % clf_name, np.std(cv_scores))
return train, test
cat_train, cat_test = cv_model(CatBoostClassifier, x_train, y_train, x_test, "cat")
pd.DataFrame(
{
'uuid': test['uuid'],
'Label': np.where(cat_test>0.5, 1, 0)
}
).to_csv('submit.csv', index=None)
五、数据分析
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
train = pd.read_excel('./dataset-new/traindata-new.xlsx')
test = pd.read_excel('./dataset-new/testdata-new.xlsx')
train = train.drop(['DC50 (nM)', 'Dmax (%)'], axis=1)
drop_cols = []
for f in test.columns:
if test[f].notnull().sum() < 300:
drop_cols.append(f)
train = train.drop(drop_cols, axis=1)
test = test.drop(drop_cols, axis=1)
train.info()
可以看到有20列数据大于300行，我们简单对其进行分析，其中uuid不需要分析。
1.Label
data = train['Label']
lable0 = sum(data == 0)
lable1 = sum(data == 1)
print(lable0, lable1)
133 218
2.Uniprot
Uniprot是 一个蛋白质数据库，提供关于蛋白质序列、结构和功能的信息。
data = train.loc[:,['Label','Uniprot']]
x1 = []
l0 = []
l1 = []
for index, row in data.iterrows():
if row['Uniprot'] not in x1:
x1.append(row['Uniprot'])
l0.append(0)
l1.append(0)
i = x1.index(row['Uniprot'])
if row['Label'] == 0:
l0[i] += 1
else:
l1[i] += 1
plt.figure(figsize=(10,6))
plt.xticks(rotation=90)
plt.plot(x1, l0)
plt.plot(x1, l1)
3.Target
Target: 目标蛋白，即PROTACs设计来降解的特定蛋白质。
# 分析label与target的关系
data1 = train.loc[:,['Label','Target']]
x1 = []
l0 = []
l1 = []
for index, row in data1.iterrows():
if row['Target'] not in x1:
x1.append(row['Target'])
l0.append(0)
l1.append(0)
i = x1.index(row['Target'])
if row['Label'] == 0:
l0[i] += 1
else:
l1[i] += 1
plt.figure(figsize=(10,6))
plt.xticks(rotation=90)
plt.plot(x1, l0)
plt.plot(x1, l1)
4.E3 ligase
E3 ligase: E3连接酶，一种在泛素化过程中扮演角色的酶，帮助标记蛋白质以供降解。
data2 = train.loc[:,['Label','E3 ligase']]
x1 = []
l0 = []
l1 = []
for index, row in data2.iterrows():
if row['E3 ligase'] not in x1:
x1.append(row['E3 ligase'])
l0.append(0)
l1.append(0)
i = x1.index(row['E3 ligase'])
if row['Label'] == 0:
l0[i] += 1
else:
l1[i] += 1
plt.figure(figsize=(10,6))
plt.plot(x1, l0)
plt.plot(x1, l1)
5.Smiles
Smiles: 一种表示分子结构的字符串格式。
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
smilesList = train['Smiles']
# SMILES字符串
smiles = smilesList[0]
# 将SMILES转换为分子对象
mol = Chem.MolFromSmiles(smiles)
# 检查分子是否有效
if mol is None:
print("The SMILES string is invalid.")
else:
print("The molecule was parsed successfully.")
# 绘制分子结构
molDrawer = Draw.MolToMPL(mol)
molDrawer.show()
这是label=1的某个分子
这是label=0的某个分子
其实很难分辨出它们的区别。
但是了解到可以使用指纹来计算分子的相似性
基础教程 — RDKit 中文教程 2020.09 文档 (chenzhaoqiang.com)
# 对tarin的值计算相似性
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit import DataStructs
smilesList = train['Smiles']
ans = []
equal = []
# 全部转换为拓扑指纹
for i in range(0, len(smilesList)):
smiles1 = smilesList[i]
mol1 = Chem.MolFromSmiles(smiles1)
m1 = Chem.RDKFingerprint(mol1)
ans.append(m1)
for i in range(0, len(ans)):
for j in range(i+1, len(ans)):
m1 = ans[i]
m2 = ans[j]
# 判断两个指纹的相似度
sm=DataStructs.FingerprintSimilarity(m1,m2)
if sm == 1:
equal.append({i, j})
equal
结果显示有三百多条数据是相似的，也就是相似度为1。当我把测试集的数据加进去后，有五百多条数据是相似的。那我们是否可以将相似分子赋值为同一个数据呢？
我采用上述方法再使用TfidfVectorizer对分子变量进行计算，可以发现结果确实有一定提高，但是差别不大。
6.Assay
Assay: 实验方法，这里可能指的是用于测定DC50、Dmax等的实验。
经统计，有近百种不同实验方法，个人认为对label的差异影响很小，此处不做分析。在取消使用该特征后，k折交叉验证的有略微变化。在baseline的基础上去除该值后，结果有所下降。
7.Article DOI
Article DOI: 文章的数字对象标识符，用于唯一标识科学文献。
DOI（Digital Object Identifier）, 中文译为“数字对象唯一标识符”，是一种用于唯一标识数字化对象的持久标识符。它由国际标准化组织（ISO）制定，并以国际标准 ISO 26324-2012 的形式发布。
DOI的格式由前缀和后缀两个部分组成，用分隔符“/”区分开。 我们以下面的doi举例说明：
10.1080/10286632.2013.832233
这里的 10 是数字对象标识符的固定代号， 1080 是doi注册机构。 而Doi的后缀是各个出版物自己来定义的，没有统一的规则。 在这里 “10286632” 是该期刊的ISSN号， 2013 代表出版年份。
8.Molecular Weight
Molecular Weight: 分子质量
# 分析label与Molecular Weight的关系
import seaborn as sns
data3 = train.loc[:,['Label','Molecular Weight']]
sns.violinplot(x='Label', y='Molecular Weight', data=data3)
可以看出分子质量差异不大。
9.Exact Mass
# 分析label与Exact Mass的关系
data4 = train.loc[:,['Label','Exact Mass']]
sns.violinplot(x='Label', y='Exact Mass', data=data4, scale='count')
同上，由于二者差异不大。这里在第7点Article DOI的基础上，删除8和9的特征。此时在得分下降了。
10.XLogP3
XLogP3: 预测的辛醇/水分配系数，用于预测脂溶性。
# 分析label与XLogP3 的关系
data5 = train.loc[:,['Label','XLogP3']]
sns.violinplot(x='Label', y='XLogP3', data=data5, scale='count')
11.Heavy Atom Count
Heavy Atom Count: 重原子数的数量。
# 分析label与Heavy Atom Count的关系
data6 = train.loc[:,['Label','Heavy Atom Count']]
sns.violinplot(x='Label', y='Heavy Atom Count', data=data6, scale='count')
12.Ring Count
Ring Count: 环的数量。
# 分析label与Ring Count 的关系
data7 = train.loc[:,['Label','Ring Count']]
sns.violinplot(x='Label', y='Ring Count', data=data7, scale='count')
Ring Count不同label的差异较大，适合做一个较好的特征。
13.Hydrogen Bond Acceptor Count
Hydrogen Bond Acceptor Count: 氢键受体计数
# 分析label与Hydrogen Bond Acceptor Count 的关系
data8 = train.loc[:,['Label','Hydrogen Bond Acceptor Count']]
sns.violinplot(x='Label', y='Hydrogen Bond Acceptor Count', data=data8, scale='count')
14.Hydrogen Bond Donor Count。
Hydrogen Bond Donor Count氢键供体计数。
# 分析label与Hydrogen Bond Donor Count 的关系
data9 = train.loc[:,['Label','Hydrogen Bond Donor Count']]
sns.violinplot(x='Label', y='Hydrogen Bond Donor Count', data=data9, scale='count')
15.Rotatable Bond Count
Rotatable Bond Count: 可旋转键的数量。
# 分析label与Rotatable Bond Count的关系
data10 = train.loc[:,['Label','Rotatable Bond Count']]
sns.violinplot(x='Label', y='Rotatable Bond Count', data=data10, scale='count')
16.Topological Polar Surface Area
Topological Polar Surface Area: 拓扑极性表面积。
# 分析label与Topological Polar Surface Area的关系
data11 = train.loc[:,['Label','Topological Polar Surface Area']]
sns.violinplot(x='Label', y='Topological Polar Surface Area', data=data11, scale='count')
17.InChi
Molecular Formula可以由InChi计算得到，故不对Molecular Formula进行分析，
InChI编码是一串由斜杠（/）隔开的有层级关系的数字组成的。每个InChI编码都是由InChI版本号开始，接着一个主层号。主层下包括含化学分子式层、原子关系层和固定氢原子子层等基于分子结构的层。此外，还有附加的层，如电荷层、立体化学层（和/或）同位素信息层等。且这些附加的层也有其各自的子层。
而InChI Key本身没有其意义，只是对 InChI 运用 SHA-256 算法处理后得到的哈希值，它的出现是为了解决 InChI 长度不定的问题。 故需要转换为字符串才可读取。
# 分析label和InChI
import re
# 函数用于解析单个InChI字符串
def parse_inchi(row):
inchi_str = row['InChI']
formula = ''
element_counts = {}
# 提取分子式
formula_match = re.search(r"InChI=1S/([^/]+)/c", inchi_str)
if formula_match:
formula = formula_match.group(1)
# 计算分子量和原子计数
for element, count in re.findall(r"([A-Z][a-z]*)([0-9]*)", formula):
count = int(count) if count else 1
element_counts[element.upper()] = count
return pd.Series({
'ElementCounts': element_counts
})
# 应用函数到DataFrame的每一行
train[['ElementCounts']] = train.apply(parse_inchi, axis=1)
# 定义存在的key
keys = ['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn']
# 创建一个空的DataFrame，列名为keys
df_expanded = pd.DataFrame({key: pd.Series(dtype='float64') for key in keys})
# 遍历数据，填充DataFrame
for index, item in enumerate(train['ElementCounts'].values):
for key in keys:
# 将字典中的值填充到相应的列中
df_expanded.at[index, key] = item.get(key, 0)
df_expanded = pd.DataFrame(df_expanded)
# 删除全为0的列
drop_cols = []
for f in df_expanded.columns:
if df_expanded[f].sum() == 0:
drop_cols.append(f)
df_expanded = df_expanded.drop(drop_cols, axis=1)
print(df_expanded.info())
# 加入Label
df_expanded = pd.concat([train['Label'], df_expanded], axis=1)
df_expanded
# 分析label与InChI包含元素H的关系
data = df_expanded.loc[:,['Label','H']]
sns.violinplot(x='Label', y='H', data=df_expanded, scale='count')
# 分析label与InChI包含元素C的关系
data = df_expanded.loc[:,['Label','C']]
sns.violinplot(x='Label', y='C', data=df_expanded, scale='count')
# 分析label与InChI包含元素N的关系
data = df_expanded.loc[:,['Label','N']]
sns.violinplot(x='Label', y='N', data=df_expanded, scale='count')
# 分析label与InChI包含元素O的关系
data = df_expanded.loc[:,['Label','O']]
sns.violinplot(x='Label', y='O', data=df_expanded, scale='count')
# 分析label与InChI包含元素F的关系
data = df_expanded.loc[:,['Label','F']]
sns.violinplot(x='Label', y='F', data=df_expanded, scale='count')
# 分析label与InChI包含元素P的关系
data = df_expanded.loc[:,['Label','P']]
sns.violinplot(x='Label', y='P', data=df_expanded, scale='count')
# 分析label与InChI包含元素S的关系
data = df_expanded.loc[:,['Label','S']]
sns.violinplot(x='Label', y='S', data=df_expanded, scale='count')
六、深度学习
提交了十几次，结果最好的居然是官方给的example.csv。于是很生气的（bushi）在网上借鉴了一个神经网络，替换了cv_model那部分。
x_train = torch.tensor(x_train.values, dtype=torch.float32)
x_test = torch.tensor(x_test.values, dtype=torch.float32)
y_train = torch.tensor(y_train.values, dtype=torch.float32)
from sklearn.model_selection import train_test_split
#拆分训练集和测试集
X_train, X_test, Y_train, Y_test = train_test_split(x_train,y_train, test_size = 0.10, random_state = 42)
print(X_train.shape,Y_train.shape)
print(X_test.shape,Y_test.shape)
class Model(nn.Module):
def __init__(self):
super().__init__()
self.liner_1 = nn.Linear(939, 256)
self.liner_2 = nn.Linear(256, 64)
self.liner_3 = nn.Linear(64, 1)
self.relu = nn.ReLU()
self.sigmoid = nn.Sigmoid()
def forward(self, input):
x = self.liner_1(input)
x = self.relu(x)
x = self.liner_2(x)
x = self.relu(x)
x = self.liner_3(x)
x = self.sigmoid(x)
return x
lr = 0.001
def get_model():
model = Model()
opt = torch.optim.Adam(model.parameters(), lr=lr)
return model, opt
model, optim = get_model()
loss_fn = nn.BCELoss()
batch = 64
no_of_batches = len(data)// batch
epochs = 300
from sklearn.metrics import precision_recall_fscore_support
def evaluate(model, X_val, Y_val, loss_fn):
model.eval()
with torch.no_grad():
y_pred = model(X_val)
Y_val = Y_val.view(-1, 1)
val_loss = loss_fn(y_pred, Y_val)
y_pred_binary = (y_pred > 0.5).float()
y_pred_binary = y_pred_binary.cpu().numpy().flatten()
Y_val = Y_val.cpu().numpy().flatten()
precision, recall, f1, _ = precision_recall_fscore_support(Y_val, y_pred_binary, average='binary')
return val_loss.item(), f1
f1_sum = 0
for epoch in range(epochs):
for i in range(no_of_batches):
start = i * batch
end = start + batch
x = X_train[start: end]
y = Y_train[start: end]
y_pred = model(x)
y = y.view(-1, 1)
loss = loss_fn(y_pred, y)
optim.zero_grad()
loss.backward()
optim.step()
with torch.no_grad():
validation_loss, f1 = evaluate(model, X_test, Y_test, loss_fn)
f1_sum += f1
print('epoch:', epoch, ' train_loss:', loss_fn(model(X_train), Y_train.view(-1, 1)).data.item(),
' test_loss:', validation_loss, ' f1:', f1)
print(f1_sum / epochs)
y_pred = model(x_test)
y_pred = (y_pred > 0.5).float()
y_pred
pd.DataFrame(
{
'uuid': test['uuid'],
'Label': torch.where(y_pred > 0.5, torch.tensor(1), torch.tensor(0)).numpy().flatten()
}
).to_csv('submit.csv', index=None)
直接刷新了最好的成绩，接下来就是枯燥的炼丹环节啦。
我们有缘再会。
