笔记 GWAS 操作流程6-2：手动计算GWAS分析中的GLM和Logistic模型

时间:2022-07-26
本文章向大家介绍笔记 GWAS 操作流程6-2：手动计算GWAS分析中的GLM和Logistic模型，主要内容包括其使用实例、应用技巧、基本知识点总结和需要注意事项，具有一定的参考价值，需要的朋友可以参考一下。

1. 名词解释

「GWAS」

❝全基因组关联分析 ❞

「手动计算」

❝使用R语言编程GLM模型和Logistic模型，提取Effect和Pvalue ❞

「GLM」

❝一般线性模型 ❞

「Logistic」

❝主要分析广义线性模型，Y变量是二分类性状 ❞

「6-2」

❝这是我的GWAS学习笔记，更新到了6-2，更多专栏内容，拉到最后，点击链接阅读，或者点击开头的专辑。 ❞

2. GLM模型

GLM的手动计算GWAS分析的主要步骤：

  • 1，将SNP的分型转化为0-1-2（0位次等位基因），数字格式（x变量）
  • 2，将性状观测值作为y变量（GLM一般分析连续性状）
  • 3，对y~x做回归分析，计算x的回归系数（Effect）和显著性（P-value）
  • 4，如果有协变量，加到x后面，进行回归分析（因子变量变为数字哑变量）

「示例：」

共有1500个个体，10000个SNP

[dengfei@ny 01_assoc]$ wc -l b.*
   10000 b.map
    1500 b.ped
   11500 总用量

表型数据为连续性状：

1061 1061 -3.190926
1062 1062 +24.290128
1063 1063 -19.403765
1064 1064 -0.815962
1065 1065 -19.073081
1066 1066 -21.106496
1067 1067 +15.020220
1068 1068 -15.985445

2.1 将plink数据转化为0-1-2

plink --file b --recodeA --out c
  • --file # 紧接plink的前缀名称b
  • --recodeA # 转化为0-1-2编码
  • --out # 输出文件名c

结果文件为c.raw

2.2 表型数据整理

表型数据如果只有一个，可以放在plink文件的ped数据的第六列，也可以单独拉出来：

1061 1061 -3.190926
1062 1062 +24.290128
1063 1063 -19.403765
1064 1064 -0.815962
1065 1065 -19.073081
1066 1066 -21.106496
1067 1067 +15.020220
1068 1068 -15.985445
1069 1069 +5.849143
1070 1070 +39.513181
  • 第一列为FID # 家系ID
  • 第二列为IID # 个体ID
  • 第三列为表型值 # 表型数据

2.3 使用R中的lm函数做回归分析

  • 1，首先载入软件包data.table
  • 2，然后读取0-1-2编码的c.raw文件
  • 3，然后读取表型数据文件phe.txt
  • 4，然后将表型数据和基因型数据合并
library(data.table)
geno = fread("c.raw",header=T)
phe = fread("../phe.txt")
dd = data.frame(phe$V3,geno[,7:17])

合并后的数据类型如下：

> summary(dd)
     phe.V3              M1_0        M2_0        M3_0        M4_0        M5_0
 Min.   :-92.6750   Min.   :0   Min.   :0   Min.   :0   Min.   :0   Min.   :0
 1st Qu.:-17.4169   1st Qu.:0   1st Qu.:0   1st Qu.:0   1st Qu.:0   1st Qu.:0
 Median : -0.7173   Median :0   Median :0   Median :0   Median :0   Median :0
 Mean   :  0.4054   Mean   :0   Mean   :0   Mean   :0   Mean   :0   Mean   :0
 3rd Qu.: 19.5060   3rd Qu.:0   3rd Qu.:0   3rd Qu.:0   3rd Qu.:0   3rd Qu.:0
 Max.   : 88.1515   Max.   :0   Max.   :0   Max.   :0   Max.   :0   Max.   :0
      M6_0        M7_1             M8_0        M9_1           M10_0
 Min.   :0   Min.   :0.0000   Min.   :0   Min.   :0.000   Min.   :0
 1st Qu.:0   1st Qu.:0.0000   1st Qu.:0   1st Qu.:0.000   1st Qu.:0
 Median :0   Median :0.0000   Median :0   Median :0.000   Median :0
 Mean   :0   Mean   :0.3387   Mean   :0   Mean   :0.256   Mean   :0
 3rd Qu.:0   3rd Qu.:1.0000   3rd Qu.:0   3rd Qu.:0.000   3rd Qu.:0
 Max.   :0   Max.   :2.0000   Max.   :0   Max.   :2.000   Max.   :0
     M11_2
 Min.   :0.000
 1st Qu.:0.000
 Median :0.000
 Mean   :0.256
 3rd Qu.:0.000
 Max.   :2.000

