Obesity on gout

# Obesity on gout
setwd("/n/home04/cdadams/urate")
exposure_dat <- extract_instruments(c('90'))
exposure_dat <- clump_data(exposure_dat)

outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('UKB-b:13251'),
  proxies = 1, rsq = 0.8, align_alleles = 1, 
  palindromes = 1, maf_threshold = 0.3)

dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
mr_results <- mr(dat)
mr_results
##   id.exposure  id.outcome
## 1    ieu-a-90 ukb-b-13251
## 2    ieu-a-90 ukb-b-13251
## 3    ieu-a-90 ukb-b-13251
## 4    ieu-a-90 ukb-b-13251
## 5    ieu-a-90 ukb-b-13251
##                                                          outcome
## 1 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 2 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 3 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 4 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 5 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
##                         exposure                    method nsnp           b
## 1 Obesity class 1 || id:ieu-a-90                  MR Egger   17 0.003252025
## 2 Obesity class 1 || id:ieu-a-90           Weighted median   17 0.002193745
## 3 Obesity class 1 || id:ieu-a-90 Inverse variance weighted   17 0.002565396
## 4 Obesity class 1 || id:ieu-a-90               Simple mode   17 0.002028544
## 5 Obesity class 1 || id:ieu-a-90             Weighted mode   17 0.002137643
##             se        pval
## 1 0.0025468401 0.221055817
## 2 0.0008787420 0.012543954
## 3 0.0008672378 0.003095181
## 4 0.0014821249 0.190012691
## 5 0.0009847333 0.045338610
# Outlier's check
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome, 
  seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                 Estimate    Std.Error  t value     Pr(>|t|)
## Effect (1st) 0.002565396 0.0008672378 2.958123 3.095181e-03
## Iterative    0.002564702 0.0008670531 2.957953 3.096895e-03
## Exact (FE)   0.002648780 0.0005986705 4.424437 9.669416e-06
## Exact (RE)   0.002610278 0.0007941275 3.286976 4.645000e-03
## 
## 
## Residual standard error: 1.462 on 16 degrees of freedom
## 
## F-statistic: 8.75 on 1 and 16 DF, p-value: 0.00925
## Q-Statistic for heterogeneity: 34.18127 on 16 DF , p-value: 0.005135533
## 
##  Outliers detected 
## Number of iterations = 2
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                 Estimate    Std.Error  t value     Pr(>|t|)
## Effect (1st) 0.002565396 0.0008672378 2.958123 3.095181e-03
## Iterative    0.002564702 0.0008670531 2.957953 3.096895e-03
## Exact (FE)   0.002648780 0.0005986705 4.424437 9.669416e-06
## Exact (RE)   0.002610278 0.0008091902 3.225790 5.283350e-03
## 
## 
## Residual standard error: 1.462 on 16 degrees of freedom
## 
## F-statistic: 8.75 on 1 and 16 DF, p-value: 0.00925
## Q-Statistic for heterogeneity: 34.18127 on 16 DF , p-value: 0.005135533
## 
##  Outliers detected 
## Number of iterations = 2
dim(ivwrad$outliers)[1] 
## [1] 4
outliers=ivwrad$outliers[1]
outliers
##         SNP
## 1 rs4929923
## 2  rs523288
## 3  rs887912
## 4 rs9816226
myvars=outliers$SNP
dat2 <- dat[ ! dat$SNP %in% myvars, ]
mr_results <- mr(dat2, parameters = default_parameters())
mr_results
##   id.exposure  id.outcome
## 1    ieu-a-90 ukb-b-13251
## 2    ieu-a-90 ukb-b-13251
## 3    ieu-a-90 ukb-b-13251
## 4    ieu-a-90 ukb-b-13251
## 5    ieu-a-90 ukb-b-13251
##                                                          outcome
## 1 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 2 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 3 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 4 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 5 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
##                         exposure                    method nsnp           b
## 1 Obesity class 1 || id:ieu-a-90                  MR Egger   13 0.002520496
## 2 Obesity class 1 || id:ieu-a-90           Weighted median   13 0.002148266
## 3 Obesity class 1 || id:ieu-a-90 Inverse variance weighted   13 0.002500919
## 4 Obesity class 1 || id:ieu-a-90               Simple mode   13 0.001204967
## 5 Obesity class 1 || id:ieu-a-90             Weighted mode   13 0.001898962
##             se         pval
## 1 0.0017649599 0.1810412034
## 2 0.0008573558 0.0122213654
## 3 0.0006546913 0.0001334533
## 4 0.0014833664 0.4324284789
## 5 0.0011062912 0.1117448790
singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
  single_method = "mr_wald_ratio", 
  all_method = c("mr_ivw",
  "mr_egger_regression", 
  "mr_weighted_mode",
  "mr_weighted_median"))
mr_forest_plot(singlesnp_results, exponentiate = TRUE) #saved as PDF, dimensions 8x8
## $`ieu-a-90.ukb-b-13251`
## Warning: Removed 1 rows containing missing values (geom_errorbarh).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure  id.outcome
## 1    ieu-a-90 ukb-b-13251
mr_scatter_plot(mr_results, dat2)
## $`ieu-a-90.ukb-b-13251`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure  id.outcome
## 1    ieu-a-90 ukb-b-13251
# Remove snps not used in MR
singlesnp_results.2=singlesnp_results[1:13,]
dim(singlesnp_results.2)
## [1] 13  9
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dat=dat2.1

# For the multivarible MR, save the list of kept obesity SNPs
obesity_SNPs=dat$SNP
obesity_SNPs=as.character(obesity_SNPs)

# r2 and Fstat
dat$samplesize.outcome
##  [1] 462933 462933 462933 462933 462933 462933 462933 462933 462933 462933
## [11] 462933 462933 462933
n=mean(dat$samplesize.exposure)
k=13
p=dat$eaf.exposure
r2.1<-(2*(dat$beta.exposure^2)*p*(1-p))/1
r2=sum(r2.1)
r2 
## [1] 0.0655269
Fstat <- r2*(n-1-k)/((1-r2)*k)
Fstat
## [1] 524.8631
# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)

write.csv(mr_results,'res_obesity_gout.csv')
write.csv(singlesnp_results,'singlesnps_obesity_gout.csv')
write.csv(plt,'plt_obesity_gout.csv')
write.csv(het,'het_obesity_gout.csv')
write.csv(dat, 'SIMEX_dat_obesity_gout.csv')
write.csv(outliers, 'outliers_obesity_gout.csv')
write.csv(obesity_SNPs, 'kept_obesity_SNPs.csv') #for MVMR

Obesity on urate

# Obesity on urate
exposure_dat <- extract_instruments(c('90'))
exposure_dat <- clump_data(exposure_dat)
outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('1055'),
  proxies = 1, rsq = 0.8, align_alleles = 1, 
  palindromes = 1, maf_threshold = 0.3)
dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
mr_results <- mr(dat)
mr_results
##   id.exposure id.outcome                outcome                       exposure
## 1    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
## 2    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
## 3    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
## 4    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
## 5    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
##                      method nsnp         b         se         pval
## 1                  MR Egger   14 0.1494135 0.04356155 4.986211e-03
## 2           Weighted median   14 0.1147244 0.02138012 8.052916e-08
## 3 Inverse variance weighted   14 0.1338264 0.01583943 2.939564e-17
## 4               Simple mode   14 0.1352328 0.03796511 3.475437e-03
## 5             Weighted mode   14 0.1229910 0.02601782 3.953487e-04
# Outlier's check
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome, 
  seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##               Estimate  Std.Error  t value     Pr(>|t|)
