diff --git a/results/simulation/simu_mscaviar_inflated_small.R b/results/simulation/simu_mscaviar_inflated_small.R new file mode 100644 index 0000000..f94c417 --- /dev/null +++ b/results/simulation/simu_mscaviar_inflated_small.R @@ -0,0 +1,103 @@ +#!/usr/bin/env Rscript +args = commandArgs(trailingOnly=TRUE) +########################################################### pop1 ########################################################### +library(mvtnorm) +library(data.table) +set.seed(1) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +# system(paste0("export MKL_NUM_THREADS=",nthreads)) +# system(paste0("export NUMEXPR_NUM_THREADS=",nthreads)) +# system(paste0("export OMP_NUM_THREADS=",nthreads)) +# export MKL_NUM_THREADS=10 +# export NUMEXPR_NUM_THREADS=10 +# export OMP_NUM_THREADS=10 + +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- as.numeric(args[1]) #5000 +n2 <- as.numeric(args[2])#20000 +K_true <- as.numeric(args[3]) +K <- K_true + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- as.numeric(args[4])#0.2 +pc_var2 <- as.numeric(args[5])#0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +LD1_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1") +LD2_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2") +LD_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_LD_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +write.table(rbind(LD1_file, LD2_file), + file = LD_input_file, + col.names = F, row.names = F, sep = " ", quote = F) + +ofile <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/output/mscaviar_output_inflate_pc1_", pc_var1,"_",pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z1_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs1_inflate_pc1_", pc_var1, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z2_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs2_inflate_pc1_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Z_input_inflate_pc1_", pc_var1,"_",pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) + + +for (i in 1:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_inflate_pc1_", pc_var1, "_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_inflate_pc1_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + zs1 <- z_scores1[loci$idx_left[i_loci]:loci$idx_right[i_loci],][idx_include,] + zs2 <- z_scores2[loci$idx_left[i_loci]:loci$idx_right[i_loci],][idx_include,] + + write.table(zs1, + file = z1_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(zs2, + file = z2_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(rbind(z1_file[i], z2_file[i]), + file = z_input_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + + mscaviar <- paste0("/home/mcaiad/MsCAVIAR/MsCAVIAR -l ", LD_input_file, " -z ", z_input_file[i], " -n ", n1, ",", n2, " -c ",K_true," -o ", ofile[i]) + system(mscaviar) + cat(i, "-th rep finished.\n") +} + diff --git a/results/simulation/simu_mscaviar_inflated_small_corrected.R b/results/simulation/simu_mscaviar_inflated_small_corrected.R new file mode 100644 index 0000000..73e2bcd --- /dev/null +++ b/results/simulation/simu_mscaviar_inflated_small_corrected.R @@ -0,0 +1,103 @@ +#!/usr/bin/env Rscript +args = commandArgs(trailingOnly=TRUE) +########################################################### pop1 ########################################################### +library(mvtnorm) +library(data.table) +set.seed(1) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +# system(paste0("export MKL_NUM_THREADS=",nthreads)) +# system(paste0("export NUMEXPR_NUM_THREADS=",nthreads)) +# system(paste0("export OMP_NUM_THREADS=",nthreads)) +# export MKL_NUM_THREADS=10 +# export NUMEXPR_NUM_THREADS=10 +# export OMP_NUM_THREADS=10 + +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- as.numeric(args[1]) #5000 +n2 <- as.numeric(args[2])#20000 +K_true <- as.numeric(args[3]) +K <- K_true + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- as.numeric(args[4])#0.2 +pc_var2 <- as.numeric(args[5])#0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +LD1_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1") +LD2_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2") +LD_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_LD_input_inflate_pc1_", pc_var1, "_", pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +write.table(rbind(LD1_file, LD2_file), + file = LD_input_file, + col.names = F, row.names = F, sep = " ", quote = F) + +ofile <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/output/mscaviar_output_inflate_pc1_", pc_var1,"_",pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z1_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs1_inflate_pc1_", pc_var1, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z2_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs2_inflate_pc1_", pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Z_input_inflate_pc1_", pc_var1,"_",pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) + + +for (i in 1:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_inflate_pc1_", pc_var1, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_inflate_pc1_", pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + zs1 <- z_scores1[loci$idx_left[i_loci]:loci$idx_right[i_loci],][idx_include,] + zs2 <- z_scores2[loci$idx_left[i_loci]:loci$idx_right[i_loci],][idx_include,] + + write.table(zs1, + file = z1_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(zs2, + file = z2_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(rbind(z1_file[i], z2_file[i]), + file = z_input_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + + mscaviar <- paste0("/home/mcaiad/MsCAVIAR/MsCAVIAR -l ", LD_input_file, " -z ", z_input_file[i], " -n ", n1, ",", n2, " -c ",K_true," -o ", ofile[i]) + system(mscaviar) + cat(i, "-th rep finished.\n") +} + diff --git a/results/simulation/simu_mscaviar_small_chr22.R b/results/simulation/simu_mscaviar_small_chr22.R new file mode 100644 index 0000000..834a9c0 --- /dev/null +++ b/results/simulation/simu_mscaviar_small_chr22.R @@ -0,0 +1,98 @@ +#!/usr/bin/env Rscript +args = commandArgs(trailingOnly=TRUE) +########################################################### pop1 ########################################################### +library(mvtnorm) +library(data.table) +set.seed(1) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +# system(paste0("export MKL_NUM_THREADS=",nthreads)) +# system(paste0("export NUMEXPR_NUM_THREADS=",nthreads)) +# system(paste0("export OMP_NUM_THREADS=",nthreads)) +# export MKL_NUM_THREADS=10 +# export NUMEXPR_NUM_THREADS=10 +# export OMP_NUM_THREADS=10 + +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- as.numeric(args[1]) #5000 +n2 <- 20000 +K_true <- as.numeric(args[2]) +K <- K_true + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +LD1_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1") +LD2_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2") +LD_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_wholeCHR_LD_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +write.table(rbind(LD1_file, LD2_file), + file = LD_input_file, + col.names = F, row.names = F, sep = " ", quote = F) + +ofile <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/output/mscaviar_output_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z1_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs1_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z2_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs2_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Z_input_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) + + +for (i in 2:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + + write.table(data.frame(SNP = z_scores1[[1]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include], Z_afr = z_scores1[[2]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), + file = z1_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(data.frame(SNP = z_scores2[[1]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include], Z_afr = z_scores2[[2]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), + file = z2_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(rbind(z1_file[i], z2_file[i]), + file = z_input_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + + mscaviar <- paste0("/home/mcaiad/MsCAVIAR/MsCAVIAR -l ", LD_input_file, " -z ", z_input_file[i], " -n ", n1, ",", n2, " -c ",K_true," -o ", ofile[i]) + system(mscaviar) + cat(i, "-th rep finished.\n") +} + diff --git a/results/simulation/simu_mscaviar_small_nobg.R b/results/simulation/simu_mscaviar_small_nobg.R new file mode 100644 index 0000000..58f0a21 --- /dev/null +++ b/results/simulation/simu_mscaviar_small_nobg.R @@ -0,0 +1,99 @@ +#!/usr/bin/env Rscript +args = commandArgs(trailingOnly=TRUE) +########################################################### pop1 ########################################################### +library(mvtnorm) +library(data.table) +set.seed(1) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +# system(paste0("export MKL_NUM_THREADS=",nthreads)) +# system(paste0("export NUMEXPR_NUM_THREADS=",nthreads)) +# system(paste0("export OMP_NUM_THREADS=",nthreads)) +# export MKL_NUM_THREADS=10 +# export NUMEXPR_NUM_THREADS=10 +# export OMP_NUM_THREADS=10 + +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- as.numeric(args[1]) #5000 +n2 <- 20000 +K_true <- as.numeric(args[2]) +K <- K_true + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +enrich2 <- enrich * 10 +rho_sigma <- 0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +LD1_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1") +LD2_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2") +LD_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_LD_input_nobg_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep) +write.table(rbind(LD1_file, LD2_file), + file = LD_input_file, + col.names = F, row.names = F, sep = " ", quote = F) + +ofile <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/output/mscaviar_output_nobg_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z1_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs1_nobg_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z2_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs2_nobg_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Z_input_nobg_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) + + +for (i in 1:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_nobg_loci_",i_loci,"_small_p",p,"_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_nobg_loci_",i_loci,"_small_p",p,"_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + + write.table(data.frame(SNP = z_scores1[[1]], Z_afr = z_scores1[[2]]), + file = z1_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(data.frame(SNP = z_scores2[[1]], Z_afr = z_scores2[[2]]), + file = z2_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(rbind(z1_file[i], z2_file[i]), + file = z_input_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + + mscaviar <- paste0("/home/mcaiad/MsCAVIAR/MsCAVIAR -l ", LD_input_file, " -z ", z_input_file[i], " -n ", n1, ",", n2, " -c ",K_true," -o ", ofile[i]) + system(mscaviar) + cat(i, "-th rep finished.\n") +} + diff --git a/results/simulation/simu_mscaviar_small_t_chr22.R b/results/simulation/simu_mscaviar_small_t_chr22.R new file mode 100644 index 0000000..7c8a3e4 --- /dev/null +++ b/results/simulation/simu_mscaviar_small_t_chr22.R @@ -0,0 +1,100 @@ +#!/usr/bin/env Rscript +args = commandArgs(trailingOnly=TRUE) +########################################################### pop1 ########################################################### +library(mvtnorm) +library(data.table) +set.seed(1) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +# system(paste0("export MKL_NUM_THREADS=",nthreads)) +# system(paste0("export NUMEXPR_NUM_THREADS=",nthreads)) +# system(paste0("export OMP_NUM_THREADS=",nthreads)) +# export MKL_NUM_THREADS=10 +# export NUMEXPR_NUM_THREADS=10 +# export OMP_NUM_THREADS=10 + +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- 20000#as.numeric(args[1]) #5000 +n2 <- 20000 +K_true <- 3#as.numeric(args[2]) +K <- K_true + + + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0#0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +LD1_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_t",df,"_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1") +LD2_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_t",df,"_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2") +LD_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_t",df,"_wholeCHR_LD_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +write.table(rbind(LD1_file, LD2_file), + file = LD_input_file, + col.names = F, row.names = F, sep = " ", quote = F) + +ofile <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/output/mscaviar_t",df,"_output_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z1_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs1_t",df,"_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z2_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Zs2_t",df,"_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) +z_input_file <- paste0("/home/share/mingxuan/fine_mapping/simulation/data/mscaviar/input/mscaviar_Z_t",df,"_input_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", 1:nrep) + + +for (i in 2:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_t",df,"_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_t",df,"_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + + write.table(data.frame(SNP = z_scores1[[1]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include], Z_afr = z_scores1[[2]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), + file = z1_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(data.frame(SNP = z_scores2[[1]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include], Z_afr = z_scores2[[2]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), + file = z2_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + write.table(rbind(z1_file[i], z2_file[i]), + file = z_input_file[i], + col.names = F, row.names = F, sep = " ", quote = F) + + mscaviar <- paste0("/home/mcaiad/MsCAVIAR/MsCAVIAR -l ", LD_input_file, " -z ", z_input_file[i], " -n ", n1, ",", n2, " -c ",K_true," -o ", ofile[i]) + system(mscaviar) + cat(i, "-th rep finished.\n") +} + diff --git a/results/simulation/simu_paintor_inflated_small.R b/results/simulation/simu_paintor_inflated_small.R new file mode 100644 index 0000000..7d222c5 --- /dev/null +++ b/results/simulation/simu_paintor_inflated_small.R @@ -0,0 +1,100 @@ +#!/usr/bin/env Rscript +args = commandArgs(trailingOnly = TRUE) +########################################################### pop1 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +blas_set_num_threads(10) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- as.numeric(args[1]) #5000 +n2 <- as.numeric(args[2])#20000 +K_true <- as.numeric(args[3]) +K <- K_true + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- as.numeric(args[4]) #0.2 +pc_var2 <- as.numeric(args[5])#0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +R1 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt"))) +fwrite(R1[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1"), + sep = " ", col.names = F) +R2 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt"))) +fwrite(R2[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2"), + sep = " ", col.names = F) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1")) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2")) + +annot <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".annotations") +dummy <- data.table(dummy = rep(1, p)) +fwrite(dummy, file = annot, sep = " ", col.names = TRUE) + +input_folder <- "/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/" +zfile <- paste0("paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +input_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/files_apaintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +for (i in 1:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_inflate_pc1_", pc_var1, "_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_inflate_pc1_", pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + zs1 <- z_scores1[loci$idx_left[i_loci]:loci$idx_right[i_loci],][idx_include,] + zs2 <- z_scores2[loci$idx_left[i_loci]:loci$idx_right[i_loci],][idx_include,] + + z_paintor <- data.table(RSID = zs1[[1]], ZSCORE.P1 = zs1[[2]], ZSCORE.P2 = zs2[[2]]) + fwrite(z_paintor, paste0(input_folder, zfile), sep = " ", col.names = TRUE) + + + write.table(zfile, file = input_file, sep = " ", col.names = FALSE, row.names = F, quote = F) + + paintor <- paste0("/home/mcaiad/PAINTOR_V3.0/PAINTOR -input ", input_file, " -in ", input_folder, " -out /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/output/ -Zhead ZSCORE.P1,ZSCORE.P2 -LDname LD1,LD2 -enumerate ", K, " -annotations dummy -RESname ", i, "_exact_results") + system(paintor) + cat(i, "-th rep finished.\n") +} + diff --git a/results/simulation/simu_paintor_inflated_small_corrected.R b/results/simulation/simu_paintor_inflated_small_corrected.R new file mode 100644 index 0000000..65db4f0 --- /dev/null +++ b/results/simulation/simu_paintor_inflated_small_corrected.R @@ -0,0 +1,99 @@ +#!/usr/bin/env Rscript +args = commandArgs(trailingOnly = TRUE) +########################################################### pop1 ########################################################### +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +blas_set_num_threads(10) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- as.numeric(args[1]) #5000 +n2 <- as.numeric(args[2])#20000 +K_true <- as.numeric(args[3]) +K <- K_true + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- as.numeric(args[4]) #0.2 +pc_var2 <- as.numeric(args[5])#0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +R1 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt"))) +fwrite(R1[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1"), + sep = " ", col.names = F) +R2 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt"))) +fwrite(R2[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2"), + sep = " ", col.names = F) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1")) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2")) + +annot <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".annotations") +dummy <- data.table(dummy = rep(1, p)) +fwrite(dummy, file = annot, sep = " ", col.names = TRUE) + +input_folder <- "/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/" +zfile <- paste0("paintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +input_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/files_apaintor_mcmc_input_inflate_pc1_", pc_var1, "_", pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +for (i in 1:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_inflate_pc1_", pc_var1, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_inflate_pc1_", pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + zs1 <- z_scores1[loci$idx_left[i_loci]:loci$idx_right[i_loci],][idx_include,] + zs2 <- z_scores2[loci$idx_left[i_loci]:loci$idx_right[i_loci],][idx_include,] + + z_paintor <- data.table(RSID = zs1[[1]], ZSCORE.P1 = zs1[[2]], ZSCORE.P2 = zs2[[2]]) + fwrite(z_paintor, paste0(input_folder, zfile), sep = " ", col.names = TRUE) + + + write.table(zfile, file = input_file, sep = " ", col.names = FALSE, row.names = F, quote = F) + + paintor <- paste0("/home/mcaiad/PAINTOR_V3.0/PAINTOR -input ", input_file, " -in ", input_folder, " -out /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/output/ -Zhead ZSCORE.P1,ZSCORE.P2 -LDname LD1,LD2 -enumerate ", K, " -annotations dummy -RESname ", i, "_exact_results") + system(paintor) + cat(i, "-th rep finished.\n") +} + diff --git a/results/simulation/simu_paintor_small_chr22.R b/results/simulation/simu_paintor_small_chr22.R new file mode 100644 index 0000000..a1ee5af --- /dev/null +++ b/results/simulation/simu_paintor_small_chr22.R @@ -0,0 +1,94 @@ +########################################################### pop1 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +threads <- 20 +blas_set_num_threads(threads) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- 5000 +n2 <- 20000 +K_true <- 2 +K <- K_true + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +R1 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt"))) +fwrite(R1[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_wholeCHR_input_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1"), + sep = " ", col.names = F) +R2 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt"))) +fwrite(R2[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_wholeCHR_input_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2"), + sep = " ", col.names = F) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1")) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2")) + +annot <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_wholeCHR_input_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".annotations") +dummy <- data.table(dummy = rep(1, p)) +fwrite(dummy, file = annot, sep = " ", col.names = TRUE) + +input_folder <- "/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/" +zfile <- paste0("paintor_mcmc_wholeCHR_input_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +input_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/files_paintor_mcmc_wholeCHR_input_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +for (i in 1:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_wholeCHR_loci_",i_loci,"_small_p",p,"_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_wholeCHR_loci_",i_loci,"_small_p",p,"_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + + z_paintor <- data.table(RSID = z_scores1[[1]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include], ZSCORE.P1 = z_scores1[[2]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include], ZSCORE.P2 = z_scores2[[2]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]) + fwrite(z_paintor, paste0(input_folder, zfile), sep = " ", col.names = TRUE) + + + write.table(zfile, file = input_file, sep = " ", col.names = FALSE, row.names = F, quote = F) + + paintor <- paste0("export NUMEXPR_NUM_THREADS=",threads,"; export OMP_NUM_THREADS=",threads,";/home/mcaiad/PAINTOR_V3.0/PAINTOR -input ", input_file, " -in ", input_folder, " -out /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/output/ -Zhead ZSCORE.P1,ZSCORE.P2 -LDname LD1,LD2 -enumerate ",K," -annotations dummy -RESname ", i, "_exact_results") + system(paintor) + cat(i,"-th rep finished.\n") +} + diff --git a/results/simulation/simu_paintor_small_nobg.R b/results/simulation/simu_paintor_small_nobg.R new file mode 100644 index 0000000..8a17f15 --- /dev/null +++ b/results/simulation/simu_paintor_small_nobg.R @@ -0,0 +1,94 @@ +########################################################### pop1 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +blas_set_num_threads(30) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- 15000 +n2 <- 20000 +K_true <- 3 +K <- K_true + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +enrich2 <- enrich * 10 +rho_sigma <- 0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +R1 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt"))) +fwrite(R1[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_nobg_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1"), + sep = " ", col.names = F) +R2 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt"))) +fwrite(R2[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_nobg_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2"), + sep = " ", col.names = F) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1")) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2")) + +annot <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_nobg_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, ".annotations") +dummy <- data.table(dummy = rep(1, p)) +fwrite(dummy, file = annot, sep = " ", col.names = TRUE) + +input_folder <- "/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/" +zfile <- paste0("paintor_mcmc_input_nobg_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep) +input_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/files_apaintor_mcmc_input_nobg_loci_",i_loci,"_small_p",p,"_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep) +for (i in 1:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_nobg_loci_",i_loci,"_small_p",p,"_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_nobg_loci_",i_loci,"_small_p",p,"_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + + z_paintor <- data.table(RSID = z_scores1[[1]], ZSCORE.P1 = z_scores1[[2]], ZSCORE.P2 = z_scores2[[2]]) + fwrite(z_paintor, paste0(input_folder, zfile), sep = " ", col.names = TRUE) + + + write.table(zfile, file = input_file, sep = " ", col.names = FALSE, row.names = F, quote = F) + + paintor <- paste0("/home/mcaiad/PAINTOR_V3.0/PAINTOR -input ", input_file, " -in ", input_folder, " -out /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/output/ -Zhead ZSCORE.P1,ZSCORE.P2 -LDname LD1,LD2 -enumerate ",K," -annotations dummy -RESname ", i, "_exact_results") + system(paintor) + cat(i,"-th rep finished.\n") +} + diff --git a/results/simulation/simu_paintor_small_t_chr22.R b/results/simulation/simu_paintor_small_t_chr22.R new file mode 100644 index 0000000..e43233a --- /dev/null +++ b/results/simulation/simu_paintor_small_t_chr22.R @@ -0,0 +1,96 @@ +########################################################### pop1 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +threads <- 20 +blas_set_num_threads(threads) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +# n2_all <- c(20000) +n1 <- 20000 +n2 <- 20000 +K_true <- 3 +K <- K_true + +df <- 4 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0 #0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/home/share/mingxuan/fine_mapping/loci_chr22_1mb.txt")) +i_loci <- 29 +idx_include <- (501:1000) + loci$idx_start[i_loci] - 1 +p <- length(idx_include) + + +# per-snp variance for bg +omega1 <- h1_omega / p +omega2 <- h2_omega / p + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +R1 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt"))) +fwrite(R1[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_t",df,"_new_wholeCHR_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1"), + sep = " ", col.names = F) +R2 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt"))) +fwrite(R2[idx_include, idx_include], paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_t",df,"_new_wholeCHR_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2"), + sep = " ", col.names = F) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R1_wg_in_sample_ref_n", n1, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD1")) +# system(paste0("cp /import/home/share/mingxuan/fine_mapping/simulation/data/simu_R2_ukb_in_sample_ref_n", n2, "_chr22_loci", i_loci, ".txt /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_input_loci_34_chr22_paintor_in_sample_n", n1, "_", n2, +# "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, +# "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".LD2")) + +annot <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/paintor_mcmc_t",df,"_new_wholeCHR_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".annotations") +dummy <- data.table(dummy = rep(1, p)) +fwrite(dummy, file = annot, sep = " ", col.names = TRUE) + +input_folder <- "/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/" +zfile <- paste0("paintor_mcmc_t",df,"_new_wholeCHR_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +input_file <- paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/paintor/working_input/files_paintor_mcmc_t",df,"_new_wholeCHR_input_loci_", i_loci, "_small_p", p, "_chr22_paintor_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep) +for (i in 1:nrep) { + z_scores1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs1_t",df,"_new_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + z_scores2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_zs2_t",df,"_new_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i, ".txt")) + + + z_paintor <- data.table(RSID = z_scores1[[1]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include], ZSCORE.P1 = z_scores1[[2]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include], ZSCORE.P2 = z_scores2[[2]][loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]) + fwrite(z_paintor, paste0(input_folder, zfile), sep = " ", col.names = TRUE) + + + write.table(zfile, file = input_file, sep = " ", col.names = FALSE, row.names = F, quote = F) + + paintor <- paste0("export NUMEXPR_NUM_THREADS=", threads, "; export OMP_NUM_THREADS=", threads, ";/home/mcaiad/PAINTOR_V3.0/PAINTOR -input ", input_file, " -in ", input_folder, " -out /import/home/share/mingxuan/fine_mapping/simulation/data/paintor/output/ -Zhead ZSCORE.P1,ZSCORE.P2 -LDname LD1,LD2 -enumerate ", K, " -annotations dummy -RESname ", i, "_exact_results") + system(paintor) + cat(i, "-th rep finished.\n") +} +