用户工具

侧边栏

用R计算多基因病遗传力

Falconer方法

根据Falconer方法, 多基因病的遗传力计算方法如下:


一些教材上使用查表来获取Xg和ag以及Xr的值,然而现在电脑已经普及,使用电脑可以非常容易的计算这些值,无需再去查表。这里介绍怎么使用R来计算,可以将公式写成函数,简单带入一般人群发病率、患者亲属发病率以及亲缘系数即可一步出结果。

正态分布

关于正态分布,基本的分布曲线以及累计概率曲线如下图,可以很直观地看到我们所需要求的值使用哪个函数(图中红色曲线为标准正态分布):

这里不要被这些函数吓到,了解点理论,然后知道怎么用电脑算即可。

分步计算

从上边的图中可以看出,将易患性分布标准化之后,X实际上是发病率为q时对应的易患性的值,也就是累计概率为1-q时的X值,求X的值,使用的是累计概率曲线函数的反函数,φ-1, 对应R里的qnorm函数。
比如,qg为0.1%,那么我们可以这样来求Xg:

q_g <- 0.1/100
X_g <- qnorm(1 - q_g)
X_g

运行示例如下:


同样的办法求Xr即可。

对于ag的计算,需要先算出Xg,根据X算Z使用的是dnorm函数,示例如下:

q_g <- 0.1/100
X_g <- qnorm(1 - q_g)
z_g <- dnorm(X_g)
a_g <- z_g/q_g

运行示例如下:

现在,我们可以把这些连起来,对于指定的qg、qr以及r来计算h2:

q_g <- 0.1/100
q_r <- 22/669    # 669个一级亲属中有22个患病
r   <- 1/2       # 一级亲属的亲缘系数是1/2
X_g <- qnorm(1 - q_g)
z_g <- dnorm(X_g)
a_g <- z_g/q_g
X_r <- qnorm(1 - q_r)
b   <- (X_g - X_r)/a_g
h2  <- b/r   
h2 

这里仅展示一下运行结果:


如果要多次计算,直接把计算过程封装成为函数,那么可以省去很多的步骤,如下

h <- function(q_g, q_r, r=0.5){   #r=0.5是说r的默认值为0.5,这样在算一级亲属的情况的时候,可以不提供这个参数
    X_g = qnorm(1-q_g)
    X_r = qnorm(1-q_r)
    z_g = dnorm(X_g)
    a_g = z_g/q_g
    b   = (X_g - X_r)/a_g
    return(b/r)
}
 
h(0.1/100, 22/669)

看看运行的效果:

这样我们就能够很方便的计算别的情况了,比如:

自己练习

对于没有一般人群的发病率数据时,使用对照组,比如使用对照组一级亲属数据的情况如下:

试试看,写个函数出来。


自2021年9月22日以来,本网页已被访问5567次。


页面工具

Powered by DokuWiki | ©2018-8102 Hong Wang | 桂ICP备18010616号 | 桂公网安备 45010302002017号