hJAM is a hierarchical model which unifies the framework of Mendelian Randomization and Transcriptome-wide association studies.
You can install the development version from GitHub with:
if (!require("devtools")) { install.packages("devtools") } else {}
devtools::install_github("lailylajiang/hJAM")
This is a basic example of fitting hJAM model:
library(hJAM)
# Download the data for data example 2 from the package
data(Gl)
data(betas.Gy)
data(marginal_A)
If you don't have conditional A matrix, you can use get_cond_A
(if more than one X) or get_cond_alpha
(if only one X) to convert the marginal effects to conditional A matrix with the reference panel.
conditional_A = get_cond_A(marginal_A = marginal_A, Gl = Gl, N.Gx = 339224, ridgeTerm = T)
conditional_alpha = get_cond_alpha(alphas = marginal_A[, 1], Gl = Gl, N.Gx = 339224, ridgeTerm = T)
After obtained the conditional A matrix, fit hJAM model with function hJAM_lnreg
.
# fit the hJAM model
hJAM_lnreg(betas.Gy = betas.Gy, Gl = Gl, N.Gy = 459324, A = conditional_A, ridgeTerm = T)
#> ------------------------------------------------------
#> hJAM output
#> ------------------------------------------------------
#> Number of SNPs used in model: 210
#>
#> Estimate StdErr 95% CI Pvalue
#> bmi 0.322 0.061 (0.202, 0.442) 3.189944e-07
#> t2d 0.119 0.017 (0.085, 0.153) 4.124526e-11
#> ------------------------------------------------------
In the package, you could also implement hJAM with Egger regression, which is designed to detect the unmeasured pleiotropy effect. The function for hJAM with Egger regression is hJAM_egger
.
# fit the hJAM model
hJAM_egger(betas.Gy = betas.Gy, Gl = Gl, N.Gy = 459324, A = conditional_A, ridgeTerm = T)
#> ------------------------------------------------------
#> hJAM egger output
#> ------------------------------------------------------
#> Number of SNPs used in model: 210
#>
#> Estimate StdErr 95% CI Pvalue
#> bmi 0.302 0.070 (0.163, 0.44) 2.817908e-05
#> t2d 0.107 0.027 (0.055, 0.16) 8.175202e-05
#>
#> Intercept
#> Est.Int StdErr.Int 95% CI.Int Pvalue.Int
#> [1,] "0.453" "0.787" "(-1.099, 2.005)" "0.565"
#> ------------------------------------------------------
The user could use SNPs_heatmap
and SNPs_scatter_plot
to display the correlation and pattern of the SNPs that the user used in the analysis.
scatter_plot_p = SNPs_scatter_plot(A = conditional_A, betas.Gy = betas.Gy, num_X = 2)
scatter_plot_p
heatmap_p = SNPs_heatmap(Gl)
heatmap_p
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.