GWAS power的计算

2023-12-07 20:37
文章标签 计算 power gwas

本文主要是介绍GWAS power的计算,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

import math
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.special import chdtri
from collections import defaultdict
%matplotlib inline
  1. 对于GWAS中power值(statistical power)的计算,用自己的话来说,Power值就是我们在Pvalue<α时(α:通常设为5e-8)显著水平下,原假设(H0)为假,接受备择假设(H1)的概率。
def pchisq(q,df,ncp=0):"""Calculates the cumulative of the chi-square distribution"""from scipy.stats import chi2,ncx2if ncp==0:result=chi2.cdf(x=q,df=df,loc=0,scale=1)else:result=ncx2.cdf(x=q,df=df,nc=ncp,loc=0,scale=1)return resultdef qchisq(p,df,ncp=0):"""Calculates the quantile function of the chi-square distribution"""from scipy.stats import chi2,ncx2if ncp==0:result=chi2.ppf(q=p,df=df,loc=0,scale=1)else:result=ncx2.ppf(q=p,df=df,nc=ncp,loc=0,scale=1)return result
def get_power(MAF, beta_alt, count, pvalue=5e-8):
#     MAF = 0.5
#     beta_alt = 0.2
#     count = 2000sigma = math.sqrt(1 - 2*MAF*(1-MAF)*beta_alt**2) # error sd after SNP effect is accounted for (see next part for explanation)ses = sigma/math.sqrt(count*2*MAF*(1-MAF)) # q_thresh = scipy.stats.chi2.ppf(q= 5e-8, df = 1)q_thresh = qchisq(1-pvalue, 1)pwr = 1- pchisq(q_thresh, 1, (beta_alt/ses)**2) return pwr
(1)给定效应值(Effect size,beta)为 1.2,计算在不同 MAF 和样本量下它所对应的 Power,并画出图。这里的 MAF 分为五个档位,突变频率从小到大分别是:0.01, 0.05, 0.1,0.15 和 0.2。把结果画在同一个图中,X 轴为样本量,Y 轴为 Power,画图时每一个 MAF 对应一条线。alpha 的阈值都定为 5e-8
# 计算在指定sample数据下,每个beta值对应的power值,并用字典保存
beta_value = np.arange(0, 1.2, 0.01)
count = 2000
pwr_dict=defaultdict(dict)
for maf in [0.01, 0.05, 0.1, 0.15 ,0.2]:pwr_dict[maf]=defaultdict(list)for beta in beta_value:pwr_dict[maf][beta]=get_power(maf, beta, count)
# 字典用pandas转为dataframe,再画图
df = pd.DataFrame(pwr_dict)
df.plot()
plt.title("Sample = 2000",fontsize=15)
plt.xlabel("Beta",fontsize=13)
plt.ylabel("Statistical Power",fontsize=13)
plt.legend(title = 'MAF')
<matplotlib.legend.Legend at 0x7f95af583670>

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

(2)给定样本量(sample size)为 2000,计算在不同 MAF 和效应值(Effect size, beta)下它所对应的 Power,并画出图。这里的 MAF 分为五个档位,突变频率从小到大分别是:0.01, 0.05, 0.1,0.15 和 0.2。把结果画在同一个图中,X 轴为效应值,Y 轴为 Power,画图时每一个 MAF 对应一条线。alpha 的阈值都定为 5e-8。
# 计算在指定sample数据下,每个beta值对应的power值,并用字典保存
beta_value = 0.2
count = np.arange(1, 4000, 10)
pwr_dict=defaultdict(dict)
for maf in [0.01, 0.05, 0.1, 0.15 ,0.2]:pwr_dict[maf]=defaultdict(list)for cnt in count:pwr_dict[maf][cnt]=get_power(maf, beta_value, cnt)
# 字典用pandas转为dataframe,再画图
df = pd.DataFrame(pwr_dict)
df.plot()
plt.title("Beta = 0.2",fontsize=15)
plt.xlabel("Sample Count",fontsize=13)
plt.ylabel("Statistical Power",fontsize=13)
plt.legend(title = 'MAF')
<matplotlib.legend.Legend at 0x7f95b7304b50>

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传


