一、基础配置与参数设定

1. 导入必要的库 本节代码引入了深度学习模型构建、数据处理、文件操作以及性能评估所需的Python库，为DNA序列分类任务提供全面的支持。

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. 定义模型超参数 以下参数配置了模型的结构、训练过程以及数据处理的细节，它们是影响模型性能和训练稳定性的关键因素。

TEMP_MODEL_SAVE_PATH = "temporary_model.pth"
GLOBAL_RANDOM_SEED = 5
torch.manual_seed(GLOBAL_RANDOM_SEED) # 设置PyTorch的随机种子
TRAINING_EPOCHS = 20
MODEL_LEARNING_RATE = 1e-4
SEQUENCE_LENGTH = 200 # DNA序列的统一长度
POOLING_KERNEL_SIZE = 4 # 最大池化层的核大小
CONVOLUTION_KERNEL_SIZE = 3 # 卷积层的核大小
CONV_FILTERS_BLOCK1 = 32 # 第一个卷积块的输出通道数
CONV_FILTERS_BLOCK2 = 64 # 第二个卷积块的输出通道数
# CONV_FILTERS_BLOCK3 = 128 # 原始代码中定义但未在模型中使用
DENSE_LAYER_DIM = 256 # 全连接层的神经元数量

OUTPUT_MODEL_DIR = 'model_output_seed' + str(GLOBAL_RANDOM_SEED)
if not os.path.exists(OUTPUT_MODEL_DIR):
    os.makedirs(OUTPUT_MODEL_DIR) # 创建模型保存目录

二、数据处理辅助函数与数据集类 本节定义了处理DNA序列数据和K-mer特征提取的辅助函数及自定义数据集类。这些工具旨在将原始的生物序列信息转化为深度学习模型可理解的数值特征。

def get_one_hot_encoding(idx, dimension):
    """
    生成指定维度下的独热编码向量。
    :param idx: 需要进行独热编码的索引。
    :param dimension: 独热编码向量的总维度。
    :return: 独热编码的NumPy数组。
    """
    encoding_vec = np.zeros((dimension))
    encoding_vec[idx] = 1
    return encoding_vec

def read_sequence_file(filepath):
    """
    从文本文件中加载DNA序列数据。
    文件格式假定为FASTA样式的，序列在奇数行（从1开始计数）。
    :param filepath: 文本文件路径。
    :return: 包含所有DNA序列的字符串列表。
    """
    with open(filepath) as f:
        lines = f.readlines()
        sequences = [line.strip().upper() for line in lines[1::2]] # 跳过FASTA头信息，读取序列行
        return sequences

NUCLEOTIDE_MAP = {'A': 0, 'C': 1, 'G': 2, 'T': 3,
                  'a': 0, 'c': 1, 'g': 2, 't': 3} # 核苷酸到整数的映射

NucleotideSequenceDataset 类是 PyTorch 的 Dataset 子类，专门用于处理DNA序列数据及其对应的标签。其核心功能是将原始DNA序列转换为模型可接受的数值特征，包括独热编码和多种K-mer频率特征。

• __init__: 构造函数，接收DNA序列列表 (X_data) 和标签列表 (Y_labels)。
• __getitem__: 这是PyTorch Dataset 的核心方法，用于按索引获取单个样本。
• 从数据集中获取单个DNA序列及其对应的标签。
• 将DNA序列进行独热编码，生成一个维度为 (4, SEQUENCE_LENGTH) 的特征矩阵。
• 进一步提取1-mer、2-mer和3-mer的全局频率，并以特定方式将其映射到与序列长度相关的特征通道中。
• 将独热编码特征与所有K-mer频率特征在通道维度上拼接，形成一个总维度为 (8, SEQUENCE_LENGTH) 的复合特征矩阵。
• 将该复合特征矩阵转换为PyTorch张量并返回，同时返回对应的标签。

