diff --git a/results/simulation/simu_susieinf_inflated_small.R b/results/simulation/simu_susieinf_inflated_small.R new file mode 100644 index 0000000..8041281 --- /dev/null +++ b/results/simulation/simu_susieinf_inflated_small.R @@ -0,0 +1,155 @@ +#!/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(20) +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) +n1 <- 15000#as.numeric(args[1]) #10000 +K_true <- 3#as.numeric(args[2]) #1 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- 0.05 +pc_var2 <- 0.2 + +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) + +bim_loci <- bim[bim$V4 >= loci$left[i_loci] & bim$V4 <= loci$right[i_loci],] +bp_start <- bim_loci$V4[idx_include[1]] +bp_end <- bim_loci$V4[idx_include[p]] +snps <- bim_loci$V2[idx_include] + +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"))) +R1 <- R1[idx_include, idx_include] +fwrite(R1, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") +Beta1 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_loci_", i_loci, "_small_p", p, "_Beta1_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".txt"))) + + +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")) + + zs1 <- z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] + + fwrite(data.frame(Z = zs1), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_inflate_susieinf_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")) + + susieinf1 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf1 <- paste0(susieinf1, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_inflate_susieinf_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") + susieinf1 <- paste0(susieinf1, " --num-sparse-effects ", K) + susieinf1 <- paste0(susieinf1, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz") + susieinf1 <- paste0(susieinf1, " --n ", n1) + susieinf1 <- paste0(susieinf1, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf1_inflate_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf1) + + + cat(i, "-th rep finished.\n") +} + + + + + +########################################################### pop1 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +blas_set_num_threads(20) +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 + +# n2_all <- c(10000, 15000, 20000) +n2 <- 20000 +K_true <- 2#as.numeric(args[2]) #1 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- 0.05 +pc_var2 <- 0.2 + +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) + +bim_loci <- bim[bim$V4 >= loci$left[i_loci] & bim$V4 <= loci$right[i_loci],] +bp_start <- bim_loci$V4[idx_include[1]] +bp_end <- bim_loci$V4[idx_include[p]] +snps <- bim_loci$V2[idx_include] + + + +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"))) +R2 <- R2[idx_include, idx_include] +fwrite(R2, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") +Beta2 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_loci_", i_loci, "_small_p", p, "_Beta2_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".txt"))) + +for (i in 1:nrep) { + + 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")) + + zs2 <- z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] + + + fwrite(data.frame(Z = zs2), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_inflate_susieinf_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")) + + + susieinf2 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf2 <- paste0(susieinf2, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_inflate_susieinf_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") + susieinf2 <- paste0(susieinf2, " --num-sparse-effects ", K) + susieinf2 <- paste0(susieinf2, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz") + susieinf2 <- paste0(susieinf2, " --n ", n2) + susieinf2 <- paste0(susieinf2, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf2_inflate_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf2) + + + cat(i, "-th rep finished.\n") +} + + diff --git a/results/simulation/simu_susieinf_inflated_small_corrected.R b/results/simulation/simu_susieinf_inflated_small_corrected.R new file mode 100644 index 0000000..ec56c01 --- /dev/null +++ b/results/simulation/simu_susieinf_inflated_small_corrected.R @@ -0,0 +1,155 @@ +#!/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(15) +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) +n1 <- 10000#as.numeric(args[1]) #10000 +K_true <- 1#as.numeric(args[2]) #1 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- 0.05 +pc_var2 <- 0.2 + +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) + +bim_loci <- bim[bim$V4 >= loci$left[i_loci] & bim$V4 <= loci$right[i_loci],] +bp_start <- bim_loci$V4[idx_include[1]] +bp_end <- bim_loci$V4[idx_include[p]] +snps <- bim_loci$V2[idx_include] + +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"))) +R1 <- R1[idx_include, idx_include] +fwrite(R1, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") +Beta1 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_loci_", i_loci, "_small_p", p, "_Beta1_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".txt"))) + + +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")) + + zs1 <- z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] + + fwrite(data.frame(Z = zs1), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_inflate_susieinf_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")) + + susieinf1 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf1 <- paste0(susieinf1, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_inflate_susieinf_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") + susieinf1 <- paste0(susieinf1, " --num-sparse-effects ", K) + susieinf1 <- paste0(susieinf1, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz") + susieinf1 <- paste0(susieinf1, " --n ", n1) + susieinf1 <- paste0(susieinf1, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf1_inflate_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_corrected_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf1) + + + cat(i, "-th rep finished.\n") +} + + + + + +########################################################### pop1 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +blas_set_num_threads(20) +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 + +# n2_all <- c(10000, 15000, 20000) +n2 <- 20000 +K_true <- 1#as.numeric(args[2]) #1 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- 0.05 +pc_var2 <- 0.2 + +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) + +bim_loci <- bim[bim$V4 >= loci$left[i_loci] & bim$V4 <= loci$right[i_loci],] +bp_start <- bim_loci$V4[idx_include[1]] +bp_end <- bim_loci$V4[idx_include[p]] +snps <- bim_loci$V2[idx_include] + + + +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"))) +R2 <- R2[idx_include, idx_include] +fwrite(R2, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") +Beta2 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_loci_", i_loci, "_small_p", p, "_Beta2_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".txt"))) + +for (i in 1:nrep) { + + 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")) + + zs2 <- z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] + + + fwrite(data.frame(Z = zs2), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_inflate_susieinf_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")) + + + susieinf2 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf2 <- paste0(susieinf2, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_inflate_susieinf_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") + susieinf2 <- paste0(susieinf2, " --num-sparse-effects ", K) + susieinf2 <- paste0(susieinf2, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz") + susieinf2 <- paste0(susieinf2, " --n ", n2) + susieinf2 <- paste0(susieinf2, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf2_inflate_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_corrected_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf2) + + + cat(i, "-th rep finished.\n") +} + + diff --git a/results/simulation/simu_susieinf_nobg.R b/results/simulation/simu_susieinf_nobg.R new file mode 100644 index 0000000..1f5ab04 --- /dev/null +++ b/results/simulation/simu_susieinf_nobg.R @@ -0,0 +1,151 @@ +########################################################### pop1 ########################################################### +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +blas_set_num_threads(20) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +bim <- fread("/import/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +n1 <- 20000 +K_true <- 2 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +enrich2 <- enrich*20 +rho_sigma <- 0 + +loci <- fread(paste0("/import/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"))) +R1 <- R1[idx_include, idx_include] +fwrite(R1, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") +Beta1 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_loci_", i_loci, "_small_p", p, "_Beta1_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".txt"))) + + +for (i in 3) { + 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")) + + zs1 <- z_scores1$V2 + + fwrite(data.frame(Z = zs1), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_nobg_susieinf_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")) + + susieinf1 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf1 <- paste0(susieinf1, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_nobg_susieinf_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") + susieinf1 <- paste0(susieinf1, " --num-sparse-effects ", K) + susieinf1 <- paste0(susieinf1, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz") + susieinf1 <- paste0(susieinf1, " --n ", n1) + susieinf1 <- paste0(susieinf1, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf1_nobg_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf1) + + + cat(i, "-th rep finished.\n") +} + + + + + +########################################################### pop1 ########################################################### +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +blas_set_num_threads(20) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +bim <- fread("/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n2_all <- c(10000, 15000, 20000) +n2 <- 20000 +K_true <- 2 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +enrich2 <- enrich*20 +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 + + +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"))) +R2 <- R2[idx_include, idx_include] +fwrite(R2, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") +Beta2 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_loci_", i_loci, "_small_p", p, "_Beta2_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".txt"))) + +for (i in 1:nrep) { + + 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")) + + zs2 <- z_scores2$V2 + + + fwrite(data.frame(Z = zs2), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_nobg_susieinf_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")) + + + susieinf2 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf2 <- paste0(susieinf2, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_nobg_susieinf_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") + susieinf2 <- paste0(susieinf2, " --num-sparse-effects ", K) + susieinf2 <- paste0(susieinf2, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz") + susieinf2 <- paste0(susieinf2, " --n ", n2) + susieinf2 <- paste0(susieinf2, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf2_nobg_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf2) + + + cat(i, "-th rep finished.\n") +} + + diff --git a/results/simulation/simu_susieinf_small_chr22.R b/results/simulation/simu_susieinf_small_chr22.R new file mode 100644 index 0000000..14efd65 --- /dev/null +++ b/results/simulation/simu_susieinf_small_chr22.R @@ -0,0 +1,182 @@ +########################################################### pop1 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +threads <- 30 +blas_set_num_threads(threads) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/import/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/import/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +n1 <- 15000 +K_true <- 1 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/import/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 + +Omega <- matrix(c(omega1, rho_omega * sqrt(omega1 * omega2), rho_omega * sqrt(omega1 * omega2), omega2), 2, 2) + +# 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"))) +R1 <- R1[idx_include, idx_include] +fwrite(R1, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") +Beta1 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_loci_", i_loci, "_small_p", p, "_Beta1_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".txt"))) + + +ldscore <- fread(paste0("/import/home/share/mingxuan/fine_mapping/LD_mat/LDscore_UKB_wegene_chr22.txt")) + +# output +pip_inf1 <- pip_inf2 <- matrix(0, p, 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")) + + + zs1 <- z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] + + fwrite(data.frame(Z = zs1), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_wholeCHR_susieinf_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")) + + + susieinf1 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf1 <- paste0(susieinf1, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_wholeCHR_susieinf_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") + susieinf1 <- paste0(susieinf1, " --num-sparse-effects ", K) + susieinf1 <- paste0(susieinf1, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz") + susieinf1 <- paste0(susieinf1, " --n ", n1) + susieinf1 <- paste0(susieinf1, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf1_wholeCHR_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf1) + + + + + + cat(i, "-th rep finished.\n") +} + + + + +thr <- 0.5 +FDR <- power <- rep(0,nrep) +lfdr <- 1-pip_inf2 +for(i in 1:nrep){ + idx_na <- is.na(lfdr[, i]) + gfdr <- cumsum(sort(lfdr[!idx_na, i], decreasing = F)) / (1:sum(!is.na(lfdr[, i]))) + idx_i <- order(lfdr[!idx_na, i], decreasing = F)[gfdr < 1 - thr] + FDR[i] <- ifelse(length(idx_i) == 0, 0, mean(!idx_i %in% idx_causal)) + power[i] <- mean(idx_causal %in% idx_i) +} +########################################################### pop2 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +threads <- 30 +blas_set_num_threads(threads) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/import/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/import/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n2_all <- c(10000, 15000, 20000) +n2 <- 20000 +K_true <- 1 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/import/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 + +Omega <- matrix(c(omega1, rho_omega * sqrt(omega1 * omega2), rho_omega * sqrt(omega1 * omega2), omega2), 2, 2) + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +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"))) +R2 <- R2[idx_include, idx_include] +fwrite(R2, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") +Beta2 <- data.matrix(fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_loci_", i_loci, "_small_p", p, "_Beta2_Ktrue", K_true, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, ".txt"))) + +ldscore <- fread(paste0("/import/home/share/mingxuan/fine_mapping/LD_mat/LDscore_UKB_wegene_chr22.txt")) + +# output +pip_inf1 <- pip_inf2 <- matrix(0, p, nrep) +for (i in 1:nrep) { + 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")) + + + zs2 <- z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] + + fwrite(data.frame(Z = zs2), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_wholeCHR_susieinf_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")) + + + susieinf2 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf2 <- paste0(susieinf2, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_wholeCHR_susieinf_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") + susieinf2 <- paste0(susieinf2, " --num-sparse-effects ", K) + susieinf2 <- paste0(susieinf2, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz") + susieinf2 <- paste0(susieinf2, " --n ", n2) + susieinf2 <- paste0(susieinf2, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf2_wholeCHR_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf2) + + + + cat(i, "-th rep finished.\n") +} + + + diff --git a/results/simulation/simu_susieinf_small_t_chr22.R b/results/simulation/simu_susieinf_small_t_chr22.R new file mode 100644 index 0000000..ec67929 --- /dev/null +++ b/results/simulation/simu_susieinf_small_t_chr22.R @@ -0,0 +1,183 @@ +########################################################### pop1 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +threads <- 30 +blas_set_num_threads(threads) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/import/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/import/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n1_all <- c(5000, 10000, 15000, 20000) +n1 <- 20000 +K_true <- 3 +K <- 5 + +df <- 16 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0#0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/import/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 + +Omega <- matrix(c(omega1, rho_omega * sqrt(omega1 * omega2), rho_omega * sqrt(omega1 * omega2), omega2), 2, 2) + +# 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"))) +R1 <- R1[idx_include, idx_include] +fwrite(R1, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") + + + +ldscore <- fread(paste0("/import/home/share/mingxuan/fine_mapping/LD_mat/LDscore_UKB_wegene_chr22.txt")) + +# output +pip_inf1 <- pip_inf2 <- matrix(0, p, 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")) + + + zs1 <- z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] + + fwrite(data.frame(Z = zs1), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_t",df,"_new_wholeCHR_susieinf_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")) + + + susieinf1 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf1 <- paste0(susieinf1, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs1_t",df,"_new_wholeCHR_susieinf_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") + susieinf1 <- paste0(susieinf1, " --num-sparse-effects ", K) + susieinf1 <- paste0(susieinf1, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R1_ref_n", n1, "_chr22_loci", i_loci, ".ld.gz") + susieinf1 <- paste0(susieinf1, " --n ", n1) + susieinf1 <- paste0(susieinf1, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf1_t",df,"_new_wholeCHR_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n1, "_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf1) + + + + + + cat(i, "-th rep finished.\n") +} + + + + +thr <- 0.5 +FDR <- power <- rep(0,nrep) +lfdr <- 1-pip_inf2 +for(i in 1:nrep){ + idx_na <- is.na(lfdr[, i]) + gfdr <- cumsum(sort(lfdr[!idx_na, i], decreasing = F)) / (1:sum(!is.na(lfdr[, i]))) + idx_i <- order(lfdr[!idx_na, i], decreasing = F)[gfdr < 1 - thr] + FDR[i] <- ifelse(length(idx_i) == 0, 0, mean(!idx_i %in% idx_causal)) + power[i] <- mean(idx_causal %in% idx_i) +} +########################################################### pop2 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +threads <- 30 +blas_set_num_threads(threads) +source("/import/home/share/mingxuan/fine_mapping/utils.R") +source("/import/home/share/mingxuan/fine_mapping/XMAP.R") +source("/import/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/import/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22.bim") +bp_chr <- bim$V4 + +nrep <- 50 + +# n2_all <- c(10000, 15000, 20000) +n2 <- 20000 +K_true <- 3 +K <- 5 + +df <- 16 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0#0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/import/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 + +Omega <- matrix(c(omega1, rho_omega * sqrt(omega1 * omega2), rho_omega * sqrt(omega1 * omega2), omega2), 2, 2) + +# per-snp variance for causal +Sigma1 <- omega1 * enrich +Sigma2 <- omega2 * enrich + +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"))) +R2 <- R2[idx_include, idx_include] +fwrite(R2, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz"), col.names = F, quote = F, sep = "\t") + +ldscore <- fread(paste0("/import/home/share/mingxuan/fine_mapping/LD_mat/LDscore_UKB_wegene_chr22.txt")) + +# output +pip_inf1 <- pip_inf2 <- matrix(0, p, nrep) +for (i in 1:nrep) { + 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")) + + + zs2 <- z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] + + fwrite(data.frame(Z = zs2), file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_t",df,"_new_wholeCHR_susieinf_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")) + + + susieinf2 <- "python fine-mapping-inf/run_fine_mapping.py --z-col-name Z --save-tsv" + susieinf2 <- paste0(susieinf2, " --sumstats /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/simu_zs2_t",df,"_new_wholeCHR_susieinf_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") + susieinf2 <- paste0(susieinf2, " --num-sparse-effects ", K) + susieinf2 <- paste0(susieinf2, " --ld-file /import/home/share/mingxuan/fine_mapping/simulation/susieinf/RunDirectory/R2_ref_n", n2, "_chr22_loci", i_loci, ".ld.gz") + susieinf2 <- paste0(susieinf2, " --n ", n2) + susieinf2 <- paste0(susieinf2, " --output-prefix /import/home/share/mingxuan/fine_mapping/simulation/susieinf/output/susieinf2_t",df,"_new_wholeCHR_susieinf_loci_", i_loci, "_small_p", p, "_chr22_n", n2, "_loci", i_loci, "_", "_Ktrue", K_true, "_K",K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + system(susieinf2) + + + + cat(i, "-th rep finished.\n") +} + + + diff --git a/results/simulation/simu_susiex_inflated _small_corrected.R b/results/simulation/simu_susiex_inflated _small_corrected.R new file mode 100644 index 0000000..49ef2e0 --- /dev/null +++ b/results/simulation/simu_susiex_inflated _small_corrected.R @@ -0,0 +1,78 @@ +#!/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("/import/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/import/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(10000, 15000, 20000) +n1 <- 20000#as.numeric(args[1]) #10000 +n2 <- 20000 +K_true <- 1#as.numeric(args[2]) #1 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- 0.05 +pc_var2 <- 0.2 + +loci <- fread(paste0("/import/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) + +bim_loci <- bim[bim$V4 >= loci$left[i_loci] & bim$V4 <= loci$right[i_loci],] +bp_start <- bim_loci$V4[idx_include[1]] +bp_end <- bim_loci$V4[idx_include[p]] +snps <- bim_loci$V2[idx_include] + + + +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")) + ss1 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] / sqrt(n1), se = 1 / sqrt(n1), pval = 2 * pnorm(abs(z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), lower.tail = F)) + ss2 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] / sqrt(n2), se = 1 / sqrt(n2), pval = 2 * pnorm(abs(z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), lower.tail = F)) + + write.table(ss1, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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"), col.names = T, row.names = F, quote = F) + + write.table(ss2, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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"), col.names = T, row.names = F, quote = F) + + susiex <- paste0("export NUMEXPR_NUM_THREADS=10; export OMP_NUM_THREADS=10; python SuSiEx/SuSiEx.py --pval_thresh=1 --min_purity=0 --max_iter=500 --out_dir=/import/home/share/mingxuan/fine_mapping/simulation/susiex/output --chr=22 --bp=", bp_start - 1, ",", bp_end + 1, " --chr_col=1,1 --snp_col=2,2 --bp_col=3,3 --a1_col=4,4 --a2_col=5,5 --eff_col=6,6 --se_col=7,7 --pval_col=8,8 --keep-ambig=True --mult-step=False --maf=0.0001 --full_out=True") + susiex <- paste0(susiex, " --n_sig=", K) + susiex <- paste0(susiex, " --sst_file=/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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,/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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") + susiex <- paste0(susiex, " --ld_file=/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1, ",/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2) + susiex <- paste0(susiex, " --n_gwas=", n1, ",", n2) + susiex <- paste0(susiex, " --out_name=susiex_output_inflate_pc1_", pc_var1,"_",pc_var2, "_corrected_loci_", i_loci, "_small_p", p, "_chr22_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_K", K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + + system(susiex) + cat(i, "-th rep finished.\n") +} \ No newline at end of file diff --git a/results/simulation/simu_susiex_inflated_small.R b/results/simulation/simu_susiex_inflated_small.R new file mode 100644 index 0000000..0fdabf2 --- /dev/null +++ b/results/simulation/simu_susiex_inflated_small.R @@ -0,0 +1,78 @@ +#!/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(20) +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(10000, 15000, 20000) +n1 <- 20000#as.numeric(args[1]) #10000 +n2 <- 20000 +K_true <- 3#as.numeric(args[2]) #1 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +pc_var1 <- 0.05 +pc_var2 <- 0.2 + +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) + +bim_loci <- bim[bim$V4 >= loci$left[i_loci] & bim$V4 <= loci$right[i_loci],] +bp_start <- bim_loci$V4[idx_include[1]] +bp_end <- bim_loci$V4[idx_include[p]] +snps <- bim_loci$V2[idx_include] + + + +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")) + ss1 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] / sqrt(n1), se = 1 / sqrt(n1), pval = 2 * pnorm(abs(z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), lower.tail = F)) + ss2 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] / sqrt(n2), se = 1 / sqrt(n2), pval = 2 * pnorm(abs(z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), lower.tail = F)) + + write.table(ss1, file = paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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"), col.names = T, row.names = F, quote = F) + + write.table(ss2, file = paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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"), col.names = T, row.names = F, quote = F) + + susiex <- paste0("export NUMEXPR_NUM_THREADS=10; export OMP_NUM_THREADS=10; python SuSiEx/SuSiEx.py --pval_thresh=1 --min_purity=0 --max_iter=500 --out_dir=/home/share/mingxuan/fine_mapping/simulation/susiex/output --chr=22 --bp=", bp_start - 1, ",", bp_end + 1, " --chr_col=1,1 --snp_col=2,2 --bp_col=3,3 --a1_col=4,4 --a2_col=5,5 --eff_col=6,6 --se_col=7,7 --pval_col=8,8 --keep-ambig=True --mult-step=False --maf=0.0001 --full_out=True") + susiex <- paste0(susiex, " --n_sig=", K) + susiex <- paste0(susiex, " --sst_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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,/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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") + susiex <- paste0(susiex, " --ld_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1, ",/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2) + susiex <- paste0(susiex, " --n_gwas=", n1, ",", n2) + susiex <- paste0(susiex, " --out_name=susiex_output_inflate_pc1_", pc_var1,"_",pc_var2, "_loci_", i_loci, "_small_p", p, "_chr22_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_K", K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + + system(susiex) + cat(i, "-th rep finished.\n") +} \ No newline at end of file diff --git a/results/simulation/simu_susiex_small_chr22.R b/results/simulation/simu_susiex_small_chr22.R new file mode 100644 index 0000000..c01d551 --- /dev/null +++ b/results/simulation/simu_susiex_small_chr22.R @@ -0,0 +1,115 @@ +########################################################### 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("/import/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/import/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(10000, 15000, 20000) +n1 <- 10000 +n2 <- 20000 +K_true <- 2 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/import/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) + +bim_loci <- bim[bim$V4>=loci$left[i_loci]&bim$V4<=loci$right[i_loci],] +bp_start <- bim_loci$V4[idx_include[1]] +bp_end <- bim_loci$V4[idx_include[p]] +snps <- bim_loci$V2[idx_include] + +# construct in-sample genotype for gwas and ld matrix in pop1 +# n1_all <- c(5000, 10000, 15000, 20000) +# for (n1 in n1_all) { +# +# ind_n1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_indiv_wg_in_sample_ref_n", n1, ".txt"))$x +# iid_fid <- cbind(ind_n1, ind_n1) +# write.table(iid_fid, paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_wg_in_sample_ref_n", n1, ".txt"), col.names = F, row.names = F, quote = F) +# +# susiex_ld <- paste0("python SuSiEx/SuSiEx_LD.py --ref_file=/home/share/mingxuan/wegene_qc2/height_merge_qc2_chr22 --chr=22 --bp=", bp_start, ",", bp_end, " --plink=/import/home2/mcaiad/plink/plink --maf=0.0001") +# susiex_ld <- paste0(susiex_ld, " --keep=/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_wg_in_sample_ref_n", n1, ".txt") +# susiex_ld <- paste0(susiex_ld, " --ld_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1) +# +# system(susiex_ld) +# +# system(paste0("mv /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n",n1,".ld.gz /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n",n1,".ld.gz.plink")) +# 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"))) +# R1 <- R1[idx_include, idx_include] +# fwrite(R1,file=paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1,".ld.gz"),col.names = F,row.names = F,sep="\t") +# +# } +# +# +# n2_all <- c(10000, 15000, 20000) +# for (n2 in n2_all) { +# ind_n2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_indiv_ukb_in_sample_ref_n", n2, ".txt"))$x +# iid_fid <- cbind(ind_n2, ind_n2) +# write.table(iid_fid, paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_ukb_in_sample_ref_n", n2, ".txt"), col.names = F, row.names = F, quote = F) +# +# susiex_ld <- paste0("python SuSiEx/SuSiEx_LD.py --ref_file=/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22 --chr=22 --bp=", bp_start, ",", bp_end, " --plink=/import/home2/mcaiad/plink/plink --maf=0.0001") +# susiex_ld <- paste0(susiex_ld, " --keep=/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_ukb_in_sample_ref_n", n2, ".txt") +# susiex_ld <- paste0(susiex_ld, " --ld_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2) +# +# system(susiex_ld) +# +# system(paste0("mv /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n",n2,".ld.gz /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n",n2,".ld.gz.plink")) +# 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"))) +# R2 <- R2[idx_include, idx_include] +# fwrite(R2,file=paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2,".ld.gz"),col.names = F,row.names = F,sep="\t") +# } + + + +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")) + ss1 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] / sqrt(n1), + se = 1 / sqrt(n1), pval = 2 * pnorm(abs(z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), lower.tail = F)) + ss2 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] / sqrt(n2), + se = 1 / sqrt(n2), pval = 2 * pnorm(abs(z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), lower.tail = F)) + + write.table(ss1, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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"), col.names = T, row.names = F, quote = F) + + write.table(ss2, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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"), col.names = T, row.names = F, quote = F) + + susiex <- paste0("export NUMEXPR_NUM_THREADS=10; export OMP_NUM_THREADS=10; python SuSiEx/SuSiEx.py --pval_thresh=1 --min_purity=0 --max_iter=500 --out_dir=/import/home/share/mingxuan/fine_mapping/simulation/susiex/output --chr=22 --bp=",bp_start,",",bp_end," --chr_col=1,1 --snp_col=2,2 --bp_col=3,3 --a1_col=4,4 --a2_col=5,5 --eff_col=6,6 --se_col=7,7 --pval_col=8,8 --keep-ambig=True --mult-step=False --maf=0.0001 --full_out=True") + susiex <- paste0(susiex," --n_sig=",K) + susiex <- paste0(susiex," --sst_file=/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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,/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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") + susiex <- paste0(susiex," --ld_file=/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n",n1,",/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2) + susiex <- paste0(susiex," --n_gwas=",n1,",",n2) + susiex <- paste0(susiex," --out_name=susiex_output_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_K",K,"_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_",i) + + system(susiex) + cat(i, "-th rep finished.\n") +} \ No newline at end of file diff --git a/results/simulation/simu_susiex_small_nobg.R b/results/simulation/simu_susiex_small_nobg.R new file mode 100644 index 0000000..ea30999 --- /dev/null +++ b/results/simulation/simu_susiex_small_nobg.R @@ -0,0 +1,113 @@ +########################################################### pop1 ########################################################### +library(VCM) +library(XPASS) +library(mvtnorm) +library(susieR) +library(data.table) +library(RhpcBLASctl) +set.seed(1) +blas_set_num_threads(15) +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(10000, 15000, 20000) +n1 <- 5000 +n2 <- 20000 +K_true <- 1 +K <- 5 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0.8 + +enrich <- 50 +enrich2 <- enrich * 20 +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) + +bim_loci <- bim[bim$V4 >= loci$left[i_loci] & bim$V4 <= loci$right[i_loci],] +bp_start <- bim_loci$V4[idx_include[1]] +bp_end <- bim_loci$V4[idx_include[p]] +snps <- bim_loci$V2[idx_include] + +# construct in-sample genotype for gwas and ld matrix in pop1 +# n1_all <- c(5000, 10000, 15000, 20000) +# for (n1 in n1_all) { +# +# ind_n1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_indiv_wg_in_sample_ref_n", n1, ".txt"))$x +# iid_fid <- cbind(ind_n1, ind_n1) +# write.table(iid_fid, paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_wg_in_sample_ref_n", n1, ".txt"), col.names = F, row.names = F, quote = F) +# +# susiex_ld <- paste0("python SuSiEx/SuSiEx_LD.py --ref_file=/home/share/mingxuan/wegene_qc2/height_merge_qc2_chr22 --chr=22 --bp=", bp_start, ",", bp_end, " --plink=/import/home2/mcaiad/plink/plink --maf=0.0001") +# susiex_ld <- paste0(susiex_ld, " --keep=/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_wg_in_sample_ref_n", n1, ".txt") +# susiex_ld <- paste0(susiex_ld, " --ld_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1) +# +# system(susiex_ld) +# +# system(paste0("mv /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n",n1,".ld.gz /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n",n1,".ld.gz.plink")) +# 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"))) +# R1 <- R1[idx_include, idx_include] +# fwrite(R1,file=paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1,".ld.gz"),col.names = F,row.names = F,sep="\t") +# +# } +# +# +# n2_all <- c(10000, 15000, 20000) +# for (n2 in n2_all) { +# ind_n2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_indiv_ukb_in_sample_ref_n", n2, ".txt"))$x +# iid_fid <- cbind(ind_n2, ind_n2) +# write.table(iid_fid, paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_ukb_in_sample_ref_n", n2, ".txt"), col.names = F, row.names = F, quote = F) +# +# susiex_ld <- paste0("python SuSiEx/SuSiEx_LD.py --ref_file=/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22 --chr=22 --bp=", bp_start, ",", bp_end, " --plink=/import/home2/mcaiad/plink/plink --maf=0.0001") +# susiex_ld <- paste0(susiex_ld, " --keep=/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_ukb_in_sample_ref_n", n2, ".txt") +# susiex_ld <- paste0(susiex_ld, " --ld_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2) +# +# system(susiex_ld) +# +# system(paste0("mv /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n",n2,".ld.gz /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n",n2,".ld.gz.plink")) +# 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"))) +# R2 <- R2[idx_include, idx_include] +# fwrite(R2,file=paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2,".ld.gz"),col.names = F,row.names = F,sep="\t") +# } + + +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")) + ss1 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores1$V2 / sqrt(n1), se = 1 / sqrt(n1), pval = 2 * pnorm(abs(z_scores1$V2), lower.tail = F)) + ss2 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores2$V2 / sqrt(n2), se = 1 / sqrt(n2), pval = 2 * pnorm(abs(z_scores2$V2), lower.tail = F)) + + write.table(ss1, file = paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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"), col.names = T, row.names = F, quote = F) + + write.table(ss2, file = paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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"), col.names = T, row.names = F, quote = F) + + susiex <- paste0("export NUMEXPR_NUM_THREADS=10; export OMP_NUM_THREADS=10; python SuSiEx/SuSiEx.py --pval_thresh=1 --min_purity=0 --max_iter=500 --out_dir=/home/share/mingxuan/fine_mapping/simulation/susiex/output --chr=22 --bp=", bp_start - 1, ",", bp_end + 1, " --chr_col=1,1 --snp_col=2,2 --bp_col=3,3 --a1_col=4,4 --a2_col=5,5 --eff_col=6,6 --se_col=7,7 --pval_col=8,8 --keep-ambig=True --mult-step=False --maf=0.0001 --full_out=True") + susiex <- paste0(susiex, " --n_sig=", K) + susiex <- paste0(susiex, " --sst_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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,/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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") + susiex <- paste0(susiex, " --ld_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1, ",/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2) + susiex <- paste0(susiex, " --n_gwas=", n1, ",", n2) + susiex <- paste0(susiex, " --out_name=susiex_output_nobg_loci_", i_loci, "_small_p", p, "_chr22_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_K", K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich2, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + + system(susiex) + cat(i, "-th rep finished.\n") +} \ No newline at end of file diff --git a/results/simulation/simu_susiex_small_t_chr22.R b/results/simulation/simu_susiex_small_t_chr22.R new file mode 100644 index 0000000..d7e2fef --- /dev/null +++ b/results/simulation/simu_susiex_small_t_chr22.R @@ -0,0 +1,115 @@ +########################################################### 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("/import/home/share/mingxuan/taam/code/sldxr.R") +bim <- fread("/import/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(10000, 15000, 20000) +n1 <- 20000 +n2 <- 20000 +K_true <- 3 +K <- 5 +df <- 16 + +h1_omega <- 0.005 +h2_omega <- 0.005 +rho_omega <- 0 #0.8 + +enrich <- 50 +rho_sigma <- 0 + +loci <- fread(paste0("/import/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) + +bim_loci <- bim[bim$V4 >= loci$left[i_loci] & bim$V4 <= loci$right[i_loci],] +bp_start <- bim_loci$V4[idx_include[1]] +bp_end <- bim_loci$V4[idx_include[p]] +snps <- bim_loci$V2[idx_include] + +# construct in-sample genotype for gwas and ld matrix in pop1 +# n1_all <- c(5000, 10000, 15000, 20000) +# for (n1 in n1_all) { +# +# ind_n1 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_indiv_wg_in_sample_ref_n", n1, ".txt"))$x +# iid_fid <- cbind(ind_n1, ind_n1) +# write.table(iid_fid, paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_wg_in_sample_ref_n", n1, ".txt"), col.names = F, row.names = F, quote = F) +# +# susiex_ld <- paste0("python SuSiEx/SuSiEx_LD.py --ref_file=/home/share/mingxuan/wegene_qc2/height_merge_qc2_chr22 --chr=22 --bp=", bp_start, ",", bp_end, " --plink=/import/home2/mcaiad/plink/plink --maf=0.0001") +# susiex_ld <- paste0(susiex_ld, " --keep=/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_wg_in_sample_ref_n", n1, ".txt") +# susiex_ld <- paste0(susiex_ld, " --ld_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1) +# +# system(susiex_ld) +# +# system(paste0("mv /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n",n1,".ld.gz /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n",n1,".ld.gz.plink")) +# 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"))) +# R1 <- R1[idx_include, idx_include] +# fwrite(R1,file=paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1,".ld.gz"),col.names = F,row.names = F,sep="\t") +# +# } +# +# +# n2_all <- c(10000, 15000, 20000) +# for (n2 in n2_all) { +# ind_n2 <- fread(paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_indiv_ukb_in_sample_ref_n", n2, ".txt"))$x +# iid_fid <- cbind(ind_n2, ind_n2) +# write.table(iid_fid, paste0("/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_ukb_in_sample_ref_n", n2, ".txt"), col.names = F, row.names = F, quote = F) +# +# susiex_ld <- paste0("python SuSiEx/SuSiEx_LD.py --ref_file=/home/share/mingxuan/wegene_qc2/height_ukb_qc2_chr22 --chr=22 --bp=", bp_start, ",", bp_end, " --plink=/import/home2/mcaiad/plink/plink --maf=0.0001") +# susiex_ld <- paste0(susiex_ld, " --keep=/import/home/share/mingxuan/fine_mapping/simulation/data/simu_iid_fid_ukb_in_sample_ref_n", n2, ".txt") +# susiex_ld <- paste0(susiex_ld, " --ld_file=/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2) +# +# system(susiex_ld) +# +# system(paste0("mv /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n",n2,".ld.gz /home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n",n2,".ld.gz.plink")) +# 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"))) +# R2 <- R2[idx_include, idx_include] +# fwrite(R2,file=paste0("/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2,".ld.gz"),col.names = F,row.names = F,sep="\t") +# } + + +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")) + ss1 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] / sqrt(n1), + se = 1 / sqrt(n1), pval = 2 * pnorm(abs(z_scores1$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), lower.tail = F)) + ss2 <- data.frame(chr = 22, snp = snps, bp = bim_loci$V4[idx_include], a1 = bim_loci$V5[idx_include], a2 = bim_loci$V6[idx_include], + BETA = z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include] / sqrt(n2), + se = 1 / sqrt(n2), pval = 2 * pnorm(abs(z_scores2$V2[loci$idx_left[i_loci]:loci$idx_right[i_loci]][idx_include]), lower.tail = F)) + + write.table(ss1, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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"), col.names = T, row.names = F, quote = F) + + write.table(ss2, file = paste0("/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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"), col.names = T, row.names = F, quote = F) + + susiex <- paste0("export NUMEXPR_NUM_THREADS=10; export OMP_NUM_THREADS=10; python SuSiEx/SuSiEx.py --pval_thresh=1 --min_purity=0 --max_iter=500 --out_dir=/import/home/share/mingxuan/fine_mapping/simulation/susiex/output --chr=22 --bp=", bp_start, ",", bp_end, " --chr_col=1,1 --snp_col=2,2 --bp_col=3,3 --a1_col=4,4 --a2_col=5,5 --eff_col=6,6 --se_col=7,7 --pval_col=8,8 --keep-ambig=True --mult-step=False --maf=0.0001 --full_out=True") + susiex <- paste0(susiex, " --n_sig=", K) + susiex <- paste0(susiex, " --sst_file=/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss1_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,/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/ss2_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") + susiex <- paste0(susiex, " --ld_file=/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD1_n", n1, ",/import/home/share/mingxuan/fine_mapping/simulation/susiex/RunDirectory/LD2_n", n2) + susiex <- paste0(susiex, " --n_gwas=", n1, ",", n2) + susiex <- paste0(susiex, " --out_name=susiex_t", df, "_new_output_wholeCHR_loci_", i_loci, "_small_p", p, "_chr22_in_sample_n", n1, "_", n2, + "_Ktrue", K_true, "_K", K, "_Omega", h1_omega, "_", h2_omega, "_", rho_omega, + "_enrich", enrich, "_rhoSigma", rho_sigma, "_nrep", nrep, "_", i) + + system(susiex) + cat(i, "-th rep finished.\n") +} \ No newline at end of file