2024年6月25日09:43:36
1.Questions
多个整合
2. 资源
2.1 IBD
(1)2017-NG-IBD GWAS 共计241个SNP
2017-ng-GWAS-IBD_sup.xlsx
(2)LD score regression. We performed LD score regression using LDSC v1.0.0
and European LD scores from the 1000 Genomes Project (downloaded from
https://data.broadinstitute.org/alkesgroup/LDSCORE/eur_w_ld_chr.tar.bz2)
on our filtered meta-analysis summary statistics for all sites with INFO >0.95.
This INFO threshold was chosen to avoid confounding due to poor imputation,
as recommended by the authors45.
2.2 数据库
2.3 教程
1.孟德尔随机化
概述
1. 孟德尔随机化正是利用遗传数据评估各种危险因素间因果关系的方法
2.`混杂因素`和反向因果
(1) 孟德尔随机化之研究背景
(2) 孟德尔随机化之遗传学概述
(3)https://zhuanlan.zhihu.com/p/240828209
(4) TwoSampleMR包实战教程之读取暴露文件; Major changes to the IEU GWAS resources for 2020 官网
library(TwoSampleMR) #加载R包
bmi <-extract_instruments(outcomes='ieu-a-2',access_token = NULL) #获取暴露数据
head(bmi) #查看暴露数据
# 读取本地文件
exp_dat <- read_exposure_data(
filename = 'MICAD_gwas.txt',
clump = FALSE,#是否去除LD
sep= "\t",
snp_col = "SNPID", #snp_col用于指定SNP所在的列名
beta_col = "log_OR", #指定beta值所在的列名
se_col = "se",
effect_allele_col ="effect_allele",# 用于指定效应等位基因列名
other_allele_col = "other_allele",
eaf_col = "effect_allele_freq",
pval_col = "Pvalue"
)
head(exp_dat)
#转化为TwoSampleMR的格式
micad <-read.table('MICAD_gwas.txt',header=T)
exp_dat <- format_data(
micad,
type='exposure',
snp_col = "SNPID",
beta_col = "log_OR",
se_col = "se",
effect_allele_col ="effect_allele",
other_allele_col = "other_allele",
eaf_col = "effect_allele_freq",
pval_col = "Pvalue"
)
(5) TwoSampleMR包实战教程之去除连锁不平衡(LD)
library(TwoSampleMR)
bmi <-extract_instruments(outcomes='ieu-a-2',access_token = NULL) #用默认参数获取IV
dim(bmi) #查看IV个数
# [1] 79 15 结果显示提取79个IV
bmi<-extract_instruments(outcomes='ieu-a-2',clump=FALSE, access_token = NULL)
dim(bmi)
#[1] 2041 15 结果显示共提取2041个SNP
exp_dat <-clump_data(bmi,clump_r2=0.01,clump_kb=5000)
dim(exp_dat)
# [1] 90 15 结果显示clump之后还剩90个IV
# 方法二
bmi <-extract_instruments(outcomes='ieu-a-2', clump=TRUE, r2=0.01, kb=5000,access_token = NULL)
(6) TwoSampleMR实战教程之提取IV在结局中的信息
# 1. 利用TwoSampleMR获取MR base提供的结局信息
library(TwoSampleMR)
bmi_exp <- extract_instruments(
outcomes='ieu-a-835',
clump=TRUE, r2=0.01,
kb=5000,access_token = NULL
)
dim(bmi_exp)
# [1] 80 15
t2d_out <- extract_outcome_data(
snps=bmi_exp$SNP,
outcomes='ieu-a-26',
proxies = FALSE,
maf_threshold = 0.01,
access_token = NULL
)
dim(t2d_out)
# [1] 80 16
# 2. 从自己的GWAS结果中提取IV在结局中的信息
#install.packages('data.table') 安装data.table包
library(data.table) # 加载R包
t2d <-fread('DIAGRAMv3.2012DEC17.txt',header=T) # 使用fread()快速读取大文件
head(t2d) # 查看数据
t2d$phenotype <- 'Type 2 diabetes' # 添加phenotype列
t2d$beta <- log(t2d$OR) # 转换OR值为beta
t2d$se <-abs(t2d$beta/qnorm(t2d$P_VALUE/2,lower.tail=F)) # 计算se
head(t2d) # 查看数据
t2d_out <- format_data(
dat=t2d,
type = "outcome",
snps = bmi_exp$SNP,
header = TRUE,
phenotype_col = "phenotype",
snp_col = "SNP",
beta_col = "beta",
se_col = "se",
effect_allele_col = "RISK_ALLELE",
other_allele_col = "OTHER_ALLELE",
pval_col = "P_VALUE",
ncase_col = "N_CASES",
ncontrol_col = "N_CONTROLS",
chr_col = "CHROMOSOME",
pos_col = "POSITION"
)
head(t2d_out)
(7)TwoSampleMR实战教程之计算并解读MR结果 & TwoSampleMR包实战教程之敏感性分析& )TwoSampleMR包实战教程之MR结果可视化
#1. 准备好了暴露和结局的信息
library(TwoSampleMR)
bmi_exp <- extract_instruments(
outcomes='ieu-a-835',
clump=TRUE, r2=0.01,
kb=5000,access_token= NULL
)
dim(bmi_exp)
# [1] 80 15
t2d_out <- extract_outcome_data(
snps=bmi_exp$SNP,
outcomes='ieu-a-26',
proxies = FALSE,
maf_threshold = 0.01,
access_token = NULL
)
dim(t2d_out)
# [1] 80 16
#2. 敏感性分析
mydata <- harmonise_data(
exposure_dat=bmi_exp,
outcome_dat=t2d_out,
action= 2
)
res <- mr(mydata)
res
# 第一部分 异质性检验
het <- mr_heterogeneity(mydata)
het
mr(mydata,method_list=c('mr_ivw_mre')) #使用随机效应模型
#第二部分 多效性检验
pleio <- mr_pleiotropy_test(mydata)
pleio
#第三部分 逐个剔除检验
single <- mr_leaveoneout(mydata)
mr_leaveoneout_plot(single)
#3.可视化
#绘制散点图
mr_scatter_plot(res,mydata)
#绘制森林图
res_single <- mr_singlesnp(mydata)
mr_forest_plot(res_single)
# 绘制漏斗图
mr_funnel_plot(res_single)
GWAS
3.概念
4.分析代码
代码
install.packages("MatrixEQTL") # 安装R包
library("MatrixEQTL") # 加载R包
base.dir = find.package("MatrixEQTL") # 获取该包的路径信息
useModel = modelLINEAR # 三种模型可选(modelANOVA or modelLINEAR or modelLINEAR_CROSS)
SNP_file_name = paste(base.dir, "/data/SNP.txt", sep="") # 获取SNP文件位置
SNP_file = data.table::fread(SNP_file_name, header=T) # 读取SNP文件，可以在R中查看
expression_file_name = paste(base.dir, "/data/GE.txt", sep="") # 获取基因表达量文件位置
expression_file = data.table::fread(expression_file_name, header=T) # 读取基因表达量文件，可以在R中查看
covariates_file_name = paste(base.dir, "/data/Covariates.txt", sep="") # 读取协变量文件
covariates_file = data.table::fread(covariates_file_name, header=T) # 读取协变量文件，可在R中查看
output_file_name = tempfile() # 将输出文件设置为临时文件
pvOutputThreshold = 1e-2 # 定义gene-SNP associations的显著性P值
errorCovariance = numeric() # 定义误差项的协方差矩阵 (用的很少)
snps = SlicedData$new() # 创建SNP文件为S4对象（S4对象是该包的默认输入方式）
snps$fileDelimiter = "\t" # 指定SNP文件的分隔符
snps$fileOmitCharacters = "NA" # 定义缺失值
snps$fileSkipRows = 1 # 跳过第一行（适用于第一行是列名的情况）
snps$fileSkipColumns = 1 # 跳过第一列（适用于第一列是SNP ID的情况
snps$fileSliceSize = 2000 # 每次读取2000条数据
snps$LoadFile( SNP_file_name ) # 载入SNP文件
# 下面的代码同snps的那部分一致
gene = SlicedData$new()
gene$fileDelimiter = "\t"
gene$fileOmitCharacters = "NA"
gene$fileSkipRows = 1
gene$fileSkipColumns = 1
gene$fileSliceSize = 2000
gene$LoadFile( expression_file_name )
cvrt = SlicedData$new()
cvrt$fileDelimiter = "\t"
cvrt$fileOmitCharacters = "NA"
cvrt$fileSkipRows = 1
cvrt$fileSkipColumns = 1
cvrt$fileSliceSize = 2000
cvrt$LoadFile( covariates_file_name )
# 文件的输入部分结束
me = Matrix_eQTL_engine( # 这是进行eQTL分析的主要函数
snps = snps, # 指定SNP 文件
gene = gene, # 指定基因表达量文件
cvrt = cvrt, # 指定协变量文件
output_file_name = output_file_name, # 指定输出文件
pvOutputThreshold = pvOutputThreshold, # 指定显著性P值
useModel = useModel, # 指定使用的计算模型
errorCovariance = errorCovariance, # 指定误差项的协方差矩阵
verbose = TRUE,
pvalue.hist = TRUE,
min.pv.by.genesnp = FALSE,
noFDRsaveMemory = FALSE)
res <- me$all$eqtls # 把eQTL的显著结果存储到变量res里
res # 查看结果
4072
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
