Deep-gami部分代码解读

import argparse
import pandas as pd
import os
import time
import random
from collections import Counter
import torch
from torch import nn
from torch.utils.data import TensorDataset, DataLoader
import numpy as np
import sklearn.metrics as skm
from sklearn.model_selection import train_test_split, KFold
import DeepGamiUtils as ut
from DeepGamiModel import DeepGami
from scipy.stats import bernoulli
import os

argparse 是一个 Python 模块，提供了一种简单的方式来解析命令行参数和选项。它是 Python 标准库的一部分，常用于为 Python 程序创建命令行界面（CLI）。

通过 argparse，你可以定义程序期望的参数和选项，指定它们的类型，提供帮助信息，并根据用户输入处理它们的值。

import argparse
# 创建 ArgumentParser 对象
parser = argparse.ArgumentParser(description='一个简单的程序')
# 添加参数
parser.add_argument('input', help='输入文件')
parser.add_argument('-o', '--output', help='输出文件')
# 解析命令行参数
args = parser.parse_args()
# 访问参数的值
input_file = args.input
output_file = args.output
# 在程序中使用参数的值
print(f'输入文件: {input_file}')
print(f'输出文件: {output_file}')

在这个示例中，程序期望一个必需的输入文件参数，并允许使用 -o 或 --output 选项来指定可选的输出文件，help 参数用于为每个参数提供描述。

当你在命令行中运行程序时，你需要提供输入文件，并且可以选择指定输出文件：

argparse 模块负责解析命令行参数并将它们的值存储在 args 对象中。然后，你可以访问这些值并根据需要在程序中使用。argparse 还提供了处理不同类型参数（字符串、整数、浮点数）、定义默认值、指定选择项等功能。它是一个多功能且强大的用于在 Python 中构建命令行界面的模块。

os.environ["CUDA_VISIBLE_DEVICES"] = "0"
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

这段代码用于设置CUDA可见设备和选择PyTorch的计算设备。os.environ["CUDA_VISIBLE_DEVICES"] = "0" 这行代码将环境变量 CUDA_VISIBLE_DEVICES 设置为 "0"，表示只使用 CUDA 设备索引为 0 的 GPU。这可以用来限制程序在多个 GPU 上的运行。

device = torch.device("cuda" if torch.cuda.is_available() else "cpu") 这行代码创建一个 PyTorch 设备对象，根据系统中是否有可用的 CUDA 设备来选择是使用 CUDA 还是 CPU。如果有可用的 CUDA 设备，将选择 "cuda"，否则选择 "cpu"。这样，你可以在程序中使用 device 变量来指定计算设备，以便将张量和模型移动到相应的设备上进行计算。

注意，要运行这段代码，你需要确保已经正确安装了 PyTorch 并且有可用的 CUDA 设备。如果没有 CUDA 设备或者 CUDA 配置有问题，代码将会选择使用 CPU 进行计算。

def preprocess(inp1, inp2, oup):
    """ Function to direct the input and ouput to CPU vs GPU"""
    return inp1.float().to(device), inp2.float().to(device), oup.int().to(device)

这是一个名为 preprocess 的函数，用于将输入和输出数据移动到 CPU 或 GPU 上进行处理。它接受三个参数 inp1、inp2 和 oup，表示输入和输出数据。

这个函数的目的是将输入和输出数据转换为特定的数据类型，并将它们移动到预先选择好的计算设备（CPU 或 GPU）上进行计算。在函数内部，首先使用 .float() 方法将 inp1 和 inp2 转换为浮点数类型，然后使用 .int() 方法将 oup 转换为整数类型。接着，使用 .to(device) 方法将转换后的张量移动到事先选择好的计算设备上。这个函数假设你已经在代码中定义了 device 变量，并根据可用的计算设备选择了合适的值（例如 "cuda" 或 "cpu"）。调用 preprocess 函数时，它将根据 device 变量的值来决定将数据移动到 CPU 还是 GPU 上进行处理，并返回处理后的结果。

class WrappedDataLoader:
    def __init__(self, dl, func):
        self.dl = dl
        self.func = func
    
    def __len__(self):
        return len(self.dl)
    
    def __iter__(self):
        batches = iter(self.dl)
        for b in batches:
            yield self.func(*b)

这是一个名为 WrappedDataLoader 的类，用于包装一个数据加载器并应用给定的函数到每个批次数据上进行处理。这个类可以用于自定义数据预处理或转换。

这个类有三个方法：__init__、__len__ 和 __iter__

