add Figure for SUP note 4 (simulation)

Figures_manuscript/scripts/Fig_SUPnote4_simulation_test_boundaries.R

+library(MASS)
+library(ggplot2)
+library(data.table)
+
+h1 = 0.4
+h2 = 0.4
+M = 1000
+
+pgm= 0.4 * min(h1,h1) # unused
+pgh = 0.5 * min(h1,h2)
+
+epsgh = matrix(c(h1, pgh, pgh, h2), ncol=2)/ M
+epsgoh = matrix(c(h1, -pgh, -pgh, h2), ncol=2)/ M
+
+# heritable trait but not correlated
+epsh = matrix(c(h1, 0, 0, h2), ncol=2)/M
+
+epshlow = matrix(c(0.1, 0, 0, 0.1), ncol=2)/M
+
+epsg1 = matrix(c(h1, pgh, pgh, h2), ncol=2)/ M
+epsg2 = matrix(c(h1, -0.2, -0.2, h2), ncol=2)/ M
+
+# draw genetic effect according to several scenarios
+BG_mat <-list()
+#set.seed(1234)
+#BG_mat <- append(BG_mat, list("G_cov_low_het_null_Gcor" = data.frame(mvrnorm(M, c(0,0),epshlow))))	# genetic effects with zero genetic correlation between traits (with low heritability)
+set.seed(1234)
+BG_mat <- append(BG_mat, list("G_cov_high"=  data.frame(mvrnorm(M, c(0,0),epsgh))))			# genetic effects with a given positive genetic covariance
+#set.seed(1234)
+#BG_mat <- append(BG_mat, list("G_cov_heritable_null_Gcor" = data.frame(mvrnorm(M, c(0,0),epsh))))	# genetic effects with zero genetic correlation between traits
+#set.seed(1234)
+#BG_mat <- append(BG_mat, list("G_cov_opp_high"=  data.frame(mvrnorm(M, c(0,0),epsgoh))))		# genetic effects with a given negative genetic covariance
+#set.seed(1234)
+#BG_mat <- append(BG_mat, list("G_bicov" = data.frame(rbind(mvrnorm(M/2, c(0,0),epsg1), mvrnorm(M/2, c(0,0), epsg2)))))	# genetic effects with a given mixed (local) genetic covariance
+
+
+sumZ <- function(Z, w){ return((as.matrix(Z) %*% matrix(w))^2 /(as.matrix(t(w))%*%sigr %*%as.matrix(w))[1,1])} # unused
+Univ <- function(Z){max(abs(Z))}
+
+hypothesis ="G_cov_high" # unused
+N= 10000 # unused
+overlap=  0.5 # unused
+sign_threshold = 5 *(10^-8) #/M
+C_env = seq(-0.9, 0.9, len=5) # Control the environmental correlation.
+Ns_frac = c(1, 0, 0.5) # Sample overlap
+N_samp = c(50000)
+
+
+ce=0.9 # unused
+
+for(hypothesis in names(BG_mat)){
+  print(hypothesis)
+
+  if(hypothesis=="null"){
+    h1 = 0
+    h2 = 0}else{
+      h1 = 0.4
+    h2 = 0.4
+  }
+
+  for(overlap in Ns_frac){
+    simu = 1
+    D_pval = data.frame(matrix(NA, nrow=length(N_samp)* length(C_env), ncol=4))
+    names(D_pval) = c("N_samp", "env_cor","Omnibus", "Univ")
+
+    for(ce in C_env){
+
+      pe = ce * min((1-h1), (1-h2))	
+      eps = matrix(c(1-h1, pe, pe, 1-h2), ncol=2)	# correlation in environmental and non-additive genetic effects between two traits
+
+      Beta = BG_mat[[hypothesis]]
+      # Genetic effects for the two trait
+      print(ce)
+      print(simu)
+
+      for(N in N_samp){
+        X1 = matrix(0, nrow=N, ncol=M)
+        X2 = matrix(0, nrow=N, ncol=M)
+
+        D_pval[simu,"N_samp"] = N
+        D_pval[simu,"env_cor"] = ce
+        # Draw SNPs :
+        MAF = runif(M,0.05, 0.95)
+        i= 1
+
+        for(p in MAF){
+          X1[,i] = rbinom(N, 2, p)
+          X2[,i] = rbinom(N, 2, p)
+          i = i +1
+        }
+
+        NS = dim(X1)[1] * overlap
+
+        X2[1:NS, ] = X1[1:NS,]
+
+        X1 = scale(X1)
+        X2 = scale(X2)
+
+	set.seed(4321)
+        env1 = mvrnorm(N, c(0,0), eps)
+        env2 = mvrnorm(N, c(0,0), eps)
+        env2[1:NS, ] = env1[1:NS,]
+
+        Y1 = as.matrix(X1) %*% Beta[,1] + env1[,1]
+        Y2 =as.matrix(X2) %*% Beta[,2] + env2[,2]
+
+
+        rho_ov= cor(Y1, Y2) * overlap # cor(Y1,Y2)=rho (rho=pe+pg); overlap=NS/N (fraction of overlapping samples out of all samples)
+
+        D = data.frame(matrix(c(Y1, Y2), ncol = 2))
+        png(paste0("/pasteur/zeus/projets/p02/GGS_JASS/5._ARTICLE_DISCOVERABILITY/Figures/outputs/simulations/Phenotype_", hypothesis,"_ov_",overlap,".png"))
+          p = ggplot(D, aes(x=Y1, y=Y2)) + geom_point()
+          print(p)
+        dev.off()
+
+        Beta_estimate1 = t((t(Y1) %*% X1) %*% solve(t(X1) %*% X1, tol=10^-25))
+        Beta_estimate2 = t((t(Y2) %*% X2) %*% solve(t(X2) %*% X2, tol=10^-25))
+
+        Beta_est = data.frame(B1 = Beta_estimate1, B2=Beta_estimate2)
+
+        png(paste0("/pasteur/zeus/projets/p02/GGS_JASS/5._ARTICLE_DISCOVERABILITY/Figures/outputs/simulations/Beta_", hypothesis,"_ov_",overlap,".png"))
+          print(ggplot(Beta_est, aes(x=B1, y=B2)) + geom_point())
+        dev.off()
+        Z = Beta_est
+
+        Z[,1] = Beta_est[,1] *sqrt(N)
+        Z[,2] = Beta_est[,2] *sqrt(N)
+        Z = as.data.frame(Z)
+        names(Z) = c("z1", "z2")
+
+        png(paste0("/pasteur/zeus/projets/p02/GGS_JASS/5._ARTICLE_DISCOVERABILITY/Figures/outputs/simulations/Z_", hypothesis,"_ov_",overlap,".png"))
+          print(ggplot(Z, aes(x=z1, y = z2)) + geom_point(alpha=0.5, size=1, color="blue"))
+        dev.off()
+        #
+        sigr = matrix(c(1,0,0,1),nrow=2,ncol=2)
+        sigr[1,2 ] = rho_ov
+        sigr[2, 1] = rho_ov
+        sigr
+        # Omnibus
+        Omni_stat= diag(as.matrix(Z[,1:2]) %*% solve(sigr) %*% as.matrix(t(Z[,1:2])))
+        pval = 1-pchisq(Omni_stat, df=2)
+
+        D_pval[simu, "Omnibus"]= mean(pval < sign_threshold)
+        Z["Omnibus"] = pval
+
+        # Stat Univariate test
+
+        stat = apply(Z[,1:2], 1, Univ)
+        pval = 1-pnorm(stat)
+        D_pval[simu,"Univ"] = mean(pval < sign_threshold)
+        Z["Univ"] = pval
+
+        Z["Significance status"] = "None"
+        Z[(Z$Univ < sign_threshold) & (Z$Omnibus < sign_threshold) , "Significance status"] = "Both"
+        Z[(Z$Univ < sign_threshold) & (Z$Omnibus > sign_threshold) , "Significance status"] = "Univariate"
+        Z[(Z$Univ > sign_threshold) & (Z$Omnibus < sign_threshold) , "Significance status"] = "Omnibus"
+
+	range_z <- max(abs(Z$z1),abs(Z$z2))
+
+        write.table(Z, paste0("/pasteur/zeus/projets/p02/GGS_JASS/5._ARTICLE_DISCOVERABILITY/Figures/outputs/simulations/Z_scores_", hypothesis, "_ov_", overlap,"_N_", N,"ce_",ce, ".csv"), sep="\t", row.names=TRUE)
+
+        png(paste0("/pasteur/zeus/projets/p02/GGS_JASS/5._ARTICLE_DISCOVERABILITY/Figures/outputs/simulations/", hypothesis, "_ov_", overlap,"_N_", N,"ce_",ce, ".png"), width=1500, height=1500, res=300)
+          p = ggplot(Z , id.var=c("z1", "z2"), aes(x=z1, y=z2, color=`Significance status`)) + geom_point(size=1.5) + scale_colour_manual(values = c("#3ba3ec","grey", "#f77189","#50b131")) + xlim(-range_z,range_z) + ylim(-range_z,range_z)
+           print(p + theme_minimal()+ theme(legend.position = "top", panel.spacing = unit(1.5, "lines"), text=element_text(size=16)) + labs(color="Significance Status"))
+        dev.off()
+          simu = simu+1
+      }
+    }
+      write.table(D_pval, paste0("/pasteur/zeus/projets/p02/GGS_JASS/5._ARTICLE_DISCOVERABILITY/Figures/outputs/simulations/", hypothesis, "_ov_", overlap,".csv"), sep="\t", row.names=TRUE)
+  }
+}