生命科学siRNA药物药效预测之baseline（Datawhale AI 夏令营）

背景

        随着mRNA疫苗在新冠预防领域取得成功，核酸类药物的研发获得了越来越多的关注。本次比赛聚焦于通过机器学习技术，利用化学修饰后的siRNA序列来预测RNA干扰（RNAi）机制下对靶基因（target gene）的沉默效率。这一指标与药物实际疗效直接相关。RNAi是生物体内天然存在的一种基因表达调控机制，通过抑制靶基因的表达来实现降低目标蛋白量的目的，这一机制一般可通过siRNA实现。

        目前开源的数据库中，以RNA主干序列（裸序列）为主，缺少相应的化学修饰数据。而本赛题特别关注了化学修饰对siRNA序列功能的影响，化学修饰对siRNA的毒性、体内稳定性、靶向效果、药效等具有重大影响，在实际药物设计中至关重要。参赛者将接触到领域独特的包含靶基因、siRNA裸序列、经过化学修饰的siRNA修饰序列以及实验室测定的沉默效率值。这些数据反映了当今siRNA设计的最新科技成果，包括化学修饰的种类和位置，以及它们如何影响siRNA对靶基因的沉默效率。

任务

       目标是构建一个能够准确预测siRNA沉默效率的模型。提供了一部分公开文献中提取的siRNA修饰序列以及相应实验条件数据，使用这些数据训练模型并提交预测结果。

数据集介绍

        数据集包括四个文件：train_data.csv、sample_submission.csv、vocab.csv、baseline.py。

• train_data.csv：每行为一条训练记录，包含siRNA裸序列、修饰后的siRNA序列、目标mRNA序列、实验条件（如siRNA浓度、细胞系、转染方式）及对应的mRNA Remaining值等19个字段。mRNA Remaining值为模型的训练目标，表示siRNA对靶基因沉默后的剩余mRNA百分比，值越低表示沉默效率越好。
• sample_submission.csv：测试集，格式与train_data.csv相同，不同之处在于mRNA_remaining_pct列的数值为空，参赛者需要填充这些空白处的预测结果后提交。
• vocab.csv：经过归一化后的siRNA化学修饰缩写表，包含各种修饰核苷酸的缩写及其对应的完整化学名称。
• baseline.py：提供了一个基础的基线方法，利用RNN来预测Remaining值。基线方法中仅使用了siRNA_sense_seq字段作为特征，除此之外尚有其他特征对Remaining结果有重大影响。

评价指标

        本次任务采用三个指标来进行评测：

1. 平均绝对误差（MAE）：衡量模型预测值与真实值之间的平均绝对误差。

   mae = np.mean(np.abs(y_true - y_pred))
   

2. 预测值在一定范围内的平均绝对误差（Range-MAE）：在特定范围内的预测值的平均绝对误差。在本次比赛中，低Remaining范围为0,300, 300,30。

   mask = (y_pred >= 0) & (y_pred <= threshold)
   range_mae = mean_absolute_error(y_true[mask], y_pred[mask]) if mask.sum() > 0 else 100
   

3. 预测值在一定范围内的F1指标（F1）：衡量在特定范围内的预测值的precision以及在相应区间预测的recall情况，综合得到F1值。在本次比赛中，低Remaining范围为0,300, 300,30。

   precision = precision_score(y_true_binary, y_pred_binary, average='binary')
   recall = recall_score(y_true_binary, y_pred_binary, average='binary')
   f1 = 2 * precision * recall / (precision + recall)
   

   最终得分基于上述三个指标计算得到，计算公式：

score = 50% × (1−MAE/100) + 50% × F1 × (1−Range-MAE/100)

相关生物背景知识

