一个基因的GSEA分析代码

以下是一个简单的基因的GSEA分析代码示例，供参考：

import pandas as pd
import numpy as np
from scipy.stats import norm

# 读取基因表达数据
gene_expression_data = pd.read_csv("gene_expression_data.csv")

# 读取基因集合数据
gene_set_data = pd.read_csv("gene_set_data.csv")

# 计算每个基因的表达值的均值
gene_expression_mean = gene_expression_data.mean(axis=1)

# 基于基因表达均值计算基因的排名
gene_ranks = gene_expression_mean.rank(method='min')

# 初始化GSEA参数
n_genes = len(gene_ranks)
n_perms = 1000
gene_set_size_min = 15
gene_set_size_max = 500
alpha = 0.05

# 计算基因集合的ES值
def calculate_es(gene_set):
    gene_ranks_in_set = gene_ranks[gene_ranks.index.isin(gene_set)]
    gene_set_size = len(gene_ranks_in_set)
    if gene_set_size < gene_set_size_min:
        return np.nan
    else:
        es = 0
        running_sum = 0
        for gene_rank in gene_ranks:
            if gene_rank in gene_ranks_in_set:
                running_sum += 1
            else:
                running_sum -= 1
            es = max(es, running_sum)
        return es / gene_set_size

# 计算基因集合的NES值
def calculate_nes(gene_set, gene_set_es, gene_set_size):
    gene_set_mean = gene_expression_mean[gene_expression_mean.index.isin(gene_set)]
    gene_set_sd = np.std(gene_set_mean, ddof=1)
    if gene_set_sd == 0:
        return np.nan
    else:
        nes = gene_set_es / gene_set_sd
        return nes

# 计算随机排列基因集合的ES值
def calculate_permutation_es(gene_set):
    permuted_ranks = gene_ranks.sample(frac=1).reset_index(drop=True)
    permuted_es = calculate_es(gene_set, permuted_ranks)
    return permuted_es

# 计算基因集合的p值和FDR校正的p值
def calculate_p_value(gene_set_es, permutation_es):
    p_value = (np.sum(permutation_es >= gene_set_es) + 1) / (len(permutation_es) + 1)
    fdr_p_value = p_value * n_perms
    return p_value, fdr_p_value

# 进行GSEA分析
results = []
for index, row in gene_set_data.iterrows():
    gene_set = row['gene_set'].split(',')
    gene_set_es = calculate_es(gene_set)
    if np.isnan(gene_set_es):
        continue
    gene_set_size = len(gene_set)
    if gene_set_size < gene_set_size_min or gene_set_size > gene_set_size_max:
        continue
    permutation_es = [calculate_permutation_es(gene_set) for i in range(n_perms)]
    p_value, fdr_p_value = calculate_p_value(gene_set_es, permutation_es)
    if fdr_p_value <= alpha:
        nes = calculate_nes(gene_set, gene_set_es, gene_set_size)
        results.append({
            'gene_set': row['gene_set'],
            'es': gene_set_es,
            'nes': nes,
            'p_value': p_value,
            'fdr_p_value': fdr_p_value,
            'gene_set_size': gene_set_size
        })
        
# 输出结果
results_df = pd.DataFrame(results)
results_df = results_df.sort_values(by=['fdr_p_value'], ascending=True)
print(results_df)

这个代码是一个简单的基因的GSEA分析代码示例，具体实现过程如下：

1. 读取基因表达数据和基因集合数据，计算每个基因的表达值的均值。
2. 基于基因表达均值计算每个基因的排名。
3. 初始化GSEA参数，包括基因集合的大小范围、置信水平等。
4. 定义计算基因集合的ES值、NES值、随机排列基因集合的ES值、p值和FDR校正的p值的函数。
5. 对每个基因集合进行GSEA分析，计算其ES值、NES值、p值和FDR校正的p值，筛选出FDR校正的p值小于等于置信水平的基因集合，并将结果保存在结果列表中。
6. 输出结果，按照FDR校正的p值升序排列。

需要注意的是，这个代码只是一个简单的示例，实际应用中可能需要根据具体情况进行修改和优化。