iEnhancer-ENCC_train_layer1.py源码阅读

cauwj / 2023-05-19 / 原文

一、基本准备¶

1、导入包¶

In [ ]:

import numpy as np

import torch 
from torch import nn
from torch.autograd import Variable
from torch.utils.data import Dataset, DataLoader
import torch.nn.functional as F

import pickle
import time
import math
import os
import csv
import glob

from sklearn import metrics
from sklearn.model_selection import KFold, StratifiedKFold

2、设置超参数¶

In [ ]:

FILE_MODEL_TMP = "model_tmp.pkl"

MY_RANDOM_STATE = 5
torch.manual_seed(MY_RANDOM_STATE)
NUMBER_EPOCHS = 20

LEARNING_RATE = 1e-4

SAMPLE_LENGTH = 200
AVGPOOL1D_KERNEL_SIZE = 4
CONV1D_KERNEL_SIZE = 3
CONV1D_FEATURE_SIZE_BLOCK1 = 32
CONV1D_FEATURE_SIZE_BLOCK2 = 64
CONV1D_FEATURE_SIZE_BLOCK3 = 128

FULLY_CONNECTED_LAYER_SIZE = 256

MODEL_DIR = 'model_layer1_seed' + str(MY_RANDOM_STATE)
if not os.path.exists(MODEL_DIR):
    os.makedirs(MODEL_DIR)

二、定义所需使用函数¶

In [ ]:

def one_hot(index, dimension):
    data = np.zeros((dimension))
    data[index] = 1
    return data

In [ ]:

def load_text_file(file_text):
    with open(file_text) as f:
        lines = f.readlines()
        my_data = [line.strip().upper() for line in lines[1::2]]
        return my_data

In [ ]:

my_dict = {'A': 0,
        'C': 1, 
        'G': 2,
        'T':3,
        'a':0,
        'c':1,
        'g':2,
        't':3}

定义了一个名为EnhancerDataset的Python类，它是PyTorch中Dataset类的子类。EnhancerDataset类通过定义，主要用于处理DNA序列数据和相关的标签数据，用于训练深度学习模型。具体实现方式如下：

构造函数__init__接受两个参数X和Y，分别表示DNA序列数据和相关的标签数据
__getitem__：
- 从X中取出一个样本（DNA序列），同时从Y中取出相应的标签
- 将DNA序列转换为一个one-hot编码的矩阵，其中每一列对应于DNA序列中的一个字符，每一行对应于字符的one-hot编码的矩阵，其中每一列对应于DNA序列中的一个字符，每一行对应于字符的one-hot编码
- 从DNA序列中提取出1-mer、2-mer和3-mer，分别对应于单个字符、相邻两个字符和相邻三个字符的频率
- 最后，它将one-hot编码矩阵和提取出的1-mer、2-mer和3-mer矩阵进行拼接，得到一个维度为(7, 200)的特征矩阵。这个特征矩阵被转换为PyTorch张量，并和标签一起返回
三个成员函数extract_1_mer、extract_2_mer和extract_3_mer，用于从DNA序列中提取1-mer、2-mer和3-mer
- extract_1_mer函数统计DNA序列中每种字符（A、C、G、T）的出现次数，然后将出现次数除以序列长度得到该字符的频率
- extract_2_mer函数统计DNA序列中每种相邻两个字符（AA、AC、AG、AT等）的出现次数，然后将出现次数除以序列长度减1得到相邻两个字符的频率
- extract_3_mer函数统计DNA序列中每种相邻三个字符（AAA、AAC、AAG、AAT等）的出现次数，然后将出现次数除以序列长度得到相邻三个字符的频率
__len__：用于返回数据集的大小，即DNA序列的个数

In [ ]:

