基于深度学习的在线小分子Kinome选择性预测平台的Python实现详解

基于深度学习的在线小分子Kinome选择性预测平台的Python实现详解

本文将详细介绍如何使用Python实现KinomeProDL，一个基于深度学习的多任务神经网络模型，用于预测小分子对激酶的选择性谱。我们将从数据准备、模型构建、模型训练、模型评估和模型应用等方面，逐步讲解实现过程，力求每一个方面都考虑到。

目录

1. 环境准备
2. 数据准备
   • 数据收集
   • 数据清洗和去重
   • 分子表示转换
3. 模型构建
   • 消息传递神经网络（MPNN）的实现
   • 多任务学习架构
4. 模型训练
   • 损失函数定义
   • 模型训练流程
   • 训练超参数设置
5. 模型评估
   • 评估指标计算
   • 模型性能可视化
6. 模型应用
   • 新化合物的预测
   • 模型微调
7. 总结

1. 环境准备

首先，确保您的计算环境中安装了以下主要库：

• Python 3.7+
• PyTorch
• PyTorch Geometric（用于图神经网络）
• RDKit（用于化学信息学处理）
• scikit-learn（用于模型评估）

安装PyTorch和PyTorch Geometric的命令如下：

# 安装PyTorch
pip install torch torchvision

# 安装PyTorch Geometric及其依赖
pip install torch-scatter torch-sparse torch-cluster torch-spline-conv torch-geometric

安装RDKit：

# 如果使用Anaconda
conda install -c conda-forge rdkit

2. 数据准备

2.1 数据收集

我们需要收集涵盖多个激酶和大量小分子化合物的活性数据。可以从以下公开数据库获取数据：

• ChEMBL
• BindingDB
• PubChem BioAssay
• Davis Dataset
• Metz Dataset
• SARfari

由于实际数据量较大，本文以示例数据集为例。

2.2 数据清洗和去重

数据清洗的主要步骤包括：

• 去重：移除重复的化合物和实验数据。
• 数据筛选：根据实验条件、数据可靠性等筛选有效数据。
• 缺失值处理：处理缺失的实验值或化合物结构信息。

以下是数据清洗的示例代码：

import pandas as pd

# 读取数据
data = pd.read_csv('kinase_activity_data.csv')

# 去重
data.drop_duplicates(subset=['compound_id', 'kinase_id'], inplace=True)

# 移除缺失值
data.dropna(subset=['activity_value', 'smiles'], inplace=True)

# 数据筛选，例如只保留IC50值小于10μM的活性数据
data = data[data['activity_value'] <= 10000]

2.3 分子表示转换

将SMILES字符串转换为分子图，以供图神经网络使用。使用RDKit将SMILES转换为分子对象，然后构建图表示。

from rdkit import Chem
from torch_geometric.data import Data
from rdkit.Chem import AllChem

def mol_to_graph_data_obj(smiles):
    mol = Chem.MolFromSmiles(smiles)
    # 添加氢原子
    mol = Chem.AddHs(mol)
    # 计算2D坐标（可选）
    AllChem.Compute2DCoords(mol)

    # 获取原子特征
    atom_features = []
    for atom in mol.GetAtoms():
        atom_feature = get_atom_features(atom)
        atom_features.append(atom_feature)
    
    # 获取边（键）信息
    edge_index = []
    edge_attrs = []
    for bond in mol.GetBonds():
        start = bond.GetBeginAtomIdx()
        end = bond.GetEndAtomIdx()
        edge_index.append((start, end))
        edge_index.append((end, start))  # 无向图
        bond_feature = get_bond_features(bond)
        edge_attrs.append(bond_feature)
        edge_attrs.append(bond_feature)

    # 构建PyTorch Geometric的数据对象
    data = Data(x=torch.tensor(atom_features, dtype=torch.float),
                edge_index=torch.tensor(edge_index, dtype=torch.long).t().contiguous(),
                edge_attr=torch.tensor(edge_attrs, dtype=torch.float))
    return data

