siRNA药物药效预测（一）（Datawhale AI 夏令营）

最新推荐文章于 2024-09-14 16:43:45 发布

数学难

最新推荐文章于 2024-09-14 16:43:45 发布

阅读量1k

点赞数 17

文章标签：人工智能

本文链接：https://blog.csdn.net/m0_72765822/article/details/140782246

版权

赛事链接：http://competition.sais.com.cn/competitionDetail/532230/format

一、赛题解析

1、基本概念

1、RNA干扰

RNA干扰（RNAi）是一种利用双链RNA（dsRNA）特异性降解同源mRNA，从而抑制基因表达的现象。自发现以来，RNAi因其高效性和特异性在基因功能研究、疾病治疗及药物研发等领域得到广泛应用。通过化学合成或体外转录制备的短双链RNA（siRNA）被导入细胞，与靶mRNA结合并导致其降解，从而精确调控基因表达。RNAi技术不仅帮助科学家揭示基因功能和信号通路，还促进了疾病模型的构建和基因治疗的发展。尽管面临稳定性、细胞毒性和非特异性干扰等挑战，但随着技术的不断进步，RNAi有望在更多领域发挥重要作用。

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

2、化学修饰siRNA

化学修饰siRNA是一种关键的技术手段，它通过化学方法调整siRNA分子的结构，旨在增强其稳定性、特异性、生物活性，并降低潜在的毒性。这些修饰包括改变磷酸二酯键的组成（如引入硫代磷酸酯），核苷酸碱基的替代（如2'-氧甲基化），以及siRNA末端的特殊基团添加等，旨在保护siRNA免受体内核酸酶的降解，同时提高其与靶标mRNA的结合效率和特异性。化学修饰siRNA在基因功能研究、疾病治疗及药物研发中展现出广泛应用前景，特别是在提高基因沉默效果、延长药物作用时间和降低治疗副作用方面具有重要价值。然而，在选择和应用这些修饰方式时，需综合考虑其对siRNA性能的潜在影响，以确保最佳的研究效果和临床安全性。

2、赛题介绍

本次比赛旨在利用机器学习技术，预测化学修饰后的siRNA序列在RNA干扰（RNAi）机制下对靶基因的沉默效率。RNAi是一种重要的基因表达调控机制，通过干扰特定基因的表达，可以用于疾病治疗。这次比赛的目标是通过构建并优化模型，准确预测siRNA的沉默效率，从而提升药物设计的效率和效果。

比赛流程

初赛阶段：提供一部分公开文献中提取的siRNA修饰序列和实验数据，参赛者需要使用这些数据训练模型并提交预测结果。初赛的重点是评估在训练集中出现过的目标mRNA序列，不同siRNA的沉默效率预测的准确性。
复赛阶段：在初赛基础上，增加部分尚未公开的专利数据作为测试数据，评估模型在未见过的目标mRNA序列上的预测准确性。

数据集

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

3、评分体系

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

平均绝对误差（MAE）
预测值在一定范围内的平均绝对误差（Range-MAE）
预测值在一定范围内的F1指标（F1）

4、训练集结构

（官方资料）

表头字段	表头说明
id	数据唯一识别号
publication_id	公开文献号，已经过脱敏处理
gene_target_symbol_name	靶基因符号名称，靶基因即为siRNA想要沉默的目标基因
gene_target_ncbi_id	靶基因的NCBI标识
gene_target_species	靶基因参考序列的物种，一般物种应为人类
siRNA_duplex_id	siRNA双链编号
siRNA_sense_seq	siRNA的sense序列
siRNA_antisense_seq	siRNA的antisense序列
cell_line_donor	实验使用的细胞系，一般认为对沉默效率测定存在较大影响
siRNA_concentration	实验使用的siRNA浓度，一般认为对沉默效率测定存在较大影响
concentration_unit	siRNA浓度单位，比如nM
Transfection_method	转染方法，即将siRNA植入细胞的方法，一般认为对沉默效率测定存在较大影响
Duration_after_transfection_h	转染后持续时间（小时）
modified_siRNA_sense_seq	带修饰的siRNA的sense序列
modified_siRNA_antisense_seq	带修饰的siRNA的antisense序列
modified_siRNA_sense_seq_list	带修饰的siRNA的sense序列分解列表，基于标准词表将sense序列分解，以空格分隔
modified_siRNA_antisense_seq_list	带修饰的siRNA的antisense序列分解列表，基于标准词表将antisense序列分解，以空格分隔
gene_target_seq	靶基因的参考序列，siRNA一般靶向序列中的一部分片段以达到沉默的效果
mRNA_remaining_pct	实验后mRNA的剩余百分比