RNA干扰（RNAi）

        RNA干扰（RNAi）是一种天然存在的基因表达调控机制，通过小干扰RNA（siRNA）等分子来沉默特定基因的表达。这一机制在细胞中起着重要作用，能精确地抑制目标基因的表达，从而减少相应蛋白质的产生。siRNA通过与靶mRNA结合，诱导RNA诱导沉默复合物（RISC）切割mRNA，最终导致mRNA降解和基因沉默。在基因治疗和疾病治疗中，RNAi技术有望通过沉默致病基因来发挥治疗作用。

化学修饰siRNA

        化学修饰siRNA是指在siRNA分子中引入化学修饰，以增强其稳定性、靶向性和有效性。这些修饰可以增加siRNA在体内的稳定性，减少其毒性和副作用，提高其基因沉默效率。常见的化学修饰包括磷酸酯骨架修饰、核苷酸修饰和末端修饰等。这些修饰不仅能提高siRNA的药效，还能减少非特异性沉默，提升siRNA药物的临床应用潜力。

深度学习与RNN

        深度学习是一种基于人工神经网络的机器学习方法，擅长处理复杂的非线性关系和高维数据。递归神经网络（RNN）是一类深度学习模型，特别适用于处理序列数据。RNN通过在隐藏层中引入循环连接，可以有效捕捉序列中的时间依赖关系。在RNAi效率预测任务中，RNN能够通过学习siRNA序列和靶mRNA序列之间的复杂关系，准确预测其基因沉默效果。

词汇表与序列编码

        在处理基因序列数据时，通常需要将核酸序列转换为数值表示形式，以便输入到深度学习模型中。词汇表（vocab）是一种将序列中的每个元素（如核苷酸或核苷酸组合）映射到一个唯一的数值索引的结构。在本文中，使用了一个基于3-gram的词汇表，这意味着每三个连续的核苷酸组合成一个“单词”。这种方法能够捕捉序列中的局部模式，并提高模型的预测能力。

数据处理与特征选择

        数据处理是机器学习中的关键步骤，包含数据清洗、预处理和特征选择。在本次比赛中，需要对原始数据进行清洗，去除缺失值和异常值。特征选择则是选择最能代表数据特征的字段，以提高模型的性能和训练效率。这里的特征包括siRNA序列、修饰后的siRNA序列、靶mRNA序列以及实验条件（如药物浓度、细胞系、转染方式等）。

模型训练与评估

        模型训练是指通过优化算法（如Adam优化器）调整模型参数，使其在训练数据上表现良好。评估模型性能时，常用的指标包括均方误差（MSE）、平均绝对误差（MAE）、精确率（Precision）和召回率（Recall）。在本次比赛中，模型的最终得分由MAE和预测值在一定范围内的F1指标（F1）综合计算而得。训练过程中需要监控验证集上的模型表现，以防止过拟合。

PyTorch框架

        PyTorch是一个开源的深度学习框架，广泛用于研究和生产环境中。它提供了灵活的动态计算图，使得模型的定义和训练更加直观和便捷。在本次比赛的baseline代码中，使用了PyTorch构建和训练RNN模型，包括数据加载、序列编码、模型定义、训练循环和评估等步骤。PyTorch的优势在于其简洁的API和强大的功能，能够快速实现复杂的深度学习模型。

代码实现

1. 依赖库的导入

import os  # 文件操作
import torch  # 深度学习框架
import random  # 随机数生成
import numpy as np  # 数值计算
import pandas as pd  # 数据处理

import torch.nn as nn  # 神经网络模块
import torch.optim as optim  # 优化器模块

from tqdm import tqdm  # 进度条显示
from rich import print  # 美化打印输出
from collections import Counter  # 计数器工具

from torch.utils.data import Dataset, DataLoader  # 数据集和数据加载器
from sklearn.model_selection import train_test_split  # 数据集划分
from sklearn.metrics import precision_score, recall_score, mean_absolute_error  # 模型评估指标

这些库包括了文件操作、深度学习、数据处理、模型评估等必要的工具。

2. 设置随机种子

def set_random_seed(seed):
    np.random.seed(seed)
    random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True
    torch.backends.cudnn.benchmark = False