__init__(self, dl, func) 方法用于初始化 WrappedDataLoader 类的实例。它接受两个参数：
dl：一个数据加载器对象，比如 PyTorch 的 DataLoader。
func：一个函数，用于对每个批次数据进行处理的操作。
__len__(self) 方法返回包装数据加载器的长度，即批次的数量。这使得你可以像使用原始数据加载器一样使用 len() 函数来获取数据集的大小。
__iter__(self) 方法返回一个迭代器，它通过原始数据加载器 dl 迭代批次数据，并在每个批次上应用函数 func 进行处理。使用 yield 语句可以将处理后的批次数据逐个地返回。
这个类的作用是将数据加载器中的每个批次数据传递给指定的函数进行处理，并以迭代器的形式返回处理后的结果。这样，你可以在训练过程中对批次数据进行自定义的操作或转换，例如数据预处理、数据增强等。

在 Python 类的方法中，self 是一个特殊的参数，用于表示该方法所属的实例对象。在 WrappedDataLoader 类中，self 在方法定义中作为第一个参数出现，如 def __init__(self, dl, func)。

self 的作用是引用当前类的实例对象，它允许你在类的方法中访问和操作实例的属性和方法。通过使用 self，你可以在类的方法中访问类的属性、调用其他方法，以及与实例对象进行交互。

在 WrappedDataLoader 类中，self.dl 和 self.func 是类的属性，通过 self 可以在类的方法中对它们进行访问。例如，self.dl 引用了传递给 WrappedDataLoader 构造函数的数据加载器对象，而 self.func 引用了传递给构造函数的处理函数。

通过在方法中使用 self，你可以在整个类的范围内共享和操作这些属性，包括在 __len__ 和 __iter__ 方法中使用它们。

总而言之，self 的作用是引用当前类的实例对象，使得你可以在类的方法中访问和操作实例的属性和方法。

def bio_dropout(data, drop_type='random'):
    """ Function to get bernaoulli dropout with biological prior"""
    new_data = data.copy()
    
    pos_len = len(np.where(data == 1)[0])
    neg_len = len(np.where(data == 0)[0])
    
    if drop_type == 'random':
        pos_rv = bernoulli(0.5)
        neg_rv = bernoulli(0.5)
    else:
        pos_rv = bernoulli(0.8)
        neg_rv = bernoulli(0.2)

    new_data[data == 1] = pos_rv.rvs(pos_len)
    new_data[data == 0] = neg_rv.rvs(neg_len)
    return new_data

这是一个名为 bio_dropout 的函数，用于实现带有生物先验的伯努利丢弃（dropout）操作。

这个函数接受两个参数：

data：一个 numpy 数组，表示输入的数据。
drop_type：一个字符串参数，表示丢弃类型，默认为 'random'。
函数的实现逻辑如下：

首先，创建一个名为 new_data 的副本，用于存储处理后的数据。
接下来，计算输入数据中值为 1 和 0 的元素的数量，分别保存在 pos_len 和 neg_len 变量中。
根据 drop_type 参数的值，选择不同的伯努利分布：
如果 drop_type 为 'random'，则使用参数为 0.5 的伯努利分布，表示正负例丢弃的概率相同。
如果 drop_type 不为 'random'，则使用参数为 0.8 和 0.2 的伯努利分布，表示正例和负例丢弃的概率不同。
对于输入数据中值为 1 的元素，将其替换为从正例伯努利分布中采样的随机变量的值。
对于输入数据中值为 0 的元素，将其替换为从负例伯努利分布中采样的随机变量的值。
返回处理后的 new_data 数组。
这个函数的目的是模拟使用伯努利丢弃操作对输入数据进行处理，并引入生物先验。具体来说，它根据给定的丢弃类型，在输入数据中的正例和负例之间进行随机丢弃，从而得到处理后的数据。

def get_binary_performance(y_true, y_score):
    """Function to return the acc, bacc, and auc"""
    y_pred = np.argmax(y_score, 1)
    y_truth = np.argmax(y_true, 1)

    #auc = skm.roc_auc_score(y_true, y_score)
    acc = skm.accuracy_score(y_truth, y_pred)
    bacc = skm.balanced_accuracy_score(y_truth, y_pred)
    return acc, bacc#, auc

这是一个名为 get_binary_performance 的函数，用于计算二分类问题的性能指标，包括准确率（accuracy）、平衡准确率（balanced accuracy）和曲线下面积（AUC）。

这个函数接受两个参数：

y_true：一个 numpy 数组，表示真实的标签（ground truth）。
y_score：一个 numpy 数组，表示分类器的预测分数。
函数的实现逻辑如下：