extract_monomer_frequencies: 计算序列中每个核苷酸（A, C, G, T）的全局频率，并将其"广播"到 (1, SEQUENCE_LENGTH) 的特征矩阵中，作为位置相关的特征。 extract_dimer_frequencies: 计算所有16种二核苷酸的全局频率。然后，通过两种不同的映射方式（基于二核苷酸的结束位置和开始位置），将其填充到 (2, SEQUENCE_LENGTH) 的特征矩阵中。 extract_trimer_frequencies: 计算所有64种三核苷酸的全局频率，并映射到 (1, SEQUENCE_LENGTH) 的特征矩阵中。 __len__: 返回数据集中样本的总数。

class NucleotideSequenceDataset(Dataset):
    """
    用于增强子预测的DNA序列数据集。
    输入X_data为DNA序列字符串列表，Y_labels为对应的二分类标签列表。
    """
    def __init__(self, X_data, Y_labels):
        self.sequences = X_data
        self.labels = Y_labels

    def __getitem__(self, index):
        current_label = self.labels[index]
        current_sequence = self.sequences[index]

        # 1. 独热编码特征 (4个通道)
        one_hot_features = np.zeros((4, SEQUENCE_LENGTH))
        for i in range(SEQUENCE_LENGTH):
            char_index = NUCLEOTIDE_MAP[current_sequence[i]]
            one_hot_features[char_index, i] = 1

        # 2. 提取并拼接k-mer频率特征
        monomer_freq_feat = self.extract_monomer_frequencies(current_sequence) # (1, SEQUENCE_LENGTH)
        dimer_freq_feat = self.extract_dimer_frequencies(current_sequence)     # (2, SEQUENCE_LENGTH)
        trimer_freq_feat = self.extract_trimer_frequencies(current_sequence)   # (1, SEQUENCE_LENGTH)

        # 最终特征矩阵维度为 (4 (独热) + 1 (1-mer) + 2 (2-mer) + 1 (3-mer), SEQUENCE_LENGTH) = (8, SEQUENCE_LENGTH)
        combined_features = np.concatenate((one_hot_features,
                                            monomer_freq_feat,
                                            dimer_freq_feat,
                                            trimer_freq_feat), axis=0)

        input_tensor = torch.from_numpy(combined_features).float() # 确保特征张量为浮点型
        return input_tensor, float(current_label) # 确保标签为浮点型

    def extract_monomer_frequencies(self, sequence):
        """
        计算单核苷酸的全局频率，并生成对应的特征矩阵。
        该特征矩阵的每个位置存储了该位置核苷酸在整个序列中的全局频率。
        """
        monomer_counts = {'A': 0.0, 'C': 0.0, 'G': 0.0, 'T': 0.0}
        for char in sequence:
            monomer_counts[char] += 1

        feature_matrix = np.zeros((1, SEQUENCE_LENGTH))
        for i in range(SEQUENCE_LENGTH):
            # 将对应位置核苷酸的全局频率作为该位置的特征
            feature_matrix[0, i] = monomer_counts[sequence[i]] / SEQUENCE_LENGTH
        return feature_matrix

    def extract_dimer_frequencies(self, sequence):
        """
        计算二核苷酸的全局频率，并生成对应的特征矩阵（两种表示）。
        feature_matrix[0, i] 存储以 i 为结束位置的二核苷酸的全局频率。
        feature_matrix[1, i] 存储以 i 为开始位置的二核苷酸的全局频率。
        """
        dimer_counts = {nt1 + nt2: 0.0 for nt1 in 'ACGT' for nt2 in 'ACGT'}
        for i in range(SEQUENCE_LENGTH - 1): # 计算序列中所有二核苷酸的全局计数
            dimer = sequence[i:i+2]
            dimer_counts[dimer] += 1

        feature_matrix = np.zeros((2, SEQUENCE_LENGTH))
        for i in range(1, SEQUENCE_LENGTH - 1): # 注意边界处理，跳过首尾位置
            left_dimer = sequence[i-1:i+1] # 结束于位置 i 的二核苷酸
            right_dimer = sequence[i:i+2]  # 开始于位置 i 的二核苷酸
            
            # 除以 (SEQUENCE_LENGTH - 1) 是因为有这么多可能的二核苷酸起始位置
            feature_matrix[0, i] = dimer_counts[left_dimer] / (SEQUENCE_LENGTH - 1)
            feature_matrix[1, i] = dimer_counts[right_dimer] / (SEQUENCE_LENGTH - 1)
        return feature_matrix

    def extract_trimer_frequencies(self, sequence):
        """
        计算三核苷酸的全局频率，并生成对应的特征矩阵。
        feature_matrix[0, i] 存储以 i-1 为开始位置的三核苷酸的全局频率。
        """
        trimer_counts = {nt1 + nt2 + nt3: 0.0 for nt1 in 'ACGT' for nt2 in 'ACGT' for nt3 in 'ACGT'}
        for i in range(SEQUENCE_LENGTH - 2): # 计算序列中所有三核苷酸的全局计数
            trimer = sequence[i:i+3]
            trimer_counts[trimer] += 1

        feature_matrix = np.zeros((1, SEQUENCE_LENGTH))
        for i in range(1, SEQUENCE_LENGTH - 2): # 注意边界处理
            trimer = sequence[i-1:i+2] # 以 i-1 为起始，i+1 为结束的三核苷酸
            feature_matrix[0, i] = trimer_counts[trimer] / (SEQUENCE_LENGTH - 2) # 可能的三核苷酸起始位置总数
        return feature_matrix

    def __len__(self):
        return len(self.sequences)