## Effect (1st) 0.1293876 0.01723304 7.508115 5.998497e-14
## Iterative    0.1293017 0.01721436 7.511271 5.855616e-14
## Exact (FE)   0.1322927 0.01394040 9.489876 2.313088e-21
## Exact (RE)   0.1312928 0.01632525 8.042318 5.178224e-07
## 
## 
## Residual standard error: 1.29 on 16 degrees of freedom
## 
## F-statistic: 56.37 on 1 and 16 DF, p-value: 1.25e-06
## Q-Statistic for heterogeneity: 26.62542 on 16 DF , p-value: 0.04583906
## 
##  Outliers detected 
## Number of iterations = 3
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##               Estimate  Std.Error  t value     Pr(>|t|)
## Effect (1st) 0.1293876 0.01723304 7.508115 5.998497e-14
## Iterative    0.1293017 0.01721436 7.511271 5.855616e-14
## Exact (FE)   0.1322927 0.01394040 9.489876 2.313088e-21
## Exact (RE)   0.1312928 0.01597545 8.218412 3.903143e-07
## 
## 
## Residual standard error: 1.29 on 16 degrees of freedom
## 
## F-statistic: 56.37 on 1 and 16 DF, p-value: 1.25e-06
## Q-Statistic for heterogeneity: 26.62542 on 16 DF , p-value: 0.04583906
## 
##  Outliers detected 
## Number of iterations = 3
dim(ivwrad$outliers)[1] 
## [1] 2
outliers=ivwrad$outliers[1]
outliers
##         SNP
## 1 rs7531118
## 2 rs9816226
myvars=outliers$SNP
dat2 <- dat[ ! dat$SNP %in% myvars, ]
mr_results <- mr(dat2)
mr_results
##   id.exposure id.outcome                outcome                       exposure
## 1    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
## 2    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
## 3    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
## 4    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
## 5    ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055 Obesity class 1 || id:ieu-a-90
##                      method nsnp         b         se         pval
## 1                  MR Egger   13 0.1625215 0.03919417 1.625851e-03
## 2           Weighted median   13 0.1140041 0.02184793 1.807928e-07
## 3 Inverse variance weighted   13 0.1274276 0.01490288 1.225204e-17
## 4               Simple mode   13 0.1302263 0.03514084 3.003549e-03
## 5             Weighted mode   13 0.1203524 0.02689369 7.588563e-04
singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
  single_method = "mr_wald_ratio", 
  all_method = c("mr_ivw",
  "mr_egger_regression", 
  "mr_weighted_mode",
  "mr_weighted_median"))
mr_forest_plot(singlesnp_results, exponentiate = FALSE)#saved as PDF, dimensions 600x600
## $`ieu-a-90.ieu-a-1055`
## Warning: Removed 1 rows containing missing values (geom_errorbarh).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1    ieu-a-90 ieu-a-1055
mr_scatter_plot(mr_results, dat2)
## $`ieu-a-90.ieu-a-1055`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1    ieu-a-90 ieu-a-1055
# Remove snps not used in MR
singlesnp_results.2=singlesnp_results[1:13,]
dim(singlesnp_results.2)
## [1] 13  9
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dat=dat2.1

# r2 and Fstat
dat$samplesize.outcome
##  [1] 109975 102860 103766 108212 110092 105467 109249 109766 109757 110147
## [11] 109975 105492 103183
n=mean(dat$samplesize.exposure)
k=13
p=dat$eaf.exposure
r2.1<-(2*(dat$beta.exposure^2)*p*(1-p))/1
r2=sum(r2.1)
r2 
## [1] 0.06418817
Fstat <- r2*(n-1-k)/((1-r2)*k)
Fstat
## [1] 513.5612
# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)

# all_res<-combine_all_mrresults(res,het,plt,sin,ao_slc=T,
#   Exp=T,split.exposure=T,split.outcome=T)
# all_res$lower_intercept=all_res$intercept-1.96*all_res$intercept_se
# all_res$upper_intercept=all_res$intercept+1.96*all_res$intercept_se
# all_res$r2=r2
# all_res$Fstat=Fstat
# head(all_res[,c("Method","outcome","exposure","nsnp","b","se","or", "or_lci95",  
#   "or_uci95","pval","intercept","intercept_se","intercept_pval",
#   "Q","Q_df","Q_pval","consortium","ncase","ncontrol","pmid","population")])
# write.csv(all_res,"combined_results_obesity_urate.csv")

write.csv(mr_results,'res_obesity_urate.csv')
write.csv(singlesnp_results,'singlesnps_obesity_urate.csv')
write.csv(plt,'plt_obesity_urate.csv')
write.csv(het,'het_obesity_urate.csv')
write.csv(dat, 'SIMEX_obesity_urate.csv')
write.csv(outliers, 'outliers_obesity_urate.csv')

Urate on gout

# Urate on gout
exposure_dat <- extract_instruments(c('1055')) #urate 
exposure_dat <- clump_data(exposure_dat)
outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('UKB-b:13251'),
  proxies = 1, rsq = 0.8, align_alleles = 1, 
  palindromes = 1, maf_threshold = 0.3)

dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
mr_results <- mr(dat)
mr_results
##   id.exposure  id.outcome
## 1  ieu-a-1055 ukb-b-13251
## 2  ieu-a-1055 ukb-b-13251
## 3  ieu-a-1055 ukb-b-13251
## 4  ieu-a-1055 ukb-b-13251
## 5  ieu-a-1055 ukb-b-13251
##                                                          outcome
## 1 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 2 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 3 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 4 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 5 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
##                 exposure                    method nsnp          b          se
## 1 Urate || id:ieu-a-1055                  MR Egger   25 0.02947976 0.006161929
## 2 Urate || id:ieu-a-1055           Weighted median   25 0.01612214 0.001931061
## 3 Urate || id:ieu-a-1055 Inverse variance weighted   25 0.02835044 0.003292101
## 4 Urate || id:ieu-a-1055               Simple mode   25 0.03066390 0.004642865
## 5 Urate || id:ieu-a-1055             Weighted mode   25 0.01418291 0.001315246
##           pval
## 1 7.972866e-05
## 2 6.893160e-17
## 3 7.201195e-18
## 4 7.859271e-07
## 5 1.104429e-10
# Outlier's check
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome, 
  seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                Estimate    Std.Error   t value      Pr(>|t|)
## Effect (1st) 0.02830933 0.0031455052  8.999929  2.258628e-19
## Iterative    0.02806637 0.0031315217  8.962533  3.173028e-19
## Exact (FE)   0.03085771 0.0008447108 36.530499 3.638134e-292
## Exact (RE)   0.02848601 0.0064786507  4.396904  1.650842e-04
## 
## 
## Residual standard error: 4.544 on 26 degrees of freedom
## 
## F-statistic: 81 on 1 and 26 DF, p-value: 1.81e-09
## Q-Statistic for heterogeneity: 536.8916 on 26 DF , p-value: 7.965029e-97
## 
##  Outliers detected 
## Number of iterations = 3
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                Estimate    Std.Error   t value      Pr(>|t|)
## Effect (1st) 0.02830933 0.0031455052  8.999929  2.258628e-19
## Iterative    0.02806637 0.0031315217  8.962533  3.173028e-19
## Exact (FE)   0.03085771 0.0008447108 36.530499 3.638134e-292
## Exact (RE)   0.02848601 0.0064811586  4.395203  1.658257e-04
## 
## 
## Residual standard error: 4.544 on 26 degrees of freedom
## 
## F-statistic: 81 on 1 and 26 DF, p-value: 1.81e-09
## Q-Statistic for heterogeneity: 536.8916 on 26 DF , p-value: 7.965029e-97
## 
##  Outliers detected 
## Number of iterations = 3
dim(ivwrad$outliers)[1] 
## [1] 9
outliers=ivwrad$outliers[1]
outliers
##          SNP
## 1 rs11722228
## 2  rs1178977
## 3  rs1260326
## 4 rs17050272
## 5  rs1825043
## 6  rs2231142
## 7  rs2941484
## 8   rs675209
## 9  rs7654258
myvars=outliers$SNP
dat2 <- dat[ ! dat$SNP %in% myvars, ]