class EnhancerDataset(Dataset):
    # X: list Enhancer sequence (每个序列200个字符)
    # Y: list label [0, 1]; 0: 负样本, 1: 正样本
    def __init__(self, X, Y):
        self.X = X
        self.Y = Y

    def __getitem__(self, index):
        label = self.Y[index]
        sample = self.X[index]
        
        values = np.zeros((4, SAMPLE_LENGTH))
        for i in range(SAMPLE_LENGTH):
            char_idx = my_dict[sample[i]]
            values[char_idx, i] = 1 
        
        values_one_mer = self.extract_1_mer(sample)
        #values = np.concatenate((values, values_one_mer), axis=0)
        values_two_mer = self.extract_2_mer(sample)
        #values = np.concatenate((values, values_two_mer), axis=0)
        values_three_mer = self.extract_3_mer(sample)
        #values = np.concatenate((values, values_three_mer), axis=0)
        values = np.concatenate((values, values_one_mer, values_two_mer, 
                        values_three_mer), axis=0)
        
        input = torch.from_numpy(values)
        return input, label
    
    def extract_1_mer(self, sample):
        my_count = {'A': 0.0, 'C': 0.0, 'G': 0.0, 'T': 0.0}        
        values = np.zeros((1, SAMPLE_LENGTH))
        for i in range(SAMPLE_LENGTH):
            my_count[sample[i]] += 1
        
        #for one_mer in my_count:
        #    print("one mer: ", one_mer, " : ", my_count[one_mer])
        
        for i in range(SAMPLE_LENGTH):
            values[0, i] = my_count[sample[i]] / SAMPLE_LENGTH;
        
        #print("values: ", values)    
        return values
    
    def extract_2_mer(self, sample):
        my_count = {'AA': 0.0, 'AC': 0.0, 'AG': 0.0, 'AT': 0.0,
                    'CA': 0.0, 'CC': 0.0, 'CG': 0.0, 'CT': 0.0,
                    'GA': 0.0, 'GC': 0.0, 'GG': 0.0, 'GT': 0.0,
                    'TA': 0.0, 'TC': 0.0, 'TG': 0.0, 'TT': 0.0} 
        values = np.zeros((2, SAMPLE_LENGTH))
        for i in range(SAMPLE_LENGTH - 1):
            two_mer = sample[i:i+2]
            #print("two_mer: ", two_mer)
            my_count[two_mer] += 1
        
        #for two_mer in my_count:
        #    print("two mer: ", two_mer, " : ", my_count[two_mer])
        
        values = np.zeros((2, SAMPLE_LENGTH))
        for i in range(1,SAMPLE_LENGTH-1):
            two_mer_left = sample[i-1:i+1]
            two_mer_right = sample[i:i+2]
            
            values[0, i] = my_count[two_mer_left] / (SAMPLE_LENGTH - 1);
            values[1, i] = my_count[two_mer_right] / (SAMPLE_LENGTH - 1);
        
        #print("values: ", values) 
        return values
    
    def extract_3_mer(self, sample):
        my_count = {}
                                        
        for firchCh in ['A', 'C', 'G', 'T']:
            for secondCh in ['A', 'C', 'G', 'T']:
                for thirdCh in ['A', 'C', 'G', 'T']:
                    three_mer = firchCh + secondCh + thirdCh
                    my_count[three_mer] = 0.0
        for i in range(SAMPLE_LENGTH - 2):
            three_mer = sample[i:i+3]
            #print("two_mer: ", two_mer)
            my_count[three_mer] += 1
                    
        values = np.zeros((1, SAMPLE_LENGTH))
        for i in range(1,SAMPLE_LENGTH-2):
            three_mer = sample[i-1:i+2]
            values[0, i] = my_count[three_mer] / SAMPLE_LENGTH;
                    
        return values
        
    def __len__(self):
        #return 100
        return len(self.X)

三、搭建CNN模型¶

该模型为一个卷积神经网络，通过多层卷积和池化操作，最终将特征映射转换为一个实数，并经过 Sigmoid 函数输出二元分类结果。其中，CONV1D_FEATURE_SIZE_BLOCK1、CONV1D_FEATURE_SIZE_BLOCK2、CONV1D_KERNEL_SIZE、AVGPOOL1D_KERNEL_SIZE、FULLY_CONNECTED_LAYER_SIZE、LEARNING_RATE 等超参数已定义和赋值

In [ ]:

# import torch.nn as nn
# import torch.nn.functional as F
# import torch

