使用带有 apply() 函数的 pspearman 包计算相关矩阵

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)

导致错误：

Error in pspearman::spearman.test(as.numeric(x), as.numeric(y)) : (list) object cannot be coerced to type 'double'

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

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

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

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

谢谢！

Answer 1

考虑从 row.names 和 combn 中创建基因的成对组合，然后通过定义的函数遍历成对列表。确保 return 来自 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"

使用带有 apply() 函数的 pspearman 包计算相关矩阵

Calculating a correlation matrix with pspearman package with apply() function

r

apply

correlation

sapply

mapply