#Datawhale AI 夏令营第三期—siRNA药物药效预测学习

JeasonLi ——Wenrui

已于 2024-08-03 17:26:18 修改

阅读量1.3k

点赞数 30

文章标签：人工智能学习 python 深度学习

于 2024-07-26 18:22:25 首次发布

本文链接：https://blog.csdn.net/2301_80270213/article/details/140697806

版权

赛事概要

一、赛题背景

随着mRNA疫苗在新冠预防领域取得成功，核酸类药物的研发获得了越来越多的关注。本次比赛聚焦于通过机器学习技术，利用化学修饰后的siRNA序列来预测RNA干扰（RNAi）机制下对靶基因（target gene）的沉默效率，这一指标与药物实际疗效直接相关。RNAi是生物体内天然存在的一种基因表达调控机制，通过抑制靶基因的表达来实现降低目标蛋白量的目的，这一机制一般可通过siRNA实现。目前开源的数据库中，以RNA主干序列（裸序列）为主，缺少相应的化学修饰数据。而本赛题特别关注了化学修饰对siRNA序列功能的影响，化学修饰对siRNA的毒性、体内稳定性、靶向效果、药效等具有重大影响，在实际药物设计中至关重要。参赛者将接触到领域独特的包含靶基因、siRNA裸序列、经过化学修饰的siRNA修饰序列以及实验室测定的沉默效率值。这些数据反映了当今siRNA设计的最新科技成果，包括化学修饰的种类和位置，以及它们如何影响siRNA对靶基因的沉默效率。

本赛题的研究成果将可以直接用于优化siRNA药物在医学研究和临床应用中的效率和安全性。优化siRNA的设计可以提高其沉默特定基因的能力，降低非特异性作用，从而增强疗效和降低副作用。这一技术的提升在疾病治疗和基因治疗领域具有极高的应用潜力和商业价值。此外通过建立对应预测模型，可以加速新型siRNA药物的虚拟筛选，促进个性化医疗的实现。

二、赛题任务

初赛

提供一部分公开文献中提取的siRNA修饰序列以及相应实验条件数据（例如药物浓度、细胞系、转染方式等），随机打散后使用约85%数据训练，剩余约15%数据用于leaderboard submission评分，主要用于测试针对在训练集中出现过的目标mRNA序列，不同siRNA的沉默效率预测的准确率。

复赛

在初赛基础上，将额外增加一部分公开文献提及的以及RONA Therapeutics（大睿生物）提供的尚未公开专利中的siRNA修饰序列数据作为测试数据，主要评估未在训练集中出现过的目标mRNA序列的siRNA的沉默效率的预测准确性。

三、评审规则

1.数据说明

train_data.csv

train_data.csv的每行为一条训练记录，包含数据记录的id、siRNA裸序列、相应的siRNA修饰序列、目标mRNA序列、siRNA浓度、细胞系、转染方法等实验室条件以及对应的实验室测量的mRNA Remaining值等总计19个字段。其中mRNA Remaining值为我们模型的训练目标，其余18个字段的全部或部分可以作为模型的输入特征。Remaining值代表了经过siRNA的沉默之后，mRNA的剩余百分比（相对于对照组）。Remaining值越低，siRNA的沉默效率越好，药效就越好。Remaining值一般位于0-100的区间内，100表示完全没有沉默效果，0表示该mRNA被彻底沉默，但是由于实验室测量的误差，可能存在少量训练记录的mRNA Remaining值在这个范围之外，这是正常的数据。

打开文件如下图所示：

表头解释：

sample_submission.csv

sample_submission.csv为初赛的leaderboard submission测试集，格式与初赛训练集train_data.csv相同。不同之处在于mRNA_remaining_pct列的数值为空，参赛者需要填充这些空白处的预测结果后提交。

vocab.csv

vocab.csv为经过归一化后的siRNA化学修饰缩写表

第一列（Abbreviation）：列出了各种修饰核苷酸的缩写。
第二列（Chemical Name）：展示了与第一列中缩写对应的完整化学名称。这些化学名称详细描述了核苷酸的结构，包括对其进行的化学修饰。这些修饰可能涉及改变核苷酸的糖部分、磷酸基团或碱基，以增强其性能，如提高稳定性或特定的生物活性。

baseline.py

本次比赛提供一个基础的基线方法，利用 RNN 来预测Remaining值。基线方法中我们仅使用了siRNA_sense_seq 字段作为特征，除此之外尚有其他特征对Remaining结果有重大影响。

2.评分机制

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

平均绝对误差（MAE）：平均绝对误差是对每个样本的绝对误差取平均值，主要用于衡量模型在整体上的预测精度。 $y \displaystyle _{_{i}}$ 为实验室的真实测量值， $\hat{y}_{i}$ 为模型的预测值，计算公式 MAE= $\frac{1}{n}\sum _{i=1}^{n}\left | y_{i} -\hat{y_{i}}\right |$ ，其中 $n$ 是样本总数，MAE∈[0,100]MAE∈[0,100]。
预测值在一定范围内的平均绝对误差（Range-MAE）：在实际场景我们会尤其关注在低Remaining区间模型的预测准确率，此指标衡量在特定范围内的预测值的平均绝对误差。在本次比赛中，低Remaining范围为[0,30]，计算公式 Range-MAE= $\frac{1}{m}\sum _{i=1}^{m}\left | y_{i} -\hat{y_{i}}\right |$ ，其中 $m$ 是预测值在[0,30][0,30]范围内的样本数量，Range-MAE∈[0,100]。
预测值在一定范围内的F1指标（F1）：此指标衡量在特定范围内的预测值的precision以及在相应区间预测的recall情况，综合得到F1值。在本次比赛中，低Remaining范围为[0,30][0,30]。

初赛

初赛评测：将公开文献中的siRNA序列数据随机打乱后，提供选手85%数据训练评估，保留15%数据用于leaderboard评分。最终得分基于上述3个指标计算得到，计算公式

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

score计算的第一项关注模型在整体区间上的预测精度情况，score计算的第二项同时关注模型在低Remaining区间的粗分类结果以及更进一步在低Remaining区间的预测精度情况。

复赛

复赛评测：额外增加一部分公开文献以及大睿生物提供的尚未公开的siRNA修饰序列实验数据作为测试数据。最终得分计算方式与初赛一致。

代码实现解释：

# score = 50% × (1−MAE/100) + 50% × F1 × (1−Range-MAE/100)
def calculate_metrics(y_true, y_pred, threshold=30):
    # 预测值和真实值之间绝对误差的平均值，衡量了模型预测的总体准确性。
    mae = np.mean(np.abs(y_true - y_pred))

    # 这里将真实值和预测值进行二值化处理：如果值小于阈值（30），则为1，否则为0。
    y_true_binary = (y_true < threshold).astype(int)
    y_pred_binary = (y_pred < threshold).astype(int)

    # 这里计算在特定区间（0到阈值30）内的平均绝对误差。如果预测值在这个区间内，才计算其误差，否则设为100。这个指标评估了模型在重要预测区间内的表现。
    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

Task 1：小白入门跑通Baseline（RNN模型）

本次比赛提供的baseline代码涵盖了数据处理、词汇表构建、序列编码和RNN模型训练等内容，帮助参赛者快速入门。通过本次比赛，参赛者将深入了解化学修饰对siRNA功能的影响，学习深度学习模型和序列处理技术。

明确本次比赛的目标是通过机器学习模型预测化学修饰siRNA对靶基因的沉默效率（即预测某类药物基因对某类疾病基因的治疗效果），这一目标决定了数据处理、特征选择和模型设计的方向。

必知概念入门

RNA干扰（RNAi）

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

化学修饰siRNA

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

siRNA修饰的目的

1.提高稳定性：防止siRNA在生物体内被核酸酶降解，从而延长其作用时间。

2.增强靶向性：确保siRNA能够特异性地结合并抑制靶基因的表达，减少脱靶效应。

3.降低毒性：减少siRNA对细胞的毒性作用，提高其安全性。

4.提高生物利用度：通过优化递送系统，提高siRNA 在体内的吸收和分布效率。

深度学习与RNN

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

补充：RNN详解(Recurrent Neural Network)-CSDN博客

