Quick estimation for empirical permutation and correlation test in R

Question

I have two matrices (a, b), each including 3 genes and 6 observations. I have estimated the Spearman correlation for gene pairs (rcorr.observed). Then I want to do an empirical permutation (1000 times) and estimate the correlation for each time, then achieve a Pval_estim file. Since I have a very huge datasets (23000*23000), the script is very slow. I wonder if you could help me to speed it up. I have also used some code from (+).

a<-structure(c(7.284640193, 8.386102403, 10.27187091, 6.56612737, 
8.969982518, 10.03221978, 7.193522121, 7.358283395, 10.92117745, 
7.162802801, 8.297228578, 11.36980277, 7.865611714, 7.999185693, 
10.33028086, 6.831671275, 8.953536984, 8.826461297), .Dim = c(3L, 
6L), .Dimnames = list(c("DDR1_MIR4640", "RFC2", "HSPA6"), c("a_1", "a_2", "a_3",     "a_4", "a_5", "a_6")))


b<-structure(c(9.1048886, 6.4114527, 5.7281808, 9.4302985, 6.2576226, 
4.871274, 9.187927, 5.9036324, 4.3635891, 8.6896685, 6.6680496, 
5.5445622, 8.9274641, 6.4394849, 5.5175364, 7.8629553, 8.4304969, 
6.1402062), .Dim = c(3L, 6L), .Dimnames = list(c("FAM174B", "SV2B", 
"RBPMS2"), c("b_1", "b_2", "b_3", "b_4", "b_5", "b_6")))

rcorr.observed<- structure(c(0.00309597514569759, -0.0526315793395042,                 -0.108359135687351, 
NA, -0.0567595474421978, -0.0299277603626251, NA, NA, -0.116615064442158
), .Dim = c(3L, 3L))


library(Hmisc)
Pval_estim<-matrix(nrow = 3, ncol = 3)
number_of_permutations = 1000
diff.random = NULL
for(i in 1:nrow(a))
  {
  for(j in 1:nrow(b))
  {
    if(i>=j)#selecting only lower part of the symmetric matrix
      {
    combined<-as.numeric(t(c(a[i,], b[j,])))
for (k in 1 : number_of_permutations)
  {
    shuffled = sample (combined, length(combined))
    a.random = shuffled[1 : 6]
    b.random = shuffled[7 : length(combined)]
    diff.random[k] = rcorr(a.random, b.random, type="spearman")$r[1,2] 
}
Pval_estim[i,j] = sum(abs(diff.random) >= abs(rcorr.observed[i,j])) /             number_of_permutations
}}
print(i)
}

Answer 1

Using apply and replicate, you can get a slight speedup (below). Setting a seed, you can check that it gives the same result as your version. In my benchmark test, it's about 10% faster. For further speedup, you could consider using clusterApply and parLapply from the parallel package.

# combining matrices into one to use apply
combinedMat <- cbind(a[rep(1:nrow(a), each=3), ], b[rep(1:nrow(b), times=3), ])
combinedMat <- cbind(combinedMat, obs = as.numeric(t(rcorr.observed)))
# using apply instead of loop
Pval_estim <- t(matrix(apply(combinedMat, 1, function(combined) {
  obs <- combined['obs']
  if (is.na(obs)) return(NA)
  combined <- as.numeric(combined[names(combined) != 'obs'])
  # using replicate instead of loop
  diff.random <- replicate(number_of_permutations,
  {
    shuffled = sample(combined)
    # you only use the correlation, so there is no need to use rcorr. 
    diff.random = cor(shuffled[1 : 6], 
                      shuffled[7 : length(combined)], 
                      method="spearman")
  })
  sum(abs(diff.random) >= abs(obs)) / number_of_permutations
}), nrow = 3, ncol = 3))