三、卷积神经网络模型定义 本节构建了一个用于DNA增强子识别的卷积神经网络 (CNN) 模型，命名为 SequenceClassifierCNN。该模型采用多层卷积、批量归一化、ReLU激活函数和最大池化操作，最终通过全连接层输出二分类预测结果。

class SequenceClassifierCNN(nn.Module):
    """
    用于DNA序列分类的卷积神经网络模型。
    输入特征通道为8，代表独热编码及K-mer频率特征。
    """
    def __init__(self):
        super(SequenceClassifierCNN, self).__init__()

        # 第一卷积块: 包含3个Conv1D层，每个层后接BatchNorm和ReLU激活。
        # 输入通道为8 (4独热 + 1一核苷酸频率 + 2二核苷酸频率 + 1三核苷酸频率)
        self.conv_block_A_layer1 = nn.Conv1d(8, CONV_FILTERS_BLOCK1, CONVOLUTION_KERNEL_SIZE, padding='same')
        self.bn_block_A_layer1 = nn.BatchNorm1d(CONV_FILTERS_BLOCK1)

        self.conv_block_A_layer2 = nn.Conv1d(CONV_FILTERS_BLOCK1, CONV_FILTERS_BLOCK1, CONVOLUTION_KERNEL_SIZE, padding='same')
        self.bn_block_A_layer2 = nn.BatchNorm1d(CONV_FILTERS_BLOCK1)

        self.conv_block_A_layer3 = nn.Conv1d(CONV_FILTERS_BLOCK1, CONV_FILTERS_BLOCK1, CONVOLUTION_KERNEL_SIZE, padding='same')
        self.bn_block_A_layer3 = nn.BatchNorm1d(CONV_FILTERS_BLOCK1)

        self.pooling_A = nn.MaxPool1d(POOLING_KERNEL_SIZE)

        # 第二卷积块: 同样包含3个Conv1D层，后接BatchNorm和ReLU激活。
        self.conv_block_B_layer1 = nn.Conv1d(CONV_FILTERS_BLOCK1, CONV_FILTERS_BLOCK2, CONVOLUTION_KERNEL_SIZE, padding='same')
        self.bn_block_B_layer1 = nn.BatchNorm1d(CONV_FILTERS_BLOCK2)

        self.conv_block_B_layer2 = nn.Conv1d(CONV_FILTERS_BLOCK2, CONV_FILTERS_BLOCK2, CONVOLUTION_KERNEL_SIZE, padding='same')
        self.bn_block_B_layer2 = nn.BatchNorm1d(CONV_FILTERS_BLOCK2)

        self.conv_block_B_layer3 = nn.Conv1d(CONV_FILTERS_BLOCK2, CONV_FILTERS_BLOCK2, CONVOLUTION_KERNEL_SIZE, padding='same')
        self.bn_block_B_layer3 = nn.BatchNorm1d(CONV_FILTERS_BLOCK2)

        self.pooling_B = nn.MaxPool1d(POOLING_KERNEL_SIZE)

        # 全连接层
        # 计算全连接层输入维度: 
        # (SEQUENCE_LENGTH / POOLING_KERNEL_SIZE / POOLING_KERNEL_SIZE) * CONV_FILTERS_BLOCK2
        # 即 (200 / 4 / 4) * 64 = 12 * 64 = 768 (使用整数除法)
        flattened_dim = int(SEQUENCE_LENGTH / POOLING_KERNEL_SIZE / POOLING_KERNEL_SIZE) * CONV_FILTERS_BLOCK2
        self.fully_connected_layer = nn.Linear(flattened_dim, DENSE_LAYER_DIM)
        
        # 输出层 (二分类，输出单个概率值)
        self.output_layer = nn.Linear(DENSE_LAYER_DIM, 1)

        # 损失函数和优化器
        self.loss_fn = nn.BCELoss() # 二元交叉熵损失
        self.optimizer = torch.optim.Adam(self.parameters(), lr=MODEL_LEARNING_RATE)

    def forward(self, input_data):
        """
        模型的前向传播逻辑。
        :param input_data: 输入的DNA序列特征张量 (batch_size, channels, sequence_length)。
        :return: 经过Sigmoid激活后的二分类预测概率。
        """
        batch_size = input_data.size(0)

        # 执行第一卷积块操作
        output = F.relu(self.bn_block_A_layer1(self.conv_block_A_layer1(input_data)))
        output = F.relu(self.bn_block_A_layer2(self.conv_block_A_layer2(output)))
        output = F.relu(self.bn_block_A_layer3(self.conv_block_A_layer3(output)))
        output = self.pooling_A(output) # 应用最大池化

        # 执行第二卷积块操作
        output = F.relu(self.bn_block_B_layer1(self.conv_block_B_layer1(output)))
        output = F.relu(self.bn_block_B_layer2(self.conv_block_B_layer2(output)))
        output = F.relu(self.bn_block_B_layer3(self.conv_block_B_layer3(output)))
        output = self.pooling_B(output) # 应用最大池化

        # 展平特征图并连接全连接层
        output = output.view(batch_size, -1) # 展平操作
        output = F.relu(self.fully_connected_layer(output))

        # 通过Sigmoid激活函数输出最终预测概率
        final_prediction = torch.sigmoid(self.output_layer(output))

        return final_prediction

