【生物信息学】PyTorch 环境下的基因富集分析
发布时间
阅读量：
阅读量
在 PyTorch 环境中实现 基因富集分析（Gene Enrichment Analysis, GEA） ，可以结合深度学习的优势，实现更灵活的分析流程，例如基因表达数据的降维、聚类、富集预测等。
1. 环境准备
在开始基因富集分析前，您需要配置适当的环境：
pip install torch torchvision torchaudio pandas numpy scipy matplotlib seaborn
pip install gseapy # GSEA 工具包，用于基因富集分析
2. 数据准备
- 输入数据格式 ：
- RNA-Seq 或微阵列数据，通常以矩阵形式表示：行代表基因，列代表样本。
- 每个基因的表达值可用于排序（GSEA）或显著性筛选（ORA）。
示例加载数据：
import pandas as pd
# 假设输入数据为差异表达分析结果
diff_expr_data = pd.read_csv('diff_expr_results.csv') # 包含基因名、log2FC、p-value 等
gene_list = diff_expr_data['Gene Name'].tolist() # 基因名称列表
以下是一个简单的模拟基因表达分析结果数据集，可以作为测试数据使用。这个数据集包含基因名称、对数2倍数变化（log2FC）和P值等列。
示例数据集内容（diff_expr_results.csv）
| Gene Name | log2FC | P-value |
|---|---|---|
| Gene1 | 2.5 | 0.0001 |
| Gene2 | -1.8 | 0.003 |
| Gene3 | 0.5 | 0.02 |
| Gene4 | -2.1 | 0.00005 |
| Gene5 | 1.3 | 0.04 |
如何生成这个文件
可以使用以下代码生成一个 diff_expr_results.csv 文件：
import pandas as pd
# 构建模拟数据
data = {
"Gene Name": ["Gene1", "Gene2", "Gene3", "Gene4", "Gene5"],
"log2FC": [2.5, -1.8, 0.5, -2.1, 1.3],
"P-value": [0.0001, 0.003, 0.02, 0.00005, 0.04]
}
# 创建 DataFrame
df = pd.DataFrame(data)
# 保存为 CSV 文件
df.to_csv("diff_expr_results.csv", index=False)
print("文件已生成：diff_expr_results.csv")
生成文件后的操作
- 上述代码后，文件
diff_expr_results.csv会保存到当前脚本的工作目录中。 - 确保新生成的文件与代码的环境路径一致，或根据需要调整路径。
3. 基于 PyTorch 的数据预处理
- 数据标准化 ： 基因表达数据通常需要标准化，以便适应深度学习模型。
import torch
# 转换为张量并归一化
expression_data = torch.tensor(diff_expr_data['log2FC'].values, dtype=torch.float32)
normalized_data = (expression_data - torch.mean(expression_data)) / torch.std(expression_data)
- 构建基因排序 ： 根据基因的 log2FC 对基因进行排序，这是 GSEA 分析的基础。
sorted_genes = [gene for _, gene in sorted(zip(normalized_data.tolist(), gene_list), reverse=True)]
4. 深度学习应用
在 PyTorch 中，可以使用神经网络方法增强传统基因富集分析，例如预测功能注释或构建新的基因集间关系。
4.1 降维分析
使用自编码器对高维基因表达数据进行降维，以捕获主要特征。
import torch.nn as nn
# 定义自编码器结构
class Autoencoder(nn.Module):
def __init__(self, input_dim, latent_dim):
super(Autoencoder, self).__init__()
self.encoder = nn.Sequential(
nn.Linear(input_dim, 128),
nn.ReLU(),
nn.Linear(128, latent_dim)
)
self.decoder = nn.Sequential(
nn.Linear(latent_dim, 128),
nn.ReLU(),
nn.Linear(128, input_dim)
)
def forward(self, x):
latent = self.encoder(x)
reconstructed = self.decoder(latent)
return latent, reconstructed
# 模型初始化
input_dim = len(diff_expr_data) # 基因数量
latent_dim = 10 # 降维目标维度
model = Autoencoder(input_dim, latent_dim)
# 生成潜在空间数据
data_tensor = normalized_data.unsqueeze(0) # 添加批次维度
latent_space, reconstructed_data = model(data_tensor)
# 将潜在空间转换为 NumPy 数组
latent_space_np = latent_space.detach().numpy()
4.2 富集特征预测
通过深度学习模型预测基因集富集特征，例如特定条件下的功能活动。
# 定义简单的预测网络
class EnrichmentPredictor(nn.Module):
def __init__(self, input_dim, output_dim):
super(EnrichmentPredictor, self).__init__()
self.network = nn.Sequential(
nn.Linear(input_dim, 128),
nn.ReLU(),
nn.Linear(128, output_dim)
)
def forward(self, x):
return self.network(x)
# 初始化预测模型
output_dim = 5 # 目标富集类别数
predictor = EnrichmentPredictor(latent_dim, output_dim)
# 生成预测分布
predicted_output = predictor(latent_space)
predicted_output_np = predicted_output.squeeze().detach().numpy()
在上述代码中，我们已经定义了一个自编码器和一个预测网络，用于基因富集分析。为了可视化分析结果，可以采用以下方法：
富集结果可视化
import matplotlib.pyplot as plt
# 富集结果柱状图
categories = [f"Category {i+1}" for i in range(output_dim)]
plt.bar(categories, predicted_output_np, color='green', alpha=0.7)
plt.xlabel("Enrichment Categories")
plt.ylabel("Predicted Scores")
plt.title("Predicted Enrichment Scores for Categories")
plt.show()
5. 基因富集分析
5.1 基于排序的 GSEA 分析
使用 GSEAPy 工具直接 GSEA 分析，并将其整合到 PyTorch 中。（要用真实的基因数据或者按照规则生成的测试数据）
from gseapy import prerank
# 保存排序好的基因列表
with open('sorted_genes.rnk', 'w') as f:
for gene, score in zip(sorted_genes, normalized_data.tolist()):
f.write(f"{gene}\t{score}\n")
# 使用 GSEAPy GSEA
prerank_res = prerank(
rnk='sorted_genes.rnk',
gene_sets='KEGG_2019_Human', # 可选：GO_BP、Reactome 等
outdir='gsea_results',
permutation_num=1000
)
5.2 富集结果可视化
PyTorch 结合 Matplotlib 绘制 GSEA 结果。
import matplotlib.pyplot as plt
# 加载富集结果
enrichment_score = prerank_res.res2d['nes'].values # 富集评分
p_values = prerank_res.res2d['pval'].values
# 绘制富集条形图
plt.barh(prerank_res.res2d.index, enrichment_score, color='skyblue', alpha=0.8)
plt.xlabel('Normalized Enrichment Score (NES)')
plt.title('GSEA Results')
plt.show()
6. 优化与扩展
- 多任务学习 ：同时预测多个功能注释类别。
- 多组学整合 ：结合转录组、蛋白质组数据进行分析。
- 自定义损失函数 ：引入基于生物学知识的正则化项。
7. 总结
在 PyTorch 中实现基因富集分析，不仅可以复用深度学习的强大功能，还能增强传统分析方法的灵活性和扩展性。通过引入神经网络、自编码器等模型，研究者可以：
- 降维高维组学数据；
- 预测新的功能注释关系；
- 提高富集分析的鲁棒性。
PyTorch 环境为基因富集分析提供了更多的可能性，有助于深入挖掘基因组数据背后的生物学意义。
全部评论 (0)
还没有任何评论哟~