首先，使用 np.argmax() 函数找到每个样本的预测类别。y_score 是一个预测分数的数组，通过 np.argmax(y_score, 1) 可以得到每个样本的预测类别。
然后，使用 np.argmax() 函数找到每个样本的真实类别。y_true 是一个真实标签的数组，通过 np.argmax(y_true, 1) 可以得到每个样本的真实类别。
接下来，使用 skm.accuracy_score() 函数计算预测准确率（accuracy），并将结果保存在 acc 变量中。
使用 skm.balanced_accuracy_score() 函数计算平衡准确率（balanced accuracy），并将结果保存在 bacc 变量中。
最后，函数返回计算得到的准确率和平衡准确率。
这个函数的目的是计算二分类问题的性能指标，其中准确率表示预测正确的样本比例，平衡准确率考虑了不同类别之间的样本不平衡问题。函数还注释掉了计算曲线下面积（AUC）的代码行，如果需要计算 AUC，可以取消注释该行代码，并将 AUC 的结果一并返回。

def get_classification_performance(y_true, y_score, task='binary'):
    """Function to return various performance metrics"""
    auc = 0.0
    y_pred = None
    if task=='binary':
        auc = skm.roc_auc_score(y_true, y_score)
        y_pred = np.where(y_score<0.5, 0, 1)
    else:
        print(y_score)
        auc = skm.roc_auc_score(y_true, y_score, average="weighted", multi_class="ovr")
        y_pred = np.argmax(y_score, 1)

    acc = skm.accuracy_score(y_true, y_pred)
    bacc = skm.balanced_accuracy_score(y_true, y_pred)
    return acc, bacc, auc

这是一个名为 get_classification_performance 的函数，用于计算分类任务的性能指标，包括准确率（accuracy）、平衡准确率（balanced accuracy）和曲线下面积（AUC）。

这个函数接受三个参数：

y_true：一个 numpy 数组，表示真实的标签（ground truth）。
y_score：一个 numpy 数组，表示分类器的预测分数。
task：一个字符串，表示任务类型，默认为 'binary' 表示二分类任务。
函数的实现逻辑如下：

首先，初始化 auc 为 0.0，y_pred 为 None。
如果任务类型为二分类（task=='binary'），则使用 skm.roc_auc_score() 函数计算二分类问题的 AUC，并将结果保存在 auc 变量中。然后，使用阈值 0.5 对预测分数进行二值化，即将大于等于 0.5 的预测分数标记为 1，小于 0.5 的预测分数标记为 0，将结果保存在 y_pred 变量中。
如果任务类型不是二分类，即多分类任务，使用 skm.roc_auc_score() 函数计算多分类问题的 AUC。通过指定 average="weighted" 和 multi_class="ovr" 参数，计算加权平均的 AUC，并将结果保存在 auc 变量中。然后，使用 np.argmax() 函数找到每个样本的预测类别，并将结果保存在 y_pred 变量中。接下来，使用 skm.accuracy_score() 函数计算预测准确率（accuracy），并将结果保存在 acc 变量中。使用 skm.balanced_accuracy_score() 函数计算平衡准确率（balanced accuracy），并将结果保存在 bacc 变量中。最后，函数返回计算得到的准确率、平衡准确率和 AUC。这个函数的目的是计算分类任务的性能指标，可以根据任务类型选择计算二分类或多分类问题的 AUC，并返回准确率和平衡准确率。

def predict(model, data_dl, estimate, task):
    """预测样本的函数"""
    predictions, truth = [], []
    for snps, gex, yb in data_dl:
        yhat = model(snps, gex, estimate)
        if task == 'binary':
            predictions.extend(yhat.detach().cpu().numpy())
            truth.extend(yb.detach().cpu().numpy())
        else:
            predictions.extend(torch.softmax(yhat, dim=1).detach().cpu().numpy())
            truth.extend(np.argmax(yb.detach().cpu().numpy(), 1))

    predictions = np.asarray(predictions)
    truth = np.asarray(truth)
    return predictions, truth

这是一个名为 predict 的函数，用于对样本进行预测。这个函数接受四个参数：

model：一个预训练的模型对象，用于进行预测。
data_dl：一个数据加载器（DataLoader），用于加载待预测的数据。
estimate：一个估计参数，用于模型的预测。
task：一个字符串，表示任务类型，默认为 'binary' 表示二分类任务。
函数的实现逻辑如下：

