Datawhale AI 夏令营——siRNA药物药效预测Task1

最新推荐文章于 2024-07-28 23:52:46 发布

五倍子的代码空间

最新推荐文章于 2024-07-28 23:52:46 发布

阅读量610

点赞数 10

分类专栏：项目实战文章标签：人工智能 python ipython scikit-learn rnn 回归

本文链接：https://blog.csdn.net/qq_62738515/article/details/140751133

版权

项目实战专栏收录该内容

6 篇文章 0 订阅

订阅专栏

一、赛题背景与要求

近年来，mRNA疫苗在新冠预防领域取得巨大成功，推动了核酸类药物的研发。本次比赛旨在利用机器学习与深度学习相关技术，通过化学修饰后的siRNA序列预测RNA干扰（RNAi）机制下对靶基因的沉默效率。RNAi是一种通过抑制靶基因表达来调控基因表达的天然机制，在基因治疗和疾病治疗中具有广泛应用。本次比赛的数据集包括公开文献中的siRNA修饰序列及实验条件数据，初赛阶段提供约85%用于训练，15%用于验证。复赛将增加未公开专利数据，评估模型在新靶基因上的预测能力。本次比赛提供的baseline代码涵盖了数据处理、词汇表构建、序列编码和RNN模型训练等内容，帮助参赛者快速入门。通过本次比赛，参赛者将深入了解化学修饰对siRNA功能的影响，学习深度学习模型和序列处理技术。

二、跑通baseline

1.数据集概览

数据集包括siRNA裸序列、经过化学修饰的siRNA序列、目标mRNA序列以及实验条件（如药物浓度、细胞系、转染方式等）。最重要的字段是mRNA_remaining_pct，这是我们模型的训练目标，表示siRNA对靶基因沉默后的剩余mRNA百分比，值越低表示沉默效率越好。

2.运行baseline

官方提供的baseline代码需要在jupyter notebook中才能运行，也可以使用魔搭平台提供的免费算力我的Notebook · 魔搭社区 (modelscope.cn)https://modelscope.cn/my/mynotebook/preset

该代码主要由以下几部分组成：

①基因组分词

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

class GenomicTokenizer:
    def __init__(self, ngram=5, stride=2):
        # 初始化分词器，设置n-gram长度和步幅
        self.ngram = ngram
        self.stride = stride
        
    def tokenize(self, t):
        # 将输入序列转换为大写
        t = t.upper()
        
        if self.ngram == 1:
            # 如果n-gram长度为1，直接将序列转换为字符列表
            toks = list(t)
        else:
            # 否则，按照步幅对序列进行n-gram分词
            toks = [t[i:i+self.ngram] for i in range(0, len(t), self.stride) if len(t[i:i+self.ngram]) == self.ngram]
        
        # 如果最后一个分词长度小于n-gram，移除最后一个分词
        if len(toks[-1]) < self.ngram:
            toks = toks[:-1]
        
        # 返回分词结果
        return toks

②创建词汇表

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

class GenomicVocab:
    def __init__(self, itos):
        # 初始化词汇表，itos是一个词汇表列表
        self.itos = itos
        # 创建从词汇到索引的映射
        self.stoi = {v: k for k, v in enumerate(self.itos)}
        
    @classmethod
    def create(cls, tokens, max_vocab, min_freq):
        # 创建词汇表类方法
        # 统计每个token出现的频率
        freq = Counter(tokens)
        # 选择出现频率大于等于min_freq的token，并且最多保留max_vocab个token
        itos = ['<pad>'] + [o for o, c in freq.most_common(max_vocab - 1) if c >= min_freq]
        # 返回包含词汇表的类实例
        return cls(itos)

③加载数据集

该类用于加载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):
        # 获取数据集中的第idx个样本
        row = self.df.iloc[idx]  # 获取第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)  # 使用分词器分词
        
        # 将token转换为索引，未知token使用0（<pad>）
        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)

④构建模型

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

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)
        # 初始化GRU层
        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)  # hidden_dim * 4 因为GRU是双向的，有n_layers层
        # 初始化Dropout层
        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)  # 传入GRU层
            x = self.dropout(x[:, -1, :])  # 取最后一个隐藏状态，并进行dropout处理
            outputs.append(x)
        
        # 将所有序列的输出拼接起来
        x = torch.cat(outputs, dim=1)
        # 传入全连接层
        x = self.fc(x)
        # 返回结果
        return x.squeeze()

⑤评估

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

def calculate_metrics(y_true, y_pred, threshold=30):
    # 计算平均绝对误差
    mae = np.mean(np.abs(y_true - y_pred))

    # 将实际值和预测值转换为二进制分类（低于阈值为1，高于或等于阈值为0）
    y_true_binary = (y_true < threshold).astype(int)
    y_pred_binary = (y_pred < threshold).astype(int)

    # 创建掩码，用于筛选预测值在0和阈值之间的样本
    mask = (y_pred >= 0) & (y_pred <= threshold)
    range_mae = mean_absolute_error(y_true[mask], y_pred[mask]) if mask.sum() > 0 else 100

    # 计算精确度、召回率和F1得分
    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

⑥测试

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

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)
            # 将预测结果从GPU移到CPU，并转换为numpy数组，添加到predictions列表中
            predictions.extend(outputs.cpu().numpy())
            # 将目标值转换为numpy数组，添加到targets列表中
            targets.extend(target.numpy())

    # 将预测结果和目标值转换为numpy数组
    y_pred = np.array(predictions)
    y_test = np.array(targets)
    
    # 计算评估指标
    score = calculate_metrics(y_test, y_pred)
    # 打印测试得分
    print(f"Test Score: {score:.4f}")

⑦模型训练

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

def train_model(model, train_loader, val_loader, criterion, optimizer, num_epochs=50, device='cuda', output_dir: str=""):
    # 将模型移动到指定设备
    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()  # 更新参数
            #scheduler.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  # 返回最佳模型

⑧训练主程序

# 设置参数
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)

# 设置线性衰减学习率  
# lambda_lr = lambda epoch: lr - (epoch / epochs)  # 从 1 逐渐减小到 0  
# scheduler = LambdaLR(optimizer, lr_lambda=lambda_lr) 

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

⑨测试主程序

# 设置输出目录
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()  # 切换到评估模式，这对于某些模块如Dropout和BatchNorm是必需的

# 进行预测
preds = []
with torch.no_grad():
    for inputs in tqdm(test_loader):
        # import pdb;pdb.set_trace()
        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)

三、改进

1.调参

在官方提供的baseline代码中可以尝试增大epoch，改进学习率策略，以及调整batch size来提高精度。同时由于官方提供的数据集规模不是很大，可以尝试使用外部数据，扩大训练样本。需要注意在改进epoch时，要观察模型是否已经收敛，还可以使用早停策略。