备注：在实际工业制备过程中，siRNA的双链不一定完全对齐，可能存在1-3个不等的核苷酸数量差距。siRNA裸序列一般由AUGC4类核苷酸构成，但也有少量T的存在。上述都是真实实验数据。

二、官方baseline分析

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import Dataset, DataLoader
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from tqdm import tqdm
from collections import Counter
from rich import print
from sklearn.metrics import precision_score, recall_score, mean_absolute_error

# GenomicTokenizer类，它用于对输入的字符串进行分词处理。这个类通过ngram和stride两个参数来控制分词的方式
class GenomicTokenizer:
    # 用于设置类的实例变量ngram和stride
    def __init__(self, ngram=5, stride=2):
        self.ngram = ngram
        self.stride = stride

    # 用于对输入的字符串t进行分词处理
    def tokenize(self, t):
        # 字符串大写化
        t = t.upper()
        # 处理 ngram 为 1 的特殊情况：如果ngram的值为1，意味着每个分词只包含一个字符。此时，方法直接将字符串t转换为字符列表toks。
        if self.ngram == 1:
            toks = list(t)
        # 生成分词：如果ngram的值大于1，方法会通过一个列表推导式来生成分词。这个推导式遍历字符串t，从索引i开始，每次跳跃stride个字符，然后截取从i开始、长度为ngram的子字符串作为分词
        # 但是，这里有一个条件判断if len(t[i:i+self.ngram]) == self.ngram，确保截取的子字符串长度确实等于ngram
        # 这个条件在stride小于ngram时特别有用，可以防止在字符串末尾生成不完整的分词
        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]
        # 去除不完整的分词：由于列表推导式可能会生成一个长度小于ngram的分词（如果字符串t的长度不是stride的整数倍，且stride小于ngram），方法会检查toks列表的最后一个元素，如果其长度小于ngram，则将其从列表中移除。
        if len(toks[-1]) < self.ngram:
            toks = toks[:-1]
        return toks

# GenomicVocab类是一个用于处理基因组词汇的类，它提供了词汇到索引（itos）和索引到词汇（stoi）的映射。
class GenomicVocab:
    def __init__(self, itos):
        # itos：一个列表，其中包含了词汇到索引的映射。索引从0开始，并且通常会在列表的开始添加一个特殊的词汇（如<pad>），用于表示填充或空值。
        self.itos = itos
        # stoi：一个字典，其中包含了索引到词汇的映射。它是通过遍历itos列表并反转键值对来创建的。
        self.stoi = {v:k for k,v in enumerate(self.itos)}

    # tokens（一个包含词汇的列表或可迭代对象）
    # max_vocab（词汇表的最大大小）
    # min_freq（词汇的最低出现频率）
    @classmethod
    def create(cls, tokens, max_vocab, min_freq):
        # tokens（一个包含词汇的列表或可迭代对象），max_vocab（词汇表的最大大小），和min_freq（词汇的最低出现频率）
        freq = Counter(tokens)
        # 创建一个新的 itos 列表，它首先包含一个特殊的填充词汇 <pad>，然后添加出现频率最高的 max_vocab-1 个词汇（忽略出现频率低于 min_freq 的词汇）
        itos = ['<pad>'] + [o for o,c in freq.most_common(max_vocab-1) if c >= min_freq]
        return cls(itos)


