Two-Sample MR Practical
Caroline Brito Nunes
2026-07-01
With thanks to Kaitlin Wade and Chris Flatley.
rm(list=ls())Preliminary steps
If you haven’t already installed the necessary packages with libraries, please do so!
install.packages("metafor")
install.packages("plyr")
install.packages("meta")
install.packages("rmeta")library(metafor)
library(plyr)
library(meta)
library(rmeta)Set working directory to your local directory.
For example:
setwd("/Users/uqcbrito/GGWS26/module4/Practice_2_TSMR")Practical overview
In this practical we will be testing whether an exposure (here, BMI) has a causal effect on an outcome (here, coronary heart disease, CHD).
The basic logic of a two-sample MR (the design used in this practical) is:
- Find SNPs that are strongly and independently associated with the exposure (BMI) in a large GWAS (here, GIANT/Locke et al. 2015).
- Look up the association of those same SNPs with the outcome (CHD) in an independent GWAS (here, CARDIOGRAM).
- Make sure the effect estimates in both datasets refer to the same effect allele (“harmonization”).
- For each SNP, calculate a Wald ratio: the SNP-outcome effect divided by the SNP-exposure effect. This estimates the causal effect of the exposure on the outcome.
- Combine the SNP-specific ratios into a single overall causal estimate (meta-analysis), and run sensitivity analyses (e.g. MR-Egger, weighted median, weighted mode) to test the robustness of the causal inference.
Objectives
Identify independent SNPs from BMI GWAS for use as the instrument variable of MR.
Merge and harmonize with SNPs from the CHD GWAS.
Check for palindromic SNPs and for SNPs in opposing directions.
Estimate Wald Ratio and meta analyze results.
Calculate heterogeneity statistics.
Run sensitivity analyses.
PART 1: Understanding the data and defining our BMI genetic
instruments.
1. We will be using results from the Locke et al. 2015 paper
using data from the GIANT Consortium (downloaded from https://portals.broadinstitute.org/collaboration/giant/index.php/GIANT_consortium_data_files)
Note: If you want to have a look at the full GIANT data, then
download, save to your working directory and load in the full results
(this is a fairly large file) using the following code:
giant_full <- read.table(“./GIANT_raw.uniq”, header = T)
In the Locke et al. paper, the authors describe a certain number of
SNPs that are “approximately independently associated with BMI” across
all ancestries.
Let’s have a look at these SNPs.
Note: This file was generated by taking the first rows from the Supplementary Table 8 of the Locke et al. paper.
snps <- read.csv("./Data/giant_snps_all.csv", header = T)
a. How many SNPs do authors describe as being independently
associated with BMI in all ancestries?
dim(snps)## [1] 97 14head(snps)## SNP CHR BP NearestGene Effect_Allele Other_Allele Beta SE EAF Variance N P Sig_Analysis Sig_P
## 1 rs1558902 16 52,361,075 FTO A T 0.081 0.003 0.409 0.316% 336,974 1.13e-156 European sex-combined 7.51e-153
## 2 rs6567160 18 55,980,115 MC4R C T 0.056 0.004 0.236 0.114% 339,006 6.68e-59 European sex-combined 3.93e-53
## 3 rs13021737 2 622,348 TMEM18 G A 0.060 0.004 0.830 0.103% 333,169 5.44e-54 European sex-combined 1.11e-50
## 4 rs10938397 4 44,877,284 GNPDA2 G A 0.040 0.003 0.428 0.078% 337,092 1.42e-40 European sex-combined 3.21e-38
## 5 rs543874 1 176,156,103 SEC16B G A 0.050 0.004 0.195 0.077% 339,078 2.29e-40 European sex-combined 2.62e-35
## 6 rs2207139 6 50,953,449 TFAP2B G A 0.045 0.004 0.176 0.058% 339,089 8.06e-31 European sex-combined 4.13e-29
b. What information does this file contain that are needed for
Mendelian randomization analyses?
colnames(snps)## [1] "SNP" "CHR" "BP" "NearestGene" "Effect_Allele" "Other_Allele" "Beta" "SE" "EAF"
## [10] "Variance" "N" "P" "Sig_Analysis" "Sig_P"
c. How many are associated within only Europeans?
table(snps$Sig_Analysis)##
## All Ancestries European Men European Population Based European sex-combined European Women
## 10 3 4 77 3
d. Why might it be best to use the SNPs that have been
identified as being associated with BMI in Europeans only?
2. Read in the second sheet of this file to get the
estimates of the SNPs associated with BMI in Europeans.
euro_snps <- read.csv("./Data/giant_snps_euro.csv", header = T)
dim(euro_snps)## [1] 77 12head(euro_snps)## SNP CHR BP NearestGene Effect_Allele Other_Allele Beta SE EAF Variance N P
## 1 rs1558902 16 52,361,075 FTO A T 0.082 0.003 0.415 0.325% 320,073 7.51e-153
## 2 rs6567160 18 55,980,115 MC4R C T 0.056 0.004 0.236 0.111% 321,958 3.93e-53
## 3 rs13021737 2 622,348 TMEM18 G A 0.060 0.004 0.828 0.103% 318,287 1.11e-50
## 4 rs10938397 4 44,877,284 GNPDA2 G A 0.040 0.003 0.434 0.079% 320,955 3.21e-38
## 5 rs543874 1 176,156,103 SEC16B G A 0.048 0.004 0.193 0.072% 322,008 2.62e-35
## 6 rs2207139 6 50,953,449 TFAP2B G A 0.045 0.004 0.177 0.058% 322,019 4.13e-29Note on units: Beta and SE are the change
in BMI, in standard deviation (SD) units, per copy of
the effect allele — Locke et al. inverse-normal transformed BMI before
running the GWAS, which is what puts these estimates on an SD scale
rather than raw kg/m². Later, in Part 4, we’ll compare these genetic
estimates against an observational BMI-CHD association reported per 4.56
kg/m² (i.e., that study’s SD of BMI); we use the same 4.56 kg/m² figure
to translate the SD units here into kg/m², as an approximation for
interpretability.
a. Check that these are all associated with
BMI at a conventional level of genome-wide significance.
sort(euro_snps$P)## [1] 7.51e-153 3.93e-53 1.11e-50 3.21e-38 2.62e-35 4.13e-29 5.56e-28 2.66e-26 8.15e-24 8.78e-24 3.14e-23 1.89e-22 1.48e-18 4.59e-18 1.91e-17
## [16] 5.15e-17 6.19e-17 3.28e-15 4.81e-15 1.23e-14 1.74e-14 6.61e-14 7.03e-14 5.48e-13 1.09e-12 1.31e-12 1.83e-12 2.07e-12 1.11e-11 1.14e-11
## [31] 2.07e-11 2.25e-11 2.90e-11 5.94e-11 1.15e-10 1.17e-10 1.63e-10 1.75e-10 1.83e-10 1.94e-10 2.29e-10 3.55e-10 6.33e-10 7.47e-10 7.91e-10
## [46] 8.11e-10 1.33e-09 1.92e-09 2.48e-09 2.49e-09 3.99e-09 4.56e-09 7.34e-09 7.41e-09 7.76e-09 8.45e-09 1.20e-08 1.25e-08 1.28e-08 1.39e-08
## [61] 1.48e-08 1.61e-08 1.83e-08 1.89e-08 2.02e-08 2.31e-08 2.41e-08 2.55e-08 2.60e-08 2.67e-08 2.96e-08 2.97e-08 3.19e-08 3.42e-08 3.86e-08
## [76] 4.17e-08 4.89e-08length(which(euro_snps$P<=5E-8))## [1] 77
b. Are all of these SNPs “good instruments”? What else might we
want to check to see if they are strongly and independently associated
with BMI?
3. We’re going to make sure the effect allele is the
allele that increases BMI using the effect allele and beta column.
Browse the data, are all SNPs coded so that the effect allele increases
BMI?
a. Are all SNP effects in the same direction?
euro_snps[,c("SNP","Effect_Allele","Beta","SE","P")]## SNP Effect_Allele Beta SE P
## 1 rs1558902 A 0.082 0.003 7.51e-153
## 2 rs6567160 C 0.056 0.004 3.93e-53
## 3 rs13021737 G 0.060 0.004 1.11e-50
## 4 rs10938397 G 0.040 0.003 3.21e-38
## 5 rs543874 G 0.048 0.004 2.62e-35
## 6 rs2207139 G 0.045 0.004 4.13e-29
## 7 rs11030104 A 0.041 0.004 5.56e-28
## 8 rs3101336 C 0.033 0.003 2.66e-26
## 9 rs7138803 A 0.032 0.003 8.15e-24
## 10 rs10182181 G 0.031 0.003 8.78e-24
## 11 rs3888190 A 0.031 0.003 3.14e-23
## 12 rs1516725 C 0.045 0.005 1.89e-22
## 13 rs12446632 G 0.040 0.005 1.48e-18
## 14 rs2287019 C 0.036 0.004 4.59e-18
## 15 rs16951275 T 0.031 0.004 1.91e-17
## 16 rs3817334 T 0.026 0.003 5.15e-17
## 17 rs2112347 T 0.026 0.003 6.19e-17
## 18 rs12566985 G 0.024 0.003 3.28e-15
## 19 rs3810291 A 0.028 0.004 4.81e-15
## 20 rs7141420 T 0.024 0.003 1.23e-14
## 21 rs13078960 G 0.030 0.004 1.74e-14
## 22 rs10968576 G 0.025 0.003 6.61e-14
## 23 rs17024393 C 0.066 0.009 7.03e-14
## 24 rs657452 A 0.023 0.003 5.48e-13
## 25 rs12429545 A 0.033 0.005 1.09e-12
## 26 rs12286929 G 0.022 0.003 1.31e-12
## 27 rs13107325 T 0.048 0.007 1.83e-12
## 28 rs11165643 T 0.022 0.003 2.07e-12
## 29 rs7903146 C 0.023 0.003 1.11e-11
## 30 rs10132280 C 0.023 0.003 1.14e-11
## 31 rs17405819 T 0.022 0.003 2.07e-11
## 32 rs1016287 T 0.023 0.003 2.25e-11
## 33 rs4256980 G 0.021 0.003 2.90e-11
## 34 rs17094222 C 0.025 0.004 5.94e-11
## 35 rs12401738 A 0.021 0.003 1.15e-10
## 36 rs7599312 G 0.022 0.003 1.17e-10
## 37 rs2365389 C 0.020 0.003 1.63e-10
## 38 rs205262 G 0.022 0.004 1.75e-10
## 39 rs2820292 C 0.020 0.003 1.83e-10
## 40 rs12885454 C 0.021 0.003 1.94e-10
## 41 rs12016871 T 0.030 0.005 2.29e-10
## 42 rs16851483 T 0.048 0.008 3.55e-10
## 43 rs1167827 G 0.020 0.003 6.33e-10
## 44 rs758747 T 0.023 0.004 7.47e-10
## 45 rs1928295 T 0.019 0.003 7.91e-10
## 46 rs9925964 A 0.019 0.003 8.11e-10
## 47 rs11126666 A 0.021 0.003 1.33e-09
## 48 rs2650492 A 0.021 0.004 1.92e-09
## 49 rs6804842 G 0.019 0.003 2.48e-09
## 50 rs12940622 G 0.018 0.003 2.49e-09
## 51 rs11847697 T 0.049 0.008 3.99e-09
## 52 rs4740619 T 0.018 0.003 4.56e-09
## 53 rs13191362 A 0.028 0.005 7.34e-09
## 54 rs3736485 A 0.018 0.003 7.41e-09
## 55 rs17001654 G 0.031 0.005 7.76e-09
## 56 rs11191560 C 0.031 0.005 8.45e-09
## 57 rs1528435 T 0.018 0.003 1.20e-08
## 58 rs2075650 A 0.026 0.005 1.25e-08
## 59 rs1000940 G 0.019 0.003 1.28e-08
## 60 rs2033529 G 0.019 0.003 1.39e-08
## 61 rs11583200 C 0.018 0.003 1.48e-08
## 62 rs9400239 C 0.019 0.003 1.61e-08
## 63 rs10733682 A 0.017 0.003 1.83e-08
## 64 rs11688816 G 0.017 0.003 1.89e-08
## 65 rs11057405 G 0.031 0.006 2.02e-08
## 66 rs2121279 T 0.025 0.004 2.31e-08
## 67 rs29941 G 0.018 0.003 2.41e-08
## 68 rs11727676 T 0.036 0.006 2.55e-08
## 69 rs3849570 A 0.019 0.003 2.60e-08
## 70 rs6477694 C 0.017 0.003 2.67e-08
## 71 rs7899106 G 0.040 0.007 2.96e-08
## 72 rs2176598 T 0.020 0.004 2.97e-08
## 73 rs2245368 C 0.032 0.006 3.19e-08
## 74 rs17724992 A 0.019 0.004 3.42e-08
## 75 rs7243357 T 0.022 0.004 3.86e-08
## 76 rs1808579 C 0.017 0.003 4.17e-08
## 77 rs2033732 C 0.019 0.004 4.89e-08summary(euro_snps$Beta)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01700 0.02000 0.02300 0.02847 0.03200 0.08200E.g., a Beta of 0.05 means each copy of that allele raises BMI by
0.05 SD, roughly 0.05 × 4.56 ≈ 0.23 kg/m².
PART 2: SNP lookup in the GWAS of CHD.
1. These were downloaded from the CARDIOGRAM website http://www.cardiogramplusc4d.org/ and provided for you
(however, this is quite a large file so we have truncated for you using
the following code):
CARDIOGRAM <-
read.table(“./CARDIoGRAM_GWAS_RESULTS.txt”, header = T)
BMI_SNPS_in_CARDIOGRAM <- CARDIOGRAM[CARDIOGRAM$SNP %in% euro_snps$SNP,]
top_1000 <- CARDIOGRAM[c(1:1000),]
bottom_1000 <- CARDIOGRAM[c((nrow(CARDIOGRAM)-1000):(nrow(CARDIOGRAM))),]
CARDIOGRAM_TRUNC <- rbind(BMI_SNPS_in_CARDIOGRAM, top_1000, bottom_1000)
write.table(CARDIOGRAM_TRUNC, “./CARDIOGRAM_CLEANED.txt”, quote = F,
col.names = T, row.names = F, sep = “)
a. How many SNPs are in
this truncated CARDIOGRAM dataset?
CARDIOGRAM <- read.table("./Data/CARDIOGRAM_CLEANED.txt", sep="\t", header=T, colClasses = "character")
dim(CARDIOGRAM)## [1] 2078 12head(CARDIOGRAM)## SNP chr_pos_.b36. reference_allele other_allele ref_allele_frequency pvalue het_pvalue log_odds log_odds_se N_case N_control model
## 1 rs657452 chr1:49362434 G A 0.60513584 0.2262602 0.27224751 -0.0174996 0.0144619 19723 57869 FE
## 2 rs11583200 chr1:50332407 C T 0.39664739 0.1818971 0.45532089 0.0190464 0.0142676 20247 58547 FE
## 3 rs3101336 chr1:72523773 C T 0.60256842 0.6796559 0.57600141 0.005948 0.0144044 19760 58014 FE
## 4 rs12566985 chr1:74774781 G A 0.44899777 0.3564933 0.42546079 -0.0128532 0.0139396 20806 61606 FE
## 5 rs12401738 chr1:78219349 G A 0.62028006 0.1955482 0.24955897 -0.0195047 0.0150693 18196 58218 FE
## 6 rs11165643 chr1:96696685 C T 0.41074541 0.9557793 0.82853378 -0.0007751 0.0139775 20644 61437 FENote on units: log_odds is the log odds ratio
(log OR) for CHD, per copy of the reference allele — it’s on
the log scale, not the OR scale, until it’s exponentiated
(we’ll do that once we get to results in Part 4).
b. Does this
file contain everything that is required to perform a two-sample
Mendelian randomization analysis?
You need, at minimum: the SNP identifier, effect (reference) allele, other allele, effect estimate (log OR) and its standard error, and ideally the effect allele frequency (useful for checking strand alignment and for calculating instrument strength/F-statistics).
colnames(CARDIOGRAM)## [1] "SNP" "chr_pos_.b36." "reference_allele" "other_allele" "ref_allele_frequency" "pvalue"
## [7] "het_pvalue" "log_odds" "log_odds_se" "N_case" "N_control" "model"
2. How many of the BMI SNPs are included in the
CARDIOGRAM dataset?
BMI_SNPs <- euro_snps$SNP
BMI_SNPs <- as.vector(BMI_SNPs)
matches <- unique(grep(paste(BMI_SNPs, collapse="|"), CARDIOGRAM$SNP, value=TRUE))
matches## [1] "rs657452" "rs11583200" "rs3101336" "rs12566985" "rs12401738" "rs11165643" "rs17024393" "rs543874" "rs2820292" "rs13021737" "rs10182181"
## [12] "rs11126666" "rs1016287" "rs11688816" "rs2121279" "rs1528435" "rs7599312" "rs6804842" "rs2365389" "rs3849570" "rs13078960" "rs16851483"
## [23] "rs1516725" "rs10938397" "rs17001654" "rs13107325" "rs11727676" "rs2112347" "rs205262" "rs2033529" "rs2207139" "rs9400239" "rs13191362"
## [34] "rs1167827" "rs2245368" "rs17405819" "rs2033732" "rs4740619" "rs10968576" "rs6477694" "rs1928295" "rs10733682" "rs7899106" "rs17094222"
## [45] "rs11191560" "rs7903146" "rs4256980" "rs11030104" "rs2176598" "rs3817334" "rs12286929" "rs7138803" "rs11057405" "rs12016871" "rs12429545"
## [56] "rs10132280" "rs12885454" "rs11847697" "rs7141420" "rs3736485" "rs16951275" "rs758747" "rs12446632" "rs2650492" "rs3888190" "rs9925964"
## [67] "rs1558902" "rs1000940" "rs12940622" "rs1808579" "rs7243357" "rs6567160" "rs17724992" "rs29941" "rs2075650" "rs2287019" "rs3810291"
View the data from CARDIOGRAM for our BMI SNPs.
CARDIOGRAM_BMI <- CARDIOGRAM[grepl(paste(BMI_SNPs, collapse="|"), CARDIOGRAM$SNP),]
CARDIOGRAM_BMI## SNP chr_pos_.b36. reference_allele other_allele ref_allele_frequency pvalue het_pvalue log_odds log_odds_se N_case N_control model
## 1 rs657452 chr1:49362434 G A 0.60513584 0.2262602 0.27224751 -0.0174996 0.0144619 19723 57869 FE
## 2 rs11583200 chr1:50332407 C T 0.39664739 0.1818971 0.45532089 0.0190464 0.0142676 20247 58547 FE
## 3 rs3101336 chr1:72523773 C T 0.60256842 0.6796559 0.57600141 0.005948 0.0144044 19760 58014 FE
## 4 rs12566985 chr1:74774781 G A 0.44899777 0.3564933 0.42546079 -0.0128532 0.0139396 20806 61606 FE
## 5 rs12401738 chr1:78219349 G A 0.62028006 0.1955482 0.24955897 -0.0195047 0.0150693 18196 58218 FE
## 6 rs11165643 chr1:96696685 C T 0.41074541 0.9557793 0.82853378 -0.0007751 0.0139775 20644 61437 FE
## 7 rs17024393 chr1:109956211 C T 0.03280444 0.7560841 0.13856056 0.0139407 0.0448792 14934 53634 FE
## 8 rs543874 chr1:176156103 G A 0.19669364 0.7931315 0.27599476 -0.004636 0.0176779 21134 61643 FE
## 9 rs2820292 chr1:200050910 C A 0.56965298 0.0099703 0.58433751 0.0359029 0.0139328 21198 61870 FE
## 10 rs13021737 chr2:622348 G A 0.80895378 0.3153078 0.69171591 0.0186889 0.0186117 20834 61762 FE
## 11 rs10182181 chr2:25003800 G A 0.46240927 0.9498548 0.61182332 0.0009032 0.0143625 19568 56734 FE
## 12 rs11126666 chr2:26782315 G A 0.71929726 0.1809574 0.89136444 0.0209674 0.0156729 21443 59360 FE
## 13 rs1016287 chr2:59159129 C T 0.72117719 0.102763 0.12279346 -0.0254175 0.0155782 20468 61571 FE
## 14 rs11688816 chr2:62906552 G A 0.5348187 0.9710981 0.35503175 0.0005135 0.0141725 19915 58201 FE
## 15 rs2121279 chr2:142759755 C T 0.86261919 0.8737069 0.41814919 -0.0032328 0.0203383 20892 61352 FE
## 16 rs1528435 chr2:181259207 C T 0.37962868 0.7415367 0.81685878 0.0047002 0.0142509 20460 61311 FE
## 17 rs7599312 chr2:213121476 G A 0.70337496 0.9192913 0.10235164 0.0016027 0.0158173 20703 57717 FE
## 18 rs6804842 chr3:25081441 G A 0.57965224 0.9107156 0.67473071 0.0016094 0.0143525 20403 61351 FE
## 19 rs2365389 chr3:61211502 C T 0.55986741 0.3350863 0.91057566 0.0136276 0.0141377 20635 61523 FE
## 20 rs3849570 chr3:81874802 C A 0.61947841 0.126133 0.00935007 -0.0390811 0.0255511 20556 61394 RE
## 21 rs13078960 chr3:85890280 G T 0.22900043 0.1345178 0.03395648 -0.0258172 0.0172515 20549 61416 FE
## 22 rs16851483 chr3:142758126 G T 0.93429797 0.5088742 0.05493564 -0.0188109 0.0284759 19390 59683 FE
## 23 rs1516725 chr3:187306698 C T 0.8590188 0.1629709 0.22517731 -0.0279354 0.0200233 20476 61368 FE
## 24 rs10938397 chr4:44877284 G A 0.41363733 0.0278084 0.68643434 0.0349131 0.0158698 16024 55685 FE
## 25 rs17001654 chr4:77348592 C G 0.83729352 0.5035566 0.53373783 -0.0126514 0.0189136 21412 61536 FE
## 26 rs13107325 chr4:103407732 C T 0.97601496 0.9074287 0.45904911 -0.0049381 0.0424663 14269 35628 FE
## 27 rs11727676 chr4:145878514 C T 0.08996777 0.4402455 0.48987311 0.0253447 0.0328393 15668 53910 FE
## 28 rs2112347 chr5:75050998 G T 0.36318919 0.5414854 0.58598695 -0.0088627 0.0145155 20579 61408 FE
## 29 rs205262 chr6:34671142 G A 0.25943383 0.0001102 0.4801711 0.0614484 0.0158908 20409 58712 FE
## 30 rs2033529 chr6:40456631 G A 0.26951838 0.6939016 0.9186536 -0.0061959 0.015743 18908 56847 FE
## 31 rs2207139 chr6:50953449 G A 0.17342856 0.3550604 0.6093131 0.0170642 0.0184515 20539 61411 FE
## 32 rs9400239 chr6:109084356 C T 0.71596318 0.0227834 0.07735518 0.0348705 0.0153139 20460 61373 FE
## 33 rs13191362 chr6:162953340 G A 0.12900275 0.5722799 0.54836822 -0.0116552 0.0206397 21401 61861 FE
## 34 rs1167827 chr7:75001105 G A 0.59298726 0.2370525 0.84489268 0.0190735 0.0161314 15647 51329 FE
## 35 rs2245368 chr7:76446079 C T 0.16123262 0.5501828 0.14094314 0.0303632 0.0508183 3340 18930 FE
## 36 rs17405819 chr8:76969139 C T 0.32156614 0.9167122 0.22795208 -0.0015648 0.0149631 20724 61620 FE
## 37 rs2033732 chr8:85242264 C T 0.75652234 0.1335914 0.04706672 -0.0253631 0.0169077 18943 59038 FE
## 38 rs4740619 chr9:15624326 C T 0.46589474 0.2675504 0.95611066 -0.0154694 0.0139524 20779 61682 FE
## 39 rs10968576 chr9:28404339 G A 0.31119502 0.5871441 0.53064937 0.00811 0.0149361 20590 61535 FE
## 40 rs6477694 chr9:110972163 C T 0.3768822 0.0596142 0.73673365 -0.026955 0.0143101 21572 59220 FE
## 41 rs1928295 chr9:119418304 C T 0.44991607 0.970387 0.77734102 0.00052 0.0140084 21633 59309 FE
## 42 rs10733682 chr9:128500735 G A 0.52449229 0.1829671 0.5076359 -0.0212797 0.0159796 15670 51340 FE
## 43 rs7899106 chr10:87400884 G A 0.06251531 0.5857598 0.35079691 -0.0163075 0.0299225 20393 60788 FE
## 44 rs17094222 chr10:102385430 C T 0.22353262 0.72768 0.77948036 -0.0059791 0.0171708 21528 62065 FE
## 45 rs11191560 chr10:104859028 C T 0.11109357 0.0001767 0.53911165 -0.0966196 0.0257636 20098 60493 FE
## 46 rs7903146 chr10:114748339 C T 0.69817348 0.2034467 0.36081083 -0.0192779 0.0151581 19910 58247 FE
## 47 rs4256980 chr11:8630515 C G 0.34826703 0.299491 0.02913706 0.0150656 0.0145207 20711 61646 FE
## 48 rs11030104 chr11:27641093 G A 0.19185393 0.0496195 0.36961605 -0.034644 0.0176464 21638 59361 FE
## 49 rs2176598 chr11:43820854 C T 0.74580823 0.1264641 0.05596993 -0.0247371 0.0161872 20452 61406 FE
## 50 rs3817334 chr11:47607569 C T 0.58779101 0.1194239 0.81445328 -0.0219146 0.0140731 21525 61809 FE
## 51 rs12286929 chr11:114527614 G A 0.50744154 0.9892883 0.46198211 0.0001851 0.0137839 21869 62138 FE
## 52 rs7138803 chr12:48533735 G A 0.61809264 0.9656206 0.9396865 0.0006419 0.0148916 19151 56023 FE
## 53 rs11057405 chr12:121347850 G A 0.92089578 0.4796771 0.59224188 -0.0214346 0.0303253 13629 49312 FE
## 54 rs12016871 chr13:26915782 C T 0.80407573 0.3017645 0.36210717 -0.0181606 0.0175863 20330 61194 FE
## 55 rs12429545 chr13:53000207 G A 0.845308 0.0396211 0.93056094 -0.0480905 0.0233712 16413 52969 FE
## 56 rs10132280 chr14:24998019 C A 0.66888613 0.302347 0.2972794 0.0169764 0.0164594 16945 58264 FE
## 57 rs12885454 chr14:28806589 C A 0.631088 0.9719424 0.66514854 -0.0005252 0.0149317 19926 60456 FE
## 58 rs11847697 chr14:29584863 C T 0.96069173 0.0476032 0.26232231 -0.0707465 0.0357144 20155 61186 FE
## 59 rs7141420 chr14:78969207 C T 0.48985341 0.8946607 0.95298609 -0.0018271 0.0137989 21258 62144 FE
## 60 rs3736485 chr15:49535902 G A 0.55461993 0.5162199 0.10728276 0.0091099 0.0140328 20940 61375 FE
## 61 rs16951275 chr15:65864222 C T 0.235351 0.0024846 0.22809659 -0.0513109 0.0169611 20781 61718 FE
## 62 rs758747 chr16:3567359 C T 0.74188258 0.0076412 0.91893399 0.0613185 0.022987 8437 44300 FE
## 63 rs12446632 chr16:19842890 G A 0.84624483 0.6290302 0.33198375 -0.0101474 0.0210052 20232 58533 FE
## 64 rs2650492 chr16:28240912 G A 0.70112768 0.5639757 0.05832099 0.0126597 0.0219426 8628 44465 FE
## 65 rs3888190 chr16:28796987 C A 0.57882566 0.8709159 0.00489802 -0.0042697 0.026276 19770 59987 RE
## 66 rs9925964 chr16:31037396 G A 0.38141338 0.6696603 0.2066786 0.0060442 0.0141679 20484 61554 FE
## 67 rs1558902 chr16:52361075 T A 0.58450108 0.0205892 0.26380736 -0.0324866 0.0140305 20777 61645 FE
## 68 rs1000940 chr17:5223976 G A 0.31640352 0.5722296 0.79718646 -0.0091901 0.0162722 16633 54057 FE
## 69 rs12940622 chr17:76230166 G A 0.55719561 0.6813264 0.13926311 -0.005711 0.0139071 20824 61691 FE
## 70 rs1808579 chr18:19358886 C T 0.54056958 0.500074 0.17404987 0.0096342 0.0142861 19215 57225 FE
## 71 rs7243357 chr18:55034299 G T 0.19721707 0.5943844 0.73397076 -0.0098297 0.0184597 20404 61241 FE
## 72 rs6567160 chr18:55980115 C T 0.26395759 0.13325 0.33646334 0.024367 0.0162294 20445 61296 FE
## 73 rs17724992 chr19:18315825 G A 0.25712419 0.002919 0.87146506 -0.0493236 0.016573 19109 59755 FE
## 74 rs29941 chr19:39001372 G A 0.67182131 0.7661974 0.88467459 -0.0044497 0.0149645 20307 58647 FE
## 75 rs2075650 chr19:50087459 G A 0.1705933 0.1957039 0.33474452 0.0378692 0.0292677 6913 27904 FE
## 76 rs2287019 chr19:50894012 C T 0.79836678 0.0618941 0.086903 0.045534 0.0243882 8640 44554 FE
## 77 rs3810291 chr19:52260843 G A 0.33939316 0.0303086 0.59137189 -0.0422661 0.0195131 11287 50470 FE
3. Merge the GIANT and CARDIOGRAM SNP summary
associations.
First, make sure the column headings are
easy to understand (i.e., add “BMI” and “CHD” onto the respective
datasets).
colnames(euro_snps) <- paste("BMI", colnames(euro_snps), sep = "_")
head(euro_snps)## BMI_SNP BMI_CHR BMI_BP BMI_NearestGene BMI_Effect_Allele BMI_Other_Allele BMI_Beta BMI_SE BMI_EAF BMI_Variance BMI_N BMI_P
## 1 rs1558902 16 52,361,075 FTO A T 0.082 0.003 0.415 0.325% 320,073 7.51e-153
## 2 rs6567160 18 55,980,115 MC4R C T 0.056 0.004 0.236 0.111% 321,958 3.93e-53
## 3 rs13021737 2 622,348 TMEM18 G A 0.060 0.004 0.828 0.103% 318,287 1.11e-50
## 4 rs10938397 4 44,877,284 GNPDA2 G A 0.040 0.003 0.434 0.079% 320,955 3.21e-38
## 5 rs543874 1 176,156,103 SEC16B G A 0.048 0.004 0.193 0.072% 322,008 2.62e-35
## 6 rs2207139 6 50,953,449 TFAP2B G A 0.045 0.004 0.177 0.058% 322,019 4.13e-29colnames(CARDIOGRAM) <- paste("CHD", colnames(CARDIOGRAM), sep = "_")
head(CARDIOGRAM)## CHD_SNP CHD_chr_pos_.b36. CHD_reference_allele CHD_other_allele CHD_ref_allele_frequency CHD_pvalue CHD_het_pvalue CHD_log_odds CHD_log_odds_se
## 1 rs657452 chr1:49362434 G A 0.60513584 0.2262602 0.27224751 -0.0174996 0.0144619
## 2 rs11583200 chr1:50332407 C T 0.39664739 0.1818971 0.45532089 0.0190464 0.0142676
## 3 rs3101336 chr1:72523773 C T 0.60256842 0.6796559 0.57600141 0.005948 0.0144044
## 4 rs12566985 chr1:74774781 G A 0.44899777 0.3564933 0.42546079 -0.0128532 0.0139396
## 5 rs12401738 chr1:78219349 G A 0.62028006 0.1955482 0.24955897 -0.0195047 0.0150693
## 6 rs11165643 chr1:96696685 C T 0.41074541 0.9557793 0.82853378 -0.0007751 0.0139775
## CHD_N_case CHD_N_control CHD_model
## 1 19723 57869 FE
## 2 20247 58547 FE
## 3 19760 58014 FE
## 4 20806 61606 FE
## 5 18196 58218 FE
## 6 20644 61437 FEmerged <- merge(euro_snps,CARDIOGRAM, by.x="BMI_SNP", by.y="CHD_SNP")
dim(merged)## [1] 77 23head(merged)## BMI_SNP BMI_CHR BMI_BP BMI_NearestGene BMI_Effect_Allele BMI_Other_Allele BMI_Beta BMI_SE BMI_EAF BMI_Variance BMI_N BMI_P CHD_chr_pos_.b36.
## 1 rs1000940 17 5,223,976 RABEP1 G A 0.019 0.003 0.320 0.016% 321,836 1.28e-08 chr17:5223976
## 2 rs10132280 14 24,998,019 STXBP6 C A 0.023 0.003 0.682 0.023% 321,797 1.14e-11 chr14:24998019
## 3 rs1016287 2 59,159,129 LINC01122 T C 0.023 0.003 0.287 0.021% 321,969 2.25e-11 chr2:59159129
## 4 rs10182181 2 25,003,800 ADCY3 G A 0.031 0.003 0.462 0.047% 321,759 8.78e-24 chr2:25003800
## 5 rs10733682 9 128,500,735 LMX1B A G 0.017 0.003 0.478 0.015% 320,727 1.83e-08 chr9:128500735
## 6 rs10938397 4 44,877,284 GNPDA2 G A 0.040 0.003 0.434 0.079% 320,955 3.21e-38 chr4:44877284
## CHD_reference_allele CHD_other_allele CHD_ref_allele_frequency CHD_pvalue CHD_het_pvalue CHD_log_odds CHD_log_odds_se CHD_N_case CHD_N_control CHD_model
## 1 G A 0.31640352 0.5722296 0.79718646 -0.0091901 0.0162722 16633 54057 FE
## 2 C A 0.66888613 0.302347 0.2972794 0.0169764 0.0164594 16945 58264 FE
## 3 C T 0.72117719 0.102763 0.12279346 -0.0254175 0.0155782 20468 61571 FE
## 4 G A 0.46240927 0.9498548 0.61182332 0.0009032 0.0143625 19568 56734 FE
## 5 G A 0.52449229 0.1829671 0.5076359 -0.0212797 0.0159796 15670 51340 FE
## 6 G A 0.41363733 0.0278084 0.68643434 0.0349131 0.0158698 16024 55685 FEPART 3: Harmonizing the effect alleles in the BMI and CHD
datasets.
1. Make sure that the effect alleles in the CARDIOGRAM and
GIANT datasets are the same. We want the CARDIOGRAM effect allele to be
the allele that increases BMI.
But be careful of
palindromic SNPs or SNPs on different strands.
First we need to see
whether the effect alleles are the same.
Browse the data.
merged[,c("BMI_SNP", "BMI_Effect_Allele","CHD_reference_allele","BMI_Other_Allele","CHD_other_allele","BMI_EAF","CHD_ref_allele_frequency")]## BMI_SNP BMI_Effect_Allele CHD_reference_allele BMI_Other_Allele CHD_other_allele BMI_EAF CHD_ref_allele_frequency
## 1 rs1000940 G G A A 0.320 0.31640352
## 2 rs10132280 C C A A 0.682 0.66888613
## 3 rs1016287 T C C T 0.287 0.72117719
## 4 rs10182181 G G A A 0.462 0.46240927
## 5 rs10733682 A G G A 0.478 0.52449229
## 6 rs10938397 G G A A 0.434 0.41363733
## 7 rs10968576 G G A A 0.320 0.31119502
## 8 rs11030104 A G G A 0.792 0.19185393
## 9 rs11057405 G G A A 0.901 0.92089578
## 10 rs11126666 A G G A 0.283 0.71929726
## 11 rs11165643 T C C T 0.583 0.41074541
## 12 rs11191560 C C T T 0.089 0.11109357
## 13 rs11583200 C C T T 0.396 0.39664739
## 14 rs1167827 G G A A 0.553 0.59298726
## 15 rs11688816 G G A A 0.525 0.5348187
## 16 rs11727676 T C C T 0.910 0.08996777
## 17 rs11847697 T C C T 0.042 0.96069173
## 18 rs12016871 T C C T 0.203 0.80407573
## 19 rs12286929 G G A A 0.523 0.50744154
## 20 rs12401738 A G G A 0.352 0.62028006
## 21 rs12429545 A G G A 0.133 0.845308
## 22 rs12446632 G G A A 0.865 0.84624483
## 23 rs12566985 G G A A 0.446 0.44899777
## 24 rs12885454 C C A A 0.642 0.631088
## 25 rs12940622 G G A A 0.575 0.55719561
## 26 rs13021737 G G A A 0.828 0.80895378
## 27 rs13078960 G G T T 0.196 0.22900043
## 28 rs13107325 T C C T 0.072 0.97601496
## 29 rs13191362 A G G A 0.879 0.12900275
## 30 rs1516725 C C T T 0.872 0.8590188
## 31 rs1528435 T C C T 0.631 0.37962868
## 32 rs1558902 A T T A 0.415 0.58450108
## 33 rs16851483 T G G T 0.066 0.93429797
## 34 rs16951275 T C C T 0.784 0.235351
## 35 rs17001654 G C C G 0.153 0.83729352
## 36 rs17024393 C C T T 0.040 0.03280444
## 37 rs17094222 C C T T 0.211 0.22353262
## 38 rs17405819 T C C T 0.700 0.32156614
## 39 rs17724992 A G G A 0.746 0.25712419
## 40 rs1808579 C C T T 0.534 0.54056958
## 41 rs1928295 T C C T 0.548 0.44991607
## 42 rs2033529 G G A A 0.293 0.26951838
## 43 rs2033732 C C T T 0.747 0.75652234
## 44 rs205262 G G A A 0.273 0.25943383
## 45 rs2075650 A G G A 0.848 0.1705933
## 46 rs2112347 T G G T 0.629 0.36318919
## 47 rs2121279 T C C T 0.152 0.86261919
## 48 rs2176598 T C C T 0.251 0.74580823
## 49 rs2207139 G G A A 0.177 0.17342856
## 50 rs2245368 C C T T 0.180 0.16123262
## 51 rs2287019 C C T T 0.804 0.79836678
## 52 rs2365389 C C T T 0.582 0.55986741
## 53 rs2650492 A G G A 0.303 0.70112768
## 54 rs2820292 C C A A 0.555 0.56965298
## 55 rs29941 G G A A 0.669 0.67182131
## 56 rs3101336 C C T T 0.613 0.60256842
## 57 rs3736485 A G G A 0.454 0.55461993
## 58 rs3810291 A G G A 0.666 0.33939316
## 59 rs3817334 T C C T 0.407 0.58779101
## 60 rs3849570 A C C A 0.359 0.61947841
## 61 rs3888190 A C C A 0.403 0.57882566
## 62 rs4256980 G C C G 0.646 0.34826703
## 63 rs4740619 T C C T 0.542 0.46589474
## 64 rs543874 G G A A 0.193 0.19669364
## 65 rs6477694 C C T T 0.365 0.3768822
## 66 rs6567160 C C T T 0.236 0.26395759
## 67 rs657452 A G G A 0.394 0.60513584
## 68 rs6804842 G G A A 0.575 0.57965224
## 69 rs7138803 A G G A 0.384 0.61809264
## 70 rs7141420 T C C T 0.527 0.48985341
## 71 rs7243357 T G G T 0.812 0.19721707
## 72 rs758747 T C C T 0.265 0.74188258
## 73 rs7599312 G G A A 0.724 0.70337496
## 74 rs7899106 G G A A 0.052 0.06251531
## 75 rs7903146 C C T T 0.713 0.69817348
## 76 rs9400239 C C T T 0.688 0.71596318
## 77 rs9925964 A G G A 0.620 0.38141338
a. How can we tell if the CARDIOGRAM and GIANT SNPs are coded
using the same reference strand?
Compare the effect and other allele pairs at each SNP across the
two datasets. If, for a non-palindromic SNP, the two “effect alleles”
don’t match and neither do the complementary bases, one dataset may be
reporting alleles on the opposite DNA strand.
b. Are
CARDIOGRAM and the GIANT SNPs coded using the same reference strand?
c. Are there any palindromic SNPs?
Palindromic SNPs (A/T or G/C) are ambiguous because the same two letters are used to describe the alleles regardless of which DNA strand a dataset happens to report them on — so you cannot tell just from the letters whether two datasets agree on which allele is which.
palindromic_at<-subset(merged,BMI_Effect_Allele %in% "A" & BMI_Other_Allele %in% "T")
palindromic_ta<-subset(merged,BMI_Effect_Allele %in% "T" & BMI_Other_Allele %in% "A")
palindromic_gc<-subset(merged,BMI_Effect_Allele %in% "G" & BMI_Other_Allele %in% "C")
palindromic_cg<-subset(merged,BMI_Effect_Allele %in% "C" & BMI_Other_Allele %in% "G")
dim(palindromic_at)## [1] 1 23dim(palindromic_ta)## [1] 0 23dim(palindromic_gc)## [1] 2 23dim(palindromic_cg)## [1] 0 23
d. How can we tell whether the effect alleles are the same in
both datasets for palindromic SNPs (i.e., the allele that increases BMI
is the same as the reference allele in CARDIOGRAM)?
For palindromic SNPs, compare the effect allele frequency (EAF)
between the two datasets. If a SNP’s EAF is close to 0.5 in either
dataset it remains ambiguous and is usually best excluded; otherwise, a
similar EAF in both datasets (e.g., both close to 0.2, or both close to
0.8) implies the alleles are aligned on the same strand, while a “mirror
image” EAF (e.g., 0.2 vs 0.8) implies they are on opposite strands and
need flipping.
2. Make sure the CARDIOGRAM log
odds ratio reflects the allele that increases BMI in the GIANT
data.
First, find the positions of SNPs with different
effect alleles.
head(merged)## BMI_SNP BMI_CHR BMI_BP BMI_NearestGene BMI_Effect_Allele BMI_Other_Allele BMI_Beta BMI_SE BMI_EAF BMI_Variance BMI_N BMI_P CHD_chr_pos_.b36.
## 1 rs1000940 17 5,223,976 RABEP1 G A 0.019 0.003 0.320 0.016% 321,836 1.28e-08 chr17:5223976
## 2 rs10132280 14 24,998,019 STXBP6 C A 0.023 0.003 0.682 0.023% 321,797 1.14e-11 chr14:24998019
## 3 rs1016287 2 59,159,129 LINC01122 T C 0.023 0.003 0.287 0.021% 321,969 2.25e-11 chr2:59159129
## 4 rs10182181 2 25,003,800 ADCY3 G A 0.031 0.003 0.462 0.047% 321,759 8.78e-24 chr2:25003800
## 5 rs10733682 9 128,500,735 LMX1B A G 0.017 0.003 0.478 0.015% 320,727 1.83e-08 chr9:128500735
## 6 rs10938397 4 44,877,284 GNPDA2 G A 0.040 0.003 0.434 0.079% 320,955 3.21e-38 chr4:44877284
## CHD_reference_allele CHD_other_allele CHD_ref_allele_frequency CHD_pvalue CHD_het_pvalue CHD_log_odds CHD_log_odds_se CHD_N_case CHD_N_control CHD_model
## 1 G A 0.31640352 0.5722296 0.79718646 -0.0091901 0.0162722 16633 54057 FE
## 2 C A 0.66888613 0.302347 0.2972794 0.0169764 0.0164594 16945 58264 FE
## 3 C T 0.72117719 0.102763 0.12279346 -0.0254175 0.0155782 20468 61571 FE
## 4 G A 0.46240927 0.9498548 0.61182332 0.0009032 0.0143625 19568 56734 FE
## 5 G A 0.52449229 0.1829671 0.5076359 -0.0212797 0.0159796 15670 51340 FE
## 6 G A 0.41363733 0.0278084 0.68643434 0.0349131 0.0158698 16024 55685 FEeffect_diff <- which(merged$BMI_Effect_Allele != merged$CHD_reference_allele) # The position of SNPs where effect alleles are different
merged[effect_diff,c("BMI_SNP", "BMI_Effect_Allele","CHD_reference_allele","BMI_Other_Allele","CHD_other_allele","BMI_EAF","CHD_ref_allele_frequency")]## BMI_SNP BMI_Effect_Allele CHD_reference_allele BMI_Other_Allele CHD_other_allele BMI_EAF CHD_ref_allele_frequency
## 3 rs1016287 T C C T 0.287 0.72117719
## 5 rs10733682 A G G A 0.478 0.52449229
## 8 rs11030104 A G G A 0.792 0.19185393
## 10 rs11126666 A G G A 0.283 0.71929726
## 11 rs11165643 T C C T 0.583 0.41074541
## 16 rs11727676 T C C T 0.910 0.08996777
## 17 rs11847697 T C C T 0.042 0.96069173
## 18 rs12016871 T C C T 0.203 0.80407573
## 20 rs12401738 A G G A 0.352 0.62028006
## 21 rs12429545 A G G A 0.133 0.845308
## 28 rs13107325 T C C T 0.072 0.97601496
## 29 rs13191362 A G G A 0.879 0.12900275
## 31 rs1528435 T C C T 0.631 0.37962868
## 32 rs1558902 A T T A 0.415 0.58450108
## 33 rs16851483 T G G T 0.066 0.93429797
## 34 rs16951275 T C C T 0.784 0.235351
## 35 rs17001654 G C C G 0.153 0.83729352
## 38 rs17405819 T C C T 0.700 0.32156614
## 39 rs17724992 A G G A 0.746 0.25712419
## 41 rs1928295 T C C T 0.548 0.44991607
## 45 rs2075650 A G G A 0.848 0.1705933
## 46 rs2112347 T G G T 0.629 0.36318919
## 47 rs2121279 T C C T 0.152 0.86261919
## 48 rs2176598 T C C T 0.251 0.74580823
## 53 rs2650492 A G G A 0.303 0.70112768
## 57 rs3736485 A G G A 0.454 0.55461993
## 58 rs3810291 A G G A 0.666 0.33939316
## 59 rs3817334 T C C T 0.407 0.58779101
## 60 rs3849570 A C C A 0.359 0.61947841
## 61 rs3888190 A C C A 0.403 0.57882566
## 62 rs4256980 G C C G 0.646 0.34826703
## 63 rs4740619 T C C T 0.542 0.46589474
## 67 rs657452 A G G A 0.394 0.60513584
## 69 rs7138803 A G G A 0.384 0.61809264
## 70 rs7141420 T C C T 0.527 0.48985341
## 71 rs7243357 T G G T 0.812 0.19721707
## 72 rs758747 T C C T 0.265 0.74188258
## 77 rs9925964 A G G A 0.620 0.38141338
a. How many SNPs have effect alleles that are coded in the
opposite direction?
b. Where the effect alleles are different,
flip the direction of the log odds ratio by multiplying it by -1.
This works because for a log OR, flipping which allele is treated
as the “effect” allele simply reverses the sign of the effect (a
doubling in risk for allele X is equivalent to a halving of risk for the
other allele). After this step, CHD_flip_log_odds is the
log OR for CHD per copy of the BMI-increasing allele,
matching the direction of BMI_Beta.
If you
want, you can also generate new columns that reflect the effect allele
change but this isn’t used in the causal estimate.
merged$CHD_flip_log_odds <- as.numeric(merged$CHD_log_odds) # Make log odds ratio numeric
merged$CHD_log_odds_se <- as.numeric(merged$CHD_log_odds_se) # Make standard error numeric
merged$CHD_flip_log_odds[effect_diff] <- merged$CHD_flip_log_odds[effect_diff]*(-1)
head(merged)## BMI_SNP BMI_CHR BMI_BP BMI_NearestGene BMI_Effect_Allele BMI_Other_Allele BMI_Beta BMI_SE BMI_EAF BMI_Variance BMI_N BMI_P CHD_chr_pos_.b36.
## 1 rs1000940 17 5,223,976 RABEP1 G A 0.019 0.003 0.320 0.016% 321,836 1.28e-08 chr17:5223976
## 2 rs10132280 14 24,998,019 STXBP6 C A 0.023 0.003 0.682 0.023% 321,797 1.14e-11 chr14:24998019
## 3 rs1016287 2 59,159,129 LINC01122 T C 0.023 0.003 0.287 0.021% 321,969 2.25e-11 chr2:59159129
## 4 rs10182181 2 25,003,800 ADCY3 G A 0.031 0.003 0.462 0.047% 321,759 8.78e-24 chr2:25003800
## 5 rs10733682 9 128,500,735 LMX1B A G 0.017 0.003 0.478 0.015% 320,727 1.83e-08 chr9:128500735
## 6 rs10938397 4 44,877,284 GNPDA2 G A 0.040 0.003 0.434 0.079% 320,955 3.21e-38 chr4:44877284
## CHD_reference_allele CHD_other_allele CHD_ref_allele_frequency CHD_pvalue CHD_het_pvalue CHD_log_odds CHD_log_odds_se CHD_N_case CHD_N_control CHD_model
## 1 G A 0.31640352 0.5722296 0.79718646 -0.0091901 0.0162722 16633 54057 FE
## 2 C A 0.66888613 0.302347 0.2972794 0.0169764 0.0164594 16945 58264 FE
## 3 C T 0.72117719 0.102763 0.12279346 -0.0254175 0.0155782 20468 61571 FE
## 4 G A 0.46240927 0.9498548 0.61182332 0.0009032 0.0143625 19568 56734 FE
## 5 G A 0.52449229 0.1829671 0.5076359 -0.0212797 0.0159796 15670 51340 FE
## 6 G A 0.41363733 0.0278084 0.68643434 0.0349131 0.0158698 16024 55685 FE
## CHD_flip_log_odds
## 1 -0.0091901
## 2 0.0169764
## 3 0.0254175
## 4 0.0009032
## 5 0.0212797
## 6 0.0349131dim(merged)## [1] 77 24
c. Check that all of the effect estimates have been flipped
appropriately.
merged[effect_diff,c("BMI_SNP", "BMI_Effect_Allele","CHD_reference_allele", "CHD_log_odds", "CHD_flip_log_odds")]## BMI_SNP BMI_Effect_Allele CHD_reference_allele CHD_log_odds CHD_flip_log_odds
## 3 rs1016287 T C -0.0254175 0.0254175
## 5 rs10733682 A G -0.0212797 0.0212797
## 8 rs11030104 A G -0.034644 0.0346440
## 10 rs11126666 A G 0.0209674 -0.0209674
## 11 rs11165643 T C -0.0007751 0.0007751
## 16 rs11727676 T C 0.0253447 -0.0253447
## 17 rs11847697 T C -0.0707465 0.0707465
## 18 rs12016871 T C -0.0181606 0.0181606
## 20 rs12401738 A G -0.0195047 0.0195047
## 21 rs12429545 A G -0.0480905 0.0480905
## 28 rs13107325 T C -0.0049381 0.0049381
## 29 rs13191362 A G -0.0116552 0.0116552
## 31 rs1528435 T C 0.0047002 -0.0047002
## 32 rs1558902 A T -0.0324866 0.0324866
## 33 rs16851483 T G -0.0188109 0.0188109
## 34 rs16951275 T C -0.0513109 0.0513109
## 35 rs17001654 G C -0.0126514 0.0126514
## 38 rs17405819 T C -0.0015648 0.0015648
## 39 rs17724992 A G -0.0493236 0.0493236
## 41 rs1928295 T C 0.00052 -0.0005200
## 45 rs2075650 A G 0.0378692 -0.0378692
## 46 rs2112347 T G -0.0088627 0.0088627
## 47 rs2121279 T C -0.0032328 0.0032328
## 48 rs2176598 T C -0.0247371 0.0247371
## 53 rs2650492 A G 0.0126597 -0.0126597
## 57 rs3736485 A G 0.0091099 -0.0091099
## 58 rs3810291 A G -0.0422661 0.0422661
## 59 rs3817334 T C -0.0219146 0.0219146
## 60 rs3849570 A C -0.0390811 0.0390811
## 61 rs3888190 A C -0.0042697 0.0042697
## 62 rs4256980 G C 0.0150656 -0.0150656
## 63 rs4740619 T C -0.0154694 0.0154694
## 67 rs657452 A G -0.0174996 0.0174996
## 69 rs7138803 A G 0.0006419 -0.0006419
## 70 rs7141420 T C -0.0018271 0.0018271
## 71 rs7243357 T G -0.0098297 0.0098297
## 72 rs758747 T C 0.0613185 -0.0613185
## 77 rs9925964 A G 0.0060442 -0.0060442merged[-effect_diff,c("BMI_SNP", "BMI_Effect_Allele","CHD_reference_allele", "CHD_log_odds", "CHD_flip_log_odds")]## BMI_SNP BMI_Effect_Allele CHD_reference_allele CHD_log_odds CHD_flip_log_odds
## 1 rs1000940 G G -0.0091901 -0.0091901
## 2 rs10132280 C C 0.0169764 0.0169764
## 4 rs10182181 G G 0.0009032 0.0009032
## 6 rs10938397 G G 0.0349131 0.0349131
## 7 rs10968576 G G 0.00811 0.0081100
## 9 rs11057405 G G -0.0214346 -0.0214346
## 12 rs11191560 C C -0.0966196 -0.0966196
## 13 rs11583200 C C 0.0190464 0.0190464
## 14 rs1167827 G G 0.0190735 0.0190735
## 15 rs11688816 G G 0.0005135 0.0005135
## 19 rs12286929 G G 0.0001851 0.0001851
## 22 rs12446632 G G -0.0101474 -0.0101474
## 23 rs12566985 G G -0.0128532 -0.0128532
## 24 rs12885454 C C -0.0005252 -0.0005252
## 25 rs12940622 G G -0.005711 -0.0057110
## 26 rs13021737 G G 0.0186889 0.0186889
## 27 rs13078960 G G -0.0258172 -0.0258172
## 30 rs1516725 C C -0.0279354 -0.0279354
## 36 rs17024393 C C 0.0139407 0.0139407
## 37 rs17094222 C C -0.0059791 -0.0059791
## 40 rs1808579 C C 0.0096342 0.0096342
## 42 rs2033529 G G -0.0061959 -0.0061959
## 43 rs2033732 C C -0.0253631 -0.0253631
## 44 rs205262 G G 0.0614484 0.0614484
## 49 rs2207139 G G 0.0170642 0.0170642
## 50 rs2245368 C C 0.0303632 0.0303632
## 51 rs2287019 C C 0.045534 0.0455340
## 52 rs2365389 C C 0.0136276 0.0136276
## 54 rs2820292 C C 0.0359029 0.0359029
## 55 rs29941 G G -0.0044497 -0.0044497
## 56 rs3101336 C C 0.005948 0.0059480
## 64 rs543874 G G -0.004636 -0.0046360
## 65 rs6477694 C C -0.026955 -0.0269550
## 66 rs6567160 C C 0.024367 0.0243670
## 68 rs6804842 G G 0.0016094 0.0016094
## 73 rs7599312 G G 0.0016027 0.0016027
## 74 rs7899106 G G -0.0163075 -0.0163075
## 75 rs7903146 C C -0.0192779 -0.0192779
## 76 rs9400239 C C 0.0348705 0.0348705
d. Check that the effect allele frequencies are correlated.
If the two datasets are describing the same allele at each SNP,
their allele frequencies should be strongly, positively correlated
across SNPs. A weak or even negative correlation before harmonization
(which improves after flipping) is a useful confirmation that the
harmonization step has worked correctly.
Check correlation
of effect allele frequency between BMI and CARDIOGRAM datasets before
harmonising alleles.
merged$BMI_EAF <- as.numeric(merged$BMI_EAF)
merged$CHD_ref_allele_frequency <- as.numeric(merged$CHD_ref_allele_frequency)
cor(merged$BMI_EAF,merged$CHD_ref_allele_frequency)## [1] -0.007221126
Check correlation of effect allele frequency between BMI and
CARDIOGRAM datasets after harmonising alleles
merged$CHD_ref_allele_frequency[effect_diff] <- 1-merged$CHD_ref_allele_frequency[effect_diff]
cor(merged$BMI_EAF,merged$CHD_ref_allele_frequency)## [1] 0.9980439
e. What happened and why?
The correlation should increase (move closer to +1) after
harmonization. Before flipping, some SNPs had their frequency reported
for the opposite allele, which weakens or reverses the correlation with
BMI_EAF. Once we flip the CARDIOGRAM allele frequency for
the same SNPs whose log OR we flipped (frequency of the opposite allele
= 1 − original frequency), both datasets are describing the same,
BMI-increasing allele, and the frequencies line up as expected.
PART 4: Estimate the effect of BMI on CHD.
1. Estimate the Wald ratios for each SNP and their delta approximated standard errors.
Reminder on units: gp/segp are in
SD-BMI units (1 SD ≈ 4.56 kg/m²); gd/segd are
in log OR units for CHD. The resulting wald_ratio is
therefore the log OR for CHD per 1-SD (≈4.56 kg/m²) increase in
BMI for each SNP.
gp <- merged$BMI_Beta # The effect of the SNP on BMI (SD-BMI units per effect allele; 1 SD ~ 4.56 kg/m^2)
segp <- merged$BMI_SE # The standard error of the SNP effect on BMI (SD-BMI units)
gd <- merged$CHD_flip_log_odds # The log odds ratio for CHD (that were harmonized to reflect an increase in BMI), log OR per effect allele
segd <- merged$CHD_log_odds_se # Standard error of the log odds ratio (log OR scale)
wald_ratio <- gd/gp # The log odds ratio of CHD per unit (1-SD) change in BMI
Cov <- 0 # Only required when the SNP-BMI and SNP-CHD associations are estimated in the same participants (therefore for two-sample MR with non-overlapping samples, this is set to 0)
wald_ratio_se <- sqrt((segd^2/gp^2) + (gd^2/gp^4)*segp^2 - 2*(gd/gp^3)*Cov) # Delta approximated standard error of the wald ratio; see Thomas, D. C., Lawlor, D. a, & Thompson, J. R. (2007). Re: Estimation of bias in nongenetic observational studies using Mendelian triangulation by Bautista et al. Annals of Epidemiology, 17(7), 5113. doi:10.1016/j.annepidem.2006.12.005
z <- wald_ratio/wald_ratio_se # Z statistic for the wald ratio
p <- 2*pnorm(abs(z) ,lower.tail=F) # P value for the z statistics under the null hypothesis that there is not effect
wald_ratio_var = wald_ratio_se^2 # Variance
weight <- 1/wald_ratio_var # Inverse variance weight
snps <- merged$BMI_SNP # SNPs that we will use in the estimates
2. Combine the Wald ratios by fixed effects
meta-analysis.
This produces a single, precision-weighted, overall causal estimate across all instrument SNPs, still on the log OR per 1-SD BMI scale.
meta_results <- metagen(wald_ratio,wald_ratio_se,comb.fixed=T,sm="OR") #combine the SNPs by fixed effects meta-analysis## Warning: Use argument 'common' instead of 'comb.fixed' (deprecated).meta_results## Number of studies: k = 77
##
## OR 95%-CI z p-value
## Common effect model 1.3189 [1.1558; 1.5051] 4.11 < 0.0001
## Random effects model 1.2950 [1.1081; 1.5133] 3.25 0.0012
##
## Quantifying heterogeneity (with 95%-CIs):
## tau^2 = 0.0723 [0.0447; 0.6449]; tau = 0.2689 [0.2114; 0.8031]
## I^2 = 31.5% [9.0%; 48.5%]; H = 1.21 [1.05; 1.39]
##
## Test of heterogeneity:
## Q d.f. p-value
## 110.97 76 0.0055
##
## Details of meta-analysis methods:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on Q
3. Create a table of the results, which you could export
to other programs e.g. excel, STATA etc.
mr_results <- data.frame(matrix(c(as.character(merged$BMI_SNP),round(wald_ratio,2),round(wald_ratio_se,2),round(p,3)),nrow=length(merged$BMI_SNP),ncol=4))
names(mr_results) <- c("SNP","log_odds","se","p")
mr_results_order <- mr_results[order(mr_results$log_odds),]
overall_genetic_effect <- data.frame(matrix(c("Overall genetic effect",meta_results$TE.fixed,meta_results$seTE.fixed,meta_results$pval.fixed),nrow=1,ncol=4))
names(overall_genetic_effect) <- c("SNP","log_odds","se","p")
overall_genetic_effect## SNP log_odds se p
## 1 Overall genetic effect 0.27682496389192 0.0673660975621878 3.96925063511265e-05twosampleResults <- rbind.fill(mr_results_order,overall_genetic_effect)
twosampleResults## SNP log_odds se p
## 1 rs7138803 -0.02 0.47 0.966
## 2 rs12885454 -0.03 0.71 0.972
## 3 rs1928295 -0.03 0.74 0.97
## 4 rs543874 -0.1 0.37 0.793
## 5 rs17094222 -0.24 0.69 0.728
## 6 rs12446632 -0.25 0.53 0.63
## 7 rs29941 -0.25 0.83 0.766
## 8 rs1528435 -0.26 0.79 0.742
## 9 rs12940622 -0.32 0.77 0.682
## 10 rs9925964 -0.32 0.75 0.67
## 11 rs2033529 -0.33 0.83 0.694
## 12 rs7899106 -0.41 0.75 0.587
## 13 rs1000940 -0.48 0.86 0.574
## 14 rs3736485 -0.51 0.78 0.519
## 15 rs12566985 -0.54 0.58 0.36
## 16 rs2650492 -0.6 1.05 0.566
## 17 rs1516725 -0.62 0.45 0.168
## 18 rs11057405 -0.69 0.99 0.484
## 19 rs11727676 -0.7 0.92 0.444
## 20 rs4256980 -0.72 0.7 0.305
## 21 rs7903146 -0.84 0.67 0.21
## 22 rs13078960 -0.86 0.59 0.142
## 23 rs11126666 -1 0.76 0.189
## 24 rs2033732 -1.33 0.93 0.153
## 25 rs2075650 -1.46 1.16 0.209
## 26 rs6477694 -1.59 0.89 0.074
## 27 rs758747 -2.67 1.1 0.016
## 28 rs11191560 -3.12 0.97 0.001
## 29 rs12286929 0.01 0.63 0.989
## 30 rs10182181 0.03 0.46 0.95
## 31 rs11688816 0.03 0.83 0.971
## 32 rs11165643 0.04 0.64 0.956
## 33 rs17405819 0.07 0.68 0.917
## 34 rs7599312 0.07 0.72 0.919
## 35 rs6804842 0.08 0.76 0.911
## 36 rs7141420 0.08 0.58 0.895
## 37 rs13107325 0.1 0.88 0.907
## 38 rs2121279 0.13 0.81 0.874
## 39 rs3888190 0.14 0.85 0.871
## 40 rs3101336 0.18 0.44 0.68
## 41 rs17024393 0.21 0.68 0.756
## 42 rs13021737 0.31 0.31 0.316
## 43 rs10968576 0.32 0.6 0.588
## 44 rs2112347 0.34 0.56 0.542
## 45 rs2207139 0.38 0.41 0.357
## 46 rs16851483 0.39 0.6 0.511
## 47 rs1558902 0.4 0.17 0.021
## 48 rs17001654 0.41 0.61 0.506
## 49 rs13191362 0.42 0.74 0.574
## 50 rs6567160 0.44 0.29 0.135
## 51 rs7243357 0.45 0.84 0.596
## 52 rs1808579 0.57 0.85 0.503
## 53 rs12016871 0.61 0.59 0.309
## 54 rs2365389 0.68 0.71 0.34
## 55 rs10132280 0.74 0.72 0.307
## 56 rs657452 0.76 0.64 0.232
## 57 rs11030104 0.84 0.44 0.054
## 58 rs3817334 0.84 0.55 0.125
## 59 rs4740619 0.86 0.79 0.276
## 60 rs10938397 0.87 0.4 0.03
## 61 rs12401738 0.93 0.73 0.203
## 62 rs1167827 0.95 0.82 0.244
## 63 rs2245368 0.95 1.6 0.553
## 64 rs11583200 1.06 0.81 0.193
## 65 rs1016287 1.11 0.69 0.111
## 66 rs2176598 1.24 0.85 0.144
## 67 rs10733682 1.25 0.97 0.195
## 68 rs2287019 1.26 0.69 0.068
## 69 rs11847697 1.44 0.77 0.059
## 70 rs12429545 1.46 0.74 0.049
## 71 rs3810291 1.51 0.73 0.039
## 72 rs16951275 1.66 0.59 0.005
## 73 rs2820292 1.8 0.75 0.016
## 74 rs9400239 1.84 0.86 0.032
## 75 rs3849570 2.06 1.38 0.137
## 76 rs17724992 2.6 1.03 0.012
## 77 rs205262 2.79 0.88 0.002
## 78 Overall genetic effect 0.27682496389192 0.0673660975621878 3.96925063511265e-05
All log_odds values in this table are log
odds ratios for CHD per 1-SD (≈4.56 kg/m²) genetically-predicted
increase in BMI.
4. Calculate heterogeneity p value
Heterogeneity here means: do the individual SNP-specific causal estimates (Wald ratios) agree with each other more than expected by chance alone? High heterogeneity (low p-value) can indicate that one or more SNPs are affecting CHD through a pathway other than BMI (horizontal pleiotropy), rather than every SNP giving a consistent estimate of the same true causal effect.
p_het <- pchisq(meta_results$Q,meta_results$df.Q,lower.tail=F)
twosampleResults$p_chi<-NA
Add the p value for heterogeneity to the results table.
twosampleResults$p_chi[twosampleResults$SNP=="Overall genetic effect"] <- p_het
write.table(twosampleResults,"./twosample_results_BMI_CHD.txt",sep="\t",col.names=T,row.names=F,quote=F)
twosampleResults## SNP log_odds se p p_chi
## 1 rs7138803 -0.02 0.47 0.966 NA
## 2 rs12885454 -0.03 0.71 0.972 NA
## 3 rs1928295 -0.03 0.74 0.97 NA
## 4 rs543874 -0.1 0.37 0.793 NA
## 5 rs17094222 -0.24 0.69 0.728 NA
## 6 rs12446632 -0.25 0.53 0.63 NA
## 7 rs29941 -0.25 0.83 0.766 NA
## 8 rs1528435 -0.26 0.79 0.742 NA
## 9 rs12940622 -0.32 0.77 0.682 NA
## 10 rs9925964 -0.32 0.75 0.67 NA
## 11 rs2033529 -0.33 0.83 0.694 NA
## 12 rs7899106 -0.41 0.75 0.587 NA
## 13 rs1000940 -0.48 0.86 0.574 NA
## 14 rs3736485 -0.51 0.78 0.519 NA
## 15 rs12566985 -0.54 0.58 0.36 NA
## 16 rs2650492 -0.6 1.05 0.566 NA
## 17 rs1516725 -0.62 0.45 0.168 NA
## 18 rs11057405 -0.69 0.99 0.484 NA
## 19 rs11727676 -0.7 0.92 0.444 NA
## 20 rs4256980 -0.72 0.7 0.305 NA
## 21 rs7903146 -0.84 0.67 0.21 NA
## 22 rs13078960 -0.86 0.59 0.142 NA
## 23 rs11126666 -1 0.76 0.189 NA
## 24 rs2033732 -1.33 0.93 0.153 NA
## 25 rs2075650 -1.46 1.16 0.209 NA
## 26 rs6477694 -1.59 0.89 0.074 NA
## 27 rs758747 -2.67 1.1 0.016 NA
## 28 rs11191560 -3.12 0.97 0.001 NA
## 29 rs12286929 0.01 0.63 0.989 NA
## 30 rs10182181 0.03 0.46 0.95 NA
## 31 rs11688816 0.03 0.83 0.971 NA
## 32 rs11165643 0.04 0.64 0.956 NA
## 33 rs17405819 0.07 0.68 0.917 NA
## 34 rs7599312 0.07 0.72 0.919 NA
## 35 rs6804842 0.08 0.76 0.911 NA
## 36 rs7141420 0.08 0.58 0.895 NA
## 37 rs13107325 0.1 0.88 0.907 NA
## 38 rs2121279 0.13 0.81 0.874 NA
## 39 rs3888190 0.14 0.85 0.871 NA
## 40 rs3101336 0.18 0.44 0.68 NA
## 41 rs17024393 0.21 0.68 0.756 NA
## 42 rs13021737 0.31 0.31 0.316 NA
## 43 rs10968576 0.32 0.6 0.588 NA
## 44 rs2112347 0.34 0.56 0.542 NA
## 45 rs2207139 0.38 0.41 0.357 NA
## 46 rs16851483 0.39 0.6 0.511 NA
## 47 rs1558902 0.4 0.17 0.021 NA
## 48 rs17001654 0.41 0.61 0.506 NA
## 49 rs13191362 0.42 0.74 0.574 NA
## 50 rs6567160 0.44 0.29 0.135 NA
## 51 rs7243357 0.45 0.84 0.596 NA
## 52 rs1808579 0.57 0.85 0.503 NA
## 53 rs12016871 0.61 0.59 0.309 NA
## 54 rs2365389 0.68 0.71 0.34 NA
## 55 rs10132280 0.74 0.72 0.307 NA
## 56 rs657452 0.76 0.64 0.232 NA
## 57 rs11030104 0.84 0.44 0.054 NA
## 58 rs3817334 0.84 0.55 0.125 NA
## 59 rs4740619 0.86 0.79 0.276 NA
## 60 rs10938397 0.87 0.4 0.03 NA
## 61 rs12401738 0.93 0.73 0.203 NA
## 62 rs1167827 0.95 0.82 0.244 NA
## 63 rs2245368 0.95 1.6 0.553 NA
## 64 rs11583200 1.06 0.81 0.193 NA
## 65 rs1016287 1.11 0.69 0.111 NA
## 66 rs2176598 1.24 0.85 0.144 NA
## 67 rs10733682 1.25 0.97 0.195 NA
## 68 rs2287019 1.26 0.69 0.068 NA
## 69 rs11847697 1.44 0.77 0.059 NA
## 70 rs12429545 1.46 0.74 0.049 NA
## 71 rs3810291 1.51 0.73 0.039 NA
## 72 rs16951275 1.66 0.59 0.005 NA
## 73 rs2820292 1.8 0.75 0.016 NA
## 74 rs9400239 1.84 0.86 0.032 NA
## 75 rs3849570 2.06 1.38 0.137 NA
## 76 rs17724992 2.6 1.03 0.012 NA
## 77 rs205262 2.79 0.88 0.002 NA
## 78 Overall genetic effect 0.27682496389192 0.0673660975621878 3.96925063511265e-05 0.00550122
5. Create a forest plot of the results and compare
the genetic and observational associations.
The
observational effect is 1.23 (95% CI: 1.17, 1.29) per 4.56 kg/m2 (i.e.,
per SD) increase in BMI. This is an odds ratio for CHD per 1-SD
higher (observed, not genetically-predicted) BMI, estimated in
a standard observational (e.g., cohort/case-control) study.
Formula for SE from 95% confidence interval: (log(uci)-log(lci))/(1.96*2)
effect <- c(wald_ratio[order(wald_ratio)],meta_results$TE.fixed,log(1.23))
se <- c(wald_ratio_se[order(wald_ratio)],meta_results$seTE.fixed,(log(1.29)-log(1.17))/(1.96*2))
snps <- c(as.character(merged$BMI_SNP)[order(wald_ratio)],"Overall genetic effect","Observational effect")
metaplot(effect,se,labels=snps,conf.level=0.95,logeffect=T,nn=0.1,boxsize=0.8,
xlab="Odds ratio for CHD per 1-SD (~4.56 kg/m^2) higher BMI, with 95% CI",ylab="SNP",cex=0.7)
The plotted boxes and lines are odds ratios and 95% confidence
intervals (the logeffect=T argument tells
metaplot that effect/se were
supplied on the log scale and to exponentiate them for the x-axis); an
OR of 1 (vertical reference line) means no effect of BMI on CHD risk, OR
> 1 means higher BMI genetically predicts higher CHD risk, and OR
< 1 means the opposite.
6. Interpret the
results.
a. Is the MR-derived effect similar to the observational association?
Compare the “Overall genetic effect” OR (from the MR analysis, on the per-1-SD-BMI scale) with the observational OR of 1.23 (95% CI 1.17–1.29) per SD of BMI. Are the point estimates and confidence intervals similar, or does the MR estimate suggest a smaller/larger/absent causal effect than the observational association implies? Bear in mind both are expressed in the same units (OR per 1-SD, i.e., per ~4.56 kg/m², increase in BMI), so they are directly comparable.
b. Is there evidence of heterogeneity in the genetic effects? How do you interpret this?
Look at the heterogeneity p-value (p_chi) calculated
above. A small p-value (e.g., < 0.05) suggests the SNP-specific Wald
ratios are more different from one another than chance alone would
predict — a signal that should prompt closer inspection of individual
SNPs (e.g., via the forest plot) and motivates the sensitivity analyses
in Part 5, several of which are specifically designed to still give a
valid causal estimate even when some instruments are
pleiotropic.
c. Can you think of reasons for caution?
PART 5: Sensitivity analyses.
All that is required is summary level results for each SNP (remember
gp, segp, gd, segd from PART 4).
gp <- merged$BMI_Beta (The effect of the SNP on BMI, in SD-BMI units; 1 SD ≈ 4.56 kg/m²)
segp <- merged$BMI_SE (The standard error of the SNP effect on BMI, SD-BMI units)
gd <- merged$CHD_flip_log_odd (The log odds ratio for CHD (that were harmonized to reflect an increase in BMI), log OR units)
segd <- merged$CHD_log_odds_se (Standard error of the log odds ratio, log OR units)
These sensitivity analyses make different assumptions about how
pleiotropy (if present) behaves across the instrument SNPs, and are used
together as a “triangulation” exercise: if the IVW, MR-Egger, weighted
median and weighted mode estimates broadly agree in direction and
magnitude, that gives more confidence that the overall causal estimate
is not being driven by a small number of pleiotropic SNPs.
1. These functions define the IVW, MR-Egger, weighted median
and weighted mode estimators, respectively, and a function that wraps up
the results.
set seed for replication purposes
set.seed(50)
two.sample.iv.ivw <- function(x, y, sigmax, sigmay) {
beta.ivw.fit = summary(lm(y~x-1, weights=sigmay^-2))
beta.ivw.fit.only = lm(y~x-1, weights=sigmay^-2)
beta.ivw = beta.ivw.fit$coef[1,1]
beta.se.ivw = beta.ivw.fit$coef[1,2]/min(beta.ivw.fit$sigma,1)
beta.df.ivw = length(y) - 1
beta.p.ivw = 2*(1-pt(abs(beta.ivw/beta.se.ivw),beta.df.ivw))
beta.lower.ivw = beta.ivw + (-1*qt(df=beta.df.ivw, 0.975)*beta.se.ivw)
beta.upper.ivw = beta.ivw + (1*qt(df=beta.df.ivw, 0.975)*beta.se.ivw)
return(list(beta.ivw=beta.ivw,beta.se.ivw=beta.se.ivw,beta.lower.ivw=beta.lower.ivw,beta.upper.ivw=beta.upper.ivw,beta.t.ivw=beta.ivw/beta.se.ivw,beta.p.ivw=beta.p.ivw,
beta.ivw.fit.only=beta.ivw.fit.only,beta.df.ivw=beta.df.ivw,beta.ivw.fit=beta.ivw.fit))
}
weighted.median <- function(x, w) {
N = length(x)
ord = order(x);
x = x[ord];
w = w[ord];
Sn = cumsum(w)
S_N = Sn[N]
Pn = (100/S_N)*(Sn-w/2)
if(sort(abs(Pn-50))[1] == 0){M = which(Pn==50); return(x[M])}
Q = length(Pn[sign(Pn-50)==-1])
V1 = Q; V2 = Q+1
M = x[V1] + (50 - Pn[V1])*(x[V2]-x[V1])/(Pn[V2]-Pn[V1])
return(list(beta.median=M,CumSum.median=Sn,ordX.median=x))
}
weighted.median.boot <- function(x, y, sigmax, sigmay, Nsim, alpha, W) {
med = NULL
for (i in 1:Nsim){
y_boot = rnorm(length(y), mean=y, sd=sigmay)
x_boot = rnorm(length(x), mean=x, sd=sigmax)
iv_boot = y_boot/x_boot
run = weighted.median(iv_boot,W)
med[i] = run$beta.median
}
lower = Nsim*alpha/2
upper = Nsim*(1-alpha/2)
Sort = sort(med)
lowerCI = Sort[lower]
upperCI = Sort[upper]
se = sd(med)
t = mean(med)/se
p = 2*(1-pt(abs(t),length(y)-1))
return(list(beta.se.median=se,beta.lower.median=lowerCI,beta.upper.median=upperCI,beta.t.median=t,beta.p.median=p))
}
two.sample.iv.egger <- function(x, y, sigmax, sigmay) {
egger.fit = summary(lm(y~x, weights=sigmay^-2))
df.egger = length(y) - 2
beta.egger = egger.fit$coef[2,1]
beta.se.egger = egger.fit$coef[2,2] / min(egger.fit$sigma, 1)
beta.p.egger = 2*(1-pt(abs(beta.egger/beta.se.egger),df.egger))
beta.lower.egger = beta.egger + (-1*qt(df=df.egger, 0.975)*beta.se.egger)
beta.upper.egger = beta.egger + (1*qt(df=df.egger, 0.975)*beta.se.egger)
alpha.egger = egger.fit$coef[1,1]
alpha.se.egger = egger.fit$coef[1,2] / min(egger.fit$sigma, 1)
alpha.p.egger = 2*(1-pt(abs(alpha.egger/alpha.se.egger),df.egger))
alpha.lower.egger = alpha.egger + (-1*qt(df=df.egger, 0.975)*alpha.se.egger)
alpha.upper.egger = alpha.egger + (1*qt(df=df.egger, 0.975)*alpha.se.egger)
return(list(beta.egger=beta.egger,beta.se.egger=beta.se.egger,beta.lower.egger=beta.lower.egger,beta.upper.egger=beta.upper.egger,beta.t.egger=beta.egger/beta.se.egger,beta.p.egger=beta.p.egger,
alpha.egger=alpha.egger,alpha.se.egger=alpha.se.egger,alpha.lower.egger=alpha.lower.egger,alpha.upper.egger=alpha.upper.egger,alpha.t.egger=alpha.egger/alpha.se.egger,alpha.p.egger=alpha.p.egger))
}
ModeEstimator <- function(x, y, sigmax, sigmay, phi=c(1,0.5,0.25), n_boot=1e4, alpha=0.05) {
beta <- function(BetaIV.in, seBetaIV.in) {
s <- 0.9*(min(sd(BetaIV.in), mad(BetaIV.in)))/length(BetaIV.in)^(1/5)
weights <- seBetaIV.in^-2/sum(seBetaIV.in^-2)
beta <- NULL
for(cur_phi in phi) {
h <- s*cur_phi
densityIV <- density(BetaIV.in, weights=weights, bw=h)
beta[length(beta)+1] <- densityIV$x[densityIV$y==max(densityIV$y)]
}
return(beta)
}
boot <- function(BetaIV.in, seBetaIV.in, beta_Mode.in) {
beta.boot <- matrix(nrow=n_boot, ncol=length(beta_Mode.in))
for(i in 1:n_boot) {
BetaIV.boot <- rnorm(length(BetaIV.in), mean=BetaIV.in, sd=seBetaIV.in[,1])
BetaIV.boot_NOME <- rnorm(length(BetaIV.in), mean=BetaIV.in, sd=seBetaIV.in[,2])
beta.boot[i,1:length(phi)] <- beta(BetaIV.in=BetaIV.boot, seBetaIV.in=rep(1, length(BetaIV)))
beta.boot[i,(length(phi)+1):(2*length(phi))] <- beta(BetaIV.in=BetaIV.boot, seBetaIV.in=seBetaIV.in[,1])
beta.boot[i,(2*length(phi)+1):(3*length(phi))] <- beta(BetaIV.in=BetaIV.boot_NOME, seBetaIV.in=rep(1, length(BetaIV)))
beta.boot[i,(3*length(phi)+1):(4*length(phi))] <- beta(BetaIV.in=BetaIV.boot_NOME, seBetaIV.in=seBetaIV.in[,2])
}
return(beta.boot)
}
BetaIV <- y/x
seBetaIV <- cbind(sqrt((sigmay^2)/(x^2) + ((y^2)*(sigmax^2))/(x^4)), sigmay/abs(x))
beta_SimpleMode <- beta(BetaIV.in=BetaIV, seBetaIV.in=rep(1, length(BetaIV)))
beta_WeightedMode <- beta(BetaIV.in=BetaIV, seBetaIV.in=seBetaIV[,1])
beta_WeightedMode_NOME <- beta(BetaIV.in=BetaIV, seBetaIV.in=seBetaIV[,2])
beta_Mode <- rep(c(beta_SimpleMode, beta_WeightedMode,
beta_SimpleMode, beta_WeightedMode_NOME))
beta_Mode.boot <- boot(BetaIV.in=BetaIV, seBetaIV.in=seBetaIV, beta_Mode.in=beta_Mode)
se_Mode <- apply(beta_Mode.boot, 2, mad)
CIlow_Mode <- beta_Mode-qnorm(1-alpha/2)*se_Mode
CIupp_Mode <- beta_Mode+qnorm(1-alpha/2)*se_Mode
P_Mode <- pt(abs(beta_Mode/se_Mode), df=length(x)-1, lower.tail=F)*2
Method <- rep(c('Simple', 'Weighted', 'Simple (NOME)', 'Weighted (NOME)'), each=length(phi))
Results <- data.frame(Method, phi, beta_Mode, se_Mode, CIlow_Mode, CIupp_Mode, P_Mode)
colnames(Results) <- c('Method', 'phi', 'Estimate', 'SE', 'CI_low', 'CI_upp', 'P')
return(Results)
}
MR_output <- function(ivw,egger,median, mode) {
output = data.frame(matrix(NA, nrow=5, ncol=7))
names(output) = c("test", "parameter", "estimate", "se", "lower_CI", "upper_CI","p_value")
output[1:5,1] = c("IVW","MR-Egger","MR-Egger","Weighted_median","Weighted_mode")
output[1:5,2] = c("beta","beta","alpha","beta","beta")
output[1,3:7] = c(IVW$beta.ivw,IVW$beta.se.ivw,IVW$beta.lower.ivw,IVW$beta.upper.ivw,IVW$beta.p.ivw)
output[2,3:7] = c(Egger$beta.egger,Egger$beta.se.egger,Egger$beta.lower.egger,Egger$beta.upper.egger,Egger$beta.p.egger)
output[3,3:7] = c(Egger$alpha.egger,Egger$alpha.se.egger,Egger$alpha.lower.egger,Egger$alpha.upper.egger,Egger$alpha.p.egger)
output[4,3:7] = c(Median$beta.median,MedianBoot$beta.se.median,MedianBoot$beta.lower.median,MedianBoot$beta.upper.median,MedianBoot$beta.p.median)
output[5,3:7] = c(Mode$Estimate[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00],Mode$SE[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00],Mode$CI_low[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00],Mode$CI_upp[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00],Mode$P[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00])
return(output)
}2. Use the functions to estimate the results.
As with Part 4, every “beta” below is a log OR for CHD per 1-SD (~4.56 kg/m²) increase in BMI, except for the MR-Egger alpha (intercept), which is on the log OR scale too, but represents average pleiotropic bias across the SNPs — i.e., how far the SNP-specific ratios’ regression line sits away from zero when the SNP-BMI effect is zero. An alpha significantly different from zero suggests directional horizontal pleiotropy.
IVW <- two.sample.iv.ivw(gp,gd,segp,segd)
Egger <- two.sample.iv.egger(gp,gd,segp,segd)
Median <- weighted.median(wald_ratio,weight)
MedianBoot <- weighted.median.boot(gp,gd,segp,segd,1000,0.05,weight)
Mode <- ModeEstimator(gp,gd,segp,segd)
sensitivity <- MR_output(IVW,Egger,Median,Mode)
sensitivity## test parameter estimate se lower_CI upper_CI p_value
## 1 IVW beta 0.287658553 0.086237732 0.11590122 0.459415882 0.001318512
## 2 MR-Egger beta 0.375935570 0.209803912 -0.04201525 0.793886395 0.077191677
## 3 MR-Egger alpha -0.002791481 0.006041587 -0.01482694 0.009243978 0.645387108
## 4 Weighted_median beta 0.379652982 0.117839188 0.08299413 0.554455224 0.005920437
## 5 Weighted_mode beta 0.311258179 0.129495876 0.05745093 0.565065432 0.018671609write.table(sensitivity,"./twosample_sensitivity_BMI_CHD.txt",sep="\t",col.names=T,row.names=F,quote=F)IVW - Inverse Variance Weighted - Combines the Wald
ratios using an inverse variance weighted meta-analysis, where the
weight of each ratio is the inverse of the variance of the association
between the SNP and the outcome. It assumes that either every SNP is a
valid instrument, or that any pleiotropy “balances out” across SNPs
(balanced pleiotropy). It’s the most efficient (precise) of these
methods when its assumptions hold, so it’s usually treated as the
primary/main MR estimate, with the others used as sensitivity checks.
MR-Egger - combines the Wald ratio’s together
into a meta-regression to estimate the causal effect adjusted for any
directional pleiotropy. This approach is less powered than the IVW, so
its confidence interval is typically much wider — a non-significant
MR-Egger result doesn’t necessarily contradict a significant IVW result,
it may simply reflect lower statistical power.
MR-Egger
(alpha) - the intercept of the MR-egger meta-regression.
Provides an indication of horizontal pleiotropy when it is not null
(i.e., when its 95% CI excludes zero).
Weighted_median - assigns a weight to each SNP derived
from the inverse variance of each SNP’s effect on the outcome. This
robust method requires only 50% of the variants to be valid and not
exhibit horizontal pleiotropy, unmeasured confounding, etc. — so, unlike
IVW, it can still give a consistent estimate even if up to half of the
SNPs are invalid instruments.
Weighted_mode -
also assigns a weight to each SNP derived from the inverse variance of
each SNP’s effect on the outcome. It assumes that the largest sub-set of
SNPs giving similar causal estimates (the “mode”) are valid instruments,
even if this isn’t a strict majority — making it robust as long as no
single subset of invalid instruments giving a similar (but wrong) answer
is larger than the subset of valid instruments.
Bringing it together:
If IVW, weighted median, and weighted mode all point in a similar direction and magnitude, and the MR-Egger intercept (alpha) is close to zero, this triangulation of evidence supports a genuine causal effect of BMI on CHD of roughly the size given by the IVW/overall genetic effect estimate.
If the methods disagree substantially, or the MR-Egger intercept is significantly different from zero, this points towards horizontal pleiotropy biasing at least some of the SNP-specific estimates, and the causal conclusion should be treated with more caution — the weighted median/mode estimates (which are more robust to a minority of invalid instruments) may then be more trustworthy than the IVW estimate.