# 定义模型类
class EnhancerCnnModel(nn.Module):
    def __init__(self):
        super(EnhancerCnnModel, self).__init__()
        
        # 第一种卷积+第一种BN层
        self.c1_1 = nn.Conv1d(8, CONV1D_FEATURE_SIZE_BLOCK1, CONV1D_KERNEL_SIZE, padding=1)
        self.c1_1bn = nn.BatchNorm1d(CONV1D_FEATURE_SIZE_BLOCK1)
        
        # 第一种卷积+第一种BN层
        self.c1_2 = nn.Conv1d(CONV1D_FEATURE_SIZE_BLOCK1, CONV1D_FEATURE_SIZE_BLOCK1, 
            CONV1D_KERNEL_SIZE, padding=1)
        self.c1_2bn = nn.BatchNorm1d(CONV1D_FEATURE_SIZE_BLOCK1)
        
        # 第一种卷积+第一种BN层
        self.c1_3 = nn.Conv1d(CONV1D_FEATURE_SIZE_BLOCK1, CONV1D_FEATURE_SIZE_BLOCK1, 
            CONV1D_KERNEL_SIZE, padding=1)    
        self.c1_3bn = nn.BatchNorm1d(CONV1D_FEATURE_SIZE_BLOCK1)
        
        # 池化层
        self.p1 = nn.MaxPool1d(AVGPOOL1D_KERNEL_SIZE)
        
        # 第二种卷积+第二种BN层
        self.c2_1 = nn.Conv1d(CONV1D_FEATURE_SIZE_BLOCK1, 
            CONV1D_FEATURE_SIZE_BLOCK2, 
            CONV1D_KERNEL_SIZE, padding=1)
        self.c2_1bn = nn.BatchNorm1d(CONV1D_FEATURE_SIZE_BLOCK2)
        
        # 第二种卷积+第二种BN层
        self.c2_2 = nn.Conv1d(CONV1D_FEATURE_SIZE_BLOCK2, CONV1D_FEATURE_SIZE_BLOCK2, 
            CONV1D_KERNEL_SIZE, padding=1)
        self.c2_2bn = nn.BatchNorm1d(CONV1D_FEATURE_SIZE_BLOCK2)
        
        # 第二种卷积+第二种BN层
        self.c2_3 = nn.Conv1d(CONV1D_FEATURE_SIZE_BLOCK2, CONV1D_FEATURE_SIZE_BLOCK2, 
            CONV1D_KERNEL_SIZE, padding=1)
        self.c2_3bn = nn.BatchNorm1d(CONV1D_FEATURE_SIZE_BLOCK2)
        
        # 池化层
        self.p2 = nn.MaxPool1d(AVGPOOL1D_KERNEL_SIZE)
        
        # 全连接层
        self.fc = nn.Linear(768, FULLY_CONNECTED_LAYER_SIZE)        
        # 输出层
        self.out = nn.Linear(FULLY_CONNECTED_LAYER_SIZE, 1)
        
        # 定义损失函数
        self.criterion = nn.BCELoss()
        # 定义优化器
        self.optimizer = torch.optim.Adam(self.parameters(), lr=LEARNING_RATE)
    
    # 前向传播函数
    def forward(self, inputs):
        batch_size = inputs.size(0)
        # 将输入的形状从(batch_size x seq_len) 转换为 (batch_size x input_size x seq_len)，以适应卷积层的输入
        #inputs = inputs.transpose(1,2)
        #print("inputs size: ", inputs.size())        
        output = F.relu(self.c1_1bn(self.c1_1(inputs)))
        output = F.relu(self.c1_2bn(self.c1_2(output)))
        output = F.relu(self.c1_3bn(self.c1_3(output)))
        output = self.p1(output)
        #print("After p1: ", output.shape) 
        
        output = F.relu(self.c2_1bn(self.c2_1(output)))
        output = F.relu(self.c2_2bn(self.c2_2(output)))
        output = F.relu(self.c2_3bn(self.c2_3(output)))
        output = self.p2(output)
        #print("After p2: ", output.shape)
        
        output = output.view(batch_size, -1)
        #print("Reshape : ", output.shape)
        
        output = F.relu(self.fc(output))
        #print("After FC layer: ", output.shape)  
        
        output = torch.sigmoid(self.out(output))
        #print("Final output (After sigmoid): ", output.shape)
        #print("Final output: ", output)
        
        return output 

