【Datawhale AI夏令营——生命科学赛道 Task1】

羊村的昭昭

已于 2024-07-28 15:42:08 修改

阅读量335

点赞数 16

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

于 2024-07-28 15:40:07 首次发布

本文链接：https://blog.csdn.net/m0_74457339/article/details/140747665

版权

💡PS：本次比赛旨在利用机器学习与深度学习相关技术，通过化学修饰后的siRNA序列预测RNA干扰（RNAi）机制下对靶基因的沉默效率。
mRNA疫苗在新冠预防领域取得巨大成功，推动了核酸类药物的研发。
⭐心动不如行动，一起来参加Datawhale的活动吧：https://linklearner.com/activity

✈️一、背景介绍

大致了解mRNA siRNA RNAi 定义以及任务

2.认识数据集
提示：在该网址下载数据集：http://competition.sais.com.cn/competitionDetail/532230/competitionData

表头字段	表头说明	示例
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序列 siRNA的反义序列
cell line_donor	实验使用的细胞系，一般认为对沉默效率测定存在较大影响
siRNA concentration	实验使用的siRNA浓度，一般认为对沉默效率测定存在较大影响
concentration unit	SiRNA浓度单位，比如nM	同上
Transfection method	转染方法，即将siRNA植入细胞的方法，一般认为对沉默效率测定存在较大影响
Duration after transfection h	转染后持续时间(小时)
modified siRNA sense seg	带修饰的siRNA的sense序列（P、Q和N、O一样，只是给基因段加了空格方便识别和处理，R：疾病基因的剩余量（我们所需要预测的值））

✈️二、代码实现

提示：上面都没有理解也没有关系，那我们开始正文叭😉

2.1 引入库

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

2.2 创建基因组分词器 GenomicTokenizer 类

用于将基因组序列分割成固定长度。

主要操作：
①将输入序列转化为大写
②将长序列RNA分子分割成更小的片段，有助于后续的分析和解读

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。
tokenize函数：对序列进行分词。

2.3 创建基因组词汇表 GenomicVocab类

将基因组中的特定片段与一种数据结构（索引）相关联，以便在后续的数据分析中高效地存储、检索和分析这些片段。

主要操作：
①统计每个基因片段的出现频率
②将基因片段的出现频率从高到低排序
③选择出现频率≥min_freq的基因片段，加入到词汇表
④最多保留max_vocab个，剩下的基因片段就不要了

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)

这里主要目的是创建词汇表，将标记映射到索引
create函数：根据提供的标记列表、最大词汇表大小和最小频率创建词汇表。

2.4 SiRNA数据集转换：

加载SiRNA数据，并将序列数据转换为模型可以处理的格式。
序列数据：按照一定的顺序排列的数据集合，一个重要特性是每个数据项通常都包含了与其前后相邻数据项相关的信息

文本数据处理的常见步骤：
①文本清洗：去除文本中的特殊字符、标点符号、数字等，并进行大小写统一、去停用词等处理。
②分词：将文本分割成单词或词语。
③词汇编码：构建词汇表，为每个词汇分配整数ID，并进行词嵌入。
④序列填充或截断：将文本序列填充或截断至统一长度。
⑤转换为张量：将处理后的文本序列转换为张量形式。张量：即多维数组
⑥数据加载与批处理：将张量数据加载到模型中，并进行批处理。

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

    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]
        target = torch.tensor(row['mRNA_remaining_pct'], dtype=torch.float)

        return seqs, target

    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)

__getitem__函数：获取数据集中的单个项目。
tokenize_and_encode函数：分词并编码序列。

2.5 SiRNAModel

这是一个基于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)
        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, :])  # Use last hidden state
            outputs.append(x)
        
        x = torch.cat(outputs, dim=1)
        x = self.fc(x)
        return x.squeeze()

forward函数：定义模型的前向传播。

6.评估指标计算函数

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

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

代码具体步骤解释：

平均绝对误差（mae）：计算真实值 y_true 和预测值 y_pred 之间的平均绝对误差。
二值化处理（y_true_binary、y_pred_binary）：将真实值和预测值根据给定的阈值 threshold 转换为二进制值（0 或 1）
区间内的平均绝对误差（range_mae)：满足 mask 条件的真实值和预测值计算平均绝对误差。如果满足条件的位置数量大于 0，则计算误差；否则，将 range_mae 设为 100
计算二分类的精确度(precision)
计算二分类的召回率(recall)
根据精确度和召回率计算 F1 分数(f1)
综合得分 (score)：综合考虑 MAE 和在特定范围内的 F1 分数来计算最终的得分

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