mr_results <- mr(dat2)
mr_results
##   id.exposure  id.outcome
## 1  ieu-a-1055 ukb-b-13251
## 2  ieu-a-1055 ukb-b-13251
## 3  ieu-a-1055 ukb-b-13251
## 4  ieu-a-1055 ukb-b-13251
## 5  ieu-a-1055 ukb-b-13251
##                                                          outcome
## 1 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 2 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 3 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 4 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
## 5 Non-cancer illness code, self-reported: gout || id:ukb-b-13251
##                 exposure                    method nsnp          b          se
## 1 Urate || id:ieu-a-1055                  MR Egger   16 0.03347015 0.004042479
## 2 Urate || id:ieu-a-1055           Weighted median   16 0.03011249 0.001937618
## 3 Urate || id:ieu-a-1055 Inverse variance weighted   16 0.02952073 0.001213701
## 4 Urate || id:ieu-a-1055               Simple mode   16 0.03164869 0.003014112
## 5 Urate || id:ieu-a-1055             Weighted mode   16 0.03200583 0.002771035
##            pval
## 1  9.159888e-07
## 2  1.831525e-54
## 3 1.122339e-130
## 4  2.619395e-08
## 5  7.274011e-09
dat2$id.exposure="urate"
dat2$exposure="urate"
dat2$id.outcome="gout"
dat2$outcome="gout"

singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
  single_method = "mr_wald_ratio", 
  all_method = c("mr_ivw",
  "mr_egger_regression", 
  "mr_weighted_mode",
  "mr_weighted_median"))
mr_forest_plot(singlesnp_results, exponentiate = TRUE) #saved as PDF, dimensions 8x8
## $urate.gout
## Warning: Removed 1 rows containing missing values (geom_errorbarh).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1       urate       gout
dat2 <- dat[ ! dat$SNP %in% myvars, ]
dat2$id.outcome="gout"
dat2$outcome="gout"
dat2$id.exposure="urate"
dat2$exposure="urate"
mr_results <- mr(dat2)
mr_results
##   id.exposure id.outcome outcome exposure                    method nsnp
## 1       urate       gout    gout    urate                  MR Egger   16
## 2       urate       gout    gout    urate           Weighted median   16
## 3       urate       gout    gout    urate Inverse variance weighted   16
## 4       urate       gout    gout    urate               Simple mode   16
## 5       urate       gout    gout    urate             Weighted mode   16
##            b          se          pval
## 1 0.03347015 0.004042479  9.159888e-07
## 2 0.03011249 0.001902082  1.891840e-56
## 3 0.02952073 0.001213701 1.122339e-130
## 4 0.03164869 0.003325057  9.529902e-08
## 5 0.03200583 0.002713390  5.465389e-09
mr_scatter_plot(mr_results, dat2)
## $urate.gout

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1       urate       gout
# Remove snps not used in MR
singlesnp_results.2=singlesnp_results[1:16,]
dim(singlesnp_results.2)
## [1] 16  9
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dat=dat2.1

# For the multivariable MR, save the kept urate SNPs
urate_SNPs=dat$SNP
urate_SNPs=as.character(urate_SNPs)

#r2 and Fstat
dat$samplesize.outcome
##  [1] 462933 462933 462933 462933 462933 462933 462933 462933 462933 462933
## [11] 462933 462933 462933 462933 462933 462933
n=mean(dat$samplesize.exposure)
k=16
p=dat$eaf.exposure

# No EAF
#r2.1<-(2*(dat$beta.exposure^2)*p*(1-p))/1
#r2=sum(r2.1)
#r2 
#Fstat <- r2*(n-1-k)/((1-r2)*k)
#Fstat

# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)

write.csv(mr_results,'res_urate_gout.csv')
write.csv(singlesnp_results,'singlesnps_urate_gout.csv')
write.csv(plt,'plt_urate_gout.csv')
write.csv(het,'het_urate_gout.csv')
write.csv(dat, 'SIMEX_urate_gout.csv')
write.csv(dat, 'SIMEX_dat_urate_gout.csv')
write.csv(outliers, 'outliers_urate_gout.csv')
write.csv(urate_SNPs, "kept_urate_SNPs_urate_gout.csv")

Urate on T2D

# Urate on T2D
exposure_dat <- extract_instruments(c('1055'))
exposure_dat <- clump_data(exposure_dat)
outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('24'),
  proxies = 1, rsq = 0.8, align_alleles = 1, 
  palindromes = 1, maf_threshold = 0.3)

dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
mr_results <- mr(dat)
mr_results
##   id.exposure id.outcome                        outcome               exposure
## 1  ieu-a-1055   ieu-a-24 Type 2 diabetes || id:ieu-a-24 Urate || id:ieu-a-1055
## 2  ieu-a-1055   ieu-a-24 Type 2 diabetes || id:ieu-a-24 Urate || id:ieu-a-1055
## 3  ieu-a-1055   ieu-a-24 Type 2 diabetes || id:ieu-a-24 Urate || id:ieu-a-1055
## 4  ieu-a-1055   ieu-a-24 Type 2 diabetes || id:ieu-a-24 Urate || id:ieu-a-1055
## 5  ieu-a-1055   ieu-a-24 Type 2 diabetes || id:ieu-a-24 Urate || id:ieu-a-1055
##                      method nsnp            b         se      pval
## 1                  MR Egger    8 -0.207068199 0.24539654 0.4311253
## 2           Weighted median    8 -0.063838363 0.07945920 0.4217374
## 3 Inverse variance weighted    8 -0.117233807 0.12033812 0.3299556
## 4               Simple mode    8  0.009848429 0.13344678 0.9432338
## 5             Weighted mode    8 -0.052736047 0.08008325 0.5312615
# Outlier's check
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome, 
  seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                Estimate  Std.Error    t value   Pr(>|t|)
## Effect (1st) -0.1172338 0.12033812 -0.9742034 0.32995556
## Iterative    -0.1172364 0.12034289 -0.9741867 0.32996386
## Exact (FE)   -0.1196810 0.05892977 -2.0309082 0.04226431
## Exact (RE)   -0.1178114 0.12005463 -0.9813147 0.35912918
## 
## 
## Residual standard error: 2.045 on 7 degrees of freedom
## 
## F-statistic: 0.95 on 1 and 7 DF, p-value: 0.362
## Q-Statistic for heterogeneity: 29.27694 on 7 DF , p-value: 0.0001287918
## 
##  Outliers detected 
## Number of iterations = 2
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                Estimate  Std.Error    t value   Pr(>|t|)
## Effect (1st) -0.1172338 0.12033812 -0.9742034 0.32995556
## Iterative    -0.1172364 0.12034289 -0.9741867 0.32996386
## Exact (FE)   -0.1196810 0.05892977 -2.0309082 0.04226431
## Exact (RE)   -0.1178114 0.11844791 -0.9946260 0.35304872
## 
## 
## Residual standard error: 2.045 on 7 degrees of freedom
## 
## F-statistic: 0.95 on 1 and 7 DF, p-value: 0.362
## Q-Statistic for heterogeneity: 29.27694 on 7 DF , p-value: 0.0001287918
## 
##  Outliers detected 
## Number of iterations = 2
dim(ivwrad$outliers)[1] 
## [1] 2
outliers=ivwrad$outliers[1]
outliers
##         SNP
## 1 rs1260326
## 2  rs642803
myvars=outliers$SNP
dat2 <- dat[ ! dat$SNP %in% myvars, ]

