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加上因子協變量結果：」

「這里，我們可以測試一下：」將性別由數字，變為因子，可以發現結果是一樣的：