「用M7_1这个位点做回归分析`」

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

结果如下：

> summary(mod_M7)

Call:
lm(formula = phe.V3 ~ M7_1, data = dd)

Residuals:
    Min      1Q  Median      3Q     Max
-92.608 -17.664  -0.871  18.708  86.824

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.0667     0.8423  -0.079    0.937
M7_1          1.3940     1.3165   1.059    0.290

Residual standard error: 27.68 on 1498 degrees of freedom
Multiple R-squared:  0.0007479, Adjusted R-squared:  8.079e-05
F-statistic: 1.121 on 1 and 1498 DF,  p-value: 0.2898

可以看到，M7_1的Effect为1.394，对应的P-value为0.290

对比plink的GLM结果：

可以看到，两者结果完全一致。

3. Logistic回归模型

Logistic的手动计算GWAS分析的主要步骤：

  • 1，将SNP的分型转化为0-1-2（0位次等位基因），数字格式（x变量）
  • 2，将性状观测值作为y变量（Logistic一般分析二分类性状）
  • 3，对y~x做Logistic回归分析，计算x的回归系数（Effect）和显著性（P-value）
  • 4，如果有协变量，加到x后面，进行回归分析（因子变量变为数字哑变量）

「示例：」

共有112个个体，1113612个SNP

[dengfei@ny R_logistic]$ wc -l c*
  1113612 b.map
      112 b.ped
  1113724 总用量

表型数据为二分类性状：

1328 NA06989 2
1377 NA11891 2
1349 NA11843 1
1330 NA12341 2
1328 NA06984 2
1418 NA12275 1
13291 NA06986 1
1418 NA12272 1
13292 NA07051 2
1354 NA12400 2

3.1 将plink数据转化为0-1-2

plink --file b --recodeA --out c
  • --file # 紧接plink的前缀名称b
  • --recodeA # 转化为0-1-2编码
  • --out # 输出文件名c

结果文件为c.raw

3.2 表型数据整理

表型数据如果只有一个，可以放在plink文件的ped数据的第六列，也可以单独拉出来：

1328 NA06989 2
1377 NA11891 2
1349 NA11843 1
1330 NA12341 2
1328 NA06984 2
1418 NA12275 1
13291 NA06986 1
1418 NA12272 1
13292 NA07051 2
1354 NA12400 2
  • 第一列为FID # 家系ID
  • 第二列为IID # 个体ID
  • 第三列为表型值 # 表型数据，默认是1-2编码（case-control）

3.3 使用R中的glm函数做Logistic回归分析

  • 1，首先载入软件包data.table
  • 2，然后读取0-1-2编码的c.raw文件
  • 3，然后读取表型数据文件phe.txt
  • 4，然后将表型数据和基因型数据合并
library(data.table)
geno[1:10,1:10]
phe = fread("pheno.txt")
dd = data.frame(phe$V3,geno[,7:17])

合并后的数据类型如下：

> summary(dd)
     phe.V3     rs3131972_A      rs3131969_A      rs1048488_C
 Min.   :1.0   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000
 1st Qu.:1.0   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000
 Median :1.5   Median :0.0000   Median :0.0000   Median :0.0000
 Mean   :1.5   Mean   :0.3304   Mean   :0.2679   Mean   :0.3333
 3rd Qu.:2.0   3rd Qu.:1.0000   3rd Qu.:0.2500   3rd Qu.:1.0000
 Max.   :2.0   Max.   :2.0000   Max.   :2.0000   Max.   :2.0000
                                                 NA's   :1
  rs12562034_A     rs12124819_G     rs4040617_G      rs4970383_A
 Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000
 1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000
 Median :0.0000   Median :0.0000   Median :0.0000   Median :0.0000
 Mean   :0.2054   Mean   :0.5804   Mean   :0.2589   Mean   :0.3929
 3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:1.0000
 Max.   :2.0000   Max.   :2.0000   Max.   :2.0000   Max.   :2.0000

  rs4475691_T      rs1806509_C      rs7537756_G      rs2340587_G
 Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000
 1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000
 Median :0.0000   Median :1.0000   Median :0.0000   Median :0.0000
 Mean   :0.2857   Mean   :0.7768   Mean   :0.3214   Mean   :0.3571
 3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000
 Max.   :2.0000   Max.   :2.0000   Max.   :2.0000   Max.   :2.0000

「用rs3131972_A这个位点做Logistic回归分析`」