首先，初始化 predictions 和 truth 为空列表，用于存储预测结果和真实标签。
对于数据加载器中的每个批次（snps, gex, yb），使用预训练的模型 model 对 SNP（single nucleotide polymorphism）数据和 GEX（gene expression）数据进行预测，得到预测结果 yhat。
如果任务类型为二分类（task=='binary'），则将预测结果 yhat 转换为 numpy 数组，并将其添加到 predictions 列表中。同时，将真实标签 yb 也转换为 numpy 数组，并添加到 truth 列表中。
如果任务类型不是二分类，即多分类任务，首先使用 torch.softmax() 函数对预测结果 yhat 进行 softmax 操作，以获取每个类别的概率。然后，将 softmax 后的结果转换为 numpy 数组，并添加到 predictions 列表中。同时，使用 np.argmax() 函数找到每个样本的预测类别，并将其添加到 truth 列表中。
最后，将 predictions 和 truth 转换为 numpy 数组，并作为结果返回。
这个函数的目的是使用预训练的模型对样本进行预测，并返回预测结果和真实标签。根据任务类型，可以选择返回预测的概率（针对多分类任务）或直接的预测类别（针对二分类任务）。

def train_step(model, loss_fn, data_dl, l1_reg, corr_reg, estimate, task, opt):
    """ Function to train an epoch"""
    tot_loss = 0.0
    predictions, truth = [], []
    corr_loss_fn = nn.MSELoss()

    for snps, gex, yb in data_dl:
        loss = 0.0

        yhat = model(snps, gex, estimate)
        _, _, Cg, Cg_est = model.get_intermediate_layers(snps, gex)

        # for param in modl.parameters():
        #     loss += l1_reg * torch.sum(torch.abs(param))
        pred_loss = loss_fn(yhat, yb.float())
        corr_loss = corr_loss_fn(Cg, Cg_est)
        loss += pred_loss + corr_reg*corr_loss
        
        if opt is not None:
            opt.zero_grad()
            loss.backward(retain_graph=True)
            opt.step()

        if task == 'binary':
            predictions.extend(yhat.detach().cpu().numpy())
            truth.extend(yb.detach().cpu().numpy())
        else:
            predictions.extend(torch.softmax(yhat, dim=1).detach().cpu().numpy())
            truth.extend(np.argmax(yb.detach().cpu().numpy(), 1))
        tot_loss += loss.item()

    predictions = np.asarray(predictions)
    truth = np.asarray(truth)

    return tot_loss/len(data_dl), predictions, truth

这是一个名为 train_step 的函数，用于训练一个 epoch。

这个函数接受七个参数：

model：一个模型对象，用于训练。
loss_fn：一个损失函数对象，用于计算损失。
data_dl：一个数据加载器（DataLoader），用于加载训练数据。
l1_reg：L1 正则化的系数，用于控制 L1 正则化的强度。
corr_reg：相关性损失的系数，用于控制相关性损失的权重。
estimate：一个估计参数，用于模型的训练。
task：一个字符串，表示任务类型。
opt：一个优化器对象，用于更新模型的参数。
函数的实现逻辑如下：

首先，初始化总损失 tot_loss 为 0.0，以及空的 predictions 和 truth 列表，用于存储预测结果和真实标签。另外，创建一个均方误差损失函数 corr_loss_fn。
对于数据加载器中的每个批次（snps, gex, yb），进行如下操作：
初始化损失 loss 为 0.0。
使用模型 model 对 SNP（single nucleotide polymorphism）数据和 GEX（gene expression）数据进行预测，得到预测结果 yhat，以及模型中间层的输出 Cg 和 Cg_est。
使用给定的损失函数 loss_fn 计算预测损失 pred_loss，并将其添加到总损失 loss 中。
使用均方误差损失函数 corr_loss_fn 计算相关性损失 corr_loss，并将其添加到总损失 loss 中，乘以相关性损失的系数 corr_reg。
如果优化器 opt 不为空，则进行优化步骤：
将优化器的梯度归零，使用 opt.zero_grad()。
根据总损失 loss 进行反向传播，使用 loss.backward(retain_graph=True)。
执行一步优化，使用 opt.step()。
根据任务类型，将预测结果 yhat 转换为 numpy 数组，并将其添加到 predictions 列表中。同时，将真实标签 yb 转换为 numpy 数组，并添加到 truth 列表中。
将当前批次的损失 loss.item() 添加到总损失 tot_loss 中。
将 predictions 和 truth 转换为 numpy 数组。
返回每个批次的平均损失 tot_loss / len(data_dl)，以及预测结果和真实标签。
这个函数的目的是在一个 epoch 中训练模型，并返回平均损失、预测结果和真实标签。在训练过程中，根据任务类型，可以选择返回预测的概率（针对多分类任务）或直接的预测类别（针对二分类任务）。此外，还可以使用 L1 正则化和相关性损失来约束模型的训练。如果提供了优化器对象，还会执行优化步骤来更新模型的参数。