create variable look up codes and save methods list

### exposures

dep<- "ieu-b-102"
scz<- "ieu-b-42"
#bip<- "ieu-b-41"
adhd<- "ieu-a-1183"

### outcomes
sleepdur<- "ebi-a-GCST003839"
napping<- "ebi-a-GCST011494"
chrono<- "ukb-b-4956"

methods<- c('mr_two_sample_ml', 'mr_egger_regression', 'mr_weighted_mode', 'mr_weighted_median', 'mr_ivw')

ivw_p<- .05/8
set.seed(222)

psych_sleep_MR<- data.frame()

models to run:

depression

depression_exp_dat <- extract_instruments(outcomes = dep)
## API: public: http://gwas-api.mrcieu.ac.uk/

napping

napping_out_dat <- read_outcome_data(
    snps = depression_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/napsR.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "BETA",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "P",
    samplesize_col = "N"
)
## No phenotype name specified, defaulting to 'outcome'.
napping_out_dat$outcome<- rep('No napping (napping GWAS)', 1, nrow(napping_out_dat))

depression_napping_dat <- harmonise_data(depression_exp_dat, napping_out_dat)
## Harmonising Major depression || id:ieu-b-102 (ieu-b-102) and No napping (napping GWAS) (VfA804)
## Removing the following SNPs for being palindromic with intermediate allele frequencies:
## rs2876520, rs4730387
res <- mr(depression_napping_dat, method_list = methods)
## Analysing 'ieu-b-102' on 'VfA804'
res
##   id.exposure id.outcome                   outcome
## 1   ieu-b-102     VfA804 No napping (napping GWAS)
## 2   ieu-b-102     VfA804 No napping (napping GWAS)
## 3   ieu-b-102     VfA804 No napping (napping GWAS)
## 4   ieu-b-102     VfA804 No napping (napping GWAS)
## 5   ieu-b-102     VfA804 No napping (napping GWAS)
##                           exposure                    method nsnp           b
## 1 Major depression || id:ieu-b-102        Maximum likelihood   47 -0.06999166
## 2 Major depression || id:ieu-b-102                  MR Egger   47 -0.19743729
## 3 Major depression || id:ieu-b-102             Weighted mode   47 -0.04417438
## 4 Major depression || id:ieu-b-102           Weighted median   47 -0.05030157
## 5 Major depression || id:ieu-b-102 Inverse variance weighted   47 -0.06330905
##            se         pval
## 1 0.007229323 3.609016e-22
## 2 0.080378106 1.795887e-02
## 3 0.023572640 6.729462e-02
## 4 0.011752231 1.867486e-05
## 5 0.015593425 4.907618e-05
psych_sleep_MR<- rbind(psych_sleep_MR, res)