dat2$id.exposure="urate"
dat2$exposure="urate"
dat2$id.outcome="T2D"
dat2$outcome="T2D"
mr_results <- mr(dat2)
mr_results
##   id.exposure id.outcome outcome exposure                    method nsnp
## 1       urate        T2D     T2D    urate                  MR Egger    6
## 2       urate        T2D     T2D    urate           Weighted median    6
## 3       urate        T2D     T2D    urate Inverse variance weighted    6
## 4       urate        T2D     T2D    urate               Simple mode    6
## 5       urate        T2D     T2D    urate             Weighted mode    6
##             b         se      pval
## 1 -0.10247520 0.11673886 0.4296075
## 2 -0.04886657 0.07634258 0.5221103
## 3 -0.05257357 0.06497066 0.4184062
## 4  0.07814569 0.11827307 0.5380123
## 5 -0.05253295 0.08755673 0.5746679
singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
  single_method = "mr_wald_ratio", 
  all_method = c("mr_ivw",
  "mr_egger_regression", 
  "mr_weighted_mode",
  "mr_weighted_median"))
mr_forest_plot(singlesnp_results, exponentiate = TRUE)#saved as PDF, dimensions 8x8
## $urate.T2D
## Warning: Removed 1 rows containing missing values (geom_errorbarh).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1       urate        T2D
mr_scatter_plot(mr_results, dat2)
## $urate.T2D

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1       urate        T2D
write.csv(mr_results,'res.csv')
write.csv(singlesnp_results,'singlesnps.csv')

# Remove snps not used in MR
singlesnp_results.2=singlesnp_results[1:6,]
dim(singlesnp_results.2)
## [1] 6 9
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dat=dat2.1
write.csv(dat, 'SIMEX.csv')

# r2 and Fstat
n=mean(dat$samplesize.outcome)
#n=mean(dat$samplesize.exposure)
k=6
p=dat$eaf.exposure
r2.1<-(2*(dat$beta.exposure^2)*p*(1-p))/1
r2=sum(r2.1)
r2 
## [1] NA
Fstat <- r2*(n-1-k)/((1-r2)*k)
Fstat
## [1] NA
# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)

write.csv(mr_results,'res_urate_T2D.csv')
write.csv(singlesnp_results,'singlesnps_urate_T2D.csv')
write.csv(plt,'plt_urate_T2D.csv')
write.csv(het,'het_urate_T2D.csv')
write.csv(dat, 'SIMEX_urate_T2D.csv')
write.csv(dat, 'SIMEX_dat_urate_T2D.csv')
write.csv(outliers, 'outliers_urate_T2D.csv')

Obesity on HDL

# Obesity on HDL
exposure_dat <- extract_instruments(c('90'))
exposure_dat <- clump_data(exposure_dat)
outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('299'), #,'300'
  proxies = 1, rsq = 0.8, align_alleles = 1, palindromes = 1, maf_threshold = 0.3)
dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
mr_results <- mr(dat)

# Outlier's check
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome, 
  seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                 Estimate   Std.Error   t value     Pr(>|t|)
## Effect (1st) -0.07522770 0.015808356 -4.758730 1.948149e-06
## Iterative    -0.07515250 0.015818934 -4.750794 2.026190e-06
## Exact (FE)   -0.07883826 0.008731888 -9.028776 1.736033e-19
## Exact (RE)   -0.07629542 0.015380417 -4.960556 1.417166e-04
## 
## 
## Residual standard error: 1.882 on 16 degrees of freedom
## 
## F-statistic: 22.65 on 1 and 16 DF, p-value: 0.000214
## Q-Statistic for heterogeneity: 56.64575 on 16 DF , p-value: 1.902833e-06
## 
##  Outliers detected 
## Number of iterations = 3
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                 Estimate   Std.Error   t value     Pr(>|t|)
## Effect (1st) -0.07522770 0.015808356 -4.758730 1.948149e-06
## Iterative    -0.07515250 0.015818934 -4.750794 2.026190e-06
## Exact (FE)   -0.07883826 0.008731888 -9.028776 1.736033e-19
## Exact (RE)   -0.07629542 0.015686312 -4.863821 1.724109e-04
## 
## 
## Residual standard error: 1.882 on 16 degrees of freedom
## 
## F-statistic: 22.65 on 1 and 16 DF, p-value: 0.000214
## Q-Statistic for heterogeneity: 56.64575 on 16 DF , p-value: 1.902833e-06
## 
##  Outliers detected 
## Number of iterations = 3
dim(ivwrad$outliers)[1] 
## [1] 3
outliers=ivwrad$outliers[1]
outliers
##         SNP
## 1  rs523288
## 2 rs9816226
## 3  rs987237
myvars=outliers$SNP
dat2 <- dat[ ! dat$SNP %in% myvars, ]
mr_results <- mr(dat2)
mr_results
##   id.exposure id.outcome                         outcome
## 1    ieu-a-90  ieu-a-299 HDL cholesterol || id:ieu-a-299
## 2    ieu-a-90  ieu-a-299 HDL cholesterol || id:ieu-a-299
## 3    ieu-a-90  ieu-a-299 HDL cholesterol || id:ieu-a-299
## 4    ieu-a-90  ieu-a-299 HDL cholesterol || id:ieu-a-299
## 5    ieu-a-90  ieu-a-299 HDL cholesterol || id:ieu-a-299
##                         exposure                    method nsnp           b
## 1 Obesity class 1 || id:ieu-a-90                  MR Egger   14 -0.07531494
## 2 Obesity class 1 || id:ieu-a-90           Weighted median   14 -0.07721865
## 3 Obesity class 1 || id:ieu-a-90 Inverse variance weighted   14 -0.08296343
## 4 Obesity class 1 || id:ieu-a-90               Simple mode   14 -0.07135635
## 5 Obesity class 1 || id:ieu-a-90             Weighted mode   14 -0.07658158
##            se         pval
## 1 0.024504952 9.654560e-03
## 2 0.012970290 2.624642e-09
## 3 0.009229886 2.504803e-19
## 4 0.018784278 2.212678e-03
## 5 0.014421832 1.413782e-04
singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
  single_method = "mr_wald_ratio", 
  all_method = c("mr_ivw",
  "mr_egger_regression", 
  "mr_weighted_mode",
  "mr_weighted_median"))
mr_forest_plot(singlesnp_results, exponentiate = FALSE)#saved as PDF, dimensions 8x8
## $`ieu-a-90.ieu-a-299`
## Warning: Removed 1 rows containing missing values (geom_errorbarh).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1    ieu-a-90  ieu-a-299
mr_scatter_plot(mr_results, dat2)
## $`ieu-a-90.ieu-a-299`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1    ieu-a-90  ieu-a-299
# Remove snps not used in MR
singlesnp_results.2=singlesnp_results[1:14,]
dim(singlesnp_results.2)
## [1] 14  9
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dat=dat2.1
write.csv(dat, 'SIMEX.csv')

# For the multivarible MR, save the list of kept obesity SNPs
hdl_SNPs=dat$SNP
hdl_SNPs=as.character(hdl_SNPs)

# r2 and Fstat
dat$samplesize.outcome
##  [1] 187005 184024 187119 185214 187129 187096  92820 187004 183768 187071
## [11] 187036 187039 185388 187121
n=mean(dat$samplesize.exposure)
k=14
p=dat$eaf.exposure
r2.1<-(2*(dat$beta.exposure^2)*p*(1-p))/1
r2=sum(r2.1)
r2 
## [1] 0.06768309
Fstat <- r2*(n-1-k)/((1-r2)*k)
Fstat
## [1] 505.8002
# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)