四、构建训练函数¶

下面函数用于训练一个 epoch，即遍历一遍训练集中的所有数据并更新模型参数。其中，该函数接受三个参数：model 表示要训练的模型；train_loader 表示训练集数据的 DataLoader，用于加载数据；learning_rate 表示当前训练的学习率。

在函数开始时，将模型设置为训练模式，并更新学习率。然后遍历训练集中的每一批数据，获取输入和标签，并将其转换为 float 类型。如果 GPU 可用，将输入和标签转换为 Variable 对象并移动到 GPU 上。然后清除参数的梯度，进行前向传播，计算损失函数，反向传播，更新参数，并累加损失。最后计算训练集的平均损失并返回。

In [ ]:

# 训练一个 epoch 的函数
def train_one_epoch(model, train_loader, learning_rate):
    # 设置模型为训练模式
    model.train()
    # 更新学习率
    for param_group in model.optimizer.param_groups:
        param_group['lr'] = learning_rate
    
    epoch_loss_train = 0.0
    nb_train = 0
    # 遍历训练集中的每一批数据
    for i, data in enumerate(train_loader, 0):
        # 获取输入和标签
        inputs, labels = data   
        #print("\ninputs: ", inputs.shape)
        #print("labels: ", labels.shape)
        
        inputs_length = inputs.size()[0]
        nb_train += inputs_length
        
        inputs = inputs.float()
        labels = labels.float()
        
        #######################
        #  使用 GPU 进行计算  #
        #######################
        if torch.cuda.is_available():
            inputs = Variable(inputs.cuda())
            labels = Variable(labels.cuda())            
        else:
            inputs = Variable(inputs)
            labels = Variable(labels)
        
        # 清除参数的梯度
        model.optimizer.zero_grad()
        # 进行前向传播
        outputs = model(inputs)
        
        # 计算损失函数
        loss = model.criterion(outputs, labels)
        #print("loss: ", loss)
        # 反向传播
        loss.backward()
        # 更新参数
        model.optimizer.step()
        
        # 累加损失
        epoch_loss_train = epoch_loss_train + loss.item() * inputs_length
    
    print("nb_train: ", nb_train)
    # 计算平均损失
    epoch_loss_train_avg = epoch_loss_train / nb_train
    return epoch_loss_train_avg

五、构建评估函数¶

评估给定神经网络模型在给定数据集上的性能表现

定义函数 evaluate(file_model, loader)，它有两个输入参数：file_model 表示预训练的模型文件路径，loader 表示数据集加载器
如果 CUDA 可用，则将模型 model 移动到 GPU 上，以便在 GPU 上进行计算
加载预训练的模型权重，将模型 model 的参数设置为预训练的参数
设置模型 model 为评估模式，即关闭 Dropout 和 Batch Normalization 等随机性操作
定义变量 epoch_loss 和 nb_samples，分别用于累加每个 batch 的损失和样本数
定义空列表 arr_labels、arr_labels_hyp 和 arr_prob，分别用于存储真实标签、预测标签和模型的预测概率
遍历数据集 loader，对于每个 batch：
- 获取输入数据 inputs 和标签 labels
- 计算 inputs 的长度 inputs_length，并将其累加到 nb_samples 中
- 将标签 labels 的维度从 (batch_size, 1) 转换为 (batch_size,)，并将其转换为 numpy 数组后添加到 arr_labels 列表中
- 将输入数据 inputs 和标签 labels 转换为 float 类型，以便于后续计算
- 如果 CUDA 可用，则将输入数据 inputs 和标签 labels 移动到 GPU 上，以便于在 GPU 上进行计算
- 将输入数据 inputs 作为输入，通过模型 model 进行前向传递，得到输出 outputs
- 计算模型的损失 loss，并将其累加到变量 epoch_loss 中
- 将模型的预测概率从 tensor 转换为 numpy 数组，并将其添加到 arr_prob 列表中
打印样本总数 nb_samples 和平均损失 epoch_loss_avg
计算 AUC 值，并打印输出
将模型的预测概率转换为二分类标签，并将其添加到 arr_labels_hyp 列表中
将真实标签 arr_labels 和预测标签 arr_labels_hyp 转换为整型，并分别存储到变量 arr_labels 和 arr_labels_hyp 中
调用函数 calculate_confusion_matrix(arr_labels, arr_labels_hyp)，计算准确率、混淆矩阵、灵敏度、特异度和 MCC 值，并将它们存储到变量 acc、confusion_matrix、sensitivity、specificity 和 mcc 中