mr_heterogeneity(depression_napping_dat)
##   id.exposure id.outcome                   outcome
## 1   ieu-b-102     VfA804 No napping (napping GWAS)
## 2   ieu-b-102     VfA804 No napping (napping GWAS)
##                           exposure                    method        Q Q_df
## 1 Major depression || id:ieu-b-102                  MR Egger 241.7814   45
## 2 Major depression || id:ieu-b-102 Inverse variance weighted 257.3041   46
##         Q_pval
## 1 9.469448e-29
## 2 3.665733e-31
mr_pleiotropy_test(depression_napping_dat)
##   id.exposure id.outcome                   outcome
## 1   ieu-b-102     VfA804 No napping (napping GWAS)
##                           exposure egger_intercept         se       pval
## 1 Major depression || id:ieu-b-102     0.004142462 0.00243714 0.09608645
p1 <- mr_scatter_plot(res, depression_napping_dat)
p1
## $`ieu-b-102.VfA804`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1   ieu-b-102     VfA804
res_single <- mr_singlesnp(depression_napping_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##                            exposure                   outcome id.exposure
## 1  Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 2  Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 3  Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 4  Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 5  Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 6  Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 7  Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 8  Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 9  Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 10 Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 11 Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 12 Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 13 Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
## 14 Major depression || id:ieu-b-102 No napping (napping GWAS)   ieu-b-102
##    id.outcome samplesize                             SNP           b
## 1      VfA804     452633                      rs10913112 -0.18214885
## 2      VfA804     452633                      rs17641524 -0.14454067
## 3      VfA804     452633                       rs1950829 -0.17748923
## 4      VfA804     452633                       rs2232423 -0.14813790
## 5      VfA804     452633                        rs247910  0.31002363
## 6      VfA804     452633                         rs30266 -0.34985246
## 7      VfA804     452633                       rs4936276 -0.20496295
## 8      VfA804     452633                      rs62535714 -0.16461386
## 9      VfA804     452633                      rs66511648 -0.17439899
## 10     VfA804     452633                      rs72948506 -0.16635434
## 11     VfA804     452633                       rs9536381  0.16584078
## 12     VfA804     452633        All - Maximum likelihood -0.06999166
## 13     VfA804     452633           All - Weighted median -0.05030157
## 14     VfA804     452633 All - Inverse variance weighted -0.06330905
##             se            p
## 1  0.047761069 1.368730e-04
## 2  0.049636333 3.591304e-03
## 3  0.040945791 1.459347e-05
## 4  0.030541613 1.232359e-06
## 5  0.051509705 1.757776e-09
## 6  0.035314481 3.890176e-23
## 7  0.044961151 5.147914e-06
## 8  0.048553687 6.980412e-04
## 9  0.045405724 1.225825e-04
## 10 0.049899245 8.566678e-04
## 11 0.050703529 1.072461e-03
## 12 0.007229323 3.609016e-22
## 13 0.011054451 5.355780e-06
## 14 0.015593425 4.907618e-05

sleep duration

sleep_dur_out_dat <- read_outcome_data(
    snps = depression_exp_dat$SNP,
    filename = "/Users/claire/Desktop/sleepGWAS/sleep_rr_upd/selfsleepdur.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "se",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "p",
    samplesize_col = "N"
)
## No phenotype name specified, defaulting to 'outcome'.
sleep_dur_out_dat$outcome<- rep('Sleep Duration', 1, nrow(sleep_dur_out_dat))


depression_dur_dat <- harmonise_data(depression_exp_dat, sleep_dur_out_dat)
## Harmonising Major depression || id:ieu-b-102 (ieu-b-102) and Sleep Duration (ip3Q2k)
## Removing the following SNPs for being palindromic with intermediate allele frequencies:
## rs2876520, rs4730387
res <- mr(depression_dur_dat, method_list = methods)
## Analysing 'ieu-b-102' on 'ip3Q2k'
res
##   id.exposure id.outcome        outcome                         exposure
## 1   ieu-b-102     ip3Q2k Sleep Duration Major depression || id:ieu-b-102
## 2   ieu-b-102     ip3Q2k Sleep Duration Major depression || id:ieu-b-102
## 3   ieu-b-102     ip3Q2k Sleep Duration Major depression || id:ieu-b-102
## 4   ieu-b-102     ip3Q2k Sleep Duration Major depression || id:ieu-b-102
## 5   ieu-b-102     ip3Q2k Sleep Duration Major depression || id:ieu-b-102
##                      method nsnp            b         se      pval
## 1        Maximum likelihood   44  0.008780569 0.01401710 0.5310405
## 2                  MR Egger   44  0.244041363 0.17986623 0.1820948
## 3             Weighted mode   44 -0.016617064 0.03956786 0.6766031
## 4           Weighted median   44  0.009038949 0.02383137 0.7044741
## 5 Inverse variance weighted   44  0.007803268 0.03447098 0.8209120
psych_sleep_MR<- rbind(psych_sleep_MR, res)

mr_heterogeneity(depression_dur_dat)
##   id.exposure id.outcome        outcome                         exposure
## 1   ieu-b-102     ip3Q2k Sleep Duration Major depression || id:ieu-b-102
## 2   ieu-b-102     ip3Q2k Sleep Duration Major depression || id:ieu-b-102
##                      method        Q Q_df       Q_pval
## 1                  MR Egger 305.5986   42 9.927876e-42
## 2 Inverse variance weighted 318.6199   43 9.402937e-44
mr_pleiotropy_test(depression_dur_dat)
##   id.exposure id.outcome        outcome                         exposure
## 1   ieu-b-102     ip3Q2k Sleep Duration Major depression || id:ieu-b-102
##   egger_intercept          se      pval
## 1    -0.007322928 0.005474041 0.1881751
p1 <- mr_scatter_plot(res, depression_dur_dat)
p1
## $`ieu-b-102.ip3Q2k`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1   ieu-b-102     ip3Q2k
res_single <- mr_singlesnp(depression_dur_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##                            exposure        outcome id.exposure id.outcome
## 1  Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 2  Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 3  Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 4  Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 5  Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 6  Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 7  Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 8  Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 9  Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 10 Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 11 Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 12 Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
## 13 Major depression || id:ieu-b-102 Sleep Duration   ieu-b-102     ip3Q2k
##    samplesize        SNP          b         se            p
## 1      446118 rs10235664 -0.6037852 0.09703704 4.901886e-10
## 2      446118 rs12967143 -0.4935507 0.07217710 8.027764e-12
## 3      446118 rs17641524  0.3469733 0.09265533 1.805558e-04
## 4      446118   rs198457  0.2726159 0.09271016 3.276642e-03
## 5      446118  rs2232423  0.2571452 0.05698742 6.412400e-06
## 6      446118 rs28541419  0.2563836 0.09282534 5.744854e-03
## 7      446118  rs3807865  0.2898523 0.07400226 8.973106e-05
## 8      446118  rs4936276  0.4274065 0.08385504 3.451358e-07
## 9      446118 rs61914045 -0.2963570 0.09153528 1.205288e-03
## 10     446118 rs62535714  0.5803717 0.09060147 1.496181e-10
## 11     446118 rs66511648 -0.2862175 0.08474310 7.315518e-04
## 12     446118  rs7152906  0.2989876 0.08769922 6.514444e-04
## 13     446118  rs9536381 -0.3840435 0.09457451 4.891513e-05

insomnia

insom_out_dat <- read_outcome_data(
    snps = depression_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/insomnia_withNeffR.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "se",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "p",
    samplesize_col = "Neff"
)
## No phenotype name specified, defaulting to 'outcome'.
insom_out_dat$outcome<- rep('Insomnia', 1, nrow(insom_out_dat))

depression_insom_dat <- harmonise_data(depression_exp_dat, insom_out_dat)
## Harmonising Major depression || id:ieu-b-102 (ieu-b-102) and Insomnia (088DHx)
## Removing the following SNPs for being palindromic with intermediate allele frequencies:
## rs2876520, rs4730387
res <- mr(depression_insom_dat, method_list = methods)
## Analysing 'ieu-b-102' on '088DHx'
res
##   id.exposure id.outcome  outcome                         exposure
## 1   ieu-b-102     088DHx Insomnia Major depression || id:ieu-b-102
## 2   ieu-b-102     088DHx Insomnia Major depression || id:ieu-b-102
## 3   ieu-b-102     088DHx Insomnia Major depression || id:ieu-b-102
## 4   ieu-b-102     088DHx Insomnia Major depression || id:ieu-b-102
## 5   ieu-b-102     088DHx Insomnia Major depression || id:ieu-b-102
##                      method nsnp           b         se         pval
## 1        Maximum likelihood   44 -0.08201878 0.00842885 2.230087e-22
## 2                  MR Egger   44 -0.03882477 0.08659833 6.562183e-01
## 3             Weighted mode   44 -0.10111622 0.03257142 3.367017e-03
## 4           Weighted median   44 -0.07538862 0.01377683 4.446583e-08
## 5 Inverse variance weighted   44 -0.07827438 0.01629662 1.562229e-06
psych_sleep_MR<- rbind(psych_sleep_MR, res)


mr_heterogeneity(depression_insom_dat)
##   id.exposure id.outcome  outcome                         exposure
## 1   ieu-b-102     088DHx Insomnia Major depression || id:ieu-b-102
## 2   ieu-b-102     088DHx Insomnia Major depression || id:ieu-b-102
##                      method        Q Q_df       Q_pval
## 1                  MR Egger 189.9329   42 1.054575e-20
## 2 Inverse variance weighted 190.9065   43 1.547770e-20
mr_pleiotropy_test(depression_insom_dat)
##   id.exposure id.outcome  outcome                         exposure
## 1   ieu-b-102     088DHx Insomnia Major depression || id:ieu-b-102
##   egger_intercept          se      pval
## 1    -0.001222882 0.002635564 0.6450491
p1 <- mr_scatter_plot(res, depression_insom_dat)
p1
## $`ieu-b-102.088DHx`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1   ieu-b-102     088DHx
res_single <- mr_singlesnp(depression_insom_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##                            exposure  outcome id.exposure id.outcome samplesize
## 1  Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 2  Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 3  Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 4  Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 5  Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 6  Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 7  Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 8  Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 9  Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 10 Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 11 Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 12 Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 13 Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 14 Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
## 15 Major depression || id:ieu-b-102 Insomnia   ieu-b-102     088DHx   259365.1
##                                SNP           b         se            p
## 1                       rs10235664 -0.21517333 0.05934741 2.882279e-04
## 2                       rs10913112  0.18157214 0.05445191 8.544021e-04
## 3                       rs12967143 -0.26305739 0.04410319 2.452506e-09
## 4                        rs1950829 -0.14247744 0.04669596 2.279494e-03
## 5                        rs2111592 -0.15740875 0.05669506 5.496278e-03
## 6                        rs2522831 -0.21240250 0.05765750 2.297231e-04
## 7                        rs2568958 -0.14533560 0.03685838 8.044114e-05
## 8                          rs30266 -0.24509317 0.04028388 1.171109e-09
## 9                       rs66511648 -0.15749663 0.05172121 2.325974e-03
## 10                       rs9536381 -0.36714314 0.05767961 1.950157e-10
## 11                       rs9831648 -0.22493151 0.05666438 7.201168e-05
## 12        All - Maximum likelihood -0.08201878 0.00842885 2.230087e-22
## 13             All - Weighted mode -0.10111622 0.03268873 3.471823e-03
## 14           All - Weighted median -0.07538862 0.01411045 9.154637e-08
## 15 All - Inverse variance weighted -0.07827438 0.01629662 1.562229e-06

efficiency factor

effic_out_dat <- read_outcome_data(
    snps = depression_exp_dat$SNP,
    filename = "/Users/claire/Desktop/sleepGWAS/sleep_rr_upd/efficiency_gsem.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "est",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "Pval_Estimate",
    samplesize_col = "N_eff"
)
## No phenotype name specified, defaulting to 'outcome'.
effic_out_dat$outcome<- rep('Efficiency factor', 1, nrow(effic_out_dat))

dep_effic_dat <- harmonise_data(depression_exp_dat, effic_out_dat)
## Harmonising Major depression || id:ieu-b-102 (ieu-b-102) and Efficiency factor (0e7UX0)
## Removing the following SNPs for being palindromic with intermediate allele frequencies:
## rs2876520, rs4730387
res <- mr(dep_effic_dat, method_list = methods)
## Analysing 'ieu-b-102' on '0e7UX0'
res
##   id.exposure id.outcome           outcome                         exposure
## 1   ieu-b-102     0e7UX0 Efficiency factor Major depression || id:ieu-b-102
## 2   ieu-b-102     0e7UX0 Efficiency factor Major depression || id:ieu-b-102
## 3   ieu-b-102     0e7UX0 Efficiency factor Major depression || id:ieu-b-102
## 4   ieu-b-102     0e7UX0 Efficiency factor Major depression || id:ieu-b-102
## 5   ieu-b-102     0e7UX0 Efficiency factor Major depression || id:ieu-b-102
##                      method nsnp          b         se      pval
## 1        Maximum likelihood   43 0.01979097 0.01522363 0.1935953
## 2                  MR Egger   43 0.06406791 0.10330139 0.5385558
## 3             Weighted mode   43 0.05937683 0.05138479 0.2544037
## 4           Weighted median   43 0.03256346 0.02264556 0.1504447
## 5 Inverse variance weighted   43 0.01896666 0.01960329 0.3332820
psych_sleep_MR<- rbind(psych_sleep_MR, res)


mr_heterogeneity(dep_effic_dat)
##   id.exposure id.outcome           outcome                         exposure
## 1   ieu-b-102     0e7UX0 Efficiency factor Major depression || id:ieu-b-102
## 2   ieu-b-102     0e7UX0 Efficiency factor Major depression || id:ieu-b-102
##                      method        Q Q_df      Q_pval
## 1                  MR Egger 72.49468   41 0.001746452
## 2 Inverse variance weighted 72.84457   42 0.002206043
mr_pleiotropy_test(dep_effic_dat)
##   id.exposure id.outcome           outcome                         exposure
## 1   ieu-b-102     0e7UX0 Efficiency factor Major depression || id:ieu-b-102
##   egger_intercept          se      pval
## 1    -0.001396697 0.003139769 0.6587745
p1 <- mr_scatter_plot(res, dep_effic_dat)
p1
## $`ieu-b-102.0e7UX0`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1   ieu-b-102     0e7UX0
res_single <- mr_singlesnp(dep_effic_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##                           exposure           outcome id.exposure id.outcome
## 1 Major depression || id:ieu-b-102 Efficiency factor   ieu-b-102     0e7UX0
##   samplesize       SNP          b         se           p
## 1   284894.1 rs4497414 -0.2525844 0.09091585 0.005465739

alertness factor

alert_out_dat <- read_outcome_data(
    snps = depression_exp_dat$SNP,
    filename = "/Users/claire/Desktop/sleepGWAS/sleep_rr_upd/alertness_gsem.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "est",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "Pval_Estimate",
    samplesize_col = "Neff"
)
## No phenotype name specified, defaulting to 'outcome'.
alert_out_dat$outcome<- rep("Alertness Factor", 1, nrow(alert_out_dat))


depression_alert_dat <- harmonise_data(depression_exp_dat, alert_out_dat)
## Harmonising Major depression || id:ieu-b-102 (ieu-b-102) and Alertness Factor (TUFCQf)
## Removing the following SNPs for being palindromic with intermediate allele frequencies:
## rs2876520, rs4730387
res <- mr(depression_alert_dat, method_list = methods)
## Analysing 'ieu-b-102' on 'TUFCQf'
res
##   id.exposure id.outcome          outcome                         exposure
## 1   ieu-b-102     TUFCQf Alertness Factor Major depression || id:ieu-b-102
## 2   ieu-b-102     TUFCQf Alertness Factor Major depression || id:ieu-b-102
## 3   ieu-b-102     TUFCQf Alertness Factor Major depression || id:ieu-b-102
## 4   ieu-b-102     TUFCQf Alertness Factor Major depression || id:ieu-b-102
## 5   ieu-b-102     TUFCQf Alertness Factor Major depression || id:ieu-b-102
##                      method nsnp            b          se         pval
## 1        Maximum likelihood   43 -0.018807460 0.002328194 6.576602e-16
## 2                  MR Egger   43 -0.042336408 0.024773934 9.502844e-02
## 3             Weighted mode   43 -0.006745758 0.007344112 3.635886e-01
## 4           Weighted median   43 -0.010636670 0.003861640 5.879204e-03
## 5 Inverse variance weighted   43 -0.017165874 0.004679156 2.438978e-04
psych_sleep_MR<- rbind(psych_sleep_MR, res)


mr_heterogeneity(depression_alert_dat)
##   id.exposure id.outcome          outcome                         exposure
## 1   ieu-b-102     TUFCQf Alertness Factor Major depression || id:ieu-b-102
## 2   ieu-b-102     TUFCQf Alertness Factor Major depression || id:ieu-b-102
##                      method        Q Q_df       Q_pval
## 1                  MR Egger 195.8768   41 4.484110e-22
## 2 Inverse variance weighted 200.9906   42 1.280317e-22
mr_pleiotropy_test(depression_alert_dat)
##   id.exposure id.outcome          outcome                         exposure
## 1   ieu-b-102     TUFCQf Alertness Factor Major depression || id:ieu-b-102
##   egger_intercept           se      pval
## 1    0.0007773394 0.0007513439 0.3069204
p1 <- mr_scatter_plot(res, depression_alert_dat)
p1
## $`ieu-b-102.TUFCQf`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1   ieu-b-102     TUFCQf
res_single <- mr_singlesnp(depression_alert_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##                            exposure          outcome id.exposure id.outcome
## 1  Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 2  Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 3  Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 4  Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 5  Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 6  Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 7  Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 8  Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 9  Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 10 Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 11 Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 12 Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 13 Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 14 Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
## 15 Major depression || id:ieu-b-102 Alertness Factor   ieu-b-102     TUFCQf
##    samplesize                             SNP           b          se
## 1          NA                      rs10913112 -0.05354468 0.015002810
## 2          NA                       rs1950829 -0.05362157 0.012974782
## 3          NA                       rs2232423 -0.05142279 0.009973326
## 4          NA                       rs2522831 -0.04601827 0.015915687
## 5          NA                         rs30266 -0.09385749 0.011469009
## 6          NA                       rs3807865 -0.03521991 0.012410541
## 7          NA                       rs4936276 -0.06046558 0.014105362
## 8          NA                      rs66511648 -0.05914058 0.014372261
## 9          NA                       rs7152906 -0.04152069 0.014675531
## 10         NA                      rs72948506 -0.05058509 0.015535838
## 11         NA                       rs7725715 -0.04873891 0.013453304
## 12         NA                       rs9536381  0.05361153 0.016007250
## 13         NA        All - Maximum likelihood -0.01880746 0.002328194
## 14         NA           All - Weighted median -0.01063667 0.003804592
## 15         NA All - Inverse variance weighted -0.01716587 0.004679156
##               p
## 1  3.583784e-04
## 2  3.584434e-05
## 3  2.522381e-07
## 4  3.835562e-03
## 5  2.755457e-16
## 6  4.541100e-03
## 7  1.813396e-05
## 8  3.873284e-05
## 9  4.665776e-03
## 10 1.129834e-03
## 11 2.914075e-04
## 12 8.104441e-04
## 13 6.576602e-16
## 14 5.178026e-03
## 15 2.438978e-04

anxiety

## exposure data
anx_full<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/anxiety_withNeff.txt.gz", header=T, data.table = F)
anx<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/FUMA_out/anxiety//GenomicRiskLoci.txt", header=T, data.table=F)

anx<- anx %>% select(rsID) %>%
  rename(SNP=rsID)

anx_full<- anx_full %>% 
  right_join(anx)
## Joining with `by = join_by(SNP)`
anxiety_exp_dat <- format_data(anx_full, type="exposure",
                              snp_col = "SNP", beta_col = "effect", se_col = "SE", effect_allele_col = "A1", other_allele_col = "A2", eaf_col = "MAF", pval_col = "P", samplesize_col = "Neff")
## No phenotype name specified, defaulting to 'exposure'.
anxiety_exp_dat$exposure<- rep("Anxiety", 1, nrow(anxiety_exp_dat))

anxiety_exp_dat <- anxiety_exp_dat %>% as.data.table
anxiety_exp_dat <- anxiety_exp_dat[, F_STAT := qf(
  pval.exposure,
  df1=1,
  df2=samplesize.exposure,
lower.tail=F
)] # this piece calculates the F statistic of each SNP in the IV and adds a column called F_STAT
anxiety_exp_dat$F_STAT %>% mean # to calculate the overall F stat of your exposure
## [1] 34.68837
anxiety_exp_dat$F_STAT %>% min # you can use this to check whether any individual SNPs fall under F of 10 which would make them inva
## [1] 29.81161

insomnia

insom_out_dat <- read_outcome_data(
    snps = anxiety_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/insomnia_withNeffR.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "se",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "p",
    samplesize_col = "Neff"
)
## No phenotype name specified, defaulting to 'outcome'.
insom_out_dat$outcome<- rep('Insomnia', 1, nrow(insom_out_dat))

anx_insom_dat <- harmonise_data(anxiety_exp_dat, insom_out_dat)
## Harmonising Anxiety (a4ubWB) and Insomnia (kawpsj)
res <- mr(anx_insom_dat, method_list = methods)
## Analysing 'a4ubWB' on 'kawpsj'
res
##   id.exposure id.outcome  outcome exposure                    method nsnp
## 1      a4ubWB     kawpsj Insomnia  Anxiety        Maximum likelihood    5
## 2      a4ubWB     kawpsj Insomnia  Anxiety                  MR Egger    5
## 3      a4ubWB     kawpsj Insomnia  Anxiety             Weighted mode    5
## 4      a4ubWB     kawpsj Insomnia  Anxiety           Weighted median    5
## 5      a4ubWB     kawpsj Insomnia  Anxiety Inverse variance weighted    5
##              b          se      pval
## 1  0.007699010 0.005504224 0.1618892
## 2 -0.007793231 0.123253014 0.9535609
## 3  0.012227685 0.011625789 0.3522449
## 4  0.008446989 0.007691130 0.2720837
## 5  0.007316268 0.007913153 0.3551893
psych_sleep_MR<- rbind(psych_sleep_MR, res)


psych_sleep_MR<- rbind(psych_sleep_MR, res)


mr_heterogeneity(anx_insom_dat)
##   id.exposure id.outcome  outcome exposure                    method        Q
## 1      a4ubWB     kawpsj Insomnia  Anxiety                  MR Egger 8.698467
## 2      a4ubWB     kawpsj Insomnia  Anxiety Inverse variance weighted 8.742281
##   Q_df     Q_pval
## 1    3 0.03358056
## 2    4 0.06787414
mr_pleiotropy_test(anx_insom_dat)
##   id.exposure id.outcome  outcome exposure egger_intercept         se      pval
## 1      a4ubWB     kawpsj Insomnia  Anxiety     0.001935641 0.01574641 0.9099386
p1 <- mr_scatter_plot(res, anx_insom_dat)
p1
## $a4ubWB.kawpsj

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1      a4ubWB     kawpsj
res_single <- mr_singlesnp(anx_insom_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
## [1] exposure    outcome     id.exposure id.outcome  samplesize  SNP        
## [7] b           se          p          
## <0 rows> (or 0-length row.names)

ADHD

adhd_exp_dat <- extract_instruments(outcomes = adhd)


adhd_exp_dat <- adhd_exp_dat %>% as.data.table
adhd_exp_dat <- adhd_exp_dat[, F_STAT := qf(
  pval.exposure,
  df1=1,
  df2=samplesize.exposure,
lower.tail=F
)] # this piece calculates the F statistic of each SNP in the IV and adds a column called F_STAT
adhd_exp_dat$F_STAT %>% mean # to calculate the overall F stat of your exposure
## [1] 33.99208
adhd_exp_dat$F_STAT %>% min # 
## [1] 30.67278

insomnia

insom_out_dat <- read_outcome_data(
    snps = adhd_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/insomnia_withNeffR.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "se",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "p",
    samplesize_col = "Neff"
)
## No phenotype name specified, defaulting to 'outcome'.
insom_out_dat$outcome<- rep('Insomnia', 1, nrow(insom_out_dat))

adhd_insom_dat <- harmonise_data(adhd_exp_dat, insom_out_dat)
## Harmonising ADHD || id:ieu-a-1183 (ieu-a-1183) and Insomnia (yvWsuG)
## Removing the following SNPs for being palindromic with intermediate allele frequencies:
## rs11591402
res <- mr(adhd_insom_dat, method_list = methods)
## Analysing 'ieu-a-1183' on 'yvWsuG'
res
##   id.exposure id.outcome  outcome              exposure
## 1  ieu-a-1183     yvWsuG Insomnia ADHD || id:ieu-a-1183
## 2  ieu-a-1183     yvWsuG Insomnia ADHD || id:ieu-a-1183
## 3  ieu-a-1183     yvWsuG Insomnia ADHD || id:ieu-a-1183
## 4  ieu-a-1183     yvWsuG Insomnia ADHD || id:ieu-a-1183
## 5  ieu-a-1183     yvWsuG Insomnia ADHD || id:ieu-a-1183
##                      method nsnp            b          se       pval
## 1        Maximum likelihood    9 -0.011446985 0.006679646 0.08658172
## 2                  MR Egger    9  0.063445411 0.068729960 0.38666174
## 3             Weighted mode    9  0.006047539 0.010906590 0.59440455
## 4           Weighted median    9  0.006346930 0.009328057 0.49624306
## 5 Inverse variance weighted    9 -0.009251920 0.016516142 0.57536049
psych_sleep_MR<- rbind(psych_sleep_MR, res)


mr_heterogeneity(adhd_insom_dat)
##   id.exposure id.outcome  outcome              exposure
## 1  ieu-a-1183     yvWsuG Insomnia ADHD || id:ieu-a-1183
## 2  ieu-a-1183     yvWsuG Insomnia ADHD || id:ieu-a-1183
##                      method        Q Q_df       Q_pval
## 1                  MR Egger 52.53333    7 4.582935e-09
## 2 Inverse variance weighted 61.43170    8 2.439353e-10
mr_pleiotropy_test(adhd_insom_dat)
##   id.exposure id.outcome  outcome              exposure egger_intercept
## 1  ieu-a-1183     yvWsuG Insomnia ADHD || id:ieu-a-1183    -0.006737961
##            se      pval
## 1 0.006187877 0.3122611
p1 <- mr_scatter_plot(res, adhd_insom_dat)
p1
## $`ieu-a-1183.yvWsuG`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1  ieu-a-1183     yvWsuG
res_single <- mr_singlesnp(adhd_insom_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##                exposure  outcome id.exposure id.outcome samplesize        SNP
## 1 ADHD || id:ieu-a-1183 Insomnia  ieu-a-1183     yvWsuG   259365.1 rs10262192
## 2 ADHD || id:ieu-a-1183 Insomnia  ieu-a-1183     yvWsuG   259365.1  rs1427829
## 3 ADHD || id:ieu-a-1183 Insomnia  ieu-a-1183     yvWsuG   259365.1  rs4916723
##             b         se            p
## 1 -0.10543973 0.01906440 3.189288e-08
## 2 -0.06611766 0.01738009 1.422527e-04
## 3  0.06813146 0.01846585 2.246146e-04

ASB

## exposure data
asb_full<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/broadABC2022_Final_CombinedSex.TBL.gz", header=T, data.table = F)
asb<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/FUMA_out/ASB/GenomicRiskLoci.txt", header=T, data.table=F)

asb<- asb %>% select(rsID) %>%
  rename(SNP=rsID)

asb_full<- asb_full %>% 
  rename(SNP=rsid) %>%
  right_join(asb)
## Joining with `by = join_by(SNP)`
asb_exp_dat <- format_data(asb_full, type="exposure",
                              snp_col = "SNP", beta_col = "beta", se_col = "standard_error", effect_allele_col = "effect_allele", other_allele_col = "other_allele", eaf_col = "effect_allele_frequency", pval_col = "p_value", samplesize_col = "n")
## No phenotype name specified, defaulting to 'exposure'.
asb_exp_dat$exposure<- rep("Anti Social Behavior", 1, nrow(asb_exp_dat))

asb_exp_dat <- asb_exp_dat %>% as.data.table
asb_exp_dat <- asb_exp_dat[, F_STAT := qf(
  pval.exposure,
  df1=1,
  df2=samplesize.exposure,
lower.tail=F
)] # this piece calculates the F statistic of each SNP in the IV and adds a column called F_STAT
asb_exp_dat$F_STAT %>% mean # to calculate the overall F stat of your exposure
## [1] 38.23419
asb_exp_dat$F_STAT %>% min # you can use this to check whether any individual SNPs fall under F of 10 which would make them inva
## [1] 38.23419

actigraphy duration

actidur_out_dat <- read_outcome_data(
    snps = asb_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/actisleepdur.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "p",
    samplesize_col = "N"
)
## No phenotype name specified, defaulting to 'outcome'.
actidur_out_dat$outcome<- rep('Actigraphy Sleep Duration', 1, nrow(actidur_out_dat))

asb_actidur_dat <- harmonise_data(asb_exp_dat, actidur_out_dat)
## Harmonising Anti Social Behavior (yxUmHB) and Actigraphy Sleep Duration (LCeoUO)
res <- mr(asb_actidur_dat, method_list = methods)
## Analysing 'yxUmHB' on 'LCeoUO'
res
## [1] id.exposure id.outcome  outcome     exposure    method      nsnp       
## [7] b           se          pval       
## <0 rows> (or 0-length row.names)
psych_sleep_MR<- rbind(psych_sleep_MR, res)


mr_heterogeneity(asb_actidur_dat)
## Not enough SNPs available for Heterogeneity analysis of 'yxUmHB' on 'LCeoUO'
## data frame with 0 columns and 0 rows
mr_pleiotropy_test(asb_actidur_dat)
## Not enough SNPs available for pleiotropy analysis of 'yxUmHB' on 'LCeoUO'
## data frame with 0 columns and 0 rows
p1 <- mr_scatter_plot(res, asb_actidur_dat)
p1
## $yxUmHB.LCeoUO

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1      yxUmHB     LCeoUO
res_single <- mr_singlesnp(asb_actidur_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]

res_single %>% filter(p<ivw_p)
##               exposure                   outcome id.exposure id.outcome
## 1 Anti Social Behavior Actigraphy Sleep Duration      yxUmHB     LCeoUO
##   samplesize        SNP            b         se           p
## 1      85449 rs12536335 -0.003274664 0.00110498 0.003041128

efficiency indicator

effic_out_dat <- read_outcome_data(
    snps = asb_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/effic.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "p",
    samplesize_col = "N"
)
## No phenotype name specified, defaulting to 'outcome'.
effic_out_dat$outcome<- rep('Efficiency', 1, nrow(effic_out_dat))

asb_effic_dat <- harmonise_data(asb_exp_dat, effic_out_dat)
## Harmonising Anti Social Behavior (yxUmHB) and Efficiency (vqDN1l)
res <- mr(asb_effic_dat, method_list = methods)
## Analysing 'yxUmHB' on 'vqDN1l'
res
## [1] id.exposure id.outcome  outcome     exposure    method      nsnp       
## [7] b           se          pval       
## <0 rows> (or 0-length row.names)
psych_sleep_MR<- rbind(psych_sleep_MR, res)


mr_heterogeneity(asb_effic_dat)
## Not enough SNPs available for Heterogeneity analysis of 'yxUmHB' on 'vqDN1l'
## data frame with 0 columns and 0 rows
mr_pleiotropy_test(asb_effic_dat)
## Not enough SNPs available for pleiotropy analysis of 'yxUmHB' on 'vqDN1l'
## data frame with 0 columns and 0 rows
p1 <- mr_scatter_plot(res, asb_effic_dat)
p1
## $yxUmHB.vqDN1l

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1      yxUmHB     vqDN1l
res_single <- mr_singlesnp(asb_effic_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]

res_single %>% filter(p<ivw_p)
##               exposure    outcome id.exposure id.outcome samplesize        SNP
## 1 Anti Social Behavior Efficiency      yxUmHB     vqDN1l      84810 rs12536335
##              b          se           p
## 1 -0.003220459 0.001098143 0.003360912

drinks per week

## exposure data
dpw_full<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/DrinksPerWeek.txt.gz", header=T, data.table = F)
dpw<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/FUMA_out/dpw/GenomicRiskLoci.txt", header=T, data.table=F)

dpw<- dpw %>% select(rsID) %>%
  rename(SNP=rsID)

dpw_full<- dpw_full %>% 
  rename(SNP=RSID) %>%
  right_join(dpw)
## Joining with `by = join_by(SNP)`
dpw_exp_dat <- format_data(dpw_full, type="exposure",
                              snp_col = "SNP", beta_col = "BETA", se_col = "SE", effect_allele_col = "REF", other_allele_col = "ALT", eaf_col = "AF", pval_col = "PVALUE", samplesize_col = "EFFECTIVE_N")
## No phenotype name specified, defaulting to 'exposure'.
dpw_exp_dat$exposure<- rep("Drinks per Week", 1, nrow(dpw_exp_dat))


dpw_exp_dat <- dpw_exp_dat %>% as.data.table
dpw_exp_dat <- dpw_exp_dat[, F_STAT := qf(
  pval.exposure,
  df1=1,
  df2=samplesize.exposure,
lower.tail=F
)] # this piece calculates the F statistic of each SNP in the IV and adds a column called F_STAT
dpw_exp_dat$F_STAT %>% mean # to calculate the overall F stat of your exposure
## [1] 69.33978
dpw_exp_dat$F_STAT %>% min # you
## [1] 29.79592

chronotype

chrono_out_dat <- read_outcome_data(
    snps = dpw_exp_dat$SNP,
    filename = "/Users/claire/Desktop/sleepGWAS/sleep_rr_upd/circadian_pref_gsem.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "est",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "Pval_Estimate",
    samplesize_col = "Neff"
)
## No phenotype name specified, defaulting to 'outcome'.
chrono_out_dat$outcome<- rep("Circadian Preference Factor", 1, nrow(chrono_out_dat))


dpw_chrono_dat <- harmonise_data(dpw_exp_dat, chrono_out_dat)
## Harmonising Drinks per Week (moyRYB) and Circadian Preference Factor (uxmElJ)
res <- mr(dpw_chrono_dat, method_list = methods)
## Analysing 'moyRYB' on 'uxmElJ'
res
##   id.exposure id.outcome                     outcome        exposure
## 1      moyRYB     uxmElJ Circadian Preference Factor Drinks per Week
## 2      moyRYB     uxmElJ Circadian Preference Factor Drinks per Week
## 3      moyRYB     uxmElJ Circadian Preference Factor Drinks per Week
## 4      moyRYB     uxmElJ Circadian Preference Factor Drinks per Week
## 5      moyRYB     uxmElJ Circadian Preference Factor Drinks per Week
##                      method nsnp          b         se         pval
## 1        Maximum likelihood   39 0.08968357 0.01550589 7.301506e-09
## 2                  MR Egger   39 0.04235779 0.06219865 5.001065e-01
## 3             Weighted mode   39 0.04138732 0.02331407 8.387401e-02
## 4           Weighted median   39 0.02687537 0.02450505 2.727603e-01
## 5 Inverse variance weighted   39 0.08572313 0.03681727 1.989403e-02
psych_sleep_MR<- rbind(psych_sleep_MR, res)

mr_heterogeneity(dpw_chrono_dat)
##   id.exposure id.outcome                     outcome        exposure
## 1      moyRYB     uxmElJ Circadian Preference Factor Drinks per Week
## 2      moyRYB     uxmElJ Circadian Preference Factor Drinks per Week
##                      method        Q Q_df       Q_pval
## 1                  MR Egger 232.1803   37 4.084358e-30
## 2 Inverse variance weighted 236.8926   38 1.403787e-30
mr_pleiotropy_test(dpw_chrono_dat)
##   id.exposure id.outcome                     outcome        exposure
## 1      moyRYB     uxmElJ Circadian Preference Factor Drinks per Week
##   egger_intercept          se      pval
## 1     0.001013843 0.001169941 0.3917594
p1 <- mr_scatter_plot(res, dpw_chrono_dat)
p1
## $moyRYB.uxmElJ

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1      moyRYB     uxmElJ
res_single <- mr_singlesnp(dpw_chrono_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##           exposure                     outcome id.exposure id.outcome
## 1  Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 2  Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 3  Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 4  Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 5  Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 6  Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 7  Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 8  Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 9  Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 10 Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 11 Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 12 Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 13 Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 14 Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 15 Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
## 16 Drinks per Week Circadian Preference Factor      moyRYB     uxmElJ
##    samplesize                      SNP           b         se            p
## 1          NA               rs10085696  0.33373545 0.11999963 5.416881e-03
## 2          NA                rs1123285 -0.63900338 0.12928993 7.716408e-07
## 3          NA                 rs153106  0.31580338 0.11141808 4.591185e-03
## 4          NA                rs2049045  0.64868364 0.13973898 3.448702e-06
## 5          NA               rs28680958  0.42770439 0.13370323 1.379489e-03
## 6          NA               rs34121753  0.52569827 0.13647925 1.172265e-04
## 7          NA                rs3814877  0.56218195 0.13544724 3.316621e-05
## 8          NA                rs4233567 -0.42692065 0.12570543 6.832869e-04
## 9          NA                rs4916723  0.41268125 0.13579916 2.374368e-03
## 10         NA                 rs682011  0.41997089 0.13575463 1.977420e-03
## 11         NA                rs6951574  0.38377458 0.11860910 1.213766e-03
## 12         NA               rs75120545  0.80556298 0.14237734 1.531975e-08
## 13         NA               rs76640332  0.30446894 0.08581427 3.881699e-04
## 14         NA               rs79616692  0.67857671 0.12786839 1.115438e-07
## 15         NA                 rs838145  0.35302475 0.09580907 2.289996e-04
## 16         NA All - Maximum likelihood  0.08968357 0.01550589 7.301506e-09

insomnia

insom_out_dat <- read_outcome_data(
    snps = dpw_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/insomnia_withNeffR.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "se",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "p",
    samplesize_col = "Neff"
)
## No phenotype name specified, defaulting to 'outcome'.
insom_out_dat$outcome<- rep('Insomnia', 1, nrow(insom_out_dat))

dpw_insom_dat <- harmonise_data(dpw_exp_dat, insom_out_dat)
## Harmonising Drinks per Week (moyRYB) and Insomnia (ZzP0Zf)
res <- mr(dpw_insom_dat, method_list = methods)
## Analysing 'moyRYB' on 'ZzP0Zf'
res
##   id.exposure id.outcome  outcome        exposure                    method
## 1      moyRYB     ZzP0Zf Insomnia Drinks per Week        Maximum likelihood
## 2      moyRYB     ZzP0Zf Insomnia Drinks per Week                  MR Egger
## 3      moyRYB     ZzP0Zf Insomnia Drinks per Week             Weighted mode
## 4      moyRYB     ZzP0Zf Insomnia Drinks per Week           Weighted median
## 5      moyRYB     ZzP0Zf Insomnia Drinks per Week Inverse variance weighted
##   nsnp          b         se         pval
## 1   40 0.05222669 0.01401464 0.0001940941
## 2   40 0.05212663 0.04015406 0.2020544672
## 3   40 0.04476293 0.02006545 0.0315138187
## 4   40 0.04899379 0.02097207 0.0194837231
## 5   40 0.05127638 0.02325713 0.0274709748
psych_sleep_MR<- rbind(psych_sleep_MR, res)


mr_heterogeneity(dpw_insom_dat)
##   id.exposure id.outcome  outcome        exposure                    method
## 1      moyRYB     ZzP0Zf Insomnia Drinks per Week                  MR Egger
## 2      moyRYB     ZzP0Zf Insomnia Drinks per Week Inverse variance weighted
##         Q Q_df       Q_pval
## 1 112.705   38 2.504450e-09
## 2 112.707   39 4.393265e-09
mr_pleiotropy_test(dpw_insom_dat)
##   id.exposure id.outcome  outcome        exposure egger_intercept           se
## 1      moyRYB     ZzP0Zf Insomnia Drinks per Week    -1.94454e-05 0.0007436229
##        pval
## 1 0.9792749
p1 <- mr_scatter_plot(res, dpw_insom_dat)
p1
## $moyRYB.ZzP0Zf

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1      moyRYB     ZzP0Zf
res_single <- mr_singlesnp(dpw_insom_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##          exposure  outcome id.exposure id.outcome samplesize
## 1 Drinks per Week Insomnia      moyRYB     ZzP0Zf   259365.1
## 2 Drinks per Week Insomnia      moyRYB     ZzP0Zf   259365.1
## 3 Drinks per Week Insomnia      moyRYB     ZzP0Zf   259365.1
## 4 Drinks per Week Insomnia      moyRYB     ZzP0Zf   259365.1
## 5 Drinks per Week Insomnia      moyRYB     ZzP0Zf   259365.1
##                        SNP           b         se            p
## 1                rs2049045 -0.43800020 0.12895458 6.824201e-04
## 2               rs28929474  0.35793136 0.10277002 4.961278e-04
## 3                rs4752999  0.53035433 0.10154817 1.763451e-07
## 4                rs4916723  0.46304094 0.12549925 2.246146e-04
## 5 All - Maximum likelihood  0.05222669 0.01401464 1.940941e-04

bipolar

## exposure data
bip_full<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/BiP_withNeff.txt.gz", header=T, data.table = F)
bip<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/FUMA_out/BiP/GenomicRiskLoci.txt", header=T, data.table=F)

bip<- bip %>% select(rsID) %>%
  rename(SNP=rsID)

bip_full<- bip_full %>% 
  right_join(bip)
## Joining with `by = join_by(SNP)`
bip_exp_dat <- format_data(bip_full, type="exposure",
                              snp_col = "SNP", beta_col = "BETA", se_col = "SE", effect_allele_col = "A1", other_allele_col = "A2", pval_col = "PVAL", samplesize_col = "Neff")
## No phenotype name specified, defaulting to 'exposure'.
## Warning in format_data(bip_full, type = "exposure", snp_col = "SNP", beta_col = "BETA", : The following columns are not present but are helpful for harmonisation
## eaf
bip_exp_dat$exposure<- rep("Bipolar Disorder", 1, nrow(bip_exp_dat))

bip_exp_dat <- bip_exp_dat %>% as.data.table
bip_exp_dat <- bip_exp_dat[, F_STAT := qf(
  pval.exposure,
  df1=1,
  df2=samplesize.exposure,
lower.tail=F
)] # this piece calculates the F statistic of each SNP in the IV and adds a column called F_STAT
bip_exp_dat$F_STAT %>% mean # to calculate the overall F stat of your exposure
## [1] 37.91726
bip_exp_dat$F_STAT %>% min # you can use this to check whether any individual SNPs fall under F of 10 which would make them inva
## [1] 29.78602

napping

napping_out_dat <- read_outcome_data(
    snps = bip_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/napsR.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "BETA",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "P",
    samplesize_col = "N"
)
## No phenotype name specified, defaulting to 'outcome'.
napping_out_dat$outcome<- rep('No napping (napping GWAS)', 1, nrow(napping_out_dat))

bip_napping_dat <- harmonise_data(bip_exp_dat, napping_out_dat)
## Harmonising Bipolar Disorder (0aLoK9) and No napping (napping GWAS) (TEhbd1)
## Removing the following SNPs for being palindromic with intermediate allele frequencies:
## rs10455979, rs11062170, rs115694474, rs13417268, rs1998820, rs2011302, rs4331993, rs62011709, rs62489493
res <- mr(bip_napping_dat, method_list = methods)
## Analysing '0aLoK9' on 'TEhbd1'
res
##   id.exposure id.outcome                   outcome         exposure
## 1      0aLoK9     TEhbd1 No napping (napping GWAS) Bipolar Disorder
## 2      0aLoK9     TEhbd1 No napping (napping GWAS) Bipolar Disorder
## 3      0aLoK9     TEhbd1 No napping (napping GWAS) Bipolar Disorder
## 4      0aLoK9     TEhbd1 No napping (napping GWAS) Bipolar Disorder
## 5      0aLoK9     TEhbd1 No napping (napping GWAS) Bipolar Disorder
##                      method nsnp             b          se         pval
## 1        Maximum likelihood   50 -0.0142639783 0.003166870 6.664698e-06
## 2                  MR Egger   50 -0.0114258949 0.031246488 7.162181e-01
## 3             Weighted mode   50 -0.0004617008 0.012328126 9.702776e-01
## 4           Weighted median   50 -0.0030519864 0.005332907 5.671233e-01
## 5 Inverse variance weighted   50 -0.0126372237 0.006653232 5.751044e-02
psych_sleep_MR<- rbind(psych_sleep_MR, res)

mr_heterogeneity(bip_napping_dat)
##   id.exposure id.outcome                   outcome         exposure
## 1      0aLoK9     TEhbd1 No napping (napping GWAS) Bipolar Disorder
## 2      0aLoK9     TEhbd1 No napping (napping GWAS) Bipolar Disorder
##                      method        Q Q_df       Q_pval
## 1                  MR Egger 253.6939   48 9.120211e-30
## 2 Inverse variance weighted 253.7022   49 2.110533e-29
mr_pleiotropy_test(bip_napping_dat)
##   id.exposure id.outcome                   outcome         exposure
## 1      0aLoK9     TEhbd1 No napping (napping GWAS) Bipolar Disorder
##   egger_intercept          se      pval
## 1   -8.296928e-05 0.002090098 0.9684999
p1 <- mr_scatter_plot(res, bip_napping_dat)
p1
## $`0aLoK9.TEhbd1`

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1      0aLoK9     TEhbd1
res_single <- mr_singlesnp(bip_napping_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##            exposure                   outcome id.exposure id.outcome samplesize
## 1  Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 2  Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 3  Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 4  Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 5  Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 6  Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 7  Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 8  Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 9  Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 10 Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 11 Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
## 12 Bipolar Disorder No napping (napping GWAS)      0aLoK9     TEhbd1     452633
##                         SNP           b         se            p
## 1                rs10255167 -0.09186815 0.02197639 2.911114e-05
## 2                rs10761661 -0.13108142 0.02390380 4.165043e-08
## 3                rs12289486  0.08190674 0.02337102 4.572448e-04
## 4                rs12668848  0.07798441 0.02145833 2.788178e-04
## 5                rs12672003 -0.05732342 0.02066772 5.544438e-03
## 6                rs13195402 -0.06593031 0.01420933 3.485275e-06
## 7                 rs1487445  0.05716574 0.01639233 4.878393e-04
## 8                  rs174592 -0.11411084 0.01744741 6.140571e-11
## 9                  rs228768 -0.07187589 0.02062017 4.908376e-04
## 10                rs2336147 -0.05856798 0.01791257 1.076776e-03
## 11                 rs748455 -0.06825091 0.01976748 5.550330e-04
## 12 All - Maximum likelihood -0.01426398 0.00316687 6.664698e-06

schizophrenia

## exposure data
scz_full<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/scz_withNeff.txt.gz", header=T, data.table = F)
scz<- fread("/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/FUMA_out/scz/GenomicRiskLoci.txt", header=T, data.table=F)

scz<- scz %>% select(rsID) %>%
  rename(SNP=rsID)

scz_full<- scz_full %>% 
  right_join(scz)
## Joining with `by = join_by(SNP)`
scz_exp_dat <- format_data(scz_full, type="exposure",
                              snp_col = "SNP", beta_col = "OR", se_col = "SE", effect_allele_col = "A1", other_allele_col = "A2", eaf_col = "MAF", pval_col = "P", samplesize_col = "Neff")
## No phenotype name specified, defaulting to 'exposure'.
scz_exp_dat$exposure<- rep("Schizophrenia", 1, nrow(scz_exp_dat))


scz_exp_dat <- scz_exp_dat %>% as.data.table
scz_exp_dat <- scz_exp_dat[, F_STAT := qf(
  pval.exposure,
  df1=1,
  df2=samplesize.exposure,
lower.tail=F
)] # this piece calculates the F statistic of each SNP in the IV and adds a column called F_STAT
scz_exp_dat$F_STAT %>% mean # to calculate the overall F stat of your exposure
## [1] 41.00856
scz_exp_dat$F_STAT %>% min # you
## [1] 29.75531
nrow(scz_exp_dat)
## [1] 96

napping

napping_out_dat <- read_outcome_data(
    snps = scz_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/napsR.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "BETA",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "P",
    samplesize_col = "N"
)
## No phenotype name specified, defaulting to 'outcome'.
napping_out_dat$outcome<- rep('No napping (napping GWAS)', 1, nrow(napping_out_dat))

scz_napping_dat <- harmonise_data(scz_exp_dat, napping_out_dat)
## Harmonising Schizophrenia (cwkKk2) and No napping (napping GWAS) (TtzCNm)
res <- mr(scz_napping_dat, method_list = methods)
## Analysing 'cwkKk2' on 'TtzCNm'
res
##   id.exposure id.outcome                   outcome      exposure
## 1      cwkKk2     TtzCNm No napping (napping GWAS) Schizophrenia
## 2      cwkKk2     TtzCNm No napping (napping GWAS) Schizophrenia
## 3      cwkKk2     TtzCNm No napping (napping GWAS) Schizophrenia
## 4      cwkKk2     TtzCNm No napping (napping GWAS) Schizophrenia
## 5      cwkKk2     TtzCNm No napping (napping GWAS) Schizophrenia
##                      method nsnp             b           se         pval
## 1        Maximum likelihood   95 -0.0006936287 0.0001404643 7.888208e-07
## 2                  MR Egger   95 -0.0064158261 0.0038756909 1.012138e-01
## 3             Weighted mode   95 -0.0014020546 0.0006607064 3.646269e-02
## 4           Weighted median   95 -0.0010361511 0.0002397823 1.551721e-05
## 5 Inverse variance weighted   95 -0.0006937928 0.0003009317 2.113955e-02
psych_sleep_MR<- rbind(psych_sleep_MR, res)

mr_heterogeneity(scz_napping_dat)
##   id.exposure id.outcome                   outcome      exposure
## 1      cwkKk2     TtzCNm No napping (napping GWAS) Schizophrenia
## 2      cwkKk2     TtzCNm No napping (napping GWAS) Schizophrenia
##                      method        Q Q_df       Q_pval
## 1                  MR Egger 421.7722   93 2.261166e-43
## 2 Inverse variance weighted 431.7169   94 9.737992e-45
mr_pleiotropy_test(scz_napping_dat)
##   id.exposure id.outcome                   outcome      exposure
## 1      cwkKk2     TtzCNm No napping (napping GWAS) Schizophrenia
##   egger_intercept          se      pval
## 1     0.005775813 0.003900455 0.1420386
p1 <- mr_scatter_plot(res, scz_napping_dat)
p1
## $cwkKk2.TtzCNm

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1      cwkKk2     TtzCNm
res_single <- mr_singlesnp(scz_napping_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##         exposure                   outcome id.exposure id.outcome samplesize
## 1  Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 2  Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 3  Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 4  Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 5  Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 6  Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 7  Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 8  Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 9  Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 10 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 11 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 12 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 13 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 14 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 15 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 16 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 17 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 18 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 19 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 20 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
## 21 Schizophrenia No napping (napping GWAS)      cwkKk2     TtzCNm     452633
##                         SNP             b           se            p
## 1                rs10520163 -0.0040727801 0.0011463977 3.813331e-04
## 2                rs11210892 -0.0050093369 0.0013857385 3.004378e-04
## 3                rs11682175  0.0074942110 0.0013102995 1.068606e-08
## 4                rs11693094  0.0044155942 0.0013091508 7.438931e-04
## 5               rs140505938  0.0084152760 0.0017776938 2.203372e-06
## 6                 rs2514218  0.0059273276 0.0013680697 1.473447e-05
## 7                 rs2851447 -0.0052632692 0.0015468993 6.678116e-04
## 8                  rs301797 -0.0060258365 0.0011945768 4.551192e-07
## 9                rs35518360 -0.0087049468 0.0026173948 8.816334e-04
## 10                rs3798869 -0.0087312137 0.0013048066 2.207786e-11
## 11                rs4523957 -0.0040181378 0.0011854972 7.004282e-04
## 12               rs55661361  0.0050578959 0.0014047524 3.175319e-04
## 13               rs58120505  0.0040979660 0.0011284166 2.816615e-04
## 14                rs6466055 -0.0049165333 0.0011881441 3.503554e-05
## 15                rs7801375 -0.0059121964 0.0017995002 1.018139e-03
## 16                rs7951870 -0.0052977568 0.0017696813 2.756874e-03
## 17                rs9398171 -0.0039136621 0.0014271068 6.099664e-03
## 18                rs9461856 -0.0043954277 0.0011276335 9.702151e-05
## 19                 rs950169  0.0042597255 0.0014859935 4.149279e-03
## 20 All - Maximum likelihood -0.0006936287 0.0001404643 7.888208e-07
## 21    All - Weighted median -0.0010361511 0.0002515071 3.792537e-05

actigraphy duration

actidur_out_dat <- read_outcome_data(
    snps = scz_exp_dat$SNP,
    filename = "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR/sumstats/actisleepdur.txt.gz",
    sep = "\t",
    snp_col = "SNP",
    beta_col = "beta",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
  eaf_col = "MAF",
    pval_col = "p",
    samplesize_col = "N"
)
## No phenotype name specified, defaulting to 'outcome'.
actidur_out_dat$outcome<- rep('Actigraphy Sleep Duration', 1, nrow(actidur_out_dat))

scz_actidur_dat <- harmonise_data(scz_exp_dat, actidur_out_dat)
## Harmonising Schizophrenia (cwkKk2) and Actigraphy Sleep Duration (cnDjRg)
res <- mr(scz_actidur_dat, method_list = methods)
## Analysing 'cwkKk2' on 'cnDjRg'
res
##   id.exposure id.outcome                   outcome      exposure
## 1      cwkKk2     cnDjRg Actigraphy Sleep Duration Schizophrenia
## 2      cwkKk2     cnDjRg Actigraphy Sleep Duration Schizophrenia
## 3      cwkKk2     cnDjRg Actigraphy Sleep Duration Schizophrenia
## 4      cwkKk2     cnDjRg Actigraphy Sleep Duration Schizophrenia
## 5      cwkKk2     cnDjRg Actigraphy Sleep Duration Schizophrenia
##                      method nsnp             b           se       pval
## 1        Maximum likelihood   92 -0.0009283762 0.0005550074 0.09438094
## 2                  MR Egger   92  0.0008733544 0.0105679436 0.93431985
## 3             Weighted mode   92  0.0007831962 0.0026446011 0.76779101
## 4           Weighted median   92 -0.0002301451 0.0008567453 0.78821654
## 5 Inverse variance weighted   92 -0.0009286611 0.0008124957 0.25304957
psych_sleep_MR<- rbind(psych_sleep_MR, res)


mr_heterogeneity(scz_actidur_dat)
##   id.exposure id.outcome                   outcome      exposure
## 1      cwkKk2     cnDjRg Actigraphy Sleep Duration Schizophrenia
## 2      cwkKk2     cnDjRg Actigraphy Sleep Duration Schizophrenia
##                      method        Q Q_df       Q_pval
## 1                  MR Egger 195.0214   90 9.993724e-10
## 2 Inverse variance weighted 195.0848   91 1.462696e-09
mr_pleiotropy_test(scz_actidur_dat)
##   id.exposure id.outcome                   outcome      exposure
## 1      cwkKk2     cnDjRg Actigraphy Sleep Duration Schizophrenia
##   egger_intercept         se     pval
## 1    -0.001820945 0.01064701 0.864585
p1 <- mr_scatter_plot(res, scz_actidur_dat)
p1
## $cwkKk2.cnDjRg

## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   id.exposure id.outcome
## 1      cwkKk2     cnDjRg
res_single <- mr_singlesnp(scz_actidur_dat,all_method = methods)
p2 <- mr_forest_plot(res_single)
p2[[1]]
## Warning: Removed 1 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

res_single %>% filter(p<ivw_p)
##        exposure                   outcome id.exposure id.outcome samplesize
## 1 Schizophrenia Actigraphy Sleep Duration      cwkKk2     cnDjRg      85449
## 2 Schizophrenia Actigraphy Sleep Duration      cwkKk2     cnDjRg      85449
## 3 Schizophrenia Actigraphy Sleep Duration      cwkKk2     cnDjRg      85449
##         SNP           b          se            p
## 1 rs1702294 -0.04390911 0.006771662 8.917762e-11
## 2 rs2851447 -0.01772710 0.005977766 3.021851e-03
## 3 rs6434928  0.01659276 0.005409702 2.160582e-03

save output of res

fwrite(psych_sleep_MR, "/Users/claire/Desktop/dissertation/cotwin_mendelian/sleepMR//psych_exposure_MR.csv", sep='\t')