用R/Python计算同基因多探针成对皮尔逊相关系数及p值

以下分别提供R和Python的实现方案，两种方案均采用按基因分组计算的逻辑，避免生成27000×27000的超大相关矩阵，内存占用低、计算速度快，完全适配你的数据规模。

默认你已经准备好两类输入数据：

• 表达矩阵：行对应探针，列对应D1~D14共14个表达字段，行名为探针名称
• 探针注释表：包含至少两列，ProbeName（探针名称，和表达矩阵行名完全对应）、Gene（探针对应的基因名称）

需要依赖tidyverse和broom包，可通过install.packages(c("tidyverse", "broom"))安装。

# 加载依赖
library(tidyverse)
library(broom)

# 1. 把表达矩阵转为带探针名的数据框
expr_df <- as.data.frame(expr_mat) %>% 
  rownames_to_column("ProbeName")

# 2. 按基因分组计算探针对的相关系数与p值
result <- expr_df %>% 
  # 关联探针注释
  left_join(probe_anno, by = "ProbeName") %>% 
  group_by(Gene) %>% 
  # 过滤只有1个探针的基因，无需计算
  filter(n() >= 2) %>% 
  group_modify(~{
    # 生成当前基因下所有不重复的探针两两组合
    probe_pairs <- combn(.x$ProbeName, 2, simplify = FALSE)
    # 逐对计算相关系数
    map_dfr(probe_pairs, function(pair){
      x <- .x[.x$ProbeName == pair[1], 2:15] %>% as.numeric()
      y <- .x[.x$ProbeName == pair[2], 2:15] %>% as.numeric()
      cor_res <- cor.test(x, y, method = "pearson")
      tibble(
        ProbeName_1 = pair[1],
        ProbeName_2 = pair[2],
        PearsonCorrelationValue = cor_res$estimate,
        Pvalue = cor_res$p.value
      )
    })
  }) %>% 
  ungroup()

# 3. 导出结果为CSV
write_csv(result, "probe_cor_result.csv")

需要依赖pandas、scipy，可通过pip install pandas scipy安装（itertools为Python标准库无需额外安装）。

import pandas as pd
from scipy.stats import pearsonr
from itertools import combinations

# 1. 合并表达数据与探针注释
merged_df = expr_df.reset_index()\
  .rename(columns={'index':'ProbeName'})\
  .merge(probe_anno, on='ProbeName')

result_list = []
# 2. 按基因分组遍历计算
for gene, group in merged_df.groupby('Gene'):
    probe_list = group['ProbeName'].tolist()
    if len(probe_list) < 2:
        continue
    # 生成所有不重复的探针两两组合
    for probe1, probe2 in combinations(probe_list, 2):
        x = group[group['ProbeName'] == probe1].iloc[:, 1:-1].values.flatten()
        y = group[group['ProbeName'] == probe2].iloc[:, 1:-1].values.flatten()
        corr, pval = pearsonr(x, y)
        result_list.append({
            'ProbeName_1': probe1,
            'ProbeName_2': probe2,
            'Gene': gene,
            'PearsonCorrelationValue': corr,
            'Pvalue': pval
        })

# 3. 转为DataFrame并导出CSV
result_df = pd.DataFrame(result_list)
result_df.to_csv('probe_cor_result.csv', index=False)

注意事项 #

• 代码中提取表达值的列范围（R中的2:15、Python中的1:-1）需要根据你的实际输入数据结构调整
• 输出结果无自配对、无跨基因配对、无重复配对，完全符合要求的字段格式

内容的提问来源于stack exchange，提问作者Cp.Recker