四、模型训练流程 execute_single_epoch_training 函数负责在一个完整的 epoch 内训练模型。它会遍历所有训练数据，执行前向传播、损失计算、反向传播以及模型参数的更新。同时，该函数还支持在每个 epoch 开始时调整优化器的学习率。

def execute_single_epoch_training(model, data_loader, current_learning_rate):
    """
    执行单个训练 epoch。
    :param model: 要训练的PyTorch模型。
    :param data_loader: 训练数据加载器。
    :param current_learning_rate: 当前epoch的学习率。
    :return: 该epoch的平均训练损失。
    """
    model.train() # 设置模型为训练模式
    
    # 更新优化器的学习率
    for param_group in model.optimizer.param_groups:
        param_group['lr'] = current_learning_rate

    epoch_total_loss = 0.0
    num_trained_samples = 0

    for batch_idx, (inputs, labels) in enumerate(data_loader, 0):
        current_batch_size = inputs.size(0)
        num_trained_samples += current_batch_size

        # 数据类型转换，确保标签维度与BCELoss兼容
        inputs = inputs.float()
        labels = labels.float().unsqueeze(1) # 将标签从 (batch_size,) 转换为 (batch_size, 1)

        # 将数据移动到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()
        # 执行前向传播，获取模型预测
        predictions = model(inputs)
        # 计算损失
        loss = model.loss_fn(predictions, labels)
        # 执行反向传播，计算梯度
        loss.backward()
        # 更新模型参数
        model.optimizer.step()

        epoch_total_loss += loss.item() * current_batch_size # 累加批次损失

    avg_train_loss = epoch_total_loss / num_trained_samples # 计算平均训练损失
    print(f"当前epoch训练样本总数: {num_trained_samples}")
    return avg_train_loss

