使用 pspearman 包和 apply() 函数计算相关矩阵

Question

我正在尝试计算数据框的 Spearman 相关性和 p 值。为了更好地逼近 p 值，我必须坚持使用 pspearman 包。我期待与rcorr() 函数相似的结果。但是我在逐行执行pspearman:spearman.test() 时遇到问题。

我的数据框包含 5000 行（基因）和 200 列（点）。我想得到这些 5000*5000 基因-基因对的相关矩阵和 p 值矩阵。仅当两个基因都不是多于两个位点的 NA 时才计算相关性。

我可以通过循环来实现这一点，但对于我的大数据集来说太慢了。当我尝试使用apply(),sapply(),mapply() 来提高速度时遇到问题。

这是我尝试过的：

data = data.frame(matrix(rbinom(10*100000, 50, .5), ncol=200))
dim(data) #5000, 200
rownames(data) <- paste("gene", 1:5000, sep="") 
colnames(data) <- paste("spot",1:200,sep='')

library(pspearman)
spearFunc = function(x,y=data) {
  df = rbind(x,y)
  # Check the number of complete spots.There are no NAs in this set.
  complete = sum(!(is.na(x)) & !(is.na(y)))
  if (complete >=2 ) {
    pspearman::spearman.test(as.numeric(x),as.numeric(y))
    # This function returns a list containing 8 values, like pvalue,correlation
    }}

pair.all1 = mapply(spearFunc,data,data)
dim(pair.all1)
# 8 200, 200 is the number of columns 
pair.all2 = apply(data,1,spearFunc) 

这会导致错误：

pspearman::spearman.test(as.numeric(x), as.numeric(y)) 中的错误： (list) 对象不能被强制输入'double'

我希望用spearman.test对每个基因对用apply()来做

spearman.test(data[gene1],data[gene1]) 
spearman.test(data[gene1],data[gene2])
....
spearman.test(data[gene1],data[gene5000])
...
spearman.test(data[gene5000],data[gene5000])

它应该返回一个 8 行 25,000,000 列（5000*5000 个基因对）的数据框。

是否可以在 apply() 中使用 apply() 来实现我的目的？

谢谢！

Answer 1

考虑从row.names 和combn 创建基因的配对组合，然后通过定义的函数遍历配对列表。请务必从 if 逻辑返回 NA 结构，以避免在矩阵输出中出现 NULL。

但是，请注意，5,000 个基因 (choose(5000, 2)) 的成对排列结果非常高，达到 12,497,500 个元素！因此，sapply（循环本身）在性能上可能与for 没有太大区别。研究并行化迭代。

gene_combns <- combn(row.names(data), 2, simplify = FALSE)

spear_func <- function(x) {
  # EXTRACT ROWS BY ROW NAMES  
  row1 <- as.numeric(data[x[1],])
  row2 <- as.numeric(data[x[2],]) 

  # Check the number of complete spots.There are no NAs in this set.
  complete = sum(!(is.na(x)) & !(is.na(y)))

  if (complete >=2 ) {
    pspearman::spearman.test(row1, row2)        
  } else {
    c(statistic=NA, parameter=NA, p.value=NA, estimate=NA, 
      null.value=NA, alternative=NA,   method=NA, data.name=NA)
  }
}

pair.all2 <- sapply(gene_combns, spear_func)

---

测试

上面已经用cor.test 进行了测试（与spearman.test 完全相同但更准确的p-value 的输入参数和输出列表）使用小数据集样本（50 obs，20 vars）：

set.seed(82418)
data <- data.frame(matrix(rbinom(10*100000, 50, .5), ncol=200))[1:50, 1:20]
rownames(data) <- paste0("gene", 1:50) 
colnames(data) <- paste0("spot", 1:20)

gene_combns <- combn(row.names(data), 2, simplify = FALSE)
# [[1]]
# [1] "gene1" "gene2"    
# [[2]]
# [1] "gene1" "gene3"    
# [[3]]
# [1] "gene1" "gene4"    
# [[4]]
# [1] "gene1" "gene5"    
# [[5]]
# [1] "gene1" "gene6"    
# [[6]]
# [1] "gene1" "gene7"

test <- sapply(gene_combns, spear_func)  # SAME FUNC BUT WITH cor.test
test[,1:5]

#             [,1]                              [,2]                             
# statistic   885.1386                          1659.598                         
# parameter   NULL                              NULL                             
# p.value     0.1494607                         0.2921304                        
# estimate    0.3344823                         -0.2478179                       
# null.value  0                                 0                                
# alternative "two.sided"                       "two.sided"                      
# method      "Spearman's rank correlation rho" "Spearman's rank correlation rho"
# data.name   "row1 and row2"                   "row1 and row2"                  
#             [,3]                              [,4]                             
# statistic   1554.533                          1212.988                         
# parameter   NULL                              NULL                             
# p.value     0.4767667                         0.7122505                        
# estimate    -0.1688217                        0.08797877                       
# null.value  0                                 0                                
# alternative "two.sided"                       "two.sided"                      
# method      "Spearman's rank correlation rho" "Spearman's rank correlation rho"
# data.name   "row1 and row2"                   "row1 and row2"                  
#             [,5]                             
# statistic   1421.707                         
# parameter   NULL                             
# p.value     0.7726922                        
# estimate    -0.06895299                      
# null.value  0                                
# alternative "two.sided"                      
# method      "Spearman's rank correlation rho"
# data.name   "row1 and row2"