2
\$\begingroup\$

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)
}
\$\endgroup\$
2
\$\begingroup\$

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))
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.