write.csv(mr_results,'res_obesity_HDL.csv')
write.csv(singlesnp_results,'singlesnps_obesity_HDL.csv')
write.csv(plt,'plt_obesity_HDL.csv')
write.csv(het,'het_obesity_HDL.csv')
write.csv(dat, 'SIMEX_obesity_HDL.csv')
write.csv(dat, 'SIMEX_dat_obesity_HDL.csv')
write.csv(outliers, 'outliers_obesity_HDL.csv')
write.csv(hdl_SNPs, 'kept_hdl_SNPs.csv')

Obesity on LDL

# Obesity on LDL
exposure_dat <- extract_instruments(c('90'))
exposure_dat <- clump_data(exposure_dat)
outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('300'), #,'300'
  proxies = 1, rsq = 0.8, align_alleles = 1, palindromes = 1, maf_threshold = 0.3)
dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
mr_results <- mr(dat)

# Outlier's check
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome, 
  seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                 Estimate   Std.Error    t value  Pr(>|t|)
## Effect (1st) -0.01038755 0.027993989 -0.3710636 0.7105901
## Iterative    -0.01038480 0.027992572 -0.3709843 0.7106493
## Exact (FE)   -0.01224611 0.009014419 -1.3585023 0.1743043
## Exact (RE)   -0.01055663 0.022768157 -0.4636576 0.6491382
## 
## 
## Residual standard error: 3.108 on 16 degrees of freedom
## 
## F-statistic: 0.14 on 1 and 16 DF, p-value: 0.715
## Q-Statistic for heterogeneity: 154.5625 on 16 DF , p-value: 9.815379e-25
## 
##  Outliers detected 
## Number of iterations = 2
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                 Estimate   Std.Error    t value  Pr(>|t|)
## Effect (1st) -0.01038755 0.027993989 -0.3710636 0.7105901
## Iterative    -0.01038480 0.027992572 -0.3709843 0.7106493
## Exact (FE)   -0.01224611 0.009014419 -1.3585023 0.1743043
## Exact (RE)   -0.01055663 0.022292179 -0.4735576 0.6422124
## 
## 
## Residual standard error: 3.108 on 16 degrees of freedom
## 
## F-statistic: 0.14 on 1 and 16 DF, p-value: 0.715
## Q-Statistic for heterogeneity: 154.5625 on 16 DF , p-value: 9.815379e-25
## 
##  Outliers detected 
## Number of iterations = 2
dim(ivwrad$outliers)[1] 
## [1] 4
outliers=ivwrad$outliers[1]
outliers
##          SNP
## 1 rs10182181
## 2 rs13130484
## 3  rs2307111
## 4   rs987237
write.csv(outliers, 'adv_outliers.csv')
myvars=outliers$SNP
dat2 <- dat[ ! dat$SNP %in% myvars, ]
mr_results <- mr(dat2)
mr_results
##   id.exposure id.outcome                         outcome
## 1    ieu-a-90  ieu-a-300 LDL cholesterol || id:ieu-a-300
## 2    ieu-a-90  ieu-a-300 LDL cholesterol || id:ieu-a-300
## 3    ieu-a-90  ieu-a-300 LDL cholesterol || id:ieu-a-300
## 4    ieu-a-90  ieu-a-300 LDL cholesterol || id:ieu-a-300
## 5    ieu-a-90  ieu-a-300 LDL cholesterol || id:ieu-a-300
##                         exposure                    method nsnp             b
## 1 Obesity class 1 || id:ieu-a-90                  MR Egger   13 -0.0201584749
## 2 Obesity class 1 || id:ieu-a-90           Weighted median   13  0.0023109710
## 3 Obesity class 1 || id:ieu-a-90 Inverse variance weighted   13 -0.0112370333
## 4 Obesity class 1 || id:ieu-a-90               Simple mode   13 -0.0043224426
## 5 Obesity class 1 || id:ieu-a-90             Weighted mode   13  0.0007933091
##           se      pval
## 1 0.02750049 0.4788751
## 2 0.01381051 0.8671070
## 3 0.01000657 0.2614519
## 4 0.02422167 0.8613441
## 5 0.01560220 0.9602848
singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
  single_method = "mr_wald_ratio", 
  all_method = c("mr_ivw",
  "mr_egger_regression", 
  "mr_weighted_mode",
  "mr_weighted_median"))
mr_forest_plot(singlesnp_results, exponentiate = FALSE)#saved as PDF, dimensions 8x8
## $`ieu-a-90.ieu-a-300`
## Warning: Removed 1 rows containing missing values (geom_errorbarh).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1    ieu-a-90  ieu-a-300
mr_scatter_plot(mr_results, dat2)
## $`ieu-a-90.ieu-a-300`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1    ieu-a-90  ieu-a-300
# Remove snps not used in MR
singlesnp_results.2=singlesnp_results[1:13,]
dim(singlesnp_results.2)
## [1] 13  9
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dat=dat2.1
write.csv(dat, 'SIMEX.csv')

# For the multivarible MR, save the list of kept obesity SNPs
ldl_SNPs=dat$SNP
ldl_SNPs=as.character(ldl_SNPs)

# r2 and Fstat
dat$samplesize.outcome
##  [1] 170005 171157 173053  88433 172945 164665 169739 172999 172959 172983
## [11] 171387 173045 172813
n=mean(dat$samplesize.exposure)
k=13
p=dat$eaf.exposure
r2.1<-(2*(dat$beta.exposure^2)*p*(1-p))/1
r2=sum(r2.1)
r2 
## [1] 0.06678086
Fstat <- r2*(n-1-k)/((1-r2)*k)
Fstat
## [1] 534.3247
# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)

write.csv(mr_results,'res_obesity_LDL.csv')
write.csv(singlesnp_results,'singlesnps_obesity_LDL.csv')
write.csv(plt,'plt_obesity_LDL.csv')
write.csv(het,'het_obesity_LDL.csv')
write.csv(dat, 'SIMEX_obesity_LDL.csv')
write.csv(dat, 'SIMEX_dat_obesity_LDL.csv')
write.csv(outliers, 'outliers_obesity_LDL.csv')
write.csv(ldl_SNPs, 'kept_ldl_SNPs.csv')

Obesity on triglycerides

# Obesity on triglycerides
exposure_dat <- extract_instruments(c('90'))
exposure_dat <- clump_data(exposure_dat)
outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('302'), #,'300'
  proxies = 1, rsq = 0.8, align_alleles = 1, palindromes = 1, maf_threshold = 0.3)
dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
mr_results <- mr(dat)

# Outlier's check
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome, 
  seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                Estimate   Std.Error  t value     Pr(>|t|)