def get_atom_features(atom):
    # 自定义原子特征提取函数
    features = []
    features.append(atom.GetAtomicNum())
    features.append(atom.GetDegree())
    features.append(atom.GetFormalCharge())
    features.append(atom.GetHybridization())
    features.append(atom.GetIsAromatic())
    return features

def get_bond_features(bond):
    # 自定义键特征提取函数
    features = []
    features.append(bond.GetBondTypeAsDouble())
    features.append(bond.IsInRing())
    return features

3. 模型构建

3.1 消息传递神经网络（MPNN）的实现

使用PyTorch Geometric实现MPNN。下面是一个基本的MPNN层的实现：

import torch
import torch.nn as nn
from torch_geometric.nn import MessagePassing
from torch_geometric.utils import add_self_loops

class MPNNLayer(MessagePassing):
    def __init__(self, in_channels, out_channels):
        super(MPNNLayer, self).__init__(aggr='add')  # 使用加法聚合
        self.linear = nn.Linear(in_channels, out_channels)

    def forward(self, x, edge_index, edge_attr):
        # 添加自环
        edge_index, edge_attr = add_self_loops(edge_index, num_nodes=x.size(0))
        return self.propagate(edge_index, x=x, edge_attr=edge_attr)

    def message(self, x_i, x_j, edge_attr):
        # 消息函数，x_i为目标节点特征，x_j为源节点特征
        return self.linear(x_j)

3.2 多任务学习架构

为实现多任务学习，我们在模型的最后一层为每个激酶任务设置一个输出节点。

class KinomeProDLModel(nn.Module):
    def __init__(self, num_features, num_tasks, hidden_dim=128, num_layers=3):
        super(KinomeProDLModel, self).__init__()
        self.convs = nn.ModuleList()
        self.convs.append(MPNNLayer(num_features, hidden_dim))
        for _ in range(num_layers - 1):
            self.convs.append(MPNNLayer(hidden_dim, hidden_dim))
        self.global_pool = nn.GlobalAveragePooling1D()
        self.fc = nn.Linear(hidden_dim, num_tasks)

    def forward(self, data):
        x, edge_index, edge_attr = data.x, data.edge_index, data.edge_attr
        for conv in self.convs:
            x = conv(x, edge_index, edge_attr)
            x = nn.ReLU()(x)
        x = self.global_pool(x)
        out = self.fc(x)
        return out

注意，这里的num_tasks是激酶的数量，模型的输出是一个大小为num_tasks的向量，每个元素对应一个激酶的预测值。

4. 模型训练

4.1 损失函数定义

由于我们要同时预测多个激酶的活性，可以使用多任务的二分类交叉熵损失。

criterion = nn.BCEWithLogitsLoss()

4.2 模型训练流程

from torch_geometric.loader import DataLoader

# 将数据转换为PyTorch Geometric的数据对象列表
graph_data_list = []
for idx, row in data.iterrows():
    graph_data = mol_to_graph_data_obj(row['smiles'])
    graph_data.y = torch.tensor(row[target_columns].values, dtype=torch.float)
    graph_data_list.append(graph_data)

# 划分训练集和测试集
from sklearn.model_selection import train_test_split

train_data, test_data = train_test_split(graph_data_list, test_size=0.2, random_state=42)

# 创建数据加载器
batch_size = 32
train_loader = DataLoader(train_data, batch_size=batch_size, shuffle=True)
test_loader = DataLoader(test_data, batch_size=batch_size)

# 模型实例化
num_features = graph_data_list[0].x.shape[1]
num_tasks = len(target_columns)
model = KinomeProDLModel(num_features=num_features, num_tasks=num_tasks)

# 优化器
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)

# 训练循环
num_epochs = 100

for epoch in range(num_epochs):
    model.train()
    total_loss = 0
    for batch in train_loader:
        optimizer.zero_grad()
        out = model(batch)
        loss = criterion(out, batch.y)
        loss.backward()
        optimizer.step()
        total_loss += loss.item() * batch.num_graphs
    avg_loss = total_loss / len(train_loader.dataset)
    print(f'Epoch {epoch+1}/{num_epochs}, Loss: {avg_loss:.4f}')

