Cor.test inside for loop is very slow for large matrix processing
If you had profiled your code, you would have seen that unusually long time is taken by line:
corr_list[i - 1,] &lt;- c(rNames[1], rNames[i], r$estimate, r$p.value)
It is probably because you are combining multiple type values into one vector.
So this should be mush faster:
f3 &lt;- function(dat) {
  require(data.table)
  mainVar &lt;- dat[1,]
  rNames &lt;- rownames(dat)
  rez &lt;- lapply(2:nrow(dat), function(i) {
    r &lt;- cor.test(x = mainVar, y = dat[i,])
    c(r$estimate, r$p.value)
  })
  fin &lt;- data.table(rNames[1], rNames[-1])
  fin &lt;- cbind(fin, data.table(do.call(rbind, rez)))
  setnames(fin, c(&quot;top&quot;, &quot;correlated&quot;, &quot;cor&quot;, &quot;p.value&quot;))
  return(rez)
}

I &lt;- 60
N &lt;- 1000
dat &lt;- MASS::mvrnorm(N, mu = rep(0, I), diag(I))
rownames(dat) &lt;- paste0(&quot;G&quot;, 1:N)

# Comparison  
res1 &lt;- microbenchmark(
  f1(dat),
  f2(dat),
  f3(dat),
  times = 10
)
print(res1, unit = &quot;s&quot;)
# Unit: seconds
# expr       min        lq      mean    median        uq       max neval cld
# f1(dat) 0.4333447 0.4609346 0.4880110 0.4863479 0.5050751 0.5782526    20   c
# f2(dat) 0.3585121 0.3729532 0.4021378 0.3979831 0.4269367 0.4899676    20  b 
# f3(dat) 0.1483598 0.1610634 0.1782496 0.1787396 0.1934256 0.2273650    20 a  

autoplot(res1)


Update
If we calculate correlation with base function and p.value - ourselves, we can get rid of for loop and do all of this much faster. Only problems my arise if your data contains missing values:
pcor &lt;- function(r, n){
  t &lt;- r * sqrt(n - 2) / sqrt(1 - r^2)
  p &lt;- 2 * (1 - pt(abs(t), (n - 2)))
  p[p &gt; 1] &lt;- 1
  p
}

f4 &lt;- function(dat){
  require(data.table)
  
  rNames &lt;- rownames(dat)
  d2 &lt;- t(dat)
  cors &lt;- cor(d2[, 1], d2[, -1])
  cors &lt;- t(cors)
  cp &lt;- pcor(cors, ncol(dat))
  cp &lt;- cbind(cors, cp)
  fin &lt;- data.table(rNames[1], rNames[-1])
  fin &lt;- cbind(fin, cp)
  setnames(fin, c(&quot;top&quot;, &quot;correlated&quot;, &quot;cor&quot;, &quot;p.value&quot;))
  return(fin)
}

I &lt;- 60
N &lt;- 1000
dat &lt;- MASS::mvrnorm(N, mu = rep(0, I), diag(I))
rownames(dat) &lt;- paste0(&quot;G&quot;, 1:N)


res1 &lt;- microbenchmark(
  f1(dat),
  f2(dat),
  f3(dat),
  f4(dat),
  times = 10
)
print(res1, unit = &quot;s&quot;)
# Unit: seconds
# expr         min          lq        mean      median          uq         max neval  cld
# f1(dat) 0.408167816 0.410614574 0.428363371 0.429430744 0.442042920 0.456887399    10    d
# f2(dat) 0.334446197 0.345779175 0.366422216 0.362465218 0.378226388 0.426838502    10   c 
# f3(dat) 0.139088268 0.145095289 0.156048277 0.153321298 0.162591518 0.183544069    10  b  
# f4(dat) 0.002363351 0.002428473 0.002615854 0.002483812 0.002797716 0.003117556    10 a