循环神经网络（Recurrent Neural Network, RNN）是一类具有内部环状连接的人工神经网络，用于处理序列数据。其最大特点是网络中存在着环，使得信息能在网络中进行循环，实现对序列信息的存储和处理。

RNN的基本结构：
# 一个简单的RNN结构示例
class SimpleRNN(nn.Module):
    def __init__(self, input_size, hidden_size):
        super(SimpleRNN, self).__init__()
        self.rnn = nn.RNN(input_size, hidden_size, batch_first=True)
        
    def forward(self, x):
        out, _ = self.rnn(x)
        return out
类定义:

class SimpleRNN(nn.Module):

这里定义了一个名为 SimpleRNN 的类，它继承自 nn.Module。nn.Module 是PyTorch中所有神经网络模块的基类，提供了一些基础功能。

构造函数:

def __init__(self, input_size, hidden_size): super(SimpleRNN, self).__init__() self.rnn = nn.RNN(input_size, hidden_size, batch_first=True)

input_size: 输入特征的维度。
hidden_size: RNN隐藏层的维度。
super(SimpleRNN, self).__init__(): 调用基类的构造函数，这是Python中继承机制的一部分。
self.rnn: 创建一个 nn.RNN 实例。nn.RNN 是PyTorch中实现基本RNN的类。
第一个参数 input_size 指定了输入特征的维度。
第二个参数 hidden_size 指定了隐藏层的维度。
batch_first=True 指示输入和输出张量的第一个维度是批次大小（batch size），这使得模型的输入输出格式更加灵活。

前向传播函数:

def forward(self, x): out, _ = self.rnn(x) return out

x: 输入数据，通常是一个三维张量，其形状为 (batch_size, seq_length, input_size)，其中：
batch_size 是批次中的样本数。
seq_length 是序列的长度。
input_size 是每个时间步的输入特征数。
self.rnn(x): 调用RNN层进行前向传播。RNN层的输出是一个元组，包含：
out: RNN的输出，形状为 (batch_size, seq_length, hidden_size)。
_: 这里使用下划线（_）来忽略梯度，因为通常我们只关心输出而不关心隐藏状态的梯度。
return out: 返回RNN的输出。

工作原理

输入层：RNN能够接受一个输入序列（例如文字、股票价格、语音信号等）并将其传递到隐藏层。
隐藏层：隐藏层之间存在循环连接，使得网络能够维护一个“记忆”状态，这一状态包含了过去的信息。这使得RNN能够理解序列中的上下文信息。
输出层：RNN可以有一个或多个输出，例如在序列生成任务中，每个时间步都会有一个输出。

RNN的优缺点

优点：

能够处理不同长度的序列数据。
能够捕捉序列中的时间依赖关系。

缺点：

对长序列的记忆能力较弱，可能出现梯度消失或梯度爆炸问题。
训练可能相对复杂和时间消耗大。

RNN的应用场景

循环神经网络（RNN）因其在捕获序列数据中的时序依赖性方面的优势，在许多应用场景中都得到了广泛的使用。以下是一些主要应用领域的概述：

词汇表与序列编码

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

数据处理与特征选择

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

模型训练与评估

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

PyTorch框架

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

Baseline 分块代码如下：

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  # 模型评估指标

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

导入库：首先，代码导入了需要用到的库，包括 pandas—用于数据处理和分析，提供高效的数据结构（如DataFrame）和数据分析工具，numpy—一个科学计算库，特别是针对数组和矩阵的操作，以及torch—PyTorch的核心库，提供了张量（Tensor）的操作和自动微分系统。

pandas 和 numPy 通常一起使用，因为 pandas 的设计初衷就是为了与 numPy 无缝集成，利用 numPy 的数组结构来提高性能。

torch: 这是PyTorch的核心库，提供了张量（Tensor）的操作和自动微分系统。张量是PyTorch中的基本数据结构，可以看作是多维数组，类似于NumPy中的数组，但是PyTorch的张量可以在GPU上进行操作，从而加速计算。自动微分系统允许用户自动计算导数，这对于训练神经网络至关重要。

torch.nn: 这是PyTorch的神经网络模块，包含了构建神经网络所需的各种层（如卷积层、循环层、全连接层等）和损失函数。它还提供了一些常用的功能，比如初始化参数的方法、批量归一化、池化层等。

torch.optim: 这个模块提供了多种优化算法，用于神经网络的参数更新。这些优化算法包括但不限于：

SGD（随机梯度下降）
Adam（自适应矩估计）
RMSprop（均方根传播）
ASGD（平均随机梯度下降）
Adadelta（自适应学习率）
等等

这些优化器帮助我们在训练过程中调整网络的权重，以最小化损失函数，从而达到学习数据模式的目的。

# 该函数确保了在使用NumPy、Python内置随机数生成器和PyTorch时，所有的随机数生成都是可控的和可复现的，有助于实验结果的一致性。
def set_random_seed(seed):
    # 设置NumPy的随机种子
    np.random.seed(seed)
    # 设置Python内置的随机数生成器的种子
    random.seed(seed)
    # 设置PyTorch的随机种子
    torch.manual_seed(seed)
    # 设置CUDA的随机种子
    torch.cuda.manual_seed(seed)
    # 设置所有CUDA设备的随机种子
    torch.cuda.manual_seed_all(seed)
    # 确保每次卷积算法选择都是确定的
    torch.backends.cudnn.deterministic = True
    # 关闭CuDNN自动优化功能，确保结果可复现
    torch.backends.cudnn.benchmark = False

这段代码定义了一个名为 set_random_seed 的函数，其主要功能是设置随机种子，以确保在执行涉及随机数生成的操作时，结果的可重复性。这在进行科学实验或机器学习训练时非常重要，因为它允许实验者在不同的运行中得到相同的结果，从而验证实验的可靠性。

下面是代码中每行的详细作用解释：

np.random.seed(seed): 设置NumPy库的随机种子。NumPy是一个广泛使用的科学计算库，它提供了大量的数学函数和对多维数组的支持。通过设置随机种子，NumPy生成的随机数序列将具有可重复性。

random.seed(seed): 设置Python内置的 random 模块的随机种子。这个模块提供了生成随机数的方法，如 random.random(), random.randint() 等。设置种子后，这些方法生成的随机数也将是可重复的。

torch.manual_seed(seed): 设置PyTorch的CPU随机数生成器的种子。PyTorch是一个深度学习框架，它提供了许多用于生成随机数的操作，如初始化神经网络的权重。设置这个种子可以确保这些操作的结果是可重复的。

torch.cuda.manual_seed(seed): 如果使用CUDA（Compute Unified Device Architecture，一种由NVIDIA开发的并行计算平台和编程模型），这行代码将设置当前GPU的随机数生成器的种子。

torch.cuda.manual_seed_all(seed): 与上一行类似，但这行代码将为所有CUDA设备设置随机种子，确保在多GPU环境中结果的一致性。

torch.backends.cudnn.deterministic = True: 设置PyTorch的cudnn后端为确定性模式。cudnn是NVIDIA提供的用于深度神经网络的GPU加速库。在确定性模式下，某些操作（如卷积）将总是选择相同的算法，以确保结果的一致性。

torch.backends.cudnn.benchmark = False: 关闭cudnn的基准测试功能。cudnn的基准测试功能会尝试找到最优的算法来执行操作，但这可能会导致每次运行时选择的算法不同，从而影响结果的可重复性。关闭这个功能可以确保每次执行操作时都使用相同的算法。

通过调用这个函数并传入一个固定的种子值，可以确保在不同的实验运行中，所有涉及随机数生成的部分都能产生相同的结果，从而提高实验的可重复性。

2. 基因组分词器类

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

该类用于将基因组序列分割成固定长度的n-gram。
输入：长段基因输出：按照规则分割后的诸多短的片段基因

类初始化：`__init__` 方法接受两个参数 `ngram` 和 `stride`，用于设置分词器的 n-gram 长度和步幅。
分词方法：`tokenize` 方法将输入的序列转换为大写，并根据 `ngram` 和 `stride` 对序列进行分词。（ ngram: n-gram的长度，默认为5。n-gram是一种文本表示方法，其中文本被分割成连续的n个字符的序列。stride: 步幅，默认为2。步幅定义了在原始序列上进行n-gram分词时，每次移动的字符数）
`n-gram` 长度为 1 的处理：如果 `ngram` 为 1，直接将序列转换为字符列表。
`n-gram` 长度大于 1 的处理：按步幅进行分词，并确保每个分词的长度等于 `ngram`。
最后一个分词的处理：如果最后一个分词长度小于 `ngram`，将其移除。
返回分词结果：返回处理后的分词结果列表。