# 定义了一个名为 SiRNADataset 的类，它继承自 PyTorch 的 Dataset 类，用于处理与 siRNA（小干扰RNA）相关的数据集。
# 这个类的主要目的是将原始的、可能包含多种格式的数据（如文本序列）转换为模型训练所需的格式（即张量）
class SiRNADataset(Dataset):
    def __init__(self, df, columns, vocab, tokenizer, max_len):
        self.df = df    # 一个 pandas DataFrame，包含了数据集的所有数据
        self.columns = columns  # 一个列表，指定了 DataFrame 中哪些列（特征）将被用于模型训练。
        self.vocab = vocab  # 一个词汇表对象，通常包含了一个将单词映射到整数的字典（stoi）和一个将整数映射回单词的字典（虽然在这个类中只使用了 stoi）。
        self.tokenizer = tokenizer  # 一个分词器对象，用于将原始文本序列分割成单词或标记。
        self.max_len = max_len  # 一个整数，指定了每个序列在编码后应该被填充或截断到的最大长度。

    # 返回数据集中样本的总数，即 DataFrame 的行数。这对于模型训练过程中的迭代控制非常重要
    def __len__(self):
        return len(self.df)

    # 根据给定的索引 idx，从 DataFrame 中提取相应的行（样本），并对该样本中指定的列（特征）进行分词、编码和填充
    # 同时提取目标值（'mRNA_remaining_pct' 列的值）。最后，返回一个包含编码后的序列和目标值的元组。
    def __getitem__(self, idx):
        row = self.df.iloc[idx]
        seqs = [self.tokenize_and_encode(row[col]) for col in self.columns]
        target = torch.tensor(row['mRNA_remaining_pct'], dtype=torch.float)
        return seqs, target

    # 这个方法首先检查序列中是否包含空格（这可以视为一种特殊格式的序列，例如已经过分词的序列）
    # 如果包含空格，则直接按空格分割序列；否则，使用分词器对象 tokenizer 对序列进行分词
    # 然后，使用词汇表 vocab 将分词后的标记转换为整数索引，并将未知标记（不在词汇表中的标记）替换为0（通常用作填充标记）
    # 最后，将编码后的序列填充或截断到指定的最大长度 max_len，并转换为 PyTorch 张量返回。
    def tokenize_and_encode(self, seq):
        if ' ' in seq:  # Modified sequence
            tokens = seq.split()
        else:  # Regular sequence
            tokens = self.tokenizer.tokenize(seq)
        
        encoded = [self.vocab.stoi.get(token, 0) for token in tokens]  # Use 0 (pad) for unknown tokens
        padded = encoded + [0] * (self.max_len - len(encoded))
        return torch.tensor(padded[:self.max_len], dtype=torch.long)


#
# 模型使用了一个嵌入层（Embedding）、一个双向门控循环单元（GRU）层和一个全连接层（Linear）
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()


# 通过计算多个统计量来评估预测值（y_pred）与实际值（y_true）之间的匹配程度
# 这些统计量包括平均绝对误差（MAE）、在特定阈值下的二分类准确度（通过精确率、召回率和F1分数衡量）、以及一个自定义的加权分数（score）
def calculate_metrics(y_true, y_pred, threshold=30):
    # 计算平均绝对误差（MAE）
    mae = np.mean(np.abs(y_true - y_pred))

    # 二分类转换：将 y_true 和 y_pred 转换为二分类标签，即低于阈值的为1（正类），否则为0（负类）
    y_true_binary = (y_true < threshold).astype(int)
    y_pred_binary = (y_pred < threshold).astype(int)

    # 在预测范围内的平均绝对误差（range_mae）
    # 创建一个掩码（mask），仅选择预测值在0到阈值之间的样本。
    mask = (y_pred >= 0) & (y_pred <= threshold)
    # 如果存在这样的样本，计算这些样本上 y_true 和 y_pred 的平均绝对误差；否则，range_mae 设置为100（一个较高的默认值，表示没有有效数据）。
    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）：
    score = (1 - mae / 100) * 0.5 + (1 - range_mae / 100) * f1 * 0.5
    return score