函数在后面进行了定义

将评估指标存储到字典 result 中，并返回该字典
打印输出准确率 acc

In [ ]:

def evaluate(file_model, loader):
    # 定义一个 EnhancerCnnModel 的实例 model
    model = EnhancerCnnModel()
    
    # 如果 CUDA 可用，则将 model 移动到 GPU 上
    if torch.cuda.is_available():
        model.cuda()
    
    # 加载预训练的模型权重
    model.load_state_dict(torch.load(file_model))
    
    # 设置模型为评估模式
    model.eval()    
    
    epoch_loss = 0.0
    nb_samples = 0
    
    arr_labels = []
    arr_labels_hyp = []
    arr_prob = []
    
    # 遍历数据集
    for i, data in enumerate(loader, 0):
        # 获取输入和标签
        inputs, labels = data
        
        # inputs 的长度
        inputs_length = inputs.size()[0]
        nb_samples += inputs_length
        
        # 将标签添加到 arr_labels 列表中
        arr_labels += labels.squeeze(1).data.cpu().numpy().tolist()

        # 将 inputs 和 labels 转换为 float 类型
        inputs = inputs.float()
        labels = labels.float()
        
        # 如果 CUDA 可用，则将 inputs 和 labels 移动到 GPU 上
        if torch.cuda.is_available():
            inputs = Variable(inputs.cuda())
            labels = Variable(labels.cuda())
        else:
            inputs = Variable(inputs)
            labels = Variable(labels)
        
        # 前向传递
        outputs = model(inputs)
        
        # 计算损失使用模型 model 的损失函数来计算该模型在给定输出 outputs 和
        # 标签 labels 下的损失
        loss = model.criterion(outputs, labels)
        
        # 累加 epoch_loss
        epoch_loss = epoch_loss + loss.item() * inputs_length
        
        # 将模型的预测结果添加到 arr_prob 列表中
        arr_prob += outputs.squeeze(1).data.cpu().numpy().tolist()
    
    # 输出样本总数和平均损失
    print("nb_samples: ", nb_samples)
    epoch_loss_avg = epoch_loss / nb_samples
    print("epoch loss avg: ", epoch_loss_avg)
        
    # 输出预测概率和真实标签的列表长度
    print("arr_prob: ", len(arr_prob))
    print("arr_labels: ", len(arr_labels))
    
    # 计算 AUC 值
    auc = metrics.roc_auc_score(arr_labels, arr_prob)
    print("auc: ", auc)
    
    # 将预测概率转换为二分类标签，阈值为 0.5
    arr_labels_hyp = [int(prob > 0.5) for prob in arr_prob]
    arr_labels = [int(label) for label in arr_labels]
    
    # 计算各种分类指标
    acc, confusion_matrix, sensitivity, specificity, mcc = calculate_confusion_matrix(arr_labels, arr_labels_hyp)
    
    # 将结果存储到字典中
    result = {'epoch_loss_avg': epoch_loss_avg,
              'acc' : acc,
              'confusion_matrix' : confusion_matrix,
              'sensitivity' : sensitivity,
              'specificity' : specificity,
              'mcc' : mcc,
              'auc' : auc,
              'arr_prob': arr_prob,
              'arr_labels': arr_labels,
              'arr_labels_hyp':arr_labels_hyp
              }
    
    # 输出 accuracy
    print("acc: ", acc)
    
    # 返回结果字典
    return result

定义函数 calculate_confusion_matrix(arr_labels, arr_labels_hyp)¶

函数的作用是计算混淆矩阵和相关的分类指标，用于评估模型在给定数据集上的性能表现。函数的输入参数为 arr_labels 和 arr_labels_hyp，分别表示真实标签的列表和模型的预测标签的列表