五、模型性能评估 assess_model_performance 函数用于评估给定模型在特定数据集上的表现。它加载预训练模型权重，将其设置为评估模式，并计算包括平均损失、AUC (Area Under the Receiver Operating Characteristic Curve)、准确率 (Accuracy)、灵敏度 (Sensitivity)、特异度 (Specificity) 和马修斯相关系数 (MCC) 在内的多项关键指标。

def assess_model_performance(model_filepath, data_loader):
    """
    评估PyTorch模型在给定数据加载器上的性能。
    :param model_filepath: 预训练模型的文件路径。
    :param data_loader: 用于评估的数据加载器。
    :return: 包含各项评估指标的字典。
    """
    model = SequenceClassifierCNN() # 实例化模型
    if torch.cuda.is_available():
        model.cuda() # 将模型移动到GPU (如果可用)
    
    model.load_state_dict(torch.load(model_filepath)) # 加载模型权重
    model.eval() # 设置为评估模式 (禁用Dropout, BatchNorm等随机性操作)

    total_loss = 0.0
    total_samples = 0
    
    true_labels_list = [] # 存储真实标签
    predicted_probs_list = [] # 存储预测概率

    with torch.no_grad(): # 在评估模式下禁用梯度计算，节省内存并加速
        for batch_idx, (inputs, labels) in enumerate(data_loader, 0):
            current_batch_size = inputs.size(0)
            total_samples += current_batch_size

            inputs = inputs.float()
            labels = labels.float().unsqueeze(1) # 确保标签维度与BCELoss兼容

            # 移动数据到GPU (如果可用)
            if torch.cuda.is_available():
                inputs = Variable(inputs.cuda())
                labels = Variable(labels.cuda())
            else:
                inputs = Variable(inputs)
                labels = Variable(labels)

            predictions = model(inputs) # 执行前向传播
            loss = model.loss_fn(predictions, labels) # 计算损失

            total_loss += loss.item() * current_batch_size # 累加损失
            
            true_labels_list.extend(labels.squeeze(1).cpu().numpy().tolist()) # 收集真实标签
            predicted_probs_list.extend(predictions.squeeze(1).cpu().numpy().tolist()) # 收集预测概率

    avg_evaluation_loss = total_loss / total_samples # 计算平均评估损失
    print(f"\n评估样本总数: {total_samples}")
    print(f"平均评估损失: {avg_evaluation_loss:.4f}")
    
    # 计算AUC (ROC曲线下面积)
    auc_score = metrics.roc_auc_score(true_labels_list, predicted_probs_list)
    print(f"AUC: {auc_score:.4f}")

    # 将预测概率转换为二分类标签 (使用0.5作为阈值)
    predicted_binary_labels = [int(prob > 0.5) for prob in predicted_probs_list]
    
    # 计算其他分类指标 (准确率、混淆矩阵、灵敏度、特异度和MCC)
    accuracy, conf_matrix, sensitivity, specificity, mcc_score = \
        compute_classification_metrics(true_labels_list, predicted_binary_labels)
    
    print(f"准确率 (ACC): {accuracy:.4f}")
    
    evaluation_results = {
        'avg_loss': avg_evaluation_loss,
        'accuracy': accuracy,
        'confusion_matrix': conf_matrix,
        'sensitivity': sensitivity,
        'specificity': specificity,
        'mcc': mcc_score,
        'auc': auc_score,
        'predicted_probabilities': predicted_probs_list,
        'actual_labels': true_labels_list,
        'predicted_binary_labels': predicted_binary_labels
    }
    
    return evaluation_results

