「飞哥感言：」
❝从开始介绍plink做GWAS数据的质控，到构建模型，到定义协变量，已经灌了很多水，这篇是plink做GWAS的结尾，因为plink做GWAS只有两个模型可以用：GLM和logistic，前者分析数量性状，后者分析二分类性状。而现在GWAS更多使用LMM模型，这个模型plink没法做，以后几篇介绍GEMMA的操作方法。 ❞
1. 协变量文件整理
第一列为FID 第二列为ID 第三列以后为协变量（注意，只能是数字，不能是字符！）
这里协变量文件为：
代码语言：javascript复制[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转化为虚拟变量
代码语言：javascript复制plink --file b --covar cov2.txt --write-covar --dummy-coding
结果生成：
plink.cov
3. 使用plink获得pca结果
代码语言：javascript复制plink --file b --pca 3
「结果文件：」
4. 将pca结果和协变量结果合并
代码语言：javascript复制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模型
「代码：」
代码语言：javascript复制plink --file b --pheno phe.txt --allow-no-sex --linear --covar pca_cov.txt --out re
「日志：」
代码语言：javascript复制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
代码语言：javascript复制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系统下。」
