Note: the two steps above guarantee that the LD correlations are calculated based on the effect alleles for the SNP effects.
```{r eval=FALSE} # Estimate LD correlation matrix using R snp_coeff_id = scan("gsmr_example.xmat.gz", what="", nlines=1) snp_coeff = read.table("gsmr_example.xmat.gz", header=F, skip=2) ``` ```{r} # Match the SNP genotype data with the summary data snp_id = Reduce(intersect, list(gsmr_data$SNP, snp_coeff_id)) gsmr_data = gsmr_data[match(snp_id, gsmr_data$SNP),] snp_order = match(snp_id, snp_coeff_id) snp_coeff_id = snp_coeff_id[snp_order] snp_coeff = snp_coeff[, snp_order] # Calculate the LD correlation matrix ldrho = cor(snp_coeff) # Check the size of the correlation matrix and double-check if the order of the SNPs in the LD correlation matrix is consistent with that in the GWAS summary data colnames(ldrho) = rownames(ldrho) = snp_coeff_id ``` ```{r} dim(ldrho) # Show the first 5 rows and columns of the matrix ldrho[1:5,1:5] ``` \textcolor{red}{Note: all the analyses implemented in this R-package only require the summary data and the LD correlation matrix ("ldrho") listed above.}Note: all the analyses implemented in this R-package only require the summary data (e.g. "gsmr_data") and the LD correlation matrix (e.g. "ldrho") listed above.
### 2. Standardization This is an optional process. If the risk factor was not standardised in GWAS, the effect sizes can be scaled using the method below. Note that this process requires allele frequencies, z-statistics and sample size. After scaling, bzx is interpreted as the per-allele effect of a SNP on the exposure in standard deviation units. ```{r message=FALSE, warning=FALSE} snpfreq = gsmr_data$a1_freq # allele frequencies of the SNPs bzx = gsmr_data$bzx # effects of the instruments on risk factor bzx_se = gsmr_data$bzx_se # standard errors of bzx bzx_n = gsmr_data$bzx_n # GWAS sample size for the risk factor std_zx = std_effect(snpfreq, bzx, bzx_se, bzx_n) # perform standardisation gsmr_data$std_bzx = std_zx$b # standardized bzx gsmr_data$std_bzx_se = std_zx$se # standardized bzx_se head(gsmr_data) ``` ### 3. GSMR analysis This is the main analysis of this R-package. It uses SNPs associated with the risk factor (e.g. at p < 5e-8) as the instruments to test for putative causal effect of the risk factor on the disease. The analysis involves a step that uses the [HEIDI-outlier](#4.HEIDI-outlieranalysis) approach to remove SNPs that have effects on both the risk factor and the disease because of pleiotropy. ```{r message=FALSE, warning=FALSE} bzx = gsmr_data$std_bzx # SNP effects on the risk factor bzx_se = gsmr_data$std_bzx_se # standard errors of bzx bzx_pval = gsmr_data$bzx_pval # p-values for bzx bzy = gsmr_data$bzy # SNP effects on the disease bzy_se = gsmr_data$bzy_se # standard errors of bzy bzy_pval = gsmr_data$bzy_pval # p-values for bzy n_ref = 7703 # Sample size of the reference sample gwas_thresh = 5e-8 # GWAS threshold to select SNPs as the instruments for the GSMR analysis heidi_outlier_thresh = 0.01 # HEIDI-outlier threshold nsnps_thresh = 10 # the minimum number of instruments required for the GSMR analysis heidi_outlier_flag = T # flag for HEIDI-outlier analysis ld_r2_thresh = 0.1 # LD r2 threshold to remove SNPs in high LD ld_fdr_thresh = 0.05 # FDR threshold to remove the chance correlations between the SNP instruments gsmr_results = gsmr(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, snp_coeff_id, n_ref, heidi_outlier_flag, gwas_thresh, heidi_outlier_thresh, nsnps_thresh, ld_r2_thresh, ld_fdr_thresh) # GSMR analysis cat("The estimated effect of the exposure on outcome: ",gsmr_results$bxy) cat("Standard error of bxy: ",gsmr_results$bxy_se) cat("P-value for bxy: ", gsmr_results$bxy_pval) cat("Indexes of the SNPs used in the GSMR analysis: ", gsmr_results$used_index[1:5], "...") cat("Number of SNPs with missing estimates in the summary data: ", length(gsmr_results$na_snps)) cat("Number of non-significant SNPs: ", length(gsmr_results$weak_snps)) cat("Number of SNPs in high LD ( LD rsq >", ld_r2_thresh, "): ", length(gsmr_results$linkage_snps)) cat("Number of pleiotropic outliers: ", length(gsmr_results$pleio_snps)) ``` ### 4. HEIDI-outlier analysis The estimate of causal effect of risk factor on disease can be biased by pleiotropy ([Zhu et al. 2018 Nat. Commun](https://www.nature.com/articles/s41467-017-02317-2)). This is an analysis to detect and eliminate from the analysis instruments that show significant pleiotropic effects on both risk factor and disease. The HEIDI-outlier analysis requires bzx (effect of genetic instrument on risk factor), bzx_se (standard error of bzx), bzx_pval (p-value of bzx), bzy (effect of genetic instrument on disease), bzy_se (standard error of bzy) and ldrho (LD matrix of instruments). Note that similar to that in the GSMR analysis above, the LD matrix can be estimated from a reference sample with individual-level genotype data. **The HEIDI-outlier analysis has been integrated in the GSMR analysis above (with the heidi_outlier_flag and heidi_outlier_thresh flags).** It can also be performed separately following the example below. ```{r message=FALSE, warning=FALSE} heidi_results = heidi_outlier(bzx, bzx_se, bzx_pval, bzy, bzy_se, ldrho, snp_coeff_id, n_ref, gwas_thresh, heidi_outlier_thresh, nsnps_thresh, ld_r2_thresh, ld_fdr_thresh) # perform HEIDI-outlier analysis cat("Number of SNPs in high LD ( LD rsq >", ld_r2_thresh, "): ", length(gsmr_results$linkage_snps)) cat("Number of pleiotropic outliers: ", length(heidi_results$pleio_snps)) filtered_index = heidi_results$remain_index filtered_gsmr_data = gsmr_data[filtered_index,] # select data passed HEIDI-outlier filtering filtered_snp_id = snp_coeff_id[filtered_index] # select SNPs that passed HEIDI-outlier filtering dim(filtered_gsmr_data) # Number of SNPs in the gmsr_data with bzx_pval < 5e-8 dim(gsmr_data[gsmr_data$bzx_pval < 5e-8, ]) ``` In the example above, `r dim(gsmr_data[gsmr_data$bzx_pval <= 5e-8, ])[1] - dim(filtered_gsmr_data)[1]` SNPs are filtered out by HEIDI-outlier. ### 5. Bi-directional GSMR analysis The script below runs bi-directional GSMR analyses, i.e. a forward-GSMR analysis as described above and a reverse-GSMR analysis that uses SNPs associated with the disease (e.g. at p < 5e-8) as the instruments to test for putative causal effect of the disease on risk factor. ```{r message=FALSE, warning=FALSE} gsmr_results = bi_gsmr(bzx, bzx_se, bzx_pval, bzy, bzy_se, bzy_pval, ldrho, snp_coeff_id, n_ref, heidi_outlier_flag, gwas_thresh, heidi_outlier_thresh, nsnps_thresh, ld_r2_thresh, ld_fdr_thresh) # GSMR analysis cat("Effect of risk factor on disease: ",gsmr_results$forward_bxy) cat("Standard error of bxy in the forward-GSMR analysis: ",gsmr_results$forward_bxy_se) cat("P-value of bxy in the forward-GSMR analysis: ", gsmr_results$forward_bxy_pval) cat("Effect of disease on risk factor: ",gsmr_results$reverse_bxy) cat("Standard error of bxy in the reverse-GSMR analysis: ",gsmr_results$reverse_bxy_se) cat("P-value of bxy in the reverse-GSMR analysis: ", gsmr_results$reverse_bxy_pval) ``` ### 6. Visualization ```{r, fig.width=6, fig.height=6.5} effect_col = colors()[75] vals = c(bzx[filtered_index]-bzx_se[filtered_index], bzx[filtered_index]+bzx_se[filtered_index]) xmin = min(vals); xmax = max(vals) vals = c(bzy[filtered_index]-bzy_se[filtered_index], bzy[filtered_index]+bzy_se[filtered_index]) ymin = min(vals); ymax = max(vals) par(mar=c(5,5,4,2)) plot(bzx[filtered_index], bzy[filtered_index], pch=20, cex=0.8, bty="n", cex.axis=1.1, cex.lab=1.2, col=effect_col, xlim=c(xmin, xmax), ylim=c(ymin, ymax), xlab=expression(LDL~cholesterol~(italic(b[zx]))), ylab=expression(Coronary~artery~disease~(italic(b[zy])))) abline(0, gsmr_results$forward_bxy, lwd=1.5, lty=2, col="dim grey") nsnps = length(bzx[filtered_index]) for( i in 1:nsnps ) { # x axis xstart = bzx[filtered_index [i]] - bzx_se[filtered_index[i]]; xend = bzx[filtered_index[i]] + bzx_se[filtered_index[i]] ystart = bzy[filtered_index[i]]; yend = bzy[filtered_index[i]] segments(xstart, ystart, xend, yend, lwd=1.5, col=effect_col) # y axis xstart = bzx[filtered_index[i]]; xend = bzx[filtered_index[i]] ystart = bzy[filtered_index[i]] - bzy_se[filtered_index[i]]; yend = bzy[filtered_index[i]] + bzy_se[filtered_index[i]] segments(xstart, ystart, xend, yend, lwd=1.5, col=effect_col) } ```