## Effect (1st) 0.07417040 0.009715611 7.634147 2.273213e-14
## Iterative    0.07414486 0.009728056 7.621756 2.502481e-14
## Exact (FE)   0.07552761 0.008486643 8.899587 5.605482e-19
## Exact (RE)   0.07521264 0.008972963 8.382141 3.011434e-07
## 
## 
## Residual standard error: 1.188 on 16 degrees of freedom
## 
## F-statistic: 58.28 on 1 and 16 DF, p-value: 1.01e-06
## Q-Statistic for heterogeneity: 22.59847 on 16 DF , p-value: 0.1248971
## 
##  Outliers detected 
## Number of iterations = 3
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                Estimate   Std.Error  t value     Pr(>|t|)
## Effect (1st) 0.07417040 0.009715611 7.634147 2.273213e-14
## Iterative    0.07414486 0.009728056 7.621756 2.502481e-14
## Exact (FE)   0.07552761 0.008486643 8.899587 5.605482e-19
## Exact (RE)   0.07521264 0.008789589 8.557014 2.291083e-07
## 
## 
## Residual standard error: 1.188 on 16 degrees of freedom
## 
## F-statistic: 58.28 on 1 and 16 DF, p-value: 1.01e-06
## Q-Statistic for heterogeneity: 22.59847 on 16 DF , p-value: 0.1248971
## 
##  Outliers detected 
## Number of iterations = 3
dim(ivwrad$outliers)[1] 
## [1] 2
outliers=ivwrad$outliers[1]
outliers
##         SNP
## 1 rs4929923
## 2  rs987237
write.csv(outliers, 'adv_outliers.csv')
myvars=outliers$SNP
dat2 <- dat[ ! dat$SNP %in% myvars, ]
mr_results <- mr(dat2)
mr_results
##   id.exposure id.outcome                       outcome
## 1    ieu-a-90  ieu-a-302 Triglycerides || id:ieu-a-302
## 2    ieu-a-90  ieu-a-302 Triglycerides || id:ieu-a-302
## 3    ieu-a-90  ieu-a-302 Triglycerides || id:ieu-a-302
## 4    ieu-a-90  ieu-a-302 Triglycerides || id:ieu-a-302
## 5    ieu-a-90  ieu-a-302 Triglycerides || id:ieu-a-302
##                         exposure                    method nsnp          b
## 1 Obesity class 1 || id:ieu-a-90                  MR Egger   15 0.06773089
## 2 Obesity class 1 || id:ieu-a-90           Weighted median   15 0.08050144
## 3 Obesity class 1 || id:ieu-a-90 Inverse variance weighted   15 0.08204119
## 4 Obesity class 1 || id:ieu-a-90               Simple mode   15 0.10364179
## 5 Obesity class 1 || id:ieu-a-90             Weighted mode   15 0.07777793
##            se         pval
## 1 0.024067825 1.462754e-02
## 2 0.012740957 2.644223e-10
## 3 0.008586805 1.243831e-21
## 4 0.019774397 1.248306e-04
## 5 0.014328458 8.897672e-05
singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
  single_method = "mr_wald_ratio", 
  all_method = c("mr_ivw",
  "mr_egger_regression", 
  "mr_weighted_mode",
  "mr_weighted_median"))
mr_forest_plot(singlesnp_results, exponentiate = FALSE)#saved as PDF, dimensions 8x8
## $`ieu-a-90.ieu-a-302`
## Warning: Removed 1 rows containing missing values (geom_errorbarh).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1    ieu-a-90  ieu-a-302
mr_scatter_plot(mr_results, dat2)
## $`ieu-a-90.ieu-a-302`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1    ieu-a-90  ieu-a-302
# Remove snps not used in MR
singlesnp_results.2=singlesnp_results[1:15,]
dim(singlesnp_results.2)
## [1] 15  9
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dat=dat2.1

# r2 and Fstat
dat$samplesize.outcome
##  [1] 177702 174684 177815 175913 177823 177795  89485 169079 174409 177769
## [11] 177735 177753 176075 177818 177571
n=mean(dat$samplesize.exposure)
n
## [1] 97286.53
k=13
p=dat$eaf.exposure
r2.1<-(2*(dat$beta.exposure^2)*p*(1-p))/1
r2=sum(r2.1)
r2 
## [1] 0.07545998
Fstat <- r2*(n-1-k)/((1-r2)*k)
Fstat
## [1] 610.7139
# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)

write.csv(mr_results,'res_obesity_TRI.csv')
write.csv(singlesnp_results,'singlesnps_obesity_TRI.csv')
write.csv(plt,'plt_obesity_TRI.csv')
write.csv(het,'het_obesity_TRI.csv')
write.csv(dat, 'SIMEX_obesity_TRI.csv')
write.csv(dat, 'SIMEX_dat_obesity_TRI.csv')
write.csv(outliers, 'outliers_obesity_TRI.csv')

HDL on urate

# HDL on urate
exposure_dat <- extract_instruments(c('299'))
exposure_dat <- clump_data(exposure_dat)
outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('1055'), 
  proxies = 1, rsq = 0.8, align_alleles = 1, palindromes = 1, maf_threshold = 0.3)
dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
mr_results <- mr(dat)

# Outlier's check
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome, 
  seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                Estimate  Std.Error   t value     Pr(>|t|)
## Effect (1st) -0.1208561 0.03322906 -3.637062 2.757658e-04
## Iterative    -0.1208135 0.03324382 -3.634164 2.788831e-04
## Exact (FE)   -0.1252933 0.01530088 -8.188630 2.642175e-16
## Exact (RE)   -0.1217867 0.03018021 -4.035318 1.178026e-04
## 
## 
## Residual standard error: 2.179 on 86 degrees of freedom
## 
## F-statistic: 13.23 on 1 and 86 DF, p-value: 0.000469
## Q-Statistic for heterogeneity: 408.4026 on 86 DF , p-value: 1.952549e-43
## 
##  Outliers detected 
## Number of iterations = 3
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##                Estimate  Std.Error   t value     Pr(>|t|)
## Effect (1st) -0.1208561 0.03322906 -3.637062 2.757658e-04
## Iterative    -0.1208135 0.03324382 -3.634164 2.788831e-04
## Exact (FE)   -0.1252933 0.01530088 -8.188630 2.642175e-16
## Exact (RE)   -0.1217867 0.02980522 -4.086087 9.817034e-05
## 
## 
## Residual standard error: 2.179 on 86 degrees of freedom
## 
## F-statistic: 13.23 on 1 and 86 DF, p-value: 0.000469
## Q-Statistic for heterogeneity: 408.4026 on 86 DF , p-value: 1.952549e-43
## 
##  Outliers detected 
## Number of iterations = 3
dim(ivwrad$outliers)[1] 
## [1] 21
outliers=ivwrad$outliers[1]
outliers
##           SNP
## 1    rs102275
## 2   rs1047891
## 3  rs11065987
## 4  rs12133576
## 5  rs12748152
## 6  rs12801636
## 7  rs13099479
## 8  rs13107325
## 9  rs17145738
## 10  rs1980493
## 11  rs2075650
## 12  rs2288912
## 13  rs2606736
## 14  rs3741414
## 15  rs3861397
## 16  rs6031587
## 17  rs6567160
## 18   rs731839
## 19  rs7607980
## 20   rs970548
## 21  rs9989419
myvars=outliers$SNP
dat2 <- dat[ ! dat$SNP %in% myvars, ]
mr_results <- mr(dat2)
mr_results
##   id.exposure id.outcome                outcome                        exposure
## 1   ieu-a-299 ieu-a-1055 Urate || id:ieu-a-1055 HDL cholesterol || id:ieu-a-299
## 2   ieu-a-299 ieu-a-1055 Urate || id:ieu-a-1055 HDL cholesterol || id:ieu-a-299
## 3   ieu-a-299 ieu-a-1055 Urate || id:ieu-a-1055 HDL cholesterol || id:ieu-a-299
## 4   ieu-a-299 ieu-a-1055 Urate || id:ieu-a-1055 HDL cholesterol || id:ieu-a-299
## 5   ieu-a-299 ieu-a-1055 Urate || id:ieu-a-1055 HDL cholesterol || id:ieu-a-299
##                      method nsnp           b         se         pval
## 1                  MR Egger   64 -0.07119904 0.03924520 7.448440e-02
## 2           Weighted median   64 -0.12198397 0.02991625 4.551942e-05
## 3 Inverse variance weighted   64 -0.10905221 0.01962219 2.735137e-08
## 4               Simple mode   64 -0.07732532 0.05685162 1.786376e-01
## 5             Weighted mode   64 -0.11510603 0.03503832 1.667452e-03
singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
  single_method = "mr_wald_ratio", 
  all_method = c("mr_ivw",
  "mr_egger_regression", 
  "mr_weighted_mode",
  "mr_weighted_median"))
