用R进行gwas meta分析,原来如此简单。

转载自品略图书馆 http://www.pinlue.com/article/2020/01/0912/569886542010.html

 

在生物信息数据分析中,R语言是必备技能,简洁的语法,丰富的生态,美观的可视化,种种优势使得其成为该领域中使用最广泛的编程语言之一。用R进行meta分析当然也是可以的,本文要介绍的R包rmeta, 就是其中之一,可以用于gwas meta分析,支持随机效应模型和固定效应模型,官方文档如下:

 

作为CRAN的一员,其安装方式如下

install.packages(“rmeta”)

首先来看下其输入文件格式,其输入文件并不是常规的GWAS分析结果, 比如pvalue, beta等统计指标,而是原始的频数分布,而且是二分类数据,比如两种Allele的频数分布,对于每个study, 统计如下所示的表格数据

AlleleAa

Caseab

Controlcd

也就是说,对于每个study,我们需要有4个数值。软件内置的示例数据如下

这个数据列数很多,真正进行分析时,只需要其中的前4列数据。在该R包中,两种模型对应的函数如下

以固定效应模型为例,进行meta分析的代码如下

从代码可以看出,分析时只使用了前4列,对应每个study的2X2表格数据,因为部分study中有缺失值NA, 有缺失值的数据是没法分析的,所以通过subset函数选取了其中的部分行。

该R包最大的特点是其可视化功能,对于meta分析的结果,提供了以下两种可视化的函数1. metaplot

该函数用于展示每个study的名称和对应OR值的分布,以及meta分析后最终计算出的OR值,用法如下

输出结果如下

2. forestplot

该函数用于自定义的展示很多额外的文字信息,用法如下

输出结果如下

通过rmeta可以方便个性化的绘制森林图,更多用法请参考官方的帮助文档。

·end·

 

GWAS meta分析是一种将多个基因组关联研究结果进行综合分析的方法,以下是一个Python脚本的例子,用于进行GWAS meta分析: ```python import pandas as pd import numpy as np import statsmodels.api as sm # 读取所有研究的GWAS结果文件 study1 = pd.read_csv("study1.csv") study2 = pd.read_csv("study2.csv") study3 = pd.read_csv("study3.csv") # 将每个研究的p进行变换,转化为z study1["z"] = np.sqrt(2) * sm.stats.proportion.proportions_ztest(study1["n_cases"], study1["n_total"], value=study1["OR"])[0] study2["z"] = np.sqrt(2) * sm.stats.proportion.proportions_ztest(study2["n_cases"], study2["n_total"], value=study2["OR"])[0] study3["z"] = np.sqrt(2) * sm.stats.proportion.proportions_ztest(study3["n_cases"], study3["n_total"], value=study3["OR"])[0] # 将每个研究的z和样本量进行合并 meta_data = pd.concat([study1[["SNP", "z", "n_cases", "n_total"]], study2[["SNP", "z", "n_cases", "n_total"]], study3[["SNP", "z", "n_cases", "n_total"]]]) # 计算每个SNP的z和加权样本量 meta_data["wz"] = meta_data["z"] * np.sqrt(meta_data["n_cases"] + meta_data["n_controls"]) meta_data["w"] = np.sqrt(meta_data["n_cases"] + meta_data["n_controls"]) # 计算meta分析的z和p meta_z = meta_data["wz"].sum() / meta_data["w"].sum() meta_p = 2 * (1 - sm.stats.norm.cdf(abs(meta_z))) # 输出meta分析结果 print("Meta-analysis result:") print("z value: ", meta_z) print("p value: ", meta_p) ``` 这个脚本首先读取每个研究的GWAS结果文件,将每个研究的p转化为z,然后将每个研究的z和样本量进行合并,计算每个SNP的加权z和加权样本量,最后计算meta分析的z和p,输出结果。需要注意的是,不同研究的样本量、OR等参数可能存在差异,需要根据具体情况进行调整。
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值