如何优化R代码快速计算两个大数据框的显著相关性（FDR<0.05）

需要找出两个大数据框间的显著相关性（FDR < 0.05），但当前代码中的cor.mtest(cor)运行耗时数小时，小数据集上运行正常，求代码优化方案。

原代码：

t.bval <- t(bval)
t.expr <- t(expr)
cor <- cor(t.bval, t.expr)
colnames(cor) <- rownames(bval)
mtest <- cor.mtest(cor)
adjusted.pval <- p.adjust(mtest$p, method = "fdr")
sig.cor <- cor * (adjusted.pval < 0.05)

数据集维度：

> dim(bval)
[1] 9844  174
> dim(expr)
[1] 9844  174

t(bval)前5行5列样例：

structure(c(0.651202886519379, 0.546932468275249, 0.77990670335116, 
0.570367878904668, 0.490859703406345, 0.819127800066518, 0.896778785060273, 
0.910419995766152, 0.862881750710962, 0.895738407255165, 0.409754021776346, 
0.437562413240227, 0.404271790048482, 0.407834330819103, 0.411767445153759, 
0.111173024032316, 0.0732396067281633, 0.892561968208864, 0.0924103940922992, 
0.0745261480840073, 0.727502771503466, 0.745126721066948, 0.72446467629277, 
0.740169892735761, 0.730258358871385), dim = c(5L, 5L), dimnames = list(
    c("TCGA.2K.A9WE.01", "TCGA.2Z.A9J1.01", "TCGA.2Z.A9J3.01", 
    "TCGA.2Z.A9J6.01", "TCGA.2Z.A9J7.01"), c("A1BG", "A2M", "A4GALT", 
    "AAAS", "AACS")))

t(expr)前5行5列样例：

structure(c(3.11327062985085, 2.60736439071998, 2.51138760850467, 
2.61757871821215, 2.99116809570192, 3.94027812075928, 3.7456302068644, 
3.48141956715802, 3.70450896566585, 3.58083202637968, 3.50763791284198, 
3.60784204549275, 3.45528021497798, 3.56318932299471, 3.53594747333992, 
3.47014111307111, 3.4842944131813, 3.44268580438446, 3.5016009268885, 
3.51798634683037, 3.38519953102602, 3.45448967582629, 3.42698278688777, 
3.44786587002075, 3.38755257891847), dim = c(5L, 5L), dimnames = list(
    c("TCGA.2K.A9WE.01", "TCGA.2Z.A9J1.01", "TCGA.2Z.A9J3.01", 
    "TCGA.2Z.A9J6.01", "TCGA.2Z.A9J7.01"), c("A1BG", "A2M", "A4GALT", 
    "AAAS", "AACS")))

核心问题是cor.mtest内部多采用循环计算，面对9844×9844的相关系数矩阵（约9700万次运算）时效率极低。可通过直接利用相关系数与p值的数学关系，结合R的向量化运算大幅提速，具体步骤如下：

1. 替代cor.mtest，向量化计算p值 #

Pearson相关系数对应的p值可通过t分布推导，公式为：
[ t = r \times \sqrt{\frac{n-2}{1-r^2}} ]
其中(r)为相关系数，(n)为样本量（此处为174），p值取双侧检验结果。代码实现：

t.bval <- t(bval)
t.expr <- t(expr)
cor_mat <- cor(t.bval, t.expr)
colnames(cor_mat) <- rownames(bval)

n <- nrow(t.bval)  # 样本量174
# 计算t统计量，加极小值避免r=±1时除以0
t_stat <- cor_mat * sqrt((n - 2) / (1 - cor_mat^2 + .Machine$double.eps))
# 计算双侧p值
p_val <- 2 * pt(abs(t_stat), df = n - 2, lower.tail = FALSE)

2. FDR校正与显著相关性筛选 #

后续步骤与原代码逻辑一致，注意保持矩阵维度匹配：

adjusted_pval <- p.adjust(p_val, method = "fdr")
# 将校正后的p值转为与cor_mat同维度的矩阵
adjusted_pval_mat <- matrix(adjusted_pval, nrow = nrow(cor_mat), ncol = ncol(cor_mat))
sig_cor <- cor_mat * (adjusted_pval_mat < 0.05)

3. 额外优化建议 #

• 若仅需保留显著相关的特征对，而非完整矩阵，可直接提取索引节省内存：

  sig_indices <- which(adjusted_pval_mat < 0.05, arr.ind = TRUE)
  sig_results <- data.frame(
    bval_feature = rownames(cor_mat)[sig_indices[,1]],
    expr_feature = colnames(cor_mat)[sig_indices[,2]],
    correlation = cor_mat[sig_indices],
    fdr_pval = adjusted_pval_mat[sig_indices]
  )
  

• 若内存紧张，可分块计算（将bval/expr拆分为小批次逐块运算），但向量化方法已足够处理当前规模数据。

内容的提问来源于stack exchange，提问作者Anon