这篇关于GWAS power的计算的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/467285

相关文章

Python如何计算两个不同类型列表的相似度

《Python如何计算两个不同类型列表的相似度》在编程中,经常需要比较两个列表的相似度,尤其是当这两个列表包含不同类型的元素时,下面小编就来讲讲如何使用Python计算两个不同类型列表的相似度吧... 目录摘要引言数字类型相似度欧几里得距离曼哈顿距离字符串类型相似度Levenshtein距离Jaccard相

使用C#代码计算数学表达式实例

《使用C#代码计算数学表达式实例》这段文字主要讲述了如何使用C#语言来计算数学表达式,该程序通过使用Dictionary保存变量,定义了运算符优先级,并实现了EvaluateExpression方法来... 目录C#代码计算数学表达式该方法很长,因此我将分段描述下面的代码片段显示了下一步以下代码显示该方法如

如何用Java结合经纬度位置计算目标点的日出日落时间详解

《如何用Java结合经纬度位置计算目标点的日出日落时间详解》这篇文章主详细讲解了如何基于目标点的经纬度计算日出日落时间,提供了在线API和Java库两种计算方法,并通过实际案例展示了其应用,需要的朋友... 目录前言一、应用示例1、天安门升旗时间2、湖南省日出日落信息二、Java日出日落计算1、在线API2

poj 1113 凸包+简单几何计算

题意: 给N个平面上的点,现在要在离点外L米处建城墙,使得城墙把所有点都包含进去且城墙的长度最短。 解析: 韬哥出的某次训练赛上A出的第一道计算几何,算是大水题吧。 用convexhull算法把凸包求出来,然后加加减减就A了。 计算见下图: 好久没玩画图了啊好开心。 代码: #include <iostream>#include <cstdio>#inclu

uva 1342 欧拉定理(计算几何模板)

题意: 给几个点,把这几个点用直线连起来,求这些直线把平面分成了几个。 解析: 欧拉定理: 顶点数 + 面数 - 边数= 2。 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#inc

uva 11178 计算集合模板题

题意: 求三角形行三个角三等分点射线交出的内三角形坐标。 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#include <stack>#include <vector>#include <

XTU 1237 计算几何

题面: Magic Triangle Problem Description: Huangriq is a respectful acmer in ACM team of XTU because he brought the best place in regional contest in history of XTU. Huangriq works in a big compa

音视频入门基础:WAV专题(10)——FFmpeg源码中计算WAV音频文件每个packet的pts、dts的实现

一、引言 从文章《音视频入门基础:WAV专题(6)——通过FFprobe显示WAV音频文件每个数据包的信息》中我们可以知道,通过FFprobe命令可以打印WAV音频文件每个packet(也称为数据包或多媒体包)的信息,这些信息包含该packet的pts、dts: 打印出来的“pts”实际是AVPacket结构体中的成员变量pts,是以AVStream->time_base为单位的显

Spark MLlib模型训练—聚类算法 PIC(Power Iteration Clustering)

Spark MLlib模型训练—聚类算法 PIC(Power Iteration Clustering) Power Iteration Clustering (PIC) 是一种基于图的聚类算法,用于在大规模数据集上进行高效的社区检测。PIC 算法的核心思想是通过迭代图的幂运算来发现数据中的潜在簇。该算法适用于处理大规模图数据,特别是在社交网络分析、推荐系统和生物信息学等领域具有广泛应用。Spa

计算数组的斜率,偏移,R2

模拟Excel中的R2的计算。         public bool fnCheckRear_R2(List<double[]> lRear, int iMinRear, int iMaxRear, ref double dR2)         {             bool bResult = true;             int n = 0;             dou