data("data_matrix_final")
data("covariance_matrix_data")
data("gene_list")
data("gene_length_list")
data("selected_genotype")
data("PCA_result") ##Set basic coefficient
covariates<-PCA_result$principal.coordinates[,1:2]
n_sim=10
q=6
gamma_value<-1
phenotype_covariance<-stats::cor(covariance_matrix_data[8:13])
select_gene<-gene_list[9592]
match_char=paste(select_gene,"\\(",sep = "")
match_line=which(grepl(match_char,data_matrix_final$ANNOT))
genotype_data<-t(selected_genotype[,7:2510])
adj_geno<-stats::lm(genotype_data~covariates)$residuals
genotype_covariance<-stats::cor(adj_geno)
n_sample<-nrow(genotype_data)
geno_length<-ncol(genotype_data)
cutoff_value<-seq(from=0.1,to=1,by=0.1)
c<-as.matrix(rnorm(2504,0,1))
gamma_vec<-as.matrix(rep(gamma_value,q))
c_mat<-c%*%t(gamma_vec)
covariance_matrix<-stats::cor(covariance_matrix_data[8:13])
B_mat<-matrix(0, nrow = geno_length, ncol = q)
B_mat[6,1]=0.1
B_mat[6,3]=-0.1
B_mat[6,5]=-0.1
est_pheno_rho_mat<-covariance_matrix
total_mat<-methods::kronecker(genotype_covariance,est_pheno_rho_mat)##Generate null distribution and calculate coefficient for T3
null_distribution_T3<-ThreeWayTest::generate_null_distribution_T3(
m = q*geno_length,
n = 10000,
cov_mat = total_mat,
cutoff_value = cutoff_value)
coefficient_matrix_T3<-ThreeWayTest::approximate_distribution_coefficient_estimate_T3(
null_distribution_matrix = null_distribution_T3)##Create variables for storage
p_value_MGAS<-rep(NA,n_sim)
p_value_TWT<-rep(NA,n_sim)
p_value_metaCCA<-rep(NA,n_sim)
p_1<-rep(NA,n_sim)
p_2<-rep(NA,n_sim)
p_3<-rep(NA,n_sim)##Main simulation function
for (i in 1:n_sim) {
eplison<- MASS::mvrnorm(n=n_sample,rep(0,q),covariance_matrix)
Y_mat<-genotype_data%*%B_mat+c_mat+eplison
Z_mat<-matrix(NA,ncol(genotype_data),ncol(Y_mat))
std_beta_mat<-matrix(NA,ncol(genotype_data),ncol(Y_mat))
for (j in 1:ncol(genotype_data)) {
for (k in 1:ncol(Y_mat)) {
model<- stats::lm(Y_mat[,k]~genotype_data[,j]+covariates+c)
wald <- stats::coef(model)[2] / sqrt(diag(vcov(model)))[2]
Z_mat[j,k] <- wald
std_beta_mat[j,k]<- coef(model)[2]/sqrt(n_sample*(diag(vcov(model))[2]))
}
}
z_vector<-as.vector(t(Z_mat))
p_value_MGAS[i]<-ThreeWayTest::MGAS(z_vector,genotype_covariance,est_pheno_rho_mat)
result_TWT<-ThreeWayTest::TWT(t(Z_mat),genotype_covariance,est_pheno_rho_mat,cutoff_value,coefficient_matrix_T3)
p_value_TWT[i]<-result_TWT$p_value_final
p_1[i]<-result_TWT$p_1
p_2[i]<-result_TWT$p_2
p_3[i]<-result_TWT$p_3
p_value_metaCCA[i]<-ThreeWayTest::metaCCA(genotype_covariance,est_pheno_rho_mat,std_beta_mat,n_sample)
print(i)
}
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#> [1] 5
#> [1] 6
#> [1] 7
#> [1] 8
#> [1] 9
#> [1] 10##Display result
sum(as.numeric(p_value_MGAS<0.001),na.rm = TRUE)
#> [1] 2
sum(as.numeric(p_value_TWT<0.001),na.rm = TRUE)
#> [1] 6
sum(as.numeric(p_value_metaCCA<0.001),na.rm = TRUE)
#> [1] 4
sum(as.numeric(p_1<0.001),na.rm = TRUE)
#> [1] 0
sum(as.numeric(p_2<0.001),na.rm = TRUE)
#> [1] 6
sum(as.numeric(p_3<0.001),na.rm = TRUE)
#> [1] 4
utils::sessionInfo()
#> R Under development (unstable) (2023-03-02 r83926)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.2 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ThreeWayTest_1.0.2
#>
#> loaded via a namespace (and not attached):
#> [1] vctrs_0.5.2 cli_3.6.0 knitr_1.42 rlang_1.0.6
#> [5] xfun_0.37 stringi_1.7.12 purrr_1.0.1 textshaping_0.3.6
#> [9] data.table_1.14.8 jsonlite_1.8.4 glue_1.6.2 rprojroot_2.0.3
#> [13] htmltools_0.5.4 ragg_1.2.5 sass_0.4.5 rmarkdown_2.20
#> [17] evaluate_0.20 jquerylib_0.1.4 MASS_7.3-58.2 fastmap_1.1.1
#> [21] yaml_2.3.7 lifecycle_1.0.3 memoise_2.0.1 stringr_1.5.0
#> [25] compiler_4.3.0 fs_1.6.1 systemfonts_1.0.4 digest_0.6.31
#> [29] R6_2.5.1 magrittr_2.0.3 bslib_0.4.2 tools_4.3.0
#> [33] pkgdown_2.0.7 cachem_1.0.7 desc_1.4.2