文件:
之前的教程中,我们使用的是别人模拟的数据,数据类型是二分类数据,这里我们模拟一个数量性状的连续性状,做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表型数据:
代码解释:
–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 v3Logging to re.log.Options in effect: --allow-no-sex --assoc --file b --out re --pheno phe.txt7820 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的结果:
这里,利用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 v3Logging to re.log.Options in effect: --allow-no-sex --file b --linear --out re --pheno phe.txt7820 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是作为数值进行的回归分析,而不是因子。
最后感言:
真香!
5. 相关系列笔记 | GWAS 操作流程1:下载数据
笔记 | GWAS 操作流程2-1:缺失质控
笔记 | GWAS 操作流程2-2:性别质控
笔记 | GWAS 操作流程2-3:最小等位基因频率
笔记 | GWAS 操作流程2-4:哈温平衡检验
笔记 | GWAS 操作流程2-5:杂合率检验
笔记 | GWAS 操作流程2-6:去掉亲缘关系近的个体
笔记 | GWAS 操作流程3:plink关联分析