「注意：R中glm模型，Logistic需要Y变量为0-1分布，而我们的表型数据为1-2，所以讲表型数据减去1」

dd$phe.V3 = dd$phe.V3-1
m1 = glm(phe.V3 ~ rs3131972_A,family = "binomial",data=dd )
summary(m1)

结果如下：

> summary(m1)

Call:
glm(formula = phe.V3 ~ rs3131972_A, family = "binomial", data = dd)

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-1.43395  -1.12889  -0.09398   1.22671   1.22671

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.1152     0.2261  -0.510    0.610
rs3131972_A   0.3503     0.3773   0.928    0.353

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 155.26  on 111  degrees of freedom
Residual deviance: 154.39  on 110  degrees of freedom
AIC: 158.39

Number of Fisher Scoring iterations: 3

可以看到，rs3131972_A的Effect为0.3503，对应的P-value为0.353

对比plink的Logistic结果：

可以看到，两者结果完全一致。

「注意：」

❝plink中，默认输出的不是Effect，而是OR值，R语言中如果要输出OR值，可以用exp(coef(m1))将结果打印出来。plink中如果想要输出BETA的效应值，加上上参数--beta ❞

> exp(coef(m1))
(Intercept) rs3131972_A
  0.8911777   1.4195186

「输出OR值：」

plink --file c --pheno pheno.txt --logistic --out x1

「输出BETA值：」

plink --file c --pheno pheno.txt --logistic --out x2 --beta

4. GWAS系列相关

笔记 | GWAS 操作流程1：下载数据

笔记 | GWAS 操作流程2-1：缺失质控

笔记 | GWAS 操作流程2-2：性别质控

笔记 | GWAS 操作流程2-3：最小等位基因频率

笔记 | GWAS 操作流程2-4：哈温平衡检验

笔记 | GWAS 操作流程2-5：杂合率检验

笔记 | GWAS 操作流程2-6：去掉亲缘关系近的个体

笔记 | GWAS 操作流程3：plink关联分析

笔记 | GWAS 操作流程4-1：LM模型assoc

笔记 | GWAS 操作流程4-2：LM模型linear+数值协变量

笔记 | GWAS 操作流程4-3：LM模型+因子协变量

笔记 | GWAS 操作流程4-4：LM模型+数值+因子协变量

笔记 | GWAS 操作流程4-5：LM模型+数值+因子+PCA协变量

笔记 | GWAS 操作流程5-1：根红苗正的GWAS分析软件：GEMMA

笔记 | GWAS 操作流程5-2：利用GEMMA软件进行LMM+PCA+协变量

笔记 | GWAS 操作流程6-1：为何有些SNP效应值大却不显著？

随机文章
本站知识点必读
最近更新