该函数确保在使用NumPy、Python内置随机数生成器和PyTorch时，所有的随机数生成都是可控的和可复现的，有助于实验结果的一致性。

3. 基因组分词器类

class GenomicTokenizer:
    def __init__(self, ngram=5, stride=2):
        self.ngram = ngram
        self.stride = stride
        
    def tokenize(self, t):
        t = t.upper()
        if self.ngram == 1:
            toks = list(t)
        else:
            toks = [t[i:i+self.ngram] for i in range(0, len(t), self.stride) if len(t[i:i+self.ngram]) == self.ngram]
        if len(toks[-1]) < self.ngram:
            toks = toks[:-1]
        return toks

该类用于将基因组序列分割成固定长度的n-gram。

4. 基因组词汇类

class GenomicVocab:
    def __init__(self, itos):
        self.itos = itos
        self.stoi = {v: k for k, v in enumerate(self.itos)}
        
    @classmethod
    def create(cls, tokens, max_vocab, min_freq):
        freq = Counter(tokens)
        itos = ['<pad>'] + [o for o, c in freq.most_common(max_vocab - 1) if c >= min_freq]
        return cls(itos)

该类用于创建一个词汇表，用于将基因组片段映射为索引。

5. siRNA数据集类

class SiRNADataset(Dataset):
    def __init__(self, df, columns, vocab, tokenizer, max_len, is_test=False):
        self.df = df
        self.columns = columns
        self.vocab = vocab
        self.tokenizer = tokenizer
        self.max_len = max_len
        self.is_test = is_test

    def __len__(self):
        return len(self.df)

    def __getitem__(self, idx):
        row = self.df.iloc[idx]
        seqs = [self.tokenize_and_encode(row[col]) for col in self.columns]
        if self.is_test:
            return seqs
        else:
            target = torch.tensor(row['mRNA_remaining_pct'], dtype=torch.float)
            return seqs, target

    def tokenize_and_encode(self, seq):
        if ' ' in seq:
            tokens = seq.split()
        else:
            tokens = self.tokenizer.tokenize(seq)
        encoded = [self.vocab.stoi.get(token, 0) for token in tokens]
        padded = encoded + [0] * (self.max_len - len(encoded))
        return torch.tensor(padded[:self.max_len], dtype=torch.long)

该类用于加载siRNA数据，并将序列数据转换为模型可以处理的格式。

6. siRNA Model

class SiRNAModel(nn.Module):
    def __init__(self, vocab_size, embed_dim=200, hidden_dim=256, n_layers=3, dropout=0.5):
        super(SiRNAModel, self).__init__()
        self.embedding = nn.Embedding(vocab_size, embed_dim, padding_idx=0)
        self.gru = nn.GRU(embed_dim, hidden_dim, n_layers, bidirectional=True, batch_first=True, dropout=dropout)
        self.fc = nn.Linear(hidden_dim * 4, 1)
        self.dropout = nn.Dropout(dropout)
    
    def forward(self, x):
        embedded = [self.embedding(seq) for seq in x]
        outputs = []
        for embed in embedded:
            x, _ = self.gru(embed)
            x = self.dropout(x[:, -1, :])
            outputs.append(x)
        x = torch.cat(outputs, dim=1)
        x = self.fc(x)
        return x.squeeze()

这是一个基于GRU的神经网络模型，用于处理siRNA序列。

7. 评估指标计算函数

def calculate_metrics(y_true, y_pred, threshold=30):
    mae = np.mean(np.abs(y_true - y_pred))
    y_true_binary = (y_true < threshold).astype(int)
    y_pred_binary = (y_pred < threshold).astype(int)
    mask = (y_pred >= 0) & (y_pred <= threshold)
    range_mae = mean_absolute_error(y_true[mask], y_pred[mask]) if mask.sum() > 0 else 100
    precision = precision_score(y_true_binary, y_pred_binary, average='binary')
    recall = recall_score(y_true_binary, y_pred_binary, average='binary')
    f1 = 2 * precision * recall / (precision + recall)
    score = (1 - mae / 100) * 0.5 + (1 - range_mae / 100) * f1 * 0.5
    return score

