Two Sample MR

一、暴露数据预处理 假设研究BMI和阿尔兹海默症的关系，下载BMI的SNP数据为SNP_gwas_mc_nogc_tbl.txt。 首先将数据放于TwoSampleMR包文件夹内（可再新建文件夹（例如TwoSampleMR\xxx\SNP_gwas_mc_nogc_tbl.txt））Finn Gen，工作环境和数据位置一致。

library("Two Sample MR") ##调用R包
a <- read.table("D:/MRlearning/xxxx/xx.txt", header = T) ##读取数据，header看数据集中是否带有表头，如有即为T，否则即为F
View(a)

## 二、SNP筛选条件设置 ● 相关性：P<5e-08 ● 独立性：Clumping：LD r2<0.001, distance 1Mb ● 统计强度设置：R^2计算公式和F值计算公示，F值一般大于10 ### 1.暴露SNP的筛选

b <- subset(a,p<5e-08)                        #进行独立性筛选，即筛选a中"p"列小于5e-08的子集，注意列名一致性；
# 注：如果p值过于严格导致结果过少，可以放宽至5e-06至5e-14
View(b)
write.csv(b, file = "exposure_bmi.csv")       #输出b为一个csv文件，便于后续处理
## 从UKBB和Finn Gene来的数据这步已经完成了，从下一行代码开始看就行##
##对输出的csv修改列名，a1为effect_allele, a2为other_allele, FreqA1为beta（检验效能）, se（标准误），便于后续操作
bmi <- system.file("xxx/exposure_bmi.csv", package = "TwoSampleMR")#调用TwoSampleMR，读取exposre_bmi.csv的文件数据
bmi <- read.csv("xxxx.csv")
bmi_exp_data <- read_exposure_data(                    #数据处理，转换格式
  filename = "",        #名称同第4行
  sep = ",",            #csv文件以逗号分隔
  snp_col = "",            #标明csv文件中的列名，下同
  beta_col = "",
  se_col = "",
  effect_allele_col = "",
  other_allele_col = "",
  eaf_col = "",            #效应等位基因的频率
  pval_col = ""，
  phenotyle_col = ""            #本行可忽略
)
View(bmi_exp_data)
bmi_exp_data_clumped <- clump_data(#对上述整理出的数据进行独立性设置（LD r2<0.001, distance 1Mb）
  bmi_exp_data,            #需要进行clump的数据集
  clump_kb = 10000,        #distance设置，1Mb即为10000kb
  clump_r2 = 0.001,        #LD r2设置
  clump_p1 = 1,            #不知道什么意思
  clump_p2 = 1,            #不知道什么意思
  pop = "EUR"                            #样本人群，欧洲为EUR，东亚为EAR
)
View(bmi_exp_data_clumped)

> 2024-07-20 clump_data失效，参考clump失败的处理 ### 2.F统计量的计算

# 来源：https://github.com/HaobinZhou/Get_MR.git
# Get_MR2.0.r
get_f <- function(dat,F_value=10){
  log <- is.na(dat$eaf.exposure)
  log <- unique(log)
  if(length(log)==1)
  {if(log==TRUE){print("数据不包含eaf，无法计算F统计量")
    return(dat)}}
  if(is.null(dat$beta.exposure[1])==T || is.na(dat$beta.exposure[1])==T){print("数据不包含beta，无法计算F统计量")
    return(dat)}
  if(is.null(dat$se.exposure[1])==T || is.na(dat$se.exposure[1])==T){print("数据不包含se，无法计算F统计量")
    return(dat)}
  if(is.null(dat$samplesize.exposure[1])==T || is.na(dat$samplesize.exposure[1])==T){print("数据不包含samplesize(样本量)，无法计算F统计量")
    return(dat)}
  }
bmi_exp_data <- get_f(bmi_exp_data, F_value = 10)

## 三、导入结局数据（AD） 假设获取阿尔兹海默的结局数据IGAP_stage_1.txt ### 1.读取Outcome数据

a <- read.table("IGAP_stage_1.txt", header = T)

### 2.将暴露和结局进行Merge

data_merge <- merge(                                    #merge的作用就是将两个数据集按照SNP编号"rsxxxxxx"，将相同SNP编号整理在一起，相当于取交集，第一个在前，第二个在后
  bmi_exp_data_clumped, a,                             #参与merge的两个数据集
  by.x = "SNP", by.y = "Makername"            #by.x即第一个数据集中SNP编号的列名，by.y即第二个数据集中SNP编号的列名
)
########需要注意的是，为了满足MR的独立性假设，暴露（BMI）的工具变量（SNP）不应与结局（AD）的工具变量具有直接相关性，即上述工具变量merge后在AD数据集中的p应当≥5e-08########
View(data_merge)
write.csv(data_merge, "outcome.csv")
##在outcome.csv中可以将暴露组的相关列删掉，并对结局列重命名，只留下SNP编号、结局的effect_allele, other_allele_ beta, se, p

### 3.读取outcome.csv数据

outcome_data <- read_outcome_data(
  snps = bmi_exp_data_clumped$SNP,            #$符号意在出现并选择该数据集的列名
  filename = "outcome.csv"                            #文件名同上
  sep = ",",                                                        #对每一列的列名进行表明，同read_exposure_data
  snp_col = "",
  beta_col = "", #beta,b,eaf,FreqA1
  se_col = "",
  effect_allele_col = "",
  other_allele_col = "",
  eaf_col = "",
  pval_col = ""
)

### 4.几个问题 **#### Q1.暴露数据和结局数据的SNP对应的效应等位基因不一致怎么办？** eg：