3. 基因组词汇类

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)

该类用于创建一个词汇表，用于将基因组片段映射为索引。将基因组中的特定片段与一种数据结构（索引）相关联，以便在后续的数据分析中高效地存储、检索和分析这些片段。

1. 类初始化：`__init__` 方法接受一个参数 `itos`，它是一个词汇表列表。

`itos`：从索引到词汇的映射。
`stoi`：从词汇到索引的映射，由 `itos` 列表生成。

2. 类方法 `create`：创建词汇表的类方法，用于生成 `GenomicVocab` 类的实例。

参数：

`tokens`：所有token的列表。
`max_vocab`：词汇表的最大容量。
`min_freq`：词汇在被包含到词汇表中的最低频率。

步骤：

统计 `tokens` 中每个token出现的频率。
按照频率从高到低排序，并选择出现频率大于等于 `min_freq` 的token，最多保留 `max_vocab` 个。
在词汇表中添加一个特殊的 `<pad>` token，用于填充序列。
返回包含生成的 `itos` 列表的 `GenomicVocab` 实例。

原因：在疾病中高频率出现的基因片段，有一定概率 “是反派”，要“揪出来” 在药物中高频率出现的基因片段，有一定概率“是好人”，要“拉拢过来”

4.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)

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

序列数据：按照一定的顺序排列的数据集合，一个重要特性是每个数据项通常都包含了与其前后相邻数据项相关的信息。

1. 初始化：`__init__` 方法初始化数据集的必要参数：

`df`：包含数据的Pandas数据框。
`columns`：包含序列的列名列表。
`vocab`：词汇表对象，用于将token转换为索引。
`tokenizer`：分词器对象，用于将序列分割为token。
`max_len`：序列的最大长度，所有序列将被填充或截断到这个长度。
`is_train`：布尔值，指示数据集是否用于训练（默认为False）。

2. 获取数据集长度：`__len__` 方法返回数据集的样本数量。

3. 获取样本：__getitem__ 方法获取数据集中指定索引的样本：

获取指定行的数据。
对每个包含序列的列进行分词和编码。
获取目标值（mRNA_remaining_pct），并将其转换为张量。
返回编码后的序列和目标值。

4. 分词和编码：`tokenize_and_encode` 方法对输入序列进行分词和编码：

如果序列中包含空格，按空格分词（表示序列已经被修改）。
否则，使用分词器对序列进行分词。
将分词结果转换为索引，未知`token`用`0（<pad>）`表示。
对序列进行填充，使其长度等于最大长度 `max_len`。
返回张量格式的填充序列。

文本数据处理的常见步骤：

①文本清洗：去除文本中的特殊字符、标点符号、数字等，并进行大小写统一、去停用词等处理。

②分词：将文本分割成单词或词语。

③词汇编码：构建词汇表，为每个词汇分配整数ID，并进行词嵌入。

④序列填充或截断：将文本序列填充或截断至统一长度。

⑤转换为张量：将处理后的文本序列转换为张量形式。张量：即多维数组

⑥数据加载与批处理：将张量数据加载到模型中，并进行批处理。

5.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)
        # 初始化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()

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

1. 初始化方法 `__init__`：

`vocab_size`：词汇表大小，用于嵌入层。
`embed_dim`：嵌入维度，嵌入层将词汇映射为 `embed_dim` 维向量。
`hidden_dim`：隐藏层维度，GRU的隐藏状态维度。
`n_layers`：GRU的层数。
`dropout`：`Dropout`层的丢弃率，用于防止过拟合。

该方法初始化了模型的各层，包括嵌入层、GRU层、全连接层和Dropout层。

2. 前向传播方法 `forward`：

将输入序列传入嵌入层进行词汇嵌入。
对每个嵌入的序列进行`GRU`处理，提取最后一个隐藏状态并进行`Dropout`处理。
将所有处理后的序列输出拼接起来，并传入全连接层。
返回经过全连接层后的结果。

结合笔记开始对RNN的初步理解：

6.评估指标计算函数

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

1. 计算平均绝对误差 (MAE)：

mae = np.mean(np.abs(y_true - y_pred))：计算实际值和预测值之间的`平均绝对误差`。

2. 将实际值和预测值转换为二进制分类：

