

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

育种数据分析之放飞自我

2020-05-19

导读：神奇的是，光计划写了一篇也是没谁了……

1. GWAS笔记操作计划

之前的教程中，我们使用的是别人模拟的数据，数据类型是二分类数据，这里我们模拟一个数量性状的连续性状，做GWAS更有代表性。

我们先从没有协变量的一般线性模型（LM）模型开始，然后加入数据类型的协变量，然后加入因子类型的协变量（这里需要进行虚拟变量的转化），然后将数值协变量和因子变量放在一起作为协变量，然后将PCA的值作为协变量加进去，这样一般线性模型分析完成。

最后我们使用混合线性模型（LMM），不考虑协变量，然后加入协变量，然后将PCA加入协变量。

最最后介绍一下可视化，包括曼哈顿图，QQ图，LD衰减图等相关可视化。

最最最后，介绍一下这个数据，用这么多数据分析后，哪种模型最好，哪种模型控制假阳性较好，哪种模型控制假阴性较好。

最最最最后，介绍一下结果的有疑问的地方，比如为何有些SNP显著但是效应较小，有些SNP效应较大但是显著性较低？比如GWAS分析结果中效应值是怎么理解的？GWAS分析的效应值和rrBLUP中SNP效应值有什么关系？

计划章节：

笔记 GWAS 操作流程4-1：LM模型assoc（本篇）

❝
主要是介绍assoc和linear两个参数，不考虑协变量，然后将结果与R语言编程结果比较
❞

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

❝
主要是介绍linear参数，考虑数值协变量，然后将结果与R语言编程结果比较
❞

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

❝
主要是介绍linear参数，考虑因子协变量，然后将结果与R语言编程结果比较
❞

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

❝
主要是介绍linear参数，考虑数值+因子协变量，然后将结果与R语言编程结果比较
❞

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

❝
主要是介绍linear参数，考虑数值+因子+PCA协变量，然后将结果与R语言编程结果比较
❞

笔记 GWAS 操作流程5-1：LMM模型进行GWAS分析

❝
主要是介绍混合线性模型，不考虑协变量结果分析
❞

笔记 GWAS 操作流程5-2：LMM模型进行GWAS分析+数值协变量

❝
主要是介绍混合线性模型，考虑数值协变量结果分析
❞

笔记 GWAS 操作流程5-3：LMM模型进行GWAS分析+因子协变量

❝
主要是介绍混合线性模型，考虑因子协变量结果分析
❞

笔记 GWAS 操作流程5-4：LMM模型进行GWAS分析+数值+因子协变量

❝
主要是介绍混合线性模型，考虑数值+因子协变量结果分析
❞

笔记 GWAS 操作流程5-5：LMM模型进行GWAS分析+数值+因子+PCA协变量

❝
主要是介绍混合线性模型，考虑数值+因子+PCA协变量结果分析
❞

笔记 GWAS 操作流程6-1：GWAS分析中曼哈顿图的绘制

❝
主要是介绍GWAS分析中曼哈顿图的绘制
❞

笔记 GWAS 操作流程6-2：GWAS分析中QQ的绘制

❝
主要是介绍GWAS分析中QQ图的绘制
❞

笔记 GWAS 操作流程6-3：GWAS分析中其它图的绘制

❝
主要是介绍GWAS分析中其它图的绘制
❞

笔记 GWAS 操作流程7：GWAS分析中结果解读

❝
主要是介绍GWAS分析中结果解读
❞

2. 模拟的数据

模拟数据使用QMSim软件，模拟的数据如下：

b.map  b.ped  cov.txt  phe.txt

这里b.map和b.ped是plink格式的数据，cov.txt为协变量，phe.txt为观测值

「数据预览：」

b.map数据：

b.ped数据：

cov.txt 协变量数据：

phe.txt表型数据：

3. plink利用assoc进行LM的GWAS分析

plink --file b --pheno phe.txt --assoc --out re --allow-no-sex

「代码解释：」

--file为plink的文件名的前缀
--pheno 为plink分析GWAS的表型数据，注意前两列为FID IID，第三列为分析的性状
--assoc 利用的是LM模型
--out 输出文件名
--allow-no-sex 因为数据中没有性别信息，需要加上这个参数，允许没有性别的数据

「运行日志：」

PLINK v1.90b6.17 64-bit (28 Apr 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to re.log.
Options in effect:
  --allow-no-sex
  --assoc
  --file b
  --out re
  --pheno phe.txt

7820 MB RAM detected; reserving 3910 MB for main workspace.
.ped scan complete (for binary autoconversion).
Performing single-pass .bed write (10000 variants, 1500 people).
--file: re-temporary.bed + re-temporary.bim + re-temporary.fam written.
10000 variants loaded from .bim file.
1500 people (0 males, 0 females, 1500 ambiguous) loaded from .fam.
Ambiguous sex IDs written to re.nosex .
1500 phenotype values present after --pheno.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1500 founders and 0 nonfounders present.
Calculating allele frequencies... done.
10000 variants and 1500 people pass filters and QC.
Phenotype data is quantitative.
Writing QT --assoc report to re.qassoc ... done.

「结果文件：」re.qassoc注意，这里的SNP没有质控，有很多SNP没有分型，所以结果有很多NA，不过这里可以看出「assoc」分析GWAS的结果，其中SNP中M7和M9的结果：

M7，beta值为1.394，se为1.317，R2为0.000749，T值为1.059，P值为0.2898
M9，beta值为3.327，se为1.504，R2为0.003254，T值为2.211，P值为0.02715

4. plink利用linear进行LM的GWAS分析

这里，利用linear和上面的assoc命令类似，将linear代替assoc即可。一般来说，linear可以支持协变量，assoc不支持协变量。assoc运行速度快于linear。

plink --file b --pheno phe.txt --allow-no-sex --linear --out re

「运行日志：」

PLINK v1.90b6.17 64-bit (28 Apr 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to re.log.
Options in effect:
  --allow-no-sex
  --file b
  --linear
  --out re
  --pheno phe.txt

7820 MB RAM detected; reserving 3910 MB for main workspace.
.ped scan complete (for binary autoconversion).
Performing single-pass .bed write (10000 variants, 1500 people).
--file: re-temporary.bed + re-temporary.bim + re-temporary.fam written.
10000 variants loaded from .bim file.
1500 people (0 males, 0 females, 1500 ambiguous) loaded from .fam.
Ambiguous sex IDs written to re.nosex .
1500 phenotype values present after --pheno.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1500 founders and 0 nonfounders present.
Calculating allele frequencies... done.
10000 variants and 1500 people pass filters and QC.
Phenotype data is quantitative.
Writing linear model association results to re.assoc.linear ... done.

「结果文件：」

re.assoc.linear

结果和assoc的结果完全一致！

5. 使用R语言的一般线性回归比较结果

这里，需要将plink转化为012的形式，可以使用recodeA参数：

plink --file b --recodeA --out c

结果：c.raw

「编写R程序：」

library(data.table)
geno = fread("c.raw",header=T)
geno[1:10,1:20]
phe = fread("../phe.txt")
head(phe)
dd = data.frame(phe$V3,geno[,7:17])
head(dd)
mod_M7 = lm(phe.V3 ~ M7_1,data=dd)
summary(mod_M7)