计算混淆矩阵和分类指标函数 compute_classification_metrics 该函数根据真实标签和模型预测标签计算混淆矩阵及其衍生的分类指标，用于全面评估模型的分类性能。

def compute_classification_metrics(actual_labels, predicted_labels):
    """
    计算混淆矩阵和相关分类指标。
    混淆矩阵定义为: confusion_matrix[预测类别][真实类别]
    [[TN, FN],
     [FP, TP]]
    :param actual_labels: 真实标签列表。
    :param predicted_labels: 模型预测标签列表。
    :return: 准确率、混淆矩阵、灵敏度、特异度、MCC值。
    """
    num_correct_predictions = 0
    confusion_matrix_arr = np.zeros((2, 2)) # 初始化混淆矩阵

    for i in range(len(actual_labels)):
        pred = predicted_labels[i]
        true = actual_labels[i]
        confusion_matrix_arr[pred][true] += 1 # 填充混淆矩阵

        if true == pred:
            num_correct_predictions += 1

    total_samples = len(actual_labels)
    accuracy = num_correct_predictions / total_samples if total_samples > 0 else 0

    # 从混淆矩阵中提取True Positive (TP), True Negative (TN), False Positive (FP), False Negative (FN)
    TN = confusion_matrix_arr[0][0] # 真实为负，预测为负
    FN = confusion_matrix_arr[0][1] # 真实为正，预测为负
    FP = confusion_matrix_arr[1][0] # 真实为负，预测为正
    TP = confusion_matrix_arr[1][1] # 真实为正，预测为正

    # 计算特异度 (Specificity): 真实负例中被正确识别的比例
    specificity = TN / (TN + FP) if (TN + FP) > 0 else 0
    # 计算灵敏度 (Sensitivity / Recall): 真实正例中被正确识别的比例
    sensitivity = TP / (TP + FN) if (TP + FN) > 0 else 0

    # 计算Matthews Correlation Coefficient (MCC)
    mcc_numerator = (TP * TN) - (FP * FN)
    mcc_denominator = math.sqrt((TP + FP) * (TP + FN) * (TN + FP) * (TN + FN))
    mcc_score = mcc_numerator / mcc_denominator if mcc_denominator != 0 else 0

    return accuracy, confusion_matrix_arr, sensitivity, specificity, mcc_score

六、K折交叉验证训练与模型管理 本节包含K折交叉验证的核心逻辑。perform_fold_training_and_validation 函数负责单折的训练和验证，包括模型实例化、训练循环、最佳模型保存和性能日志记录。而 initiate_kfold_cross_validation 函数则协调整个K折交叉验证过程，加载数据并划分训练/验证集。

