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
这里第三列为性别，第四列为世代，这里，将世代作为因子，进行因子协变量的GWAS分析
2. 因子协变量
代码语言：javascript复制awk '{print $1,$2,$4}' cov.txt >cov1.txt
数据如下：
代码语言：javascript复制1061 1061 3
1062 1062 3
1063 1063 3
1064 1064 3
1065 1065 3
1066 1066 3
1067 1067 3
1068 1068 3
1069 1069 3
1070 1070 3
3. 使用plink的dummy coding转化为虚拟变量
代码语言：javascript复制plink --file b --covar cov1.txt --write-covar --dummy-coding
结果生成：
plink.cov
「注意：」这里的协变量，会减少一个水平，比如本来世代是由3,4，5三个世代，这里只有两个水平。plink文档是这样解释的：
That is, for a variable with K categories, K-1 new dummy variables are created. This new file can be used with --linear and --logistic, and a coefficient for each level would now be estimated for the first covariate (otherwise PLINK would have incorrectly treated the first covariate as an ordinal/ratio measure).
5 进行因子协变量GWAS分析LM模型
「代码：」
代码语言：javascript复制plink --file b --pheno phe.txt --allow-no-sex --linear --covar plink.cov --out re --hide-covar
「日志：」
代码语言：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 plink.cov
--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: 2 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
代码语言：javascript复制library(data.table)
geno = fread("c.raw")
geno[1:10,1:10]
phe = fread("phe.txt")
cov = fread("cov.txt")
dd = data.frame(phe$V3,cov$V4,geno[,7:20])
head(dd)
str(dd)
mod_M7 = lm(phe.V3 ~ cov.V4 M7_1,data=dd)
summary(mod_M7)
mod_M9 = lm(phe.V3 ~ cov.V4 M9_1,data=dd);summary(mod_M9)
「M7加上因子协变量结果：」
「如果是作为数值协变量的结果为：」
结果是不一样的。
5. 使用R语言进行结果比较lm plink.cov
结果和上面世代作为因子完全一样。
6. 固定即回归
「所以，怎么理解固定即回归这句话的？」
❝R语言中，所谓的因子，在进行回归分析时，也是将其转化为不通过水平的数字变量进行的分析，所以和你手动转化的虚拟变量结果是一样的。 ❞