该函数用于计算模型的各项评估指标，包括精确度、召回率、F1值和评分。

8. 模型评估函数

def evaluate_model(model, test_loader, device='cuda'):
    model.eval()
    predictions = []
    targets = []
    with torch.no_grad():
        for inputs, target in test_loader:
            inputs = [x.to(device) for x in inputs]
            outputs = model(inputs)
            predictions.extend(outputs.cpu().numpy())
            targets.extend(target.numpy())
    y_pred = np.array(predictions)
    y_test = np.array(targets)
    score = calculate_metrics(y_test, y_pred)
    print(f"Test Score: {score:.4f}")

该函数用于在测试集上评估模型性能。

9. 模型训练函数

def train_model(model, train_loader, val_loader, criterion, optimizer, num_epochs=50, device='cuda', output_dir=""):
    model.to(device)
    best_score = -float('inf')
    best_model = None
    for epoch in range(num_epochs):
        model.train()
        train_loss = 0
        for inputs, targets in tqdm(train_loader, desc=f'Epoch {epoch+1}/{num_epochs}'):
            inputs = [x.to(device) for x in inputs]
            targets = targets.to(device)
            optimizer.zero_grad()
            outputs = model(inputs)
            loss = criterion(outputs, targets)
            loss.backward()
            optimizer.step()
            train_loss += loss.item()
        model.eval()
        val_loss = 0
        val_preds = []
        val_targets = []
        with torch.no_grad():
            for inputs, targets in val_loader:
                inputs = [x.to(device) for x in inputs]
                targets = targets.to(device)
                outputs = model(inputs)
                loss = criterion(outputs, targets)
                val_loss += loss.item()
                val_preds.extend(outputs.cpu().numpy())
                val_targets.extend(targets.cpu().numpy())
        train_loss /= len(train_loader)
        val_loss /= len(val_loader)
        val_preds = np.array(val_preds)
        val_targets = np.array(val_targets)
        score = calculate_metrics(val_targets, val_preds)
        print(f'Epoch {epoch+1}/{num_epochs}')
        print(f'Train Loss: {train_loss:.4f}, Val Loss: {val_loss:.4f}')
        print(f'Learning Rate: {optimizer.param_groups[0]["lr"]:.6f}')
        print(f'Validation Score: {score:.4f}')
        if score > best_score:
            best_score = score
            best_model = model.state_dict().copy()
            torch.save(model.state_dict(), os.path.join(output_dir, "best.pt".format(epoch)))
            print(f'New best model found with score: {best_score:.4f}')
    return best_model

该函数用于训练模型，并在每个epoch后评估模型的性能，保存最佳模型。

10. 训练主程序

# 设置参数
bs = 64
epochs = 50
lr = 0.001
seed = 42
output_dir = "output/models"

# 选择设备
device = 'cuda' if torch.cuda.is_available() else 'cpu'

# 设置随机种子以确保结果可重复
set_random_seed(seed)

# 创建输出目录
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# 加载数据
train_data = pd.read_csv('train_data.csv')

# 指定需要处理的列
columns = ['siRNA_antisense_seq', 'modified_siRNA_antisense_seq_list']
# 删除包含空值的行
train_data.dropna(subset=columns + ['mRNA_remaining_pct'], inplace=True)
# 将数据分为训练集和验证集
train_data, val_data = train_test_split(train_data, test_size=0.1, random_state=42)

# 创建分词器
tokenizer = GenomicTokenizer(ngram=3, stride=3)

# 创建词汇表
all_tokens = []
for col in columns:
    for seq in train_data[col]:
        if ' ' in seq:
            all_tokens.extend(seq.split())
        else:
            all_tokens.extend(tokenizer.tokenize(seq))
vocab = GenomicVocab.create(all_tokens, max_vocab=10000, min_freq=1)