def perform_fold_training_and_validation(fold_data_split, fold_number):
    """
    执行单次K折交叉验证的训练和验证流程。
    在每个epoch结束后，模型会在验证集上进行评估，并保存具有最佳验证损失的模型。
    :param fold_data_split: 包含训练集和验证集序列及标签的字典。
    :param fold_number: 当前折的编号。
    """
    # 创建训练集和验证集Dataset对象
    train_dataset = NucleotideSequenceDataset(fold_data_split["train_sequences"], fold_data_split["train_labels"])
    val_dataset = NucleotideSequenceDataset(fold_data_split["val_sequences"], fold_data_split["val_labels"])

    # 创建DataLoader，用于批量加载数据
    train_data_loader = DataLoader(dataset=train_dataset, batch_size=32, shuffle=True, num_workers=4)
    val_data_loader = DataLoader(dataset=val_dataset, batch_size=32, shuffle=False, num_workers=4)

    # 初始化模型实例
    model_instance = SequenceClassifierCNN()
    print(f"--- 实例化折叠 {fold_number + 1} 的CNN模型 ---")
    if torch.cuda.is_available():
        model_instance.cuda() # 将模型放置到GPU上

    best_validation_loss = float('inf') # 初始化最佳验证损失为无穷大
    best_epoch_num = 0 # 记录最佳模型的epoch编号
    
    # 用于记录每个epoch的训练和验证损失
    training_losses_per_epoch = []
    validation_losses_per_epoch = []

    # 训练循环
    for epoch in range(TRAINING_EPOCHS):
        print(f"\n####### K折验证: 第 {fold_number + 1} 折 / 总 {k_folds.n_splits} 折; 迭代 {epoch + 1}/{TRAINING_EPOCHS} #######")
        # 训练一个epoch
        current_train_loss = execute_single_epoch_training(model_instance, train_data_loader, MODEL_LEARNING_RATE)
        print(f"Epoch {epoch+1} 平均训练损失: {current_train_loss:.4f}")

        # 保存当前模型权重到临时文件，以便在验证集上进行评估
        current_model_temp_path = os.path.join(OUTPUT_MODEL_DIR, TEMP_MODEL_SAVE_PATH)
        torch.save(model_instance.state_dict(), current_model_temp_path)

        # 在验证集上评估当前模型
        validation_results = assess_model_performance(current_model_temp_path, val_data_loader)
        current_val_loss = validation_results['avg_loss']
        
        # 记录本epoch的训练和验证损失
        training_losses_per_epoch.append(current_train_loss)
        validation_losses_per_epoch.append(current_val_loss)

        # 如果当前验证损失优于历史最佳，则更新最佳模型
        if current_val_loss < best_validation_loss:
            best_validation_loss = current_val_loss
            best_epoch_num = epoch
            best_model_save_path = os.path.join(OUTPUT_MODEL_DIR, "best_model_for_fold.pth")
            torch.save(model_instance.state_dict(), best_model_save_path) # 保存最佳模型权重
            print(f"模型已更新并保存，最佳验证损失: {best_validation_loss:.4f} (对应Epoch: {best_epoch_num + 1})")

    # 加载并最终保存该折的最佳模型
    final_model_name = f"enhancer_fold_{fold_number + 1}_epoch_{best_epoch_num + 1}.pth"
    final_model_path = os.path.join(OUTPUT_MODEL_DIR, final_model_name)
    
    # 从保存的最佳模型路径加载权重
    best_model_loaded_state_dict = torch.load(os.path.join(OUTPUT_MODEL_DIR, "best_model_for_fold.pth"))
    model_instance.load_state_dict(best_model_loaded_state_dict)
    torch.save(model_instance.state_dict(), final_model_path) # 保存最终命名的最佳模型

    # 记录训练和验证损失到CSV文件
    log_filename = os.path.join(OUTPUT_MODEL_DIR, f"loss_log_fold_{fold_number + 1}_epoch_{best_epoch_num + 1}.csv")
    with open(log_filename, mode='w', newline='') as log_file:
        csv_writer = csv.writer(log_file, delimiter=',')
        csv_writer.writerow(['Epoch', 'Train Loss', 'Validation Loss'])
        for i in range(TRAINING_EPOCHS):
            csv_writer.writerow([i + 1, training_losses_per_epoch[i], validation_losses_per_epoch[i]])

    print(f"\n####### 最终评估 折叠 {fold_number + 1} 的最佳模型 (对应Epoch: {best_epoch_num + 1}) #######")
    final_val_results = assess_model_performance(final_model_path, val_data_loader)
    print(f"验证集最终结果 (AUC): {final_val_results['auc']:.4f}, 准确率: {final_val_results['accuracy']:.4f}")

initiate_kfold_cross_validation 函数是整个K折交叉验证流程的入口点。它负责加载原始数据集，将其分为增强子和非增强子序列，并为每个序列分配相应的标签。然后，它使用KFold或StratifiedKFold将数据集划分为训练集和验证集，并调用 perform_fold_training_and_validation 函数对每个折进行训练和评估。

• KFold vs. StratifiedKFold:
• KFold: 简单地将数据集划分为K个互不重叠的子集，每次选取其中一个作为验证集，其余作为训练集。适用于类别分布相对均匀的数据集。
• StratifiedKFold: 在划分数据集时，会确保每个子集（折）中各类别的样本比例与原始数据集的整体比例保持一致。这在处理类别不平衡的数据集时尤为重要，可以有效避免某些折中缺少特定类别的样本，从而使模型评估更可靠。 在增强子预测这类可能存在类别不平衡问题（正负样本数量差异）的任务中，通常推荐使用StratifiedKFold。