8.模型训练函数

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

前向传播：将数据输入模型，得到预测值。
计算损失：使用损失函数计算预测值与真实值之间的差异。
反向传播：根据损失函数的梯度反向传播到模型的每一层，计算每个参数的梯度。
参数更新：使用优化器根据计算出的梯度更新模型的参数

def train_model(model, train_loader, val_loader, criterion, optimizer, num_epochs=50, device='cuda'):
    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()
            print(f'New best model found with socre: {best_score:.4f}')

    return best_model

9.主函数

if __name__ == '__main__':

    device = 'cuda' if torch.cuda.is_available() else 'cpu'

    # Load data
    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)

    # Create vocabulary
    tokenizer = GenomicTokenizer(ngram=3, stride=3)

    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 = 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
    train_dataset = SiRNADataset(train_data, columns, vocab, tokenizer, max_len)
    val_dataset = SiRNADataset(val_data, columns, vocab, tokenizer, max_len)

    # Create data loaders
    train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)
    val_loader = DataLoader(val_dataset, batch_size=32)

    # Initialize model
    model = SiRNAModel(len(vocab.itos))
    criterion = nn.MSELoss()

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

    train_model(model, train_loader, val_loader, criterion, optimizer, 50, device)

✈️三、自己的思考

🥕3.1为什么使用RNN？

循环神经网络（RNN）

时序依赖：RNN 能够记住之前处理过的元素，这使得它非常适合处理时间序列数据或任何具有时序性质的数据。
梯度消失/爆炸：RNN 在训练过程中可能会遇到梯度消失或爆炸的问题，尤其是在处理长序列时。这限制了它们处理长序列的能力。
变体：为了解决这些问题，出现了一些 RNN 的变体，如 LSTM（长短期记忆网络）和 GRU（门控循环单元）

1）长短期记忆网络（LSTM）：
通过引入门控机制解决了梯度消失问题，特别擅长捕获长期依赖关系。
2）门控循环单元（GRU）：
与LSTM类似，但结构更简单，参数更少。GRU也通过门控机制来处理长距离依赖。

🥕3.1常见的网络和特点

1.多层感知器（MLP）：
一种简单的前馈神经网络，由多个层组成，每层包含多个神经元。通常用于分类和回归任务。
2. 卷积神经网络（CNN）：
专为处理具有网格结构的数据设计的网络，如图像（2D网格）或视频（3D网格）。CNN通过卷积层自动提取特征，广泛用于图像和视频识别。
3. 循环神经网络（RNN）：
能够处理序列数据的网络，具有记忆功能，可以捕获时间序列中的动态特征。适用于文本、语音识别和时间序列预测。
4. 残差网络（ResNet）：
通过引入残差学习框架来解决深层网络训练中的退化问题。ResNet通过添加跳过连接来提高训练深层网络的能力。
5. 生成对抗网络（GAN）：
由两个网络组成：生成器和判别器。生成器生成数据，判别器评估数据的真实性。两者相互竞争，提高生成数据的质量。
6. 变分自编码器（VAE）：
一种生成模型，通过学习输入数据的潜在表示来进行数据编码和解码。VAE能够生成新的、与训练数据相似的数据点。
7. ** Transformer：**
一种基于自注意力机制的网络，特别适合处理序列数据。Transformer在自然语言处理任务中表现出色，如机器翻译和文本摘要。

🥕我要怎么提高精确度？

调整学习率：降低学习；使用学习率调度器，如torch.optim.lr_scheduler中的StepLR或ExponentialLR，可以帮助模型在训练初期快速学习，然后在训练后期逐渐降低学习率以细化权重。
增加模型复杂度：我的score一直是0.6，我判断如果模型欠拟合，可以尝试增加隐藏层的维度或层数。防止过拟合的风险，因此需要与正则化技术结合使用。
尝试不同的优化器：尝试使用不同的优化器，如 AdamW，它是带有权重衰减的 Adam 优化器变体。
损失函数：如果均方误差（MSE）损失函数不适合您的任务，您可以尝试其他损失函数，如平均绝对误差（MAE）或 Huber 损失。