混淆矩阵（confusion matrix）是一种常见的二分类模型性能评估工具，用于衡量模型的分类能力和错误分类情况。混淆矩阵的行表示真实标签，列表示模型预测的标签，每个元素表示对应的样本数。在二分类问题中，混淆矩阵通常被表示为下面的形式：

	真实情况
	正例	反例
预测情况
正例	TP	FP
反例	FN	TN

通过混淆矩阵，我们可以计算出许多有用的评价指标，例如准确率、精确率、召回率、F1 值、特异度、灵敏度、MCC 值等

In [ ]:

def calculate_confusion_matrix(arr_labels, arr_labels_hyp):
    # 初始化变量
    corrects = 0  # 正确分类的样本数
    confusion_matrix = np.zeros((2, 2))  # 混淆矩阵

    # 遍历所有样本
    for i in range(len(arr_labels)):
        # 更新混淆矩阵
        confusion_matrix[arr_labels_hyp[i]][arr_labels[i]] += 1

        # 统计正确分类的样本数
        if arr_labels[i] == arr_labels_hyp[i]:
            corrects = corrects + 1

    # 计算分类指标
    acc = corrects * 1.0 / len(arr_labels)  # 准确率
    specificity = confusion_matrix[0][0] / (confusion_matrix[0][0] + confusion_matrix[1][0])  # 特异度
    sensitivity = confusion_matrix[1][1] / (confusion_matrix[0][1] + confusion_matrix[1][1])  # 灵敏度
    tp = confusion_matrix[1][1]  # 真正例数
    tn = confusion_matrix[0][0]  # 真反例数
    fp = confusion_matrix[1][0]  # 假正例数
    fn = confusion_matrix[0][1]  # 假反例数
    mcc = (tp * tn - fp * fn ) / math.sqrt((tp + fp)*(tp + fn)*(tn + fp)*(tn + fn))  # MCC值

    # 返回分类指标
    return acc, confusion_matrix, sensitivity, specificity, mcc

六、构建K折的卷积神经网络训练评估模型¶

train_one_fold函数的输入参数为 dataset 和 fold，分别表示数据集和当前的交叉验证折数

train_set：训练集数据和标签，通过 EnhancerDataset 类从 data_train 和 label_train 中获取
val_set：验证集数据和标签，通过 EnhancerDataset 类从 data_val 和 label_val 中获取
train_loader：训练集数据加载器，通过 DataLoader 类从 train_set 中加载，每次加载32个样本，打乱顺序，使用4个 worker 进程进行数据加载
val_loader：验证集数据加载器，通过 DataLoader 类从 val_set 中加载，每次加载32个样本，不打乱顺序，使用4个 worker 进程进行数据加载
model：卷积神经网络模型，通过调用 EnhancerCnnModel 类创建，该类继承自 nn.Module
best_val_loss：最佳验证集损失，初始化为 1000
best_epoch：最佳验证集损失对应的训练轮数，初始化为 0
file_model：最佳模型的保存路径和文件名
train_loss：训练集损失列表，用于保存每个训练轮数的训练集损失
val_loss：验证集损失列表，用于保存每个训练轮数的验证集损失
for epoch in range(NUMBER_EPOCHS)：对于每个训练轮数，进行模型训练和验证
- train_one_epoch(model, train_loader, LEARNING_RATE)：训练一个训练轮数，并返回训练集损失的平均值
- torch.save(model.state_dict(), os.path.join(MODEL_DIR, FILE_MODEL_TMP))：保存当前模型的权重参数到 FILE_MODEL_TMP 文件中。
- evaluate(file_model_tmp, val_loader)：使用当前模型在验证集上进行评估，并返回验证集损失等指标
- if best_val_loss > result_val['epoch_loss_avg']：如果当前验证集损失比历史最佳验证集损失更小，则更新历史最佳验证集损失和最佳模型的权重参数，并打印保存模型和最佳验证集损失
- train_loss.append(epoch_loss_train_avg)：将当前训练轮数的训练集损失添加到 train_loss 列表中
- val_loss.append(result_val['epoch_loss_avg'])：将当前训练轮数的验证集损失添加到 val_loss 列表中。

In [ ]:

def train_one_fold(dataset, fold):
    # 创建训练集和验证集的数据集对象
    train_set = EnhancerDataset(dataset["data_train"], dataset["label_train"])
    val_set = EnhancerDataset(dataset["data_val"], dataset["label_val"])    
    
    # 创建训练集和验证集的数据加载器
    train_loader = DataLoader(dataset=train_set, batch_size=32,                              
                              shuffle=True, num_workers=4)
    val_loader = DataLoader(dataset=val_set, batch_size=32,                              
                              shuffle=False, num_workers=4)
    
    # 创建增强型卷积神经网络模型对象
    model = EnhancerCnnModel()
    print("CNN Model: ", model)
    if torch.cuda.is_available():
        model.cuda()
    
    best_val_loss = 1000
    best_epoch = 0 
    
    file_model = "best_model_saved.pkl"
    
    train_loss = []
    val_loss = []
             
    # 遍历每个训练轮次（epoch）
    for epoch in range(NUMBER_EPOCHS):
        print("\n############### EPOCH : ", epoch, " ###############")
        # 在当前轮次中训练模型并计算训练损失
        epoch_loss_train_avg = train_one_epoch(model, train_loader, LEARNING_RATE)
        print("epoch_loss_train_avg: ", epoch_loss_train_avg)
        
        # 保存当前模型参数到临时文件
        torch.save(model.state_dict(), os.path.join(MODEL_DIR, FILE_MODEL_TMP))
        
        # 使用验证集评估当前模型，并得到验证集上的损失和其他指标
        file_model_tmp = os.path.join(MODEL_DIR, FILE_MODEL_TMP)
        result_val = evaluate(file_model_tmp, val_loader)    
        
        # 如果当前轮次的验证集损失更低，则保存当前模型为最佳模型
        if best_val_loss > result_val['epoch_loss_avg']:
            torch.save(model.state_dict(), os.path.join(MODEL_DIR, file_model))
            best_val_loss = result_val['epoch_loss_avg']
            print("Save model, best_val_loss: ", best_val_loss)
            best_epoch = epoch 
        
        # 记录训练集和验证集的损失值
        train_loss.append(epoch_loss_train_avg)
        val_loss.append(result_val['epoch_loss_avg'])
            
    # 加载最佳模型参数
    model.load_state_dict(torch.load(MODEL_DIR+"/best_model_saved.pkl"))
    model_fn = "enhancer_{}_{}.pkl".format(fold, best_epoch)
    torch.save(model.state_dict(), os.path.join(MODEL_DIR, model_fn))
    
    # 保存训练集和验证集的损失到日志文件
    with open(MODEL_DIR+"/logfile_loss_model_{}_{}.csv".format(fold, best_epoch), mode='w') as lf_loss:
        lf_loss = csv.writer(lf_loss, delimiter=',')
        lf_loss.writerow(['epoch', 'train loss', 'validation loss'])
        for i in range(np.size(train_loss)):
            lf_loss.writerow([i, train_loss[i], val_loss[i]])
    
    print("\n####### EVALUATE ON VALIDATION")
    print("best_epoch: ", best_epoch)
    
    print("\n ==> VALIDATION RESULT: ") 
    file_model = os.path.join(MODEL_DIR, model_fn) 
    result_val = evaluate(file_model, val_loader)

train_kfold()的主要作用是使用 KFold 进行 k 折交叉验证，并对每一折数据进行训练和验证。具体过程如下：

加载数据集，包括强增强子序列和非增强子序列
将强增强子序列和非增强子序列的标签分别赋值为 1 和 0，并将数据和标签合并
使用 KFold 进行 k 折交叉验证

观察源码不难发现有两种k折交叉验证的方法，它们的不同之处在于使用的采样方法的不同

KFold：K 折交叉验证，将数据集分为 k 个互不重叠的折，每次使用其中一折作为验证集，其余折作为训练集。

StratifiedKFold：分层 K 折交叉验证，与 KFold 类似，但是在分折时会考虑标签的分布情况，保证每个折中各个类别的比例与原始数据集中的比例相同。

因此，如果数据集中各个类别的样本数量比例不同，或者某个类别的样本数量较少，建议使用 StratifiedKFold 进行交叉验证，以确保每个折中的样本比例与原始数据集中的样本比例相同。如果数据集中各个类别的样本数量比例相同，可以使用 KFold 进行交叉验证。

