首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >笔记 | GWAS 操作流程4-5:LM模型+数值+因子+PCA协变量

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

作者头像
邓飞
发布2020-05-29 11:00:27
2.8K0
发布2020-05-29 11:00:27
举报

「飞哥感言:」

❝从开始介绍plink做GWAS数据的质控,到构建模型,到定义协变量,已经灌了很多水,这篇是plink做GWAS的结尾,因为plink做GWAS只有两个模型可以用:GLM和logistic,前者分析数量性状,后者分析二分类性状。而现在GWAS更多使用LMM模型,这个模型plink没法做,以后几篇介绍GEMMA的操作方法。 ❞

1. 协变量文件整理

第一列为FID 第二列为ID 第三列以后为协变量(注意,只能是数字,不能是字符!)

这里协变量文件为:

[dengfei@ny 03_linear_cov]$ head cov.txt 
1061 1061 F 3
1062 1062 M 3
1063 1063 F 3
1064 1064 F 3
1065 1065 F 3
1066 1066 F 3
1067 1067 F 3
1068 1068 M 3
1069 1069 M 3
1070 1070 M 3

  • 这里,首先将F换为1,M换为2,将其转化为连续变量(数字)
  • 然后,将世代变为虚拟变量
  • 最后,将两个协变量整合到一起
sed 's/F/1/g' cov.txt >cov2.txt
sed -i 's/M/2/g' cov2.txt 

2. 使用plink的dummy coding转化为虚拟变量

plink --file b --covar cov2.txt --write-covar --dummy-coding

结果生成:

plink.cov

3. 使用plink获得pca结果

plink --file b --pca 3 

「结果文件:」

4. 将pca结果和协变量结果合并

sed '1d' plink.cov >a.txt
head a.txt 
awk '{print $3,$4,$5}' plink.eigenvec >pca.txt
head a.txt 
wc -l pca.txt a.txt 
paste a.txt pca.txt >pca_cov.txt

「合并后的协变量:」

5 进行协变量GWAS分析LM模型

「代码:」

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

「日志:」

PLINK v1.90b5.3 64-bit (21 Feb 2018)           www.cog-genomics.org/plink/1.9/
(C) 2005-2018 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to re.log.
Options in effect:
  --allow-no-sex
  --covar pca_cov.txt
  --file b
  --linear
  --out re
  --pheno phe.txt

515199 MB RAM detected; reserving 257599 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).
--covar: 6 covariates loaded.
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

「结果预览:」

4. 使用R语言进行结果比较lm+factor+pca

library(data.table)
geno = fread("c.raw")
geno[1:10,1:10]
phe = fread("phe.txt")
plink = fread("pca_cov.txt",header=F,sep=" ")
head(plink)
dd = data.frame(phe = phe$V3,cov1 = plink$V3,cov2 = plink$V4,cov3=plink$V5,pca1 = plink$V6,pca2 = plink$V7,pca3 = plink$V8,geno[,7:20])
head(dd)
mod_M7 = lm(phe ~ cov1+cov2+cov3+pca1+pca2+pca3 + M7_1,data=dd);summary(mod_M7)



「M7加上因子协变量结果:」

结果完全一样。

5. 结论

plink中一般线性模型(LM),linear可以支持数值协变量,因子协变量(经过转化),pca等等,这些过程都可以通过R语言的lm函数复现结果。

6. 一般线性模型可以用plink做,那么混合线性模型怎么做?gemma!

gemma也可以做一般线性模型,也可以做混合线性模型。plink只可以做一般线性模型,gemma可以利用plink的数据格式做一般线性模型和混合线性模型,这就很厉害了。

「遗憾:gemma只有linux版本,所以后面的分析在linux系统下。」

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2020-05-27,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 育种数据分析之放飞自我 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1. 协变量文件整理
  • 2. 使用plink的dummy coding转化为虚拟变量
  • 3. 使用plink获得pca结果
  • 4. 将pca结果和协变量结果合并
  • 5 进行协变量GWAS分析LM模型
  • 4. 使用R语言进行结果比较lm+factor+pca
  • 5. 结论
  • 6. 一般线性模型可以用plink做,那么混合线性模型怎么做?gemma!
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档