r语言怎么实现LM模型+数值+因子协变量

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

「注意：」这里的性别虽然是因子，但是其只有两个水平，也可以将作为连续的变量，计算方法是一样的。如果是三个水平的因子，就不能直接转化为变量了。

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

「代码：」

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

「日志：」

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
  --hide-covar
  --linear
  --out re
  --pheno phe.txt

Note: --hide-covar flag deprecated.  Use e.g. '--linear hide-covar'.
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: 3 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

library(data.table)
geno = fread("c.raw")
geno[1:10,1:10]
phe = fread("phe.txt")
cov = fread("cov.txt")
plink = fread("plink.cov")
dd = data.frame(phe = phe$V3,cov1 = plink$COV1,cov2 = plink$COV2_4,cov3=plink$COV2_5,geno[,7:20])
head(dd)
mod_M7 = lm(phe ~ cov1+cov2+cov3 + M7_1,data=dd);summary(mod_M7)

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

「这里，我们可以测试一下：」将性别由数字，变为因子，可以发现结果是一样的：