对每一折数据进行训练和验证：
- 根据索引得到训练集和验证集
- 将训练集和验证集存入字典中
- 检查数据集
- 对当前折数据进行训练

In [ ]:

def train_kfold():
    # 加载数据集，包括强增强子序列和非增强子序列
    print("\n ==> Loading train set")
    data_strong = load_text_file('intersect_CD.fa')
    print("data_strong: ", len(data_strong))
    # data_weak = load_text_file('data_week_enhancers.txt')
    # print("data_weak: ", len(data_weak))
    
    data_enhancers = data_strong # + data_weak
    print("data_enhancers: ", len(data_enhancers))
    data_non_enhancers = load_text_file('noenhancer_seq_CD.fa')
    print("data_non_enhancers: ", len(data_non_enhancers))    
    
    # 给增强子和非增强子标签分别赋值为 1 和 0
    label_enhancers = np.ones((len(data_enhancers),1))
    label_non_enhancers = np.zeros((len(data_non_enhancers), 1))
    
    # 将数据和标签合并
    data = np.concatenate((data_enhancers, data_non_enhancers))
    label = np.concatenate((label_enhancers, label_non_enhancers))
    
    # 使用 KFold 进行 k 折交叉验证
    kf = KFold(n_splits=5, shuffle=True, random_state=MY_RANDOM_STATE)
    
    #kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=MY_RANDOM_STATE)
    
    fold = 0
    # 对每一折数据进行训练和验证
    for train_index, val_index in kf.split(data,label):
        # 根据索引得到训练集和验证集
        data_train, data_val = data[train_index], data[val_index]
        label_train, label_val = label[train_index], label[val_index]
        
        # 将训练集和验证集存入字典中
        dataset = {"data_train": data_train, "label_train": label_train,  
            "data_val": data_val, "label_val": label_val }
        
        # 检查数据集
        #check_dataset(dataset)
        # 对当前折数据进行训练
        train_one_fold(dataset, fold)
        fold += 1
        #break

七、检查数据集中的样本¶

检查数据集中的每个样本是否符合要求，即是否由长度为 200 的序列组成

具体过程如下：

从字典 dataset 中获取训练集 data_train 和验证集 data_val
对训练集和验证集中的每个样本进行检查，如果样本长度不为 200，则将错误样本数量加 1
输出检查结果，如果训练集或验证集中有错误样本，则输出错误样本数量；否则输出“OK”

这个函数的主要作用是确保数据集中的每个样本都符合要求，以避免在训练和验证过程中出现错误。如果数据集中的样本长度不为 200，则可能会导致模型无法处理这些样本或者产生错误的预测结果。因此，在进行训练和验证之前，需要先对数据集进行检查，以确保数据集中的每个样本都符合要求。

In [ ]:

def check_dataset(dataset):
    print("\n==> Checking dataset")
    # 检查 data_train
    data_train = dataset["data_train"]
    nb_error_samples = 0
    for sample in data_train:
        # 检查样本长度是否为 200
        if len(sample) != 200:
            nb_error_samples += 1
    if nb_error_samples > 0:
        print("data_train 出现错误样本数量: ", nb_error_samples)
    else: 
        print("data_train：无错误！")

    # 检查 data_val
    data_val = dataset["data_val"]
    nb_error_samples = 0
    for sample in data_val:
        # 检查样本长度是否为 200
        if len(sample) != 200:
            nb_error_samples += 1
    if nb_error_samples > 0:
        print("data_val 出现错误样本数量: ", nb_error_samples)
    else:
        print("data_val：无错误！")

In [ ]:

# if __name__== "__main__":
#     traininglog_fn = MODEL_DIR + "/training_log.csv"    
    
#     train_kfold()   

# 这段代码的作用是运行整个程序。具体过程如下：
# 定义变量 traininglog_fn，用于存储训练日志文件路径。
# 调用函数 train_kfold() 进行 k 折交叉验证，对每一折数据进行训练和验证。
# 训练过程中会将训练日志写入到文件中，存储在 traininglog_fn 中。
# 程序运行结束后，会输出训练和验证的结果，包括准确率、精确率、召回率、F1 值等评价指标。