4.3 训练超参数设置

可以根据实际情况调整以下超参数：

• 学习率（lr）
• 批量大小（batch_size）
• 隐藏层维度（hidden_dim）
• 模型层数（num_layers）
• 训练轮数（num_epochs）

5. 模型评估

5.1 评估指标计算

使用AUROC（ROC曲线下面积）等指标评估模型性能。

from sklearn.metrics import roc_auc_score

model.eval()
all_preds = []
all_labels = []
with torch.no_grad():
    for batch in test_loader:
        out = model(batch)
        preds = torch.sigmoid(out).cpu().numpy()
        labels = batch.y.cpu().numpy()
        all_preds.append(preds)
        all_labels.append(labels)

all_preds = np.concatenate(all_preds, axis=0)
all_labels = np.concatenate(all_labels, axis=0)

# 计算每个任务（激酶）的AUROC
auroc_list = []
for i in range(num_tasks):
    auroc = roc_auc_score(all_labels[:, i], all_preds[:, i])
    auroc_list.append(auroc)
    print(f'Kinase {i+1}, AUROC: {auroc:.4f}')

5.2 模型性能可视化

可以绘制ROC曲线、Precision-Recall曲线等，以更直观地展示模型性能。

import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc

# 示例：绘制第一个激酶的ROC曲线
fpr, tpr, _ = roc_curve(all_labels[:, 0], all_preds[:, 0])
roc_auc = auc(fpr, tpr)

plt.figure()
plt.plot(fpr, tpr, label=f'ROC curve (area = {roc_auc:.2f})')
plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve for Kinase 1')
plt.legend(loc='lower right')
plt.show()

6. 模型应用

6.1 新化合物的预测

使用训练好的模型对新化合物进行预测。

# 新化合物的SMILES列表
new_smiles_list = ['CCO...', 'CCN...']

for smiles in new_smiles_list:
    graph_data = mol_to_graph_data_obj(smiles)
    graph_data = graph_data.to(device)  # 如果使用GPU
    model.eval()
    with torch.no_grad():
        out = model(graph_data)
        preds = torch.sigmoid(out).cpu().numpy()
    print(f'Predictions for {smiles}:')
    for i in range(num_tasks):
        print(f'  Kinase {i+1}: {preds[i]:.4f}')

6.2 模型微调

如果有新的私有数据集，可以对模型进行微调，以提高特定激酶的预测精度。

# 假设有新的数据new_data，格式与之前的数据相同
new_graph_data_list = []
for idx, row in new_data.iterrows():
    graph_data = mol_to_graph_data_obj(row['smiles'])
    graph_data.y = torch.tensor(row[target_columns].values, dtype=torch.float)
    new_graph_data_list.append(graph_data)

# 创建新的数据加载器
new_train_loader = DataLoader(new_graph_data_list, batch_size=batch_size, shuffle=True)

# 微调模型
for epoch in range(fine_tune_epochs):
    model.train()
    total_loss = 0
    for batch in new_train_loader:
        optimizer.zero_grad()
        out = model(batch)
        loss = criterion(out, batch.y)
        loss.backward()
        optimizer.step()
        total_loss += loss.item() * batch.num_graphs
    avg_loss = total_loss / len(new_train_loader.dataset)
    print(f'Fine-tuning Epoch {epoch+1}/{fine_tune_epochs}, Loss: {avg_loss:.4f}')

7. 总结

通过以上步骤，我们详细介绍了如何使用Python实现KinomeProDL模型，包括数据准备、模型构建、模型训练、模型评估和模型应用。该模型利用消息传递神经网络（MPNN）和多任务学习架构，能够有效地预测小分子化合物对多个激酶的选择性。

在实际应用中，可以根据具体需求和数据集，对模型进行调整和优化。例如，可以使用更复杂的原子和键特征，调整模型的层数和维度，或者使用更高级的图神经网络模型（如GAT、GraphSAGE等）。

此外，为了提高模型的泛化能力和稳定性，可以采用交叉验证、数据增强、正则化等技术。