# 训练深度学习模型
# model: 要训练的模型。
# train_loader 和 val_loader: 分别用于加载训练数据和验证数据的 PyTorch 数据加载器。
# criterion: 损失函数，用于计算预测值和真实值之间的差异。
# optimizer: 优化器，用于更新模型的权重以最小化损失函数。
# num_epochs: 训练的总轮次。
# device: 模型和数据的存放设备，默认为 'cuda'（如果可用）。
def train_model(model, train_loader, val_loader, criterion, optimizer, num_epochs=50, device='cuda'):
    model.to(device)

    # 设置最佳分数为负无穷大，用于后续比较；最佳模型初始化为None，用于保存最佳模型的状态字典。
    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()
            print(f'New best model found with socre: {best_score:.4f}')

    return best_model

def evaluate_model(model, test_loader, device='cuda'):
    # 模型设置为评估模式
    model.eval()
    # 初始化两个列表 predictions 和 targets，用于存储模型的预测结果和真实的标签
    predictions = []
    targets = []

    # 使用 torch.no_grad() 上下文管理器来禁用梯度计算，这可以节省内存并加速计算，因为在评估模型时我们不需要计算梯度
    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())
    # 将模型的输出（预测）和真实的标签从 PyTorch 张量转换为 NumPy 数组，并分别存储在 predictions 和 targets 列表中
    y_pred = np.array(predictions)
    y_test = np.array(targets)
    
    score = calculate_metrics(y_test, y_pred)
    print(f"Test Score: {score:.4f}")



if __name__ == '__main__':
    # 代码检查CUDA（NVIDIA的GPU加速技术）是否可用，以决定是在CPU还是GPU上运行模型
    device = 'cuda' if torch.cuda.is_available() else 'cpu'

    # Load data
    # 使用pd.read_csv从CSV文件中加载训练数据
    train_data = pd.read_csv('train_data.csv')
    # 接着，通过dropna方法移除包含缺失值的行，这些缺失值出现在指定的列（包括siRNA_antisense_seq、modified_siRNA_antisense_seq_list和mRNA_remaining_pct）中
    columns = ['siRNA_antisense_seq', 'modified_siRNA_antisense_seq_list']
    train_data.dropna(subset=columns + ['mRNA_remaining_pct'], inplace=True)
    # 使用train_test_split将数据集分为训练集和验证集，验证集占10%
    train_data, val_data = train_test_split(train_data, test_size=0.1, random_state=42)

    # Create vocabulary
    # 定义了一个GenomicTokenizer用于将RNA序列转换为n-gram形式的tokens。这里使用了3-gram和步长为3的配置
    tokenizer = GenomicTokenizer(ngram=3, stride=3)
    # 遍历训练集中的相关列，将每个序列（无论是否已修改）转换为tokens，并构建词汇表vocab
    # 词汇表通过GenomicVocab.create方法创建，限制最大词汇量为10000，并过滤掉频率低于1的tokens
    all_tokens = []
    for col in columns:
        for seq in train_data[col]:
            if ' ' in seq:  # Modified sequence
                all_tokens.extend(seq.split())
            else:
                all_tokens.extend(tokenizer.tokenize(seq))
    vocab = GenomicVocab.create(all_tokens, max_vocab=10000, min_freq=1)

    # Find max sequence length
    # 计算所有训练数据中序列的最大长度max_len，这是为了确保在模型训练时，所有输入序列都有相同的长度
    # 对于包含空格的序列（即已修改的序列），直接按空格分割；对于未修改的序列，使用tokenizer进行分词
    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)
    # Create datasets
    # 使用SiRNADataset类创建训练集和验证集的数据集对象。这些对象将处理数据的批处理、编码等。
    train_dataset = SiRNADataset(train_data, columns, vocab, tokenizer, max_len)
    val_dataset = SiRNADataset(val_data, columns, vocab, tokenizer, max_len)

    # Create data loaders
    # 使用DataLoader创建数据加载器，这些加载器将在训练过程中批量提供数据，并支持打乱训练数据（对于训练加载器）。
    train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)
    val_loader = DataLoader(val_dataset, batch_size=32)

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

    optimizer = optim.Adam(model.parameters())

    # 调用train_model函数来训练模型
    train_model(model, train_loader, val_loader, criterion, optimizer, 50, device)