「GWAS」
❝全基因组关联分析 ❞
「手动计算」
❝使用R语言编程GLM模型和Logistic模型,提取Effect和Pvalue ❞
「GLM」
❝一般线性模型 ❞
「Logistic」
❝主要分析广义线性模型,Y变量是二分类性状 ❞
「6-2」
❝这是我的
GWAS
学习笔记,更新到了6-2,更多专栏内容,拉到最后,点击链接阅读,或者点击开头的专辑。 ❞
GLM的手动计算GWAS分析的主要步骤:
y~x
做回归分析,计算x
的回归系数(Effect)和显著性(P-value)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
plink --file b --recodeA --out c
b
c
结果文件为c.raw
表型数据如果只有一个,可以放在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
lm
函数做回归分析data.table
c.raw
文件phe.txt
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结果:
可以看到,两者结果完全一致。
Logistic的手动计算GWAS分析的主要步骤:
y~x
做Logistic回归分析,计算x
的回归系数(Effect)和显著性(P-value)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
plink --file b --recodeA --out c
b
c
结果文件为c.raw
表型数据如果只有一个,可以放在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
glm
函数做Logistic回归分析data.table
c.raw
文件phe.txt
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
笔记 | GWAS 操作流程4-2:LM模型linear+数值协变量
笔记 | GWAS 操作流程4-4:LM模型+数值+因子协变量
笔记 | GWAS 操作流程4-5:LM模型+数值+因子+PCA协变量
笔记 | GWAS 操作流程5-1:根红苗正的GWAS分析软件:GEMMA