def initiate_kfold_cross_validation(): """ 加载数据并执行K折交叉验证训练流程。 """ print("\n ==> 正在加载训练数据...")

加载增强子和非增强子序列

strong_enhancers = read_sequence_file('intersect_CD.fa')

weak_enhancers = read_sequence_file('data_week_enhancers.txt') # 原始代码中此行被注释掉

all_enhancers = strong_enhancers # 如果有弱增强子，可在此处合并: strong_enhancers + weak_enhancers print(f"加载增强子序列数量: {len(all_enhancers)}")

non_enhancers = read_sequence_file('noenhancer_seq_CD.fa') print(f"加载非增强子序列数量: {len(non_enhancers)}")

分配标签: 增强子为1，非增强子为0

labels_enhancers = np.ones((len(all_enhancers), 1)) labels_non_enhancers = np.zeros((len(non_enhancers), 1))

合并所有序列数据和标签

all_sequences = np.concatenate((all_enhancers, non_enhancers)) all_labels = np.concatenate((labels_enhancers, labels_non_enhancers))

初始化K折交叉验证器

使用StratifiedKFold以确保每个折中类别分布均衡，特别适用于分类任务

global k_folds # 声明为全局变量，以便在perform_fold_training_and_validation中使用 k_folds = StratifiedKFold(n_splits=5, shuffle=True, random_state=GLOBAL_RANDOM_SEED)

若类别均衡，也可使用 KFold: k_folds = KFold(n_splits=5, shuffle=True, random_state=GLOBAL_RANDOM_SEED)

fold_counter = 0

遍历K折划分，获取每次训练和验证的索引

for train_indices, val_indices in k_folds.split(all_sequences, all_labels): print(f"\n--- 开始处理第 {fold_counter + 1} 折 ---")

根据索引划分训练集和验证集数据与标签

train_seqs, val_seqs = all_sequences[train_indices], all_sequences[val_indices] train_lbls, val_lbls = all_labels[train_indices], all_labels[val_indices]

将当前折的数据打包成字典

current_fold_data = { "train_sequences": train_seqs, "train_labels": train_lbls, "val_sequences": val_seqs, "val_labels": val_lbls }

(可选) 检查当前折的数据集序列长度

validate_sequence_lengths(current_fold_data)

对当前折数据进行训练和验证

perform_fold_training_and_validation(current_fold_data, fold_counter) fold_counter += 1

<h2>七、数据集有效性验证 (辅助功能)</h2>
<p>validate_sequence_lengths 函数用于检查数据集中所有DNA序列的长度是否符合预设的 SEQUENCE_LENGTH。这有助于在模型训练前发现并避免由于数据格式不一致导致的问题。</p>
def validate_sequence_lengths(dataset_split):
    """
    检查数据集中所有序列的长度是否与预定义的 SEQUENCE_LENGTH 匹配。
    :param dataset_split: 包含训练和验证序列的字典。
    """
    print("\n==> 正在检查数据集序列长度...")
    
    # 检查训练集序列
    train_sequences = dataset_split["train_sequences"]
    error_count_train = 0
    for seq in train_sequences:
        if len(seq) != SEQUENCE_LENGTH:
            error_count_train += 1
    if error_count_train > 0:
        print(f"训练集中发现 {error_count_train} 个长度不正确的序列。")
    else:
        print("训练集序列长度检查：OK！")

    # 检查验证集序列
    val_sequences = dataset_split["val_sequences"]
    error_count_val = 0
    for seq in val_sequences:
        if len(seq) != SEQUENCE_LENGTH:
            error_count_val += 1
    if error_count_val > 0:
        print(f"验证集中发现 {error_count_val} 个长度不正确的序列。")
    else:
        print("验证集序列长度检查：OK！")