博文

GWAS分析中：为何有些SNP效应值大却不显著？

已有 5230 次阅读 2020-9-16 22:17 |个人分类:农学统计|系统分类:科研笔记

之前的GWAS分析中，我查看结果时有过这个疑问，随着不断的学习理解，我已经知道了为什么，这里将我思考的结果分享一下。

1. 错误的理解

某一个SNP，效应值（Effect）越大，就越显著？？？

「事实上：」效应值和显著性是两码事！

2. 用代码说话

GWAS分析中，最简单的就是一般线性模型（GLM），而GLM模型进行的GWAS分析中，用R语言实现的代码如下：

mod_M7 = lm(phe.V3 ~ M7_1,data=dd)
summary(mod_M7)

「结果：」这个里面：

Estimate就是效应值：3.3265，这个就是SNP M9的效应值
Pr就是P值：0.0272，这个就是SNP M9的P值

可以看到，这两个是两个指标，他们之间没有必然的联系。

3. GWAS分析的思路

1，将SNP的分型转化为0, 1, 2的形式，主效纯合（major）编码为0，杂合编码为1，次等位纯合编码为2
2，x变量为0,1,2的数值，y变量为性状表型值
3，对x和y做回归分析y ~ a*x + b，a为效应值，a是否显著的P值，为SNP的P值

4. SNP效应值很大，却不显著


# SNP 分型
set.seed(666)
x = rep(c(0,1,2),3)
x

# 模拟表型值
y = 12*x + rnorm(9)*10

# 作回归分析
mod = lm(y ~ x)
summary(mod)

# 作图
dd = data.frame(x,y)
dd
ggplot(dd,aes(x=x,y=y)) + geom_point() + stat_smooth(method='lm',formula = y~x,colour='red')

看一下这个数据：

> dd
  x          y
1 0   7.533110
2 1  32.143547
3 2  20.448655
4 0  20.281678
5 1 -10.168745
6 2  31.583962
7 0 -13.061853
8 1   3.974804
9 2   6.077592

计算每种分型对应的表型平均值：

> # 平均值
> aggregate(y~x,data = dd,mean)
  x         y
1 0  4.917645
2 1  8.649869
3 2 19.370070

可以看出，分型0对应的是4.9，最小，分型1对应的是8.6，分型2对应的是19.3，趋势是比较明显的。

「看一下回归分析的结果：」

> # 作回归分析
> mod = lm(y ~ x)
> summary(mod)

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-21.148 -12.128   2.243  13.379  21.164 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    3.753      8.558   0.439    0.674
x              7.226      6.629   1.090    0.312

Residual standard error: 16.24 on 7 degrees of freedom
Multiple R-squared:  0.1451, Adjusted R-squared:  0.023 
F-statistic: 1.188 on 1 and 7 DF,  p-value: 0.3118

可以看出，效应值为7.226，但是P值却是0.312，不显著。

所以说，效应值大的SNP位点，不一定是显著的。

「为何会出现这种情况呢？」我们画一个散点图看一下分布情况：

library(ggplot2)
ggplot(dd,aes(x=x,y=y)) + geom_point() + stat_smooth(method='lm',formula = y~x,colour='red')

在这里插入图片描述可以看到，分型为1的表型值，有一个非常高，达到了31，相当于在进行T检验时，标准误se比较高，导致P值较大，不显著。

5. SNP效应值很小，却极显著

我们还模拟数据：

# SNP 分型
set.seed(666)
x = rep(c(0,1,2),3)
x

# 模拟表型值
y = 0.1*x + rnorm(9)*0.01

# 作回归分析
mod = lm(y ~ x)
summary(mod)

# 作图
dd = data.frame(x,y)
dd
library(ggplot2)
ggplot(dd,aes(x=x,y=y)) + geom_point() + stat_smooth(method='lm',formula = y~x,colour='red')


# 平均值
aggregate(y~x,data = dd,mean)

「模拟的数据如下：」

> dd
  x           y
1 0  0.00753311
2 1  0.12014355
3 2  0.19644866
4 0  0.02028168
5 1  0.07783126
6 2  0.20758396
7 0 -0.01306185
8 1  0.09197480
9 2  0.18207759

计算每种分型对应的表型平均值：

> aggregate(y~x,data = dd,mean)
  x           y
1 0 0.004917645
2 1 0.096649869
3 2 0.195370070

可以看到，三种基因型，阶梯感也很强。

「回归分析结果：」

> # 作回归分析
> mod = lm(y ~ x)
> summary(mod)

Call:
lm(formula = y ~ x)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.021148 -0.012128  0.002243  0.013379  0.021164 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.003753   0.008558   0.439    0.674    
x           0.095226   0.006629  14.365 1.89e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01624 on 7 degrees of freedom
Multiple R-squared:  0.9672, Adjusted R-squared:  0.9625 
F-statistic: 206.4 on 1 and 7 DF,  p-value: 1.886e-06