添加了去除混淆SNP和筛选F值操作
另外输出的文件mr_KEEP=TRUE
library('TwoSampleMR')
id_list = read.table('C:/Users/DELL/Desktop/cs/邱师兄越高/expo_list.txt',header = T,sep = '\t')
remove_snp = read.table('./remove_snp.txt', header = T, sep = '\t')
# 设置最大重试次数和每次重试等待的时间
max_retries <- 55
wait_time <- 5
# 定义一个函数用来连接服务器
connect_to_server <- function() {
# 尝试连接服务器
connection <- try({extract_instruments(outcomes = i)}, silent = TRUE)
# 如果连接失败,则重试
retries <- 1
while (inherits(connection, "try-error") && retries <= max_retries) {
# 等待一段时间后重新连接
Sys.sleep(wait_time)
connection <- try({extract_instruments(outcomes = i)}, silent = TRUE)
retries <- retries + 1
}
# 返回连接结果
return(connection)
}
# 循环遍历所有的 exposure_id
for (i in id_list$exposure_id) {
# 连接服务器
exposure <- connect_to_server()
try(dev.off(), silent = TRUE)
# 检查连接是否成功
if (inherits(exposure, "try-error")) {
# 如果连接失败,则输出错误信息并继续下一个循环
cat("Unable to connect to server for exposure", i, "\n")
try(dev.off(), silent = TRUE)
next
}
print(exposure$exposure[1])
# 提取 outcome 数据
exposure$RR <- 2 * exposure$eaf.exposure * (1 - exposure$eaf.exposure) * ((exposure$beta.exposure)^2)
exposure$F <- (exposure$beta.exposure)^2 / (exposure$se.exposure)^2
exposure <- exposure[exposure$F > 10, ]
print((table(exposure$id.exposure)))
print(format(sum(exposure$RR), scientific = TRUE, digits = 2))
outcome <- extract_outcome_data(snps = exposure$SNP, outcomes = 'ukb-b-19813')
# 检查 outcome 数据是否有效
if (is.null(outcome)) {
# 如果 outcome 数据无效,则输出错误信息并继续下一个循环
cat("Unable to extract outcome data for exposure", i, "\n")
try(dev.off(), silent = TRUE)
next
}
# 继续执行后续分析
dat <- harmonise_data(exposure_dat = exposure, outcome_dat = outcome)
dat$F <- (dat$beta.exposure)^2 / (dat$se.exposure)^2
dat <- dat[dat$mr_keep==TRUE,]
dat <- dat[!(dat$SNP %in% remove_snp$SNP),]
print(table(dat$mr_keep))
# filename <- paste0("./SNP/dat", i, '.txt')
# write.table(dat,file = filename, quote = F, row.names = F, sep = '\t')
# a <- generate_odds_ratios(mr_res = mr(dat))
# b <- mr_heterogeneity(dat)
# c <- mr_pleiotropy_test(dat)
# filenamea <- paste0('./mr/mr_', i, '.txt')
# filenameb <- paste0('./heter/heter_', i, '.txt')
# filenamec <- paste0('./pleio/pleio_', i, '.txt')
# write.table(a, file = filenamea, , quote = F, sep = '\t')
# write.table(b, file = filenameb, quote = F, sep = '\t')
# write.table(c, file = filenamec, quote = F, sep = '\t')
#
# filename <- paste0("./fig/outfig_", i, '.pdf')
# try(dev.off(), silent = TRUE)
# pdf(file = filename, onefile = TRUE)
# a <- mr_scatter_plot(mr_results = mr(dat, method_list = c("mr_ivw", "mr_egger_regression", "mr_weighted_median")), dat)
# b <- mr_funnel_plot(singlesnp_results = mr_singlesnp(dat))
# c <- mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))
# print(a)
# print(b)
# print(c)
# dev.off()
try(dev.off(), silent = TRUE)
Sys.sleep(wait_time)
print('over')
}