library(data.table)
library(TwoSampleMR)

#forward
exposure_dat <- extract_instruments(outcomes='ebi-a-GCST004420',p1 = 5e-06, 
                                    clump = TRUE,r2 = 0.001,
                                    kb = 10000,
                                    access_token= NULL)##opengwasID can be found in supplementary table S1. input the OPENGWASID of  each cytokine as outcome
#this is instrumental variable extraction and linkage disequilibrium removal
outcome_dat <- extract_outcome_data(snps = clump_dat$SNP, outcomes = 'ieu-b-4980')
#Extract the information of the SNP obtained in the previous step in the sepsis outcome file
dat <- harmonise_data(clump_dat, outcome_dat, action = 2)
#Perform allele alignment
dat <- subset(dat,dat$pval.outcome < 1e-5)
#Remove instrumental variables related to outcomes
mr_results <- mr(dat)
#priamry MR reasult 
mr_results
#display the result

res_single <- mr_singlesnp(dat)
mr_scatter_plot(mr_results,dat)
#Scatter plot
leaveoneout <- mr_leaveoneout(dat)
mr_leaveoneout_plot(leaveoneout)
#leave-one-out

mr_presso(BetaExposure = 'beta.exposure',BetaOutcome = 'beta.outcome',SdExposure = 'se.exposure',SdOutcome = 'se.outcome', data = dat,OUTLIERtest = TRUE,DISTORTIONtest = TRUE)
#MRPRESSO for horizontal pleiotropy

het <- mr_heterogeneity(dat)
egger <- mr_pleiotropy_test(dat)
# Egger for Heterogeneity and horizontal pleiotropy


##reverse
exposure_dat <- extract_instruments(outcomes='ieu-b-4980',p1 = 5e-06, 
                                    clump = TRUE,r2 = 0.001,
                                    kb = 10000,
                                    access_token= NULL)opengwasID can be found in supplementary table S1. input the OPENGWASID of  each cytokine as outcome
#this is instrumental variable extraction and linkage disequilibrium removal

outcome_dat <- extract_outcome_data(snps = clump_dat$SNP, outcomes = 'ebi-a-GCST004420')
##Extract the information of the SNP obtained in the previous step in the sepsis outcome file
dat <- harmonise_data(clump_dat, outcome_dat, action = 2)
#Perform allele alignment
dat <- subset(dat,dat$pval.outcome < 1e-5)
#Remove instrumental variables related to outcomes
mr_results <- mr(dat)
#priamry MR reasult 
mr_results
#display the result

res_single <- mr_singlesnp(dat)
mr_scatter_plot(mr_results,dat)
#Scatter plot
leaveoneout <- mr_leaveoneout(dat)
mr_leaveoneout_plot(leaveoneout)
#leave-one-out

mr_presso(BetaExposure = 'beta.exposure',BetaOutcome = 'beta.outcome',SdExposure = 'se.exposure',SdOutcome = 'se.outcome', data = dat,OUTLIERtest = TRUE,DISTORTIONtest = TRUE)
#MRPRESSO for horizontal pleiotropy

het <- mr_heterogeneity(dat)
egger <- mr_pleiotropy_test(dat)
# Egger for Heterogeneity and horizontal pleiotropy