# 找到最大序列长度
max_len = max(max(len(seq.split()) if ' ' in seq else len(tokenizer.tokenize(seq)) 
                    for seq in train_data[col]) for col in columns)

# 创建数据集
train_dataset = SiRNADataset(train_data, columns, vocab, tokenizer, max_len)
val_dataset = SiRNADataset(val_data, columns, vocab, tokenizer, max_len)

# 创建数据加载器
train_loader = DataLoader(train_dataset, batch_size=bs, shuffle=True)
val_loader = DataLoader(val_dataset, batch_size=bs)

# 初始化模型
model = SiRNAModel(len(vocab.itos))
criterion = nn.MSELoss()

# 初始化优化器
optimizer = optim.Adam(model.parameters(), lr=lr)

# 训练模型
best_model = train_model(model, train_loader, val_loader, criterion, optimizer, epochs, device, output_dir=output_dir)

这个部分设置了训练的参数，加载并处理数据，创建词汇表和数据集，初始化模型并进行训练。

11. 测试程序

# 设置输出目录
output_dir = "result"
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# 加载测试数据
test_data = pd.read_csv('sample_submission.csv')
columns = ['siRNA_antisense_seq', 'modified_siRNA_antisense_seq_list']
test_data.dropna(subset=columns, inplace=True)

# 创建分词器
tokenizer = GenomicTokenizer(ngram=3, stride=3)

# 创建词汇表
all_tokens = []
for col in columns:
    for seq in test_data[col]:
        if ' ' in seq:
            all_tokens.extend(seq.split())
        else:
            all_tokens.extend(tokenizer.tokenize(seq))
            
vocab = GenomicVocab.create(all_tokens, max_vocab=10000, min_freq=1)

# 找到最大序列长度
max_len = max(max(len(seq.split()) if ' ' in seq else len(tokenizer.tokenize(seq)) 
                    for seq in test_data[col]) for col in columns)

# 创建测试数据集
test_dataset = SiRNADataset(test_data, columns, vocab, tokenizer, max_len, is_test=True)

# 创建数据加载器
test_loader = DataLoader(test_dataset, batch_size=64, shuffle=False)

# 初始化模型
model = SiRNAModel(len(vocab.itos))
model.load_state_dict(best_model)  # 加载最佳模型权重
model.to(device=device)
model.eval()

# 进行预测
preds = []
with torch.no_grad():
    for inputs in tqdm(test_loader):
        inputs = [x.to(device) for x in inputs]
        outputs = model(inputs)
        preds.extend(outputs.cpu().numpy())

# 将预测结果添加到测试数据中
test_data["mRNA_remaining_pct"] = preds
df = pd.DataFrame(test_data)

# 保存预测结果
output_csv = os.path.join(output_dir, "submission.csv")
print(f"submission.csv 保存在 {output_csv}")
df.to_csv(output_csv, index=False)

结语

        通过本次比赛，我们不仅可以深入了解化学修饰对siRNA序列的沉默效率的影响，还能够掌握基于深度学习的序列数据处理和模型构建技术。在RNA干扰（RNAi）机制下，siRNA的有效性直接关系到药物的实际疗效，因此预测其沉默效率具有重要的实际应用价值。

        本次任务提供了丰富的训练数据和清晰的评分机制，为参赛者提供了一个挑战和提升自己技能的机会。通过探索不同的深度学习模型和优化算法，我们可以不断提高模型的预测精度，从而推动siRNA药物设计的进步。

        希望大家在比赛过程中能够深入理解问题的本质，积极探索和尝试新的方法，不断优化自己的模型。愿大家在挑战中获得丰硕的成果，为核酸药物研发贡献自己的力量。祝大家比赛顺利，取得优异成绩！

如果你觉得这篇博文对你有帮助，请点赞、收藏、关注我，并且可以打赏支持我！

欢迎关注我的后续博文，我将分享更多关于人工智能、自然语言处理和计算机视觉的精彩内容。

谢谢大家的支持！