mr_forest_plot(singlesnp_results, exponentiate = FALSE)#saved as PDF, dimensions 600x600
## $`ieu-a-299.ieu-a-1055`
## Warning: Removed 1 rows containing missing values (geom_errorbarh).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1   ieu-a-299 ieu-a-1055
mr_scatter_plot(mr_results, dat2)
## $`ieu-a-299.ieu-a-1055`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1   ieu-a-299 ieu-a-1055
# Remove snps not used in MR
singlesnp_results.2=singlesnp_results[1:64,]
dim(singlesnp_results.2)
## [1] 64  9
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dat=dat2.1
write.csv(dat, 'SIMEX.csv')

# For the multivarible MR, save the list of kept obesity SNPs
hdl_SNPs=dat$SNP
hdl_SNPs=as.character(hdl_SNPs)

#r2 and Fstat
dat$samplesize.outcome
##  [1] 109852 105533 108486 101977 110127 110195 102462 105556 109309 109591
## [11] 105369 110074 109981 109934 109932 109879 110047 103662 109914 110002
## [21] 108406 110212 109005 105332 110137 110142 110007 110154 110209 109717
## [31] 103817 109828 105938 110126 108363 110217 102844 109867 109684 110065
## [41] 110127 103144 102860 109984 109989 108526 110206 110048 110083 109618
## [51] 103756 109963 109325 109687 110226 109979 110219 110197 109860 107424
## [61] 105532 103012 106388  94514
n=mean(dat$samplesize.outcome)
k=64
p=dat$eaf.exposure
r2.1<-(2*(dat$beta.exposure^2)*p*(1-p))/1
r2=sum(r2.1)
r2 
## [1] 0.04586729
Fstat <- r2*(n-1-k)/((1-r2)*k)
Fstat
## [1] 81.17422
# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)

write.csv(mr_results,'res_HDL_urate.csv')
write.csv(singlesnp_results,'singlesnps_HDL_urate.csv')
write.csv(plt,'plt_HDL_urate.csv')
write.csv(het,'het_HDL_urate.csv')
write.csv(dat, 'SIMEX_HDL_urate.csv')
write.csv(dat, 'SIMEX_dat_HDL_urate.csv')
write.csv(outliers, 'outliers_HDL_urate.csv')
write.csv(hdl_SNPs, 'kept_hdl_SNPs.csv')

Triglycerides on urate

# Triglycerides on urate
exposure_dat <- extract_instruments(c('302'))
exposure_dat <- clump_data(exposure_dat)
outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('1055'), 
    proxies = 1, rsq = 0.8, align_alleles = 1, palindromes = 1, maf_threshold = 0.3)
dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
mr_results <- mr(dat)

# Outlier's check
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome, 
  seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##               Estimate  Std.Error   t value     Pr(>|t|)
## Effect (1st) 0.1961050 0.04383078  4.474139 7.671975e-06
## Iterative    0.1961577 0.04384134  4.474264 7.667483e-06
## Exact (FE)   0.2038959 0.01693124 12.042584 2.122031e-33
## Exact (RE)   0.1972147 0.07000275  2.817242 6.752584e-03
## 
## 
## Residual standard error: 2.609 on 54 degrees of freedom
## 
## F-statistic: 20.02 on 1 and 54 DF, p-value: 3.99e-05
## Q-Statistic for heterogeneity: 367.6537 on 54 DF , p-value: 3.160823e-48
## 
##  Outliers detected 
## Number of iterations = 3
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
## 
## Radial IVW
## 
##               Estimate  Std.Error   t value     Pr(>|t|)
## Effect (1st) 0.1961050 0.04383078  4.474139 7.671975e-06
## Iterative    0.1961577 0.04384134  4.474264 7.667483e-06
## Exact (FE)   0.2038959 0.01693124 12.042584 2.122031e-33
## Exact (RE)   0.1972147 0.06853726  2.877482 5.730330e-03
## 
## 
## Residual standard error: 2.609 on 54 degrees of freedom
## 
## F-statistic: 20.02 on 1 and 54 DF, p-value: 3.99e-05
## Q-Statistic for heterogeneity: 367.6537 on 54 DF , p-value: 3.160823e-48
## 
##  Outliers detected 
## Number of iterations = 3
dim(ivwrad$outliers)[1] 
## [1] 20
outliers=ivwrad$outliers[1]
outliers
##           SNP
## 1  rs10401969
## 2  rs11613352
## 3  rs11974409
## 4  rs12280753
## 5   rs1260326
## 6  rs12748152
## 7   rs2043085
## 8   rs2247056
## 9    rs247616
## 10   rs439401
## 11  rs4587594
## 12  rs4738684
## 13  rs4810479
## 14   rs588136
## 15   rs634869
## 16  rs7248104
## 17   rs731839
## 18  rs7350481
## 19  rs8077889
## 20   rs948690
write.csv(outliers, 'adv_outliers.csv')
myvars=outliers$SNP
dat2 <- dat[ ! dat$SNP %in% myvars, ]
mr_results <- mr(dat2)
mr_results
##   id.exposure id.outcome                outcome                      exposure
## 1   ieu-a-302 ieu-a-1055 Urate || id:ieu-a-1055 Triglycerides || id:ieu-a-302
## 2   ieu-a-302 ieu-a-1055 Urate || id:ieu-a-1055 Triglycerides || id:ieu-a-302
## 3   ieu-a-302 ieu-a-1055 Urate || id:ieu-a-1055 Triglycerides || id:ieu-a-302
## 4   ieu-a-302 ieu-a-1055 Urate || id:ieu-a-1055 Triglycerides || id:ieu-a-302
## 5   ieu-a-302 ieu-a-1055 Urate || id:ieu-a-1055 Triglycerides || id:ieu-a-302
##                      method nsnp         b         se         pval
## 1                  MR Egger   35 0.1141178 0.04726498 2.146612e-02
## 2           Weighted median   35 0.1742672 0.04183604 3.106987e-05
## 3 Inverse variance weighted   35 0.1983528 0.02660299 8.916419e-14
## 4               Simple mode   35 0.3576558 0.10779977 2.169696e-03
## 5             Weighted mode   35 0.1493323 0.04323404 1.498401e-03
singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
  single_method = "mr_wald_ratio", 
  all_method = c("mr_ivw",
  "mr_egger_regression", 
  "mr_weighted_mode",
  "mr_weighted_median"))
mr_forest_plot(singlesnp_results, exponentiate = FALSE)#saved as PDF, dimensions 8x8
## $`ieu-a-302.ieu-a-1055`
## Warning: Removed 1 rows containing missing values (geom_errorbarh).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1   ieu-a-302 ieu-a-1055
mr_scatter_plot(mr_results, dat2)
## $`ieu-a-302.ieu-a-1055`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1   ieu-a-302 ieu-a-1055
# Remove snps not used in MR
singlesnp_results.2=singlesnp_results[1:35,]
dim(singlesnp_results.2)
## [1] 35  9
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dat=dat2.1

