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

文件：

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参数，考虑数值+因子协变量，然后将结果与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)

可以看出R语言lm的结果中：

M7: beta为1.3940， se为1.3165，t值为1.059，p值为0.290，

M9：beta为3.3265，se为1.5042，t值为2.211，p值为0.0272

对比assoc的结果：

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

问题来了：

如果R中跑完10000个标记的回归分析，需要很长时间，而plink只需要几秒。没有比较久没有伤害。不过R适合验证性的分析，比如你在进行plink分析时，看到LM模型，说明书上写的是回归分析，但是你用R语言编写后，验证了这是lm回归分析，这种喜悦是什么？

这种喜悦就是，明知屁是臭的，偏偏还要闻一下。虽然很多人告诉你是臭的，但还是想要自己闻一下。所谓“纸上得来终觉浅，绝知此事要躬行”便是如此啊，闻过之后，可以告诉别人说：“真香！”

起码我自己的收获：

1，知道GWAS分析中，plink进行LM分析，SNP转化为012,0是major，2是minor，这顺序是不能变的，你如果要手动进行回归分析，需要将SNP的分型进行转化，而不是粗暴的将其替换为0,1，2.

2，R在进行单个SNP的分析时，0,1，2是作为数值进行的回归分析，而不是因子。

最后感言：

真香！