y_true_binary = (y_true < threshold).astype(int)：如果实际值小于阈值，设为`1`，否则设为`0`。
y_pred_binary = (y_pred < threshold).astype(int)：如果预测值小于阈值`，设为`1`，否则设为`0`。

3. 创建掩码：

mask = (y_pred >= 0) & (y_pred <= threshold)：筛选预测值在`0`和`阈值`之间的样本。
range_mae = mean_absolute_error(y_true[mask], y_pred[mask]) if mask.sum() > 0 else 100：计算这些样本的`平均绝对误差`，如果`没有`符合条件的样本，设为`100`。

4. 计算精确度、召回率和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)：计算F1得分。

5. 计算综合评分：

score = (1 - mae / 100) * 0.5 + (1 - range_mae / 100) * f1 * 0.5：综合MAE、范围内的MAE和F1得分，计算最终评分。

6. 返回评分：

return score：返回综合评分。

7.模型评估函数

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}")

该函数用于在测试集上评估模型性能。使用上一步骤中定义的函数，作为模型评价指标。

8.模型训练函数

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()  # 更新参数
            
            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后评估模型的性能，保存最佳模型。

epoch：神经网络已经遍历（或者说“学习”）了训练集中的每一个样本一次。

模型训练通常在一个或多个epoch中进行，每个epoch包含对整个训练集的完整遍历。在每个epoch中，模型会对每个批次的数据进行以下操作：

前向传播：将数据输入模型，得到预测值。

计算损失：使用损失函数计算预测值与真实值之间的差异。

反向传播：根据损失函数的梯度反向传播到模型的每一层，计算每个参数的梯度。

参数更新：使用优化器根据计算出的梯度更新模型的参数。

9.训练主程序

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

优化模型的参数，使模型能够学习数据的特征并做出准确的预测或分类。

10.测试程序并保存结果到本地

# 设置输出目录
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)

用于评估训练好的模型或软件在实际使用中的性能和表现。

思考体验

Task 1 所给出的Baseline代码略有难度，使用深度学习模型，但是只使用了部分特征['siRNA_antisense_seq', 'modified_siRNA_antisense_seq_list']进行训练，所以提交后成绩并不高，后面还会学习更多技巧来提分。首先刚上手又接触到了许多新的知识，也认识到了AI+药物的优势。

Task2：深入理解赛题，入门RNN和特征工程

本任务我们对官方的baseline进行分析解读，之后介绍RNN相关的基础知识，包括其适用范围和问题。随后，我们从特征工程构建的角度来重新分析赛题数据，并将其和lightgbm结合，最终给出一个更好的baseline。

官方baseline分析（Task 1 代码分析）

在baseline中，我们只用到了siRNA_antisense_seq和modified_siRNA_antisense_seq_list，它们都是由一串符号标记的序列，我们希望的是把这些序列特征能够输入RNN模型，因此需要对其做一定处理。在SiRNAModel类的forward方法中，展示了在得到序列特征的tensor表示后的处理步骤：

ef 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()

那么这里的输入x是什么呢？我们可以通过train_loader来查看一个batch内的输入情况，这里的inputs和上面的x是一个东西。我们首先发现inputs包含两个元素，它们分别对应的是前面提到的两个使用的特征，每个元素的尺寸都是64*25，64代表batch的大小，25代表序列的长度。这里我们可以从inputs[0][0]看到每一行数据的siRNA_antisense_seq被向量化后的情况，这个例子中我们发现前面的7位是非零数，表示其序列编码后每一位的唯一标识；而后面都是0，这是因为RNN模型的输入需要每个样本的长度一致，因此我们需要事先算出一个所有序列编码后的最大长度，然后补0。

那么我们怎么能得到这个唯一标识呢？我们首先需要把序列给进行分词，siRNA_antisense_seq的分词策略是3个一组（GenomicTokenizer的ngram和stride都取3）进行token拆分，比如AGCCGAGAU会被分为[AGC, CGA, GAU]，而modified_siRNA_antisense_seq_list会进行按照空格分词（因为它本身已经根据空格分好了）。由此我们可以从整个数据集构建出一个词汇表，他负责token到唯一标识（索引）的映射：

有了这个词汇表，我们就可以

来获得序列的最大长度

在loader获取样本的时候把token转为索引

此时，对于某一行数据，其两个特征分别为AGCCUUAGCACA和u u g g u u Cf c，假设整个数据集对应token编码后序列的最大长度为10，那么得到的特征就可能是

[25, 38, 25, 24, 0, 0, 0, 0, 0, 0]
[65, 65, 63, 63, 65, 65, 74, 50, 0, 0]

那么假设batch的大小为16，此时forword函数的x就会是两个列表，每个列表的tensor尺寸为16 * 10

RNN模型分析

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()

我们首先第一步将得到的索引进行了embedding，token的embedding是将离散的符号（如单词、字符、或基因序列片段）映射到连续的向量空间的过程。这个过程通过将高维的稀疏表示（如独热编码）转换为低维的密集向量表示，使得相似的符号在向量空间中距离更近。此时，embed的尺寸会从BatchSize * Length成为BatchSize * Length * EmbeddingSize，此处EmbeddingSize即embed_dim=200。

token：模型输入基本单元。比如中文BERT中，token可以是一个字，也可以是等标识符。
embedding：一个用来表示token的稠密的向量。token本身不可计算，需要将其映射到一个连续向量空间，才可以进行后续运算，这个映射的结果就是该token对应的embedding。

Token Embeddings 是一种将文本中的词语转化为向量表示的方法。在自然语言处理中，我们通常将文本表示为一个向量矩阵，其中每个词语对应一个向量。这些向量被称为词向量或者词嵌入。Token Embeddings 是一种词向量的扩展，它可以将不同类型的词语（如单词、字符、子词）都转化为向量表示。

一起学习大模型 - 从底层了解Token Embeddings的原理（1）-CSDN博客

RNN，全称为递归神经网络（Recurrent Neural Network），是一种人工智能模型，特别擅长处理序列数据。它和普通的神经网络不同，因为它能够记住以前的数据，并利用这些记忆来处理当前的数据。想象你在读一本书。你在阅读每一页时，不仅仅是单独理解这一页的内容，还会记住前面的情节和信息。这些记忆帮助你理解当前的情节并预测接下来的发展。这就是 RNN 的工作方式。假设你要预测一个句子中下一个单词是什么。例如，句子是：“我今天早上吃了一个”。RNN 会根据之前看到的单词（“我今天早上吃了一个”），预测下一个可能是“苹果”或“香蕉”等。它记住了之前的单词，并利用这些信息来做出预测。

RNN 在处理序列数据时具有一定的局限性：
- 长期依赖问题：RNN 难以记住和利用很久以前的信息。这是因为在长序列中，随着时间步的增加，早期的信息会逐渐被后来的信息覆盖或淡化。
- 梯度消失和爆炸问题：在反向传播过程中，RNN 的梯度可能会变得非常小（梯度消失）或非常大（梯度爆炸），这会导致训练过程变得困难。
LSTM 的改进
- LSTM 通过引入一个复杂的单元结构来解决 RNN 的局限性。LSTM 单元包含三个门（输入门、遗忘门和输出门）和一个记忆单元（细胞状态），这些门和状态共同作用，使 LSTM 能够更好地捕捉长期依赖关系。
  - 输入门：决定当前输入的信息有多少会被写入记忆单元。
  - 遗忘门：决定记忆单元中有多少信息会被遗忘。
  - 输出门：决定记忆单元的哪些部分会作为输出。
- 通过这些门的控制，LSTM 可以选择性地保留或遗忘信息，从而有效地解决长期依赖和梯度消失的问题。
GRU 的改进
- GRU 是 LSTM 的一种简化版本，它通过合并一些门来简化结构，同时仍然保留了解决 RNN 局限性的能力。GRU 仅有两个门：更新门和重置门。
  - 更新门：决定前一个时刻的状态和当前输入信息的结合程度。
  - 重置门：决定忘记多少之前的信息。
- GRU 的结构更简单，计算效率更高，同时在许多应用中表现出与 LSTM 类似的性能。

我们在pytorch的GRU文档中可以找到对应可选的参数信息，我们需要特别关注的参数如下，它们决定了模型的输入输出的张量维度

input_size（200）
hidden_size（256）
bidirectional（True）

假设输入的BatchSize为16，序列最大长度为10，即x尺寸为16 * 10 * 200，那么其输出的张量尺寸为 16 * 10 * (256 * 2)。

在从GRU模型输出后，x = self.dropout(x[:, -1, :])使得输出变为了BatchSize * (hidden_dim * 2)，此处取了序列最后一个位置的输出数据（注意RNN网络的记忆性），这里的2是因为bidirectional参数为True，随后x = torch.cat(outputs, dim=1)指定在第二个维度拼接后，通过全连接层再映射为标量，因此最后经过squeeze（去除维数为1的维度）后得到的张量尺寸为批大小，从而可以后续和target值进行loss计算，迭代模型。

进阶

数据的特征工程

在提交官方baseline后，我们会发现得分并不好，这一方面原因可能在于数据用的特征还较为简单，序列特征的构造较为粗糙，再加上对于深度学习的RNN模型而言数据量可能还不太充足的因素，这是可以预见的output。下面我们介绍一种把序列特征的问题转化为表格问题的方法，并介绍在表格数据上如何做特征工程。

处理类别型变量

如何知道一个变量是类别型的呢，只需看下其值的分布，或者唯一值的个数

df.gene_target_symbol_name.nunique()

df.gene_target_symbol_name.value_counts()

.nunique() 是Pandas库中DataFrame对象的一个方法，主要用于计数非空唯一值。当你在一个Series或DataFrame的列上调用这个函数时，它会返回该列中不同元素（即去除了NaN或缺失值）的个数。例如，在基因表达数据框（DataFrame）中，如果你有一个名为“target_gene_symbols”的列，.nunique() 可以告诉你有多少种不同的基因靶点符号。

这个函数常用于数据清洗、预处理和基本的数据探索阶段，因为它可以帮助我们快速了解某个特征的多样性和重复度。如果结果很大，可能意味着数据分布广泛；如果结果很小，那可能意味着存在大量的重复项。

.value_counts() 是Pandas DataFrame或Series对象上的一个统计函数，用于计算每个元素（在这个上下文中通常是观测值或类别）出现的次数。当应用于df.gene_target_symbol_name这样的列时，它会对列中的不同基因目标符号名称进行计数，并返回一个新的Series，其中索引是独特的基因符号名，值则是对应的频率（出现次数）。

如果相较于数据的总行数很少，那么其很可能就是类别变量了，比

如 gene_target_symbol_name。此时，我们可以使用get_dummie函数来实现one-hot特征的构造

# 如果有40个类别，那么会产生40列，如果第i行属于第j个类别，那么第j列第i行就是1，否则为0
df_gene_target_symbol_name = pd.get_dummies(df.gene_target_symbol_name)
df_gene_target_symbol_name.columns = [
    f"feat_gene_target_symbol_name_{c}" for c in df_gene_target_symbol_name.columns
]

get_dummie：pandas.get_dummies （独热编码）详解_pandas独热编码-CSDN博客

one-hot：机器学习：数据预处理之独热编码（One-Hot）详解-CSDN博客

可能的时间特征构造

在数据观察的时候发现，siRNA_duplex_id的编码方式很有意思，其格式为AD-1810676.1，我们猜测AD是某个类别，后面的.1是版本，当中的可能是按照一定顺序的序列号，因此可以构造如下特征：

siRNA_duplex_id_values = df.siRNA_duplex_id.str[3:-2].str.strip(".").astype("int")

包含某些单词

df_cell_line_donor = pd.get_dummies(df.cell_line_donor)
df_cell_line_donor.columns = [
    f"feat_cell_line_donor_{c}" for c in df_cell_line_donor.columns
]
# 包含Hepatocytes
df_cell_line_donor["feat_cell_line_donor_hepatocytes"] = (
    (df.cell_line_donor.str.contains("Hepatocytes")).fillna(False).astype("int")
)
# 包含Cells
df_cell_line_donor["feat_cell_line_donor_cells"] = (
    df.cell_line_donor.str.contains("Cells").fillna(False).astype("int")
)

根据序列模式提取特征

假设siRNA的序列为ACGCA...，此时我们可以根据上一个task中提到的rna背景知识，对碱基的模式进行特征构造：

def siRNA_feat_builder(s: pd.Series, anti: bool = False):
    name = "anti" if anti else "sense"
    df = s.to_frame()
    # 序列长度
    df[f"feat_siRNA_{name}_seq_len"] = s.str.len()
    for pos in [0, -1]:
        for c in list("AUGC"):
            # 第一个和最后一个是否是A/U/G/C
            df[f"feat_siRNA_{name}_seq_{c}_{'front' if pos == 0 else 'back'}"] = (
                s.str[pos] == c
            )
    # 是否已某一对碱基开头和某一对碱基结尾
    df[f"feat_siRNA_{name}_seq_pattern_1"] = s.str.startswith("AA") & s.str.endswith(
        "UU"
    )
    df[f"feat_siRNA_{name}_seq_pattern_2"] = s.str.startswith("GA") & s.str.endswith(
        "UU"
    )
    df[f"feat_siRNA_{name}_seq_pattern_3"] = s.str.startswith("CA") & s.str.endswith(
        "UU"
    )
    df[f"feat_siRNA_{name}_seq_pattern_4"] = s.str.startswith("UA") & s.str.endswith(
        "UU"
    )
    df[f"feat_siRNA_{name}_seq_pattern_5"] = s.str.startswith("UU") & s.str.endswith(
        "AA"
    )
    df[f"feat_siRNA_{name}_seq_pattern_6"] = s.str.startswith("UU") & s.str.endswith(
        "GA"
    )
    df[f"feat_siRNA_{name}_seq_pattern_7"] = s.str.startswith("UU") & s.str.endswith(
        "CA"
    )
    df[f"feat_siRNA_{name}_seq_pattern_8"] = s.str.startswith("UU") & s.str.endswith(
        "UA"
    )
    # 第二位和倒数第二位是否为A
    df[f"feat_siRNA_{name}_seq_pattern_9"] = s.str[1] == "A"
    df[f"feat_siRNA_{name}_seq_pattern_10"] = s.str[-2] == "A"
    # GC占整体长度的比例
    df[f"feat_siRNA_{name}_seq_pattern_GC_frac"] = (
        s.str.contains("G") + s.str.contains("C")
    ) / s.str.len()
    return df.iloc[:, 1:]

基于lightgbm的baseline

在得到了表格数据之后，我们可以使用任意适用于表格数据的机器学习回归模型来进行预测，此处我们简单使用了lightgbm模型：

train_data = lgb.Dataset(X_train, label=y_train)
test_data = lgb.Dataset(X_test, label=y_test, reference=train_data)

def print_validation_result(env):
    result = env.evaluation_result_list[-1]
    print(f"[{env.iteration}] {result[1]}'s {result[0]}: {result[2]}")

params = {
    "boosting_type": "gbdt",
    "objective": "regression",
    "metric": "root_mean_squared_error",
    "max_depth": 7,
    "learning_rate": 0.02,
    "verbose": 0,
}

gbm = lgb.train(
    params,
    train_data,
    num_boost_round=15000,
    valid_sets=[test_data],
    callbacks=[print_validation_result],
)

下面给出完整的代码并添加注释

import pandas as pd

# 读取训练数据和提交样本数据
df_original = pd.read_csv("C:/Users/32084/Desktop/train_data.csv")
n_original = df_original.shape[0]  # 获取原始数据的行数
df_submit = pd.read_csv("C:/Users/32084/Desktop/sample_submission.csv")

# 合并数据并重置索引
df = pd.concat([df_original, df_submit], axis=0).reset_index(drop=True)

# 构建siRNA特征的函数
def siRNA_feat_builder(s: pd.Series, anti: bool = False):
    name = "anti" if anti else "sense"  # 根据参数设置名称
    df = s.to_frame()  # 将系列转换为DataFrame
    df[f"feat_siRNA_{name}_seq_len"] = s.str.len()  # 序列长度特征

    # 前后位置的碱基特征
    for pos in [0, -1]:
        for c in list("AUGC"):
            df[f"feat_siRNA_{name}_seq_{c}_{'front' if pos == 0 else 'back'}"] = (
                s.str[pos] == c
            )
    
    # 不同的特定序列模式特征
    df[f"feat_siRNA_{name}_seq_pattern_1"] = s.str.startswith("AA") & s.str.endswith("UU")
    df[f"feat_siRNA_{name}_seq_pattern_2"] = s.str.startswith("GA") & s.str.endswith("UU")
    df[f"feat_siRNA_{name}_seq_pattern_3"] = s.str.startswith("CA") & s.str.endswith("UU")
    df[f"feat_siRNA_{name}_seq_pattern_4"] = s.str.startswith("UA") & s.str.endswith("UU")
    df[f"feat_siRNA_{name}_seq_pattern_5"] = s.str.startswith("UU") & s.str.endswith("AA")
    df[f"feat_siRNA_{name}_seq_pattern_6"] = s.str.startswith("UU") & s.str.endswith("GA")
    df[f"feat_siRNA_{name}_seq_pattern_7"] = s.str.startswith("UU") & s.str.endswith("CA")
    df[f"feat_siRNA_{name}_seq_pattern_8"] = s.str.startswith("UU") & s.str.endswith("UA")
    df[f"feat_siRNA_{name}_seq_pattern_9"] = s.str[1] == "A"
    df[f"feat_siRNA_{name}_seq_pattern_10"] = s.str[-2] == "A"

    # GC含量特征
    df[f"feat_siRNA_{name}_seq_pattern_GC_frac"] = (
        s.str.count("G") + s.str.count("C")
    ) / s.str.len()
    
    return df.iloc[:, 1:]  # 返回特征列

# 对publication_id进行独热编码
df_publication_id = pd.get_dummies(df.publication_id)
df_publication_id.columns = [f"feat_publication_id_{c}" for c in df_publication_id.columns]

# 对gene_target_symbol_name进行独热编码
df_gene_target_symbol_name = pd.get_dummies(df.gene_target_symbol_name)
df_gene_target_symbol_name.columns = [f"feat_gene_target_symbol_name_{c}" for c in df_gene_target_symbol_name.columns]

# 对gene_target_ncbi_id进行独热编码
df_gene_target_ncbi_id = pd.get_dummies(df.gene_target_ncbi_id)
df_gene_target_ncbi_id.columns = [f"feat_gene_target_ncbi_id_{c}" for c in df_gene_target_ncbi_id.columns]

# 对gene_target_species进行独热编码
df_gene_target_species = pd.get_dummies(df.gene_target_species)
df_gene_target_species.columns = [f"feat_gene_target_species_{c}" for c in df_gene_target_species.columns]

# 标准化siRNA_duplex_id
siRNA_duplex_id_values = df.siRNA_duplex_id.str[3:-2].str.strip(".").astype("int")
siRNA_duplex_id_values = (siRNA_duplex_id_values - siRNA_duplex_id_values.min()) / (
    siRNA_duplex_id_values.max() - siRNA_duplex_id_values.min()
)
df_siRNA_duplex_id = pd.DataFrame(siRNA_duplex_id_values)

# 对cell_line_donor进行独热编码
df_cell_line_donor = pd.get_dummies(df.cell_line_donor)
df_cell_line_donor.columns = [f"feat_cell_line_donor_{c}" for c in df_cell_line_donor.columns]

# 增加特定的细胞系特征
df_cell_line_donor["feat_cell_line_donor_hepatocytes"] = (
    (df.cell_line_donor.str.contains("Hepatocytes")).fillna(False).astype("int")
)
df_cell_line_donor["feat_cell_line_donor_cells"] = (
    df.cell_line_donor.str.contains("Cells").fillna(False).astype("int")
)

# 将siRNA浓度转换为DataFrame
df_siRNA_concentration = df.siRNA_concentration.to_frame()

# 对Transfection_method进行独热编码
df_Transfection_method = pd.get_dummies(df.Transfection_method)
df_Transfection_method.columns = [f"feat_Transfection_method_{c}" for c in df_Transfection_method.columns]

# 对Duration_after_transfection_h进行独热编码
df_Duration_after_transfection_h = pd.get_dummies(df.Duration_after_transfection_h)
df_Duration_after_transfection_h.columns = [
    f"feat_Duration_after_transfection_h_{c}" for c in df_Duration_after_transfection_h.columns
]

# 合并所有特征
feats = pd.concat(
    [
        df_publication_id,
        df_gene_target_symbol_name,
        df_gene_target_ncbi_id,
        df_gene_target_species,
        df_siRNA_duplex_id,
        df_cell_line_donor,
        df_siRNA_concentration,
        df_Transfection_method,
        df_Duration_after_transfection_h,
        siRNA_feat_builder(df.siRNA_sense_seq, False),
        siRNA_feat_builder(df.siRNA_antisense_seq, True),
        df.iloc[:, -1].to_frame(),
    ],
    axis=1,
)

import lightgbm as lgb
from sklearn.model_selection import train_test_split

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(
    feats.iloc[:n_original, :-1],
    feats.iloc[:n_original, -1],
    test_size=0.25,
    random_state=42,
)

# 构建LightGBM数据集
train_data = lgb.Dataset(X_train, label=y_train)
test_data = lgb.Dataset(X_test, label=y_test, reference=train_data)

# 打印验证结果的回调函数
def print_validation_result(env):
    result = env.evaluation_result_list[-1]
    print(f"[{env.iteration}] {result[1]}'s {result[0]}: {result[2]}")

# 设置LightGBM参数
params = {
    "boosting_type": "gbdt",
    "objective": "regression",
    "metric": "root_mean_squared_error",
    "max_depth": 7,
    "learning_rate": 0.02,
    "verbose": 0,
}

# 训练模型
gbm = lgb.train(
    params,
    train_data,
    num_boost_round=20000,
    valid_sets=[test_data],
    callbacks=[print_validation_result],
)

# 预测结果
y_pred = gbm.predict(feats.iloc[n_original:, :-1])

# 保存预测结果到CSV文件
df_submit["mRNA_remaining_pct"] = y_pred
df_submit.to_csv("submission.csv", index=False)

重新调整参数进行修改改进不断提高分数：

考虑处理缺失值和标准化以及早停法（正在尝试~）

思考体验

Task 2 详解解释了baseline的思路及其数据原理，更加深入了解数据特征工程的重要意义。同时又见到了第二期学习中的老朋友lightgbm，加深了理解。目前卡在上分瓶颈，期待进一步的学习~~

Task 3：特征工程进阶

生物学角度新特征

使用反义链与target gene序列的序列匹配结果作为特征来增强模型表现。在siRNA沉默target gene的过程中，RNA酶III家族的Dicer与TAR RNA结合蛋白和Argonaute（Ago）共同作用，产生长度为21–25个核苷酸的siRNA。siRNA做为RNA诱导沉默复合体（RISC）的组成部分，能够识别靶mRNA，从而达到沉默目标基因的作用。但是siRNA与靶mRNA的互补方式的不同意味着沉默机制的不同：siRNA与mRNA的3′非翻译区（UTR）结合，在Ago1、Ago3和Ago4的帮助下会导致翻译抑制；但如果siRNA与编码序列（CDS）完全互补，靶mRNA将被Ago2的内切核酸酶活性切割。而本比赛的预测目标是mRNA的保留水平，因此siRNA反义链与target gene的匹配程度，以及匹配位置，都会对预测目标产生影响。

  参考文献：

  Paddison PJ. RNA interference in mammalian cell systems. Curr Top Microbiol Immunol 2008; 320: 1–19. 12 Ji X.

  The mechanism of RNase III action: how dicer dices. Curr Top Microbiol Immunol 2008; 320: 99–116. 13

  Rivas FV, Tolia NH, Song J-J, Aragon JP, Liu J, Hannon GJ et al. Purified Argonaute2 and an siRNA form recombinant human RISC. Nat Struct Mol Biol 2005; 12: 340–349. 14

  Su H, Trombly MI, Chen J, Wang X. Essential and overlapping functions for mammalian Argonautes in microRNA silencing. Genes Dev 2009; 23: 304–317.
GC含量是siRNA效率中的一个重要且基本的参数，可以作为模型预测的特征。这是因为低GC含量会导致非特异性和较弱的结合，而高GC含量可能会阻碍siRNA双链在解旋酶和RISC复合体作用下的解旋。生物学研究人员曾经通过大量实验确定了siRNA的GC含量的可接受区间，有人提出过31.6%到57.9%这个区间，也有人提出过36%到52%这个区间。此外，也有研究者更深入地地探讨了这个问题，他们发现在反义链中，第2到第7个核苷酸和第8到第18个核苷酸的GC百分比分别为19%和52%是理想的。此外，功能性siRNA在第9到第14个核苷酸之间有一个不稳定区域（GC含量低于其他区域），被称为能量谷，这是选择siRNA的重要标准。这种内部的不稳定性通过在mRNA剪切过程中诱导最理想的构象，从而提高了RISC复合体的功能性。从上述生物学研究中我们可以发现，siRNA的GC含量，或者片段中特定位置的GC含量，对其沉默效率至关重要。

  参考文献：

  Reynolds A, Leake D, Boese Q, Scaringe S, Marshall WS, Khvorova A. Rational siRNA design for RNA interference. Nat Biotechnol 2004; 22: 326–330. 40.

  Ui‐Tei K, Naito Y, Takahashi F, Haraguchi T, Ohki‐Hamazaki H, Juni A et al. Guidelines for the selection of highly effective siRNA sequences for mammalian and chick RNA interference. Nucleic Acids Res 2004; 32: 936–948. 41

  Amarzguioui M, Prydz H. An algorithm for selection of functional siRNA sequences. Biochem Biophys Res Commun 2004; 316: 1050–1058.

  Khvorova A, Reynolds A, Jayasena SD. Functional siRNAs and miRNAs exhibit strand bias. Cell 2003; 115: 209–216.
将修饰过的碱基序列也进行编码，简单的编码方式是将带有修饰的核苷酸编码为和普通核苷酸不一样的输入向量，复杂的编码方式是将不同修饰在化学上的差异也加入模型中。siRNA上核苷酸的化学修饰对于siRNA发挥其功能至关重要。在siRNA治疗技术研发的早期阶段，siRNA药物都是未加修饰的，siRNA可以在体内介导基因沉默，但是可能会出现较差治疗效果和潜在的非靶向效应。化学修饰的siRNA，例如用2′-O-甲基（2′-OMe）或2′-甲氧乙基（2′-MOE）取代2′-OH，或用locked nucleic acid、unlocked nucleic acid或glycol nucleic acid取代某些核苷酸，可以有效抑制由siRNA引发的免疫刺激性内源免疫激活，增强活性和特异性，并减少非靶向诱导的毒性。根据核苷酸的自然结构，化学修饰可以施加在磷酸骨架、核糖部分或碱基上。通常，这些修饰会同时引入到siRNA中。例如，2′-OMe和磷硫酸酯（PS）修饰的结合有助于胆固醇结合的siRNA的系统性给药，并实现体内高效的基因沉默。从siRNA药物的研发经验来说，对siRNA的精确修饰可以提高其沉默效率、特异性和稳定性，并减少其毒性和免疫原性。

  参考文献：

  Ju¨rgen Soutschek, A. A. et al. Therapeutic silencing of an endogenous gene by systemic administration of modified siRNAs. nature432, 173–178 (2004).

  Khvorova, A. & Watts, J. K. The chemical evolution of oligonucleotide therapies of clinical utility. Nat. Biotechnol.35, 238–248 (2017).

代码实现

新特征引入

修饰siRNA构建特征

修饰siRNA序列进行n-gram的词频统计，同时也可对未修饰序列进行相同的操作。

class GenomicTokenizer:
    def __init__(self, ngram=5, stride=2):
        # 初始化分词器，设置n-gram长度和步幅
        self.ngram = ngram
        self.stride = stride
        
    def tokenize(self, t):

        # 字符串变list
        if isinstance(t, str):
            t = list(t)

        if self.ngram == 1:
            # 如果n-gram长度为1，直接将序列转换为字符列表
            toks = 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]

            # sub list to str
            toks = [''.join(x) for x in toks]

        # 返回分词结果
        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]
        itos = [o for o, c in freq.most_common(max_vocab - 1) if c >= min_freq]
        # 返回包含词汇表的类实例
        return cls(itos)
    
def siRNA_feat_builder_substr(se, name, patterns):
    
    # 创建一个空字典来存储特征
    features = {}

    for pattern in patterns:
        try:
            # escaped_pattern = re.escape(pattern)  # 转义模式中的特殊字符
            escaped_pattern = pattern
            features[f"feat_{name}_seq_pattern_{escaped_pattern}"] = se.str.count(escaped_pattern)
        except re.error as e:
            print(f"Error in pattern {pattern}: {e}")

    # 将字典转换为DataFrame
    feature_df = pd.DataFrame(features)

    return feature_df

siRNA序列与target序列对比

siRNA与target序列的结合能力是一个很重要的特征，其中两者的序列比对score可作为对这个特征的定量化描述（测试了500个数据）

参考：https://biopython.org/docs/1.76/api/Bio.pairwise2.html

lgm优化

低Remaining范围样本高权重

weight_ls = np.array(feats['mRNA_remaining_pct'].apply(lambda x:2 if ((x<=30)and(x>=0)) else 1))

使用官方评价指标作为损失函数

由原来的root_mean_squared_error评价指标被替换为更加复杂的官方评价分数，具体公式为:

$\text{score} = 50\% \times \left(1 - \frac{\text{MAE}}{100}\right) + 50\% \times F1 \times \left(1 - \frac{\text{Range-MAE}}{100}\right)$

# calculate_metrics函数用于计算评估指标
def calculate_metrics(preds, data, threshold=30):
    y_pred = preds
    y_true = data.get_label()
    mae = np.mean(np.abs(y_true - y_pred))
    # if mae < 0: mae = 0
    # elif mae >100: mae = 100

    y_true_binary = ((y_true <= threshold) & (y_true >= 0)).astype(int)
    y_pred_binary = ((y_pred <= threshold) & (y_pred >= 0)).astype(int)

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

    # precision = precision_score(y_true_binary, y_pred_binary, average="binary")
    # recall = recall_score(y_true_binary, y_pred_binary, average="binary")

    if np.sum(y_pred_binary) > 0:
        precision = (np.array(y_pred_binary) & y_true_binary).sum()/np.sum(y_pred_binary)
    else:
        precision = 0
    if np.sum(y_true_binary) > 0:
        recall = (np.array(y_pred_binary) & y_true_binary).sum()/np.sum(y_true_binary)
    else:
        recall = 0

    if precision + recall == 0:
        f1 = 0
    else:
        f1 = 2 * precision * recall / (precision + recall)
    score = (1 - mae / 100) * 0.5 + (1 - range_mae / 100) * f1 * 0.5
    return "custom_score", score, True  # True表示分数越高越好

自适应学习率

# 定义自适应学习率回调函数，用于在训练过程中根据性能调整学习率
def adaptive_learning_rate(decay_rate=0.8, patience=50):
    # 初始化最佳分数为负无穷大，因为我们假设分数越高代表性能越好
    best_score = float("-inf")
    # 初始化等待计数器，用于实现patience逻辑
    wait = 0

    # 定义一个嵌套的回调函数，这个函数将在训练的每一轮调用
    def callback(env):
        # 使用nonlocal关键字，使得内部函数能够访问并修改外部函数的变量
        nonlocal best_score, wait
        # 从环境对象中获取当前轮的评估结果，这里假设最后一个评估指标是我们关心的指标
        current_score = env.evaluation_result_list[-1][2]
        # 获取当前的学习率
        current_lr = env.model.params.get('learning_rate')

        # 如果当前分数高于最佳分数，则更新最佳分数
        if current_score > best_score: 
            best_score = current_score
            # 注释掉的代码可能是用于实现连续性能提升的逻辑
            # wait = 0
        else:
            # 如果当前分数不高于最佳分数，增加等待计数器
            wait += 1

        # 如果等待计数器达到patience设定的值，则降低学习率
        if wait >= patience:
            # 计算新的学习率，为当前学习率乘以衰减率
            new_lr = float(current_lr) * decay_rate
            # 重置等待计数器
            wait = 0
            # 设置模型的新学习率
            env.model.params['learning_rate'] = new_lr
            # 打印调整后的学习率
            print(f"Learning rate adjusted to {env.model.params.get('learning_rate')}")

    # 返回回调函数
    return callback

这个函数的目的是创建一个回调函数，可以在训练过程中根据模型的性能动态调整学习率。如果模型在连续的 patience 轮中性能没有提升，则学习率会按照 decay_rate 指定的衰减率进行调整。这是一种常见的技术，用于防止训练过程中的过拟合，并帮助模型在训练后期更细致地逼近最优解。

多折交叉训练

# train函数用于训练模型
def train(feats, n_original):
    # 定义k折交叉验证
    n_splits = 10
    kf = KFold(n_splits=n_splits, shuffle=True, random_state=42)
    # 开始k折交叉验证
    gbms = []
    best_scores = []
    for fold, (train_idx, val_idx) in enumerate(
        kf.split(feats.iloc[:n_original, :]), 1
    ):
        # 准备训练集和验证集
        X_train, X_val = feats.iloc[train_idx, :-1], feats.iloc[val_idx, :-1]
        y_train, y_val = feats.iloc[train_idx, -1], feats.iloc[val_idx, -1]
        w_train = weight_ls[train_idx]
        

        # 创建LightGBM数据集
        train_data = lgb.Dataset(X_train, label=y_train, weight=w_train)
        val_data = lgb.Dataset(X_val, label=y_val, reference=train_data)

        boost_round = 25000
        early_stop_rounds = int(boost_round*0.1)

        # 显示metric
        lgb_log = lgb.log_evaluation(period=200, show_stdv=True)
        lgb_stop = lgb.early_stopping(stopping_rounds=early_stop_rounds, first_metric_only=True, verbose=True, min_delta=0.00001)

        # 设置LightGBM参数
        params = {
            "boosting_type": "gbdt",
            "objective": "regression",
            "metric": "None",
            # "metric": "root_mean_squared_error",
            "max_depth": 8,
            "num_leaves": 63,
            "min_data_in_leaf": 2,
            "learning_rate": 0.05,
            "feature_fraction": 0.9,
            "lambda_l1": 0.1,
            "lambda_l2": 0.2,
            "verbose": -1, # -1时不输出
            "early_stopping_round": early_stop_rounds,
            "num_threads": 8,
        }

        # 在训练时使用自适应学习率回调函数
        adaptive_lr = adaptive_learning_rate(decay_rate=0.9, patience=1000)
        gbm = lgb.train(
            params,
            train_data,
            num_boost_round=boost_round,
            valid_sets=[val_data],
            feval=calculate_metrics,  # 将自定义指标函数作为feval参数传入
            # callbacks=[print_validation_result, adaptive_lr, lgb_log, lgb_stop],
            callbacks=[adaptive_lr, lgb_log, lgb_stop],
        )
        valid_score = gbm.best_score["valid_0"]["custom_score"]
        best_scores.append(valid_score)
        print(f"best_valid_score: {valid_score}")
        gbms.append(gbm)
    average_best_scores = sum(best_scores) / len(best_scores)

    # 选择最佳模型
    best_model_index = best_scores.index(max(best_scores))
    best_model = gbms[best_model_index]
    return best_model

利用平均最佳分数选择最优模型

average_best_scores = sum(best_scores) / len(best_scores)

K折交叉验证（K-Fold Cross-Validation）是一种在机器学习中用于评估模型泛化能力的技术。它可以帮助我们理解模型在未知数据上的表现，减少过拟合的风险，并为模型的选择和参数调整提供依据。下面是K折交叉验证的步骤：

数据分割：首先，将可用的数据集分割成K个大小相等（或尽可能相等）的子集（称为“折”或“fold”）。

循环验证：然后，模型训练和验证的过程会重复K次。在每次迭代中：

选择其中一个子集作为验证集（即用于测试模型的性能）。
剩下的K-1个子集合并起来作为训练集。

训练和评估：在每次迭代中，模型使用训练集进行训练，然后在验证集上评估其性能。这样可以得到K个性能评估结果。

结果分析：最后，通过分析这K次迭代的性能评估结果，可以得到模型的平均性能和性能波动。通常，我们计算所有K次迭代结果的平均值作为模型的最终性能估计，并且使用标准差来衡量模型性能的稳定性。

K折交叉验证的优点包括：

更准确的性能估计：相比于单次划分训练集和测试集，K折交叉验证可以更准确地估计模型在未知数据上的表现。
充分利用数据：特别是在数据量有限的情况下，K折交叉验证可以确保每条数据都被用于训练和验证，从而充分利用可用数据。
模型选择和调参：可以用来比较不同模型或参数设置的性能，帮助选择最佳模型。

K折交叉验证也有一些缺点：

计算成本：由于模型需要训练K次，因此计算成本较高。
时间消耗：对于大型数据集或复杂模型，K折交叉验证可能需要较长的时间来完成。

模型结构和参数调整的上分思路

1. 优化超参数

a. 网格搜索（Grid Search）：

参考代码如下：

from sklearn.model_selection import GridSearchCV
from lightgbm import LGBMRegressor

param_grid = {
    'max_depth': [3, 5, 7, 9],
    'learning_rate': [0.01, 0.02, 0.05, 0.1],
    'n_estimators': [100, 500, 1000, 2000],
    'min_child_samples': [20, 30, 50, 100]
}

lgbm_model = LGBMRegressor(objective='regression', metric='rmse')
grid_search = GridSearchCV(estimator=lgbm_model, param_grid=param_grid, cv=5, n_jobs=-1, verbose=2)
grid_search.fit(X_train, y_train)

print("Best parameters found: ", grid_search.best_params_)

b. 随机搜索（Random Search）：

参考代码如下：

from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import randint as sp_randint

param_dist = {
    'max_depth': sp_randint(3, 10),
    'learning_rate': [0.01, 0.02, 0.05, 0.1],
    'n_estimators': sp_randint(100, 2000),
    'min_child_samples': sp_randint(20, 100)
}

random_search = RandomizedSearchCV(estimator=lgbm_model, param_distributions=param_dist, n_iter=100, cv=5, n_jobs=-1, verbose=2)
random_search.fit(X_train, y_train)

print("Best parameters found: ", random_search.best_params_)

c. 贝叶斯优化：

使用如optuna这样的库来执行贝叶斯优化, 参考代码如下：

import optuna

def objective(trial):
    params = {
        'max_depth': trial.suggest_int('max_depth', 3, 10),
        'learning_rate': trial.suggest_loguniform('learning_rate', 1e-3, 1e-1),
        'n_estimators': trial.suggest_int('n_estimators', 100, 2000),
        'min_child_samples': trial.suggest_int('min_child_samples', 20, 100)
    }
    
    model = LGBMRegressor(**params)
    model.fit(X_train, y_train)
    return model.score(X_val, y_val)

study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=100)

print('Best trial:')
trial = study.best_trial
print('Value: {}'.format(trial.value))
print('Params: ')
for key, value in trial.params.items():
    print('    {}: {}'.format(key, value))

我尝试了贝叶斯优化，但是运行时间耗时比较长

import optuna
import lightgbm as lgb
from sklearn.model_selection import KFold

# 定义一个优化的目标函数
def objective(trial, feats, n_original):
    # 提取超参数
    max_depth = trial.suggest_int('max_depth', 3, 10),
    num_leaves = trial.suggest_int("num_leaves", 100, 1023)
    learning_rate = trial.suggest_loguniform('learning_rate', 1e-3, 1e-1)
    n_estimators = trial.suggest_int('n_estimators', 100, 2000)
    min_child_samples = trial.suggest_int('min_child_samples', 20, 100)

    # 定义k折交叉验证
    kf = KFold(n_splits=10, shuffle=True, random_state=42)
    best_scores = []

    for fold, (train_idx, val_idx) in enumerate(kf.split(feats.iloc[:n_original, :]), 1):
        # 准备训练集和验证集
        X_train, X_val = feats.iloc[train_idx, :-1], feats.iloc[val_idx, :-1]
        y_train, y_val = feats.iloc[train_idx, -1], feats.iloc[val_idx, -1]
        w_train = weight_ls[train_idx]

        # 创建LightGBM数据集
        train_data = lgb.Dataset(X_train, label=y_train, weight=w_train)
        val_data = lgb.Dataset(X_val, label=y_val, reference=train_data)
        
        boost_round = 25000
        early_stop_rounds = int(boost_round*0.1)

        # 显示metric
        lgb_log = lgb.log_evaluation(period=200, show_stdv=True)
        lgb_stop = lgb.early_stopping(stopping_rounds=early_stop_rounds, first_metric_only=True, verbose=True, min_delta=0.00001)
        
        # 设置LightGBM参数
        params = {
            "boosting_type": "gbdt",
            "objective": "regression",
            "max_depth": max_depth,
            "num_leaves": num_leaves,
            "learning_rate": learning_rate,
            "metric": "None",
            "min_data_in_leaf": 2,
            "feature_fraction": 0.9,
            "lambda_l1": 0.1,
            "lambda_l2": 0.2,
            "verbose": -1, # -1时不输出
            "early_stopping_round": early_stop_rounds,
            "num_threads": 8,
            "use_gpu": True,
        }
        
        adaptive_lr = adaptive_learning_rate(decay_rate=0.9, patience=1000)

        # 训练模型
        gbm = lgb.train(
            params,
            train_data,
            num_boost_round=boost_round,
            valid_sets=[val_data],
            feval=calculate_metrics,  # 将自定义指标函数作为feval参数传入
            # callbacks=[print_validation_result, adaptive_lr, lgb_log, lgb_stop],
            callbacks=[adaptive_lr, lgb_log, lgb_stop],
        )

        # 更新最佳分数
        valid_score = gbm.best_score["valid_0"]["custom_score"]
        best_scores.append(valid_score)
        print(f"best_valid_score: {valid_score}")
    average_best_scores = sum(best_scores) / len(best_scores)

    return average_best_scores

# 创建 Optuna 研究对象
study = optuna.create_study(direction="minimize")

# 执行研究（超参数搜索）
study.optimize(lambda trial: objective(trial, feats, n_original), n_trials=100)

# 获取最佳参数
best_params = study.best_trial.params
print("Best parameters: ", best_params)

2. 集成学习

集成学习是一种机器学习范式，通过结合多个模型来改善预测性能，通常比单个模型的表现更优。在你的项目中，可以采用以下几种常见的集成学习方法：

Bagging（自助聚合）：通过在原始数据集上进行多次重采样来创建多个子集，分别训练多个模型，最后进行平均或多数投票决策。
Boosting：训练多个模型，每个模型都尝试纠正前一个模型的错误，通常是序列处理。
Stacking：训练多个不同的模型，然后再训练一个新的模型来综合这些模型的输出。（在第二期的学习中曾尝试使用过）

3. 混合学习

在解决复杂的生物信息学问题，如siRNA效果的预测中，机器学习与深度学习的混合方法可以提供强大的工具。这种混合方法可以结合深度学习的特征学习能力和传统机器学习模型的效率与解释性。以下是一种结合机器学习和深度学习的策略，包括两个主要部分：一个用于特征提取的深度学习模型（如卷积神经网络），以及一个传统的机器学习模型（如LightGBM）来进行最终的决策。这里提供一个示例，展示如何使用PyTorch构建深度学习部分，然后将输出特征传递给LightGBM进行分类或回归。

思考体验

Task 3 整体上有些难度，下来查询了一些资料辅助理解并且尝试完成了代码实践，提高了一定的分数。在运行siRNA和target结合分析模块耗费了大量时间，后期超参数优化也做了一部分，目前的最高分数是0.8023，后期还会对代码进一步完善改进，但是今天是第三期截止日了，只能展示目前的最好成绩~~

JeasonLi ——Wenrui

关注

30
点赞
踩
10

收藏

觉得还不错? 一键收藏
0
评论
#Datawhale AI 夏令营第三期—siRNA药物药效预测学习

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

扫一扫