# r2 and Fstat
dat$samplesize.outcome
##  [1]  98285 109975 110131 109491 103080 109855 109464 109491 110044 110127
## [11] 110008 102458 105528 109976 109717 110060 110040 110028 109963 105326
## [21] 110230 109938 110212 110024 105509 109952 110126 110226  93976 109511
## [31] 109418  96432 110068 109242  94514
dat$samplesize.exposure
##  [1] 174886 177680 177823 174454 177732 177750 177758 177783 163328 177773
## [11] 174742 177504 177712 151047 174734 172850 177813 177750 174704 175934
## [21] 176201 175846 177825 177798 177775 176257 177779 177782 177495 177778
## [31] 177486 168319 176205 177050 174573
n=mean(dat$samplesize.outcome)
k=35
p=dat$eaf.exposure
r2.1<-(2*(dat$beta.exposure^2)*p*(1-p))/1
r2=sum(r2.1)
r2 
## [1] 0.02142342
Fstat <- r2*(n-1-k)/((1-r2)*k)
Fstat
## [1] 67.21719
# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)

write.csv(mr_results,'res_TRI_urate.csv')
write.csv(singlesnp_results,'singlesnps_TRI_urate.csv')
write.csv(plt,'plt_TRI_urate.csv')
write.csv(het,'het_TRI_urate.csv')
write.csv(dat, 'SIMEX_TRI_urate.csv')
write.csv(dat, 'SIMEX_dat_TRI_urate.csv')
write.csv(outliers, 'outliers_TRI_urate.csv')

Multivariable MR: obesity, HDL, & TRI* on urate

*triglycerides

# MVMR of obesity, hdl, and triglycerides on urate
id_exposure=c('ieu-a-90','299', '302')
id_outcome=c('ieu-a-1055')

exposure_dat <- mv_extract_exposures(id_exposure)

# Remove the obesity on urate obesity outliers
remove_vars_obesity=c('rs7531118','rs9816226')

# Remove the HDL on urate outliers
remove_vars_HDL=c('rs102275', 'rs1047891', 'rs11065987', 'rs12133576', 'rs12748152',
                  'rs12801636','rs13099479','rs13107325', 'rs17145738', 'rs1980493', 'rs2075650', 
                  'rs2288912','rs2606736', 'rs3741414','rs3861397','rs6031587','rs6567160',
                  'rs731839', 'rs7607980','rs970548','rs9989419')

remove_vars_tri=c('rs10401969','rs11613352', 'rs11974409','rs12280753','rs1260326',
             'rs12748152','rs2043085','rs2247056','rs247616','rs439401','rs4587594',
             'rs4738684','rs4810479','rs588136','rs634869','rs7248104',
             'rs731839','rs7350481','rs8077889','rs948690')

exposure_dat_minus <- exposure_dat[ !exposure_dat$SNP %in% remove_vars_obesity, ]
exposure_dat_minus <- exposure_dat_minus[ !exposure_dat_minus$SNP %in% remove_vars_HDL, ]
exposure_dat_minus <- exposure_dat_minus[ !exposure_dat_minus$SNP %in% remove_vars_tri, ]
exposure_dat <- clump_data(exposure_dat_minus)

# Extract the outcome data
outcome_dat <- extract_outcome_data(exposure_dat$SNP, id_outcome)

# Harmonize
mvdat <- mv_harmonise_data(exposure_dat, outcome_dat)

# Perform MVMR
res <- mv_multiple(mvdat)
res
## $result
##   id.exposure                        exposure id.outcome                outcome
## 1   ieu-a-299 HDL cholesterol || id:ieu-a-299 ieu-a-1055 Urate || id:ieu-a-1055
## 2   ieu-a-302   Triglycerides || id:ieu-a-302 ieu-a-1055 Urate || id:ieu-a-1055
## 3    ieu-a-90  Obesity class 1 || id:ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055
##   nsnp           b         se         pval
## 1   62 -0.08421779 0.02316824 2.779230e-04
## 2   27  0.06835349 0.03252732 3.560428e-02
## 3   10  0.10375762 0.01659162 4.010601e-10
write.csv(res,"mvmr_obesity_HDL_TRI_urate.csv")
write.csv(exposure_dat, 'mvmr_dat_obesity_HDL_TRI_urate.csv')

Multivariable MR: obesity & urate on gout

# Multivariable MR of obesity and urate on gout
save.image()
id_exposure=c('ieu-a-90','ieu-a-1055')
id_outcome=c('UKB-b:13251')
exposure_dat <- mv_extract_exposures(id_exposure)

obesity_SNPs=read.csv('kept_obesity_SNPs.csv')
urate_SNPs=read.csv('kept_urate_SNPs_urate_gout.csv')

# Remove the outlier obesity and urate snps
exposure_dat2 <- exposure_dat[ exposure_dat$SNP %in% obesity_SNPs$x | exposure_dat$SNP %in% urate_SNPs$x, ]
exposure_dat2 <- clump_data(exposure_dat2)

# Extract the outcome data
outcome_dat <- extract_outcome_data(exposure_dat2$SNP, id_outcome)

# Harmonize
mvdat <- mv_harmonise_data(exposure_dat2, outcome_dat)

# Perform MVMR
res <- mv_multiple(mvdat)
res
## $result
##   id.exposure                       exposure  id.outcome
## 1  ieu-a-1055         Urate || id:ieu-a-1055 ukb-b-13251
## 2    ieu-a-90 Obesity class 1 || id:ieu-a-90 ukb-b-13251
##                                                          outcome nsnp
## 1 Non-cancer illness code, self-reported: gout || id:ukb-b-13251   15
## 2 Non-cancer illness code, self-reported: gout || id:ukb-b-13251   11
##              b           se          pval
## 1  0.029628063 0.0011481188 7.643542e-147
## 2 -0.001406938 0.0006667176  3.483717e-02
write.csv(res,"mvmr_urate_obesity_on_gout.csv")
write.csv(exposure_dat2, 'mvmr_dat_obesity_urate_gout.csv')

Multivariable MR: obesity & HDL on urate

# MVMR obesity and HDL on urate
id_exposure=c('ieu-a-90','299')
id_outcome=c('ieu-a-1055')
exposure_dat <- mv_extract_exposures(id_exposure)
dim(exposure_dat)
## [1] 196   9
# Remove the obesity on urate obesity outliers
remove_vars_obesity=c('rs7531118','rs9816226')

# Remove the HDL on urate outliers
remove_vars_HDL=c('rs102275', 'rs1047891', 'rs11065987', 'rs12133576', 'rs12748152',
'rs12801636','rs13099479','rs13107325', 'rs17145738', 'rs1980493', 'rs2075650', 
'rs2288912','rs2606736', 'rs3741414','rs3861397','rs6031587','rs6567160',
'rs731839', 'rs7607980','rs970548','rs9989419')

exposure_dat_minus <- exposure_dat[ !exposure_dat$SNP %in% remove_vars_obesity, ]
exposure_dat_minus <- exposure_dat_minus[ !exposure_dat_minus$SNP %in% remove_vars_HDL, ]
exposure_dat <- clump_data(exposure_dat_minus)

# Extract the outcome data
outcome_dat <- extract_outcome_data(exposure_dat$SNP, id_outcome)

# Harmonize
mvdat <- mv_harmonise_data(exposure_dat, outcome_dat)

# Perform MVMR
res <- mv_multiple(mvdat)
res
## $result
##   id.exposure                        exposure id.outcome                outcome
## 1   ieu-a-299 HDL cholesterol || id:ieu-a-299 ieu-a-1055 Urate || id:ieu-a-1055
## 2    ieu-a-90  Obesity class 1 || id:ieu-a-90 ieu-a-1055 Urate || id:ieu-a-1055
##   nsnp           b         se         pval
## 1   63 -0.09491173 0.02028945 2.898497e-06
## 2   10  0.11225804 0.01669030 1.744578e-11
write.csv(res,"mvmr_obesity_hdl_on_urate.csv")
write.csv(exposure_dat, 'mvmr_dat_obesity_hdl_on_urate.csv')