forked from XiaojuanJiang/LBP-analysis
-
Notifications
You must be signed in to change notification settings - Fork 0
/
code for the wole study.gitignore
190 lines (128 loc) · 6.38 KB
/
code for the wole study.gitignore
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
####Primary analysis
####library packages need for analysis
install.packages("remots")
remotes::install_github("MRCIEU/TwoSampleMR")
install.packages("RadialMR")
library(TwoSampleMR)
library(RadialMR)
####read the exposure and check
exposure_online<-extract_instruments(outcomes= "ebi-a-GCST90000616")
View(exposure_online)
####set working dictionary
setwd("your work place")
####outcome extract,the file from FinnGen should be download from FinnGen website, we also provided the
#####information extract from corresponding files,name as "finngen_R7_M13_LOWBACKPAIN1.csv"
outcome_dat <- read_outcome_data( snps = exposure_online$SNP,filename = "finngen_R7_M13_LOWBACKPAIN1.csv",sep = ",",snp_col = "SNP",beta_col = "beta",
se_col = "se",effect_allele_col = "effect_allele",other_allele_col = "other_allele",pval_col = "pval")
####协同文件输出
dat<-harmonise_data(exposure_dat = exposure_online,outcome_dat = outcome_dat)
write.csv(dat,file = "datprimaryanlysis.csv")
####MR分析
mr(dat)
generate_odds_ratios(mr_res=mr(dat))
#####mr plots,export the plots as pdf and then organize them in adobe illustration
res_single <- mr_singlesnp(dat)
res_single <- mr_singlesnp(dat, single_method="mr_meta_fixed")
res_loo <- mr_leaveoneout(dat)
res <- mr(dat, method_list=c("mr_egger_regression", "mr_ivw","mr_weighted_median"))
####scatter plot
p1 <- mr_scatter_plot(res, dat)
p1[[1]]
res_single <- mr_singlesnp(dat)
####forest plot
p2 <- mr_forest_plot(res_single)
p2[[1]]
####leave one out plot
res_loo <- mr_leaveoneout(dat)
p3 <- mr_leaveoneout_plot(res_loo)
p3[[1]]
####funnel plot
res_single <- mr_singlesnp(dat)
p4 <- mr_funnel_plot(res_single)
p4[[1]]
####灵敏度检验
mr_heterogeneity(dat)
mr_funnel_plot(singlesnp_results = mr_singlesnp(dat))
mr_pleiotropy_test(dat)
mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))
####MR-PRESSO分析
run_mr_presso(dat,NbDistribution = 3000)
####identify the outliers and then discard the outliers in the "datprimaryanlysis.csv", save and then performed other analyses with reloaded "datprimaryanlysis.csv"
####Radial MR analysis
dat <- read.csv("datprimaryanlysis.csv",header = T)
radial_data <- format_radial(dat$beta.exposure,dat$beta.outcome,dat$se.exposure,
dat$se.outcome,dat$SNP)
ivw.model <- ivw_radial(r_input = radial_data, alpha = 0.05,
weights = 1, tol = 0.0001, summary = TRUE)
ivw.model$outliers
####Radial MR plot,export the plots as pdf and then organize them in adobe illustration
plot_radial(ivw.model)
######Multiple replicated IVW analysis after removing outliers identified with radial MR
dat <- read.csv("datprimaryanlysis.csv",header = T)
mr_ivw_mre(dat$beta.exposure,dat$beta.outcome,dat$se.exposure,dat$se.outcome,parameters = default_parameters())
#####Replicate analysis
####exposure were retrieve from previous literature, and was uploaded as "exposure_circulating VD.csv"
exposure<- read_exposure_data( filename = "exposure_circulating VD.csv", sep = ",", snp_col = "SNP",beta_col = "beta",se_col = "se",effect_allele_col = "effect_allele",other_allele_col = "other_allele", eaf_col = "eaf", pval_col = "pval",
phenotype_col = "exposure")
View(exposure)
###结局读取与协同
outcome_dat <- extract_outcome_data(snps = exposure$SNP,outcomes="ukb-b-9838")
View(outcome_dat)
write.csv(outcome_dat,file = "outcome_dat_replicate.csv")
####协同文件输出
dat<-harmonise_data(exposure_dat = exposure,outcome_dat = outcome_dat)
write.csv(dat,file = "datreplicate.csv")
####MR分析
mr(dat)
generate_odds_ratios(mr_res=mr(dat))
mr_method_list()
mr(dat,method_list=c("mr_ivw","mr_weighted_median","mr_egger_regression"))
mr_scatter_plot(mr_results = mr(dat,method_list = c("mr_ivw","mr_weighted_median","mr_egger_regression")),dat)
####灵敏度检验
mr_heterogeneity(dat)
mr_funnel_plot(singlesnp_results = mr_singlesnp(dat))
mr_pleiotropy_test(dat)
mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))
##########positive control analysis
####exposure----serum 25(OH)D
exposure_online<-extract_instruments(outcomes= "ebi-a-GCST90000616")
exposure<- read_exposure_data( filename = "exposure_circulating VD.csv", sep = ",", snp_col = "SNP",beta_col = "beta",se_col = "se",effect_allele_col = "effect_allele",other_allele_col = "other_allele", eaf_col = "eaf", pval_col = "pval",
phenotype_col = "exposure")
####outcome were obtained from the FinnGen,we also provide information as "outcome_exposure_online_VDD.csv" and "outcome_exposure_VDD.csv"
a<-read.table("finngen_R7_E4_VIT_D_DEF.txt",header = T,sep="\t")
View(a)
c<-merge(exposure_online,a,by.x ="SNP",by.y = "rsids")
View(c)
write.csv(c,file = "outcome_exposure_online_VDD.csv")
d<-merge(exposure,a,by.x ="SNP",by.y = "rsids")
View(d)
write.csv(d,file = "outcome_exposure_VDD.csv")
###outcome
outcome_dat1 <- read_outcome_data( snps = exposure_online$SNP,filename = "outcome_exposure_online_VDD.csv",sep = ",",snp_col = "SNP",beta_col = "beta",se_col = "sebeta",effect_allele_col = "alt",other_allele_col = "ref",pval_col = "pval")
outcome_dat2 <- read_outcome_data( snps = exposure$SNP,filename = "outcome_exposure_VDD.csv",sep = ",",snp_col = "SNP",beta_col = "beta",se_col = "sebeta",effect_allele_col = "alt",other_allele_col = "ref",pval_col = "pval")
####协同文件输出
dat1<-harmonise_data(exposure_dat = exposure_online,outcome_dat = outcome_dat1)
dat2<-harmonise_data(exposure_dat = exposure,outcome_dat = outcome_dat2)
write.csv(dat,file = "dat.csv")
####MR分析
mr(dat1)
generate_odds_ratios(mr_res=mr(dat1))
mr(dat2)
generate_odds_ratios(mr_res=mr(dat2))
#########potential risk factors analysis
####for alcohol intake frequency
exposure_online<-extract_instruments(outcomes= "ebi-a-GCST90000616")
outcome_dat <- extract_outcome_data(snps = exposure$SNP,outcomes="ukb-b-5779")
dat<-harmonise_data(exposure_dat = exposure,outcome_dat = outcome_dat)
mr(dat)
####for BMI
exposure_online<-extract_instruments(outcomes= "ebi-a-GCST90000616")
outcome_dat <- extract_outcome_data(snps = exposure$SNP,outcomes="ukb-b-19953")
dat<-harmonise_data(exposure_dat = exposure,outcome_dat = outcome_dat)
mr(dat)
####for obesity
exposure_online<-extract_instruments(outcomes= "ebi-a-GCST90000616")
outcome_dat <- extract_outcome_data(snps = exposure$SNP,outcomes="ukb-b-15541")
dat<-harmonise_data(exposure_dat = exposure,outcome_dat = outcome_dat)
mr(dat)
generate_odds_ratios(mr_res=mr(dat))