#@ fisher's z transformation
#https://stats.stackexchange.com/questions/64152/are-two-pearson-correlation-coefficients-different
cor.diff.test = function(x1, x2, y1, y2, method="pearson") {
  cor1 = cor.test(x1, x2, method=method)
  cor2 = cor.test(y1, y2, method=method)

  r1 = cor1$estimate
  r2 = cor2$estimate
  n1 = sum(complete.cases(x1, x2))
  n2 = sum(complete.cases(y1, y2))
  fisher = ((0.5*log((1+r1)/(1-r1)))-(0.5*log((1+r2)/(1-r2))))/((1/(n1-3))+(1/(n2-3)))^0.5

  p.value = (2*(1-pnorm(abs(fisher))))

  result= list(
    "cor1" = list(
      "estimate" = as.numeric(cor1$estimate),
      "p.value" = cor1$p.value,
      "n" = n1
    ),
    "cor2" = list(
      "estimate" = as.numeric(cor2$estimate),
      "p.value" = cor2$p.value,
      "n" = n2
    ),
    "p.value.twosided" = as.numeric(p.value),
    "p.value.onesided" = as.numeric(p.value) / 2
  )
  cat(paste(sep="",
          "cor1: r=", format(result$cor1$estimate, digits=3), ", p=", format(result$cor1$p.value, digits=3), ", n=", result$cor1$n, "\n",
          "cor2: r=", format(result$cor2$estimate, digits=3), ", p=", format(result$cor2$p.value, digits=3), ", n=", result$cor2$n, "\n",
          "diffence: p(one-sided)=", format(result$p.value.onesided, digits=3), ", p(two-sided)=", format(result$p.value.twosided, digits=3), "\n"
  ))
  return(result);
}


c1=cor.test(log2(all_m2[,17]),log2(all_m2[,33]))
c2=cor.test(log2(all_m2$RPKM_KO_m2.x),log2(all_m2[,32]))
tt=paired.r(c1$estimate,c2$estimate,NULL,dim(all_m2)[1])
print (c(c1$estimate,c2$estimate,tt$z,tt$p))
0.6961544  0.6024532 17.4810979  0.0000000  # sig. higher cor in WT

cor.diff.test(log2(all_m2$RPKM_WT_m2.x),log2(all_m2$RPKM_WT_m2.y),log2(all_m2$RPKM_KO_m2.x),log2(all_m2$RPKM_KO_m2.y))
cor1: r=0.696, p=0, n=14386
cor2: r=0.602, p=0, n=14386
diffence: p(one-sided)=0, p(two-sided)=0   # sig. higher cor in WT

cor.diff.test(log2(all_m1$RPKM_WT_m1.x),log2(all_m1$RPKM_WT_m1.y),log2(all_m1$RPKM_KO_m1.x),log2(all_m1$RPKM_KO_m1.y))
cor1: r=0.64, p=0, n=14164
cor2: r=0.635, p=0, n=14164
diffence: p(one-sided)=0.242, p(two-sided)=0.484

cor.diff.test(log2(all_m0$RPKM_WT_0h.x),log2(all_m0$RPKM_WT_0h.y),log2(all_m0$RPKM_KO_0h.x),log2(all_m0$RPKM_KO_0h.y))
cor1: r=0.551, p=0, n=14788
cor2: r=0.664, p=0, n=14788
diffence: p(one-sided)=0, p(two-sided)=0 # sig. higher cor in KO


a=all_m0[(!is.na(all_m0[,32])&!is.na(all_m0[,16]))&(!is.na(all_m0[,33])&!is.na(all_m0[,17])),] # is.na is handlesd in cor.diff.test
cor.diff.test(log2(a$RPKM_WT_0h.x),log2(a$RPKM_WT_0h.y),log2(a$RPKM_KO_0h.x),log2(a$RPKM_KO_0h.y))
cor1: r=0.551, p=0, n=14788
cor2: r=0.664, p=0, n=14788
diffence: p(one-sided)=0, p(two-sided)=0

a=all_m0[(all_m0[,32]>4)&(all_m0[,16]>4)&(all_m0[,33]>4)&(all_m0[,17]>4),]
cor.diff.test(log2(a$RPKM_WT_0h.x),log2(a$RPKM_WT_0h.y),log2(a$RPKM_KO_0h.x),log2(a$RPKM_KO_0h.y))
cor1: r=0.524, p=0, n=6215
cor2: r=0.59, p=0, n=6215
diffence: p(one-sided)=5.43e-08, p(two-sided)=1.09e-07

a=all_m2[(all_m2[,32]>4)&(all_m2[,16]>4)&(all_m2[,33]>4)&(all_m2[,17]>4),]
cor.diff.test(log2(a$RPKM_WT_m2.x),log2(a$RPKM_WT_m2.y),log2(a$RPKM_KO_m2.x),log2(a$RPKM_KO_m2.y))
cor1: r=0.651, p=0, n=6808
cor2: r=0.567, p=0, n=6808
diffence: p(one-sided)=4.77e-15, p(two-sided)=9.55e-15

a=all_m2[(!is.na(all_m2[,32])&!is.na(all_m2[,16]))&(!is.na(all_m2[,33])&!is.na(all_m2[,17])),] # is.na is handlesd in cor.diff.test
cor.diff.test(log2(a$RPKM_WT_m2.x),log2(a$RPKM_WT_m2.y),log2(a$RPKM_KO_m2.x),log2(a$RPKM_KO_m2.y))
cor1: r=0.696, p=0, n=14386
cor2: r=0.602, p=0, n=14386
diffence: p(one-sided)=0, p(two-sided)=0

a=all_m1[(!is.na(all_m1[,32])&!is.na(all_m1[,16]))&(!is.na(all_m1[,33])&!is.na(all_m1[,17])),] # is.na is handlesd in cor.diff.test
cor.diff.test(log2(a$RPKM_WT_m1.x),log2(a$RPKM_WT_m1.y),log2(a$RPKM_KO_m1.x),log2(a$RPKM_KO_m1.y))
cor1: r=0.64, p=0, n=14164
cor2: r=0.635, p=0, n=14164
diffence: p(one-sided)=0.242, p(two-sided)=0.484

a=all_2h[(!is.na(all_2h[,32])&!is.na(all_2h[,16]))&(!is.na(all_2h[,33])&!is.na(all_2h[,17])),] # is.na is handlesd in cor.diff.test
cor.diff.test(log2(a$RPKM_WT_2h.x),log2(a$RPKM_WT_2h.y),log2(a$RPKM_KO_2h.x),log2(a$RPKM_KO_2h.y))
cor1: r=0.531, p=0, n=14630
cor2: r=0.6, p=0, n=14630
diffence: p(one-sided)=0, p(two-sided)=0

a=all_6h[(!is.na(all_6h[,32])&!is.na(all_6h[,16]))&(!is.na(all_6h[,33])&!is.na(all_6h[,17])),] # is.na is handlesd in cor.diff.test
cor.diff.test(log2(a$RPKM_WT_6h.x),log2(a$RPKM_WT_6h.y),log2(a$RPKM_KO_6h.x),log2(a$RPKM_KO_6h.y))
cor1: r=0.65, p=0, n=14621
cor2: r=0.594, p=0, n=14621
diffence: p(one-sided)=2.22e-15, p(two-sided)=4.44e-15