利用PLINK进行GWAS分析（五）

这一篇正式介绍全基因组关联分析的PLINK的commands。

#case-control 试验设计对绵羊的有角和无角这一对性状进行研究。假定有角为case性状，则无角为control性状。在一个绵羊群体中，存在一部分个体为有角表型，另一部分个体为无角表型。现在对这一群绵羊共m只进行全基因组的SNP检测，得到全基因组的SNP数据。 整理绵羊的SNP数据与个体表型记录数据为PLINK的输入文件格式，如sheep.ped和sheep.fam文件。 1、最简单的关联分析

#对数据进行质控，得到质控后的数据sheep_qc。改写为二进制格式为了提高读取速度。
plink --file sheep --maf 0.01 --geno 0.1 --mind 0.1 --hwe 1e-5 --make-bed --out sheep_qc
#关联分析过程用最简单的assoc。输出result1.assoc文件。
plink --bfile sheep_qc --assoc --out result1

result1.assoc文件可以用文本打开，每一行是一个位点的信息，每一列的信息如下。

column1: CHR, chromosome
column2: SNP, SNP ID
column3: BP, Physical position (base-pair)
column4: A1, Minor allele name (based on whole sample)
column5: F_A, Frequency of this allele in cases
column6: F_U, Frequency of this allele in controls
column7: A2, Major allele name
column8: CHISQ, Basic allelic test chi-square (1df)
column9: P, Asympototic p-value for this test
column10: Estimated odds ratio (for A1, i.e. A2 is reference)

第八列是SNP位点的allelic test的卡方值，第9列是卡方检验的P值。一般关注第八列的数值，对其做log变换，绘图。选取P值大于基因组或染色体显著水平的SNP，这些SNPs即是与形状关联的SNP位点。

2、Fisher exact test

plink --bfile sheep_qc --fisher --out result2

与上一个分析的差别仅在于对于统计量的检验换做了Fisher检验。 result1.assoc文件可以用文本打开，每一行是一个位点的信息，每一列的信息如下。

column1: CHR, chromosome
column2: SNP, SNP ID
column3: BP, Physical position (base-pair)
column4: A1, Minor allele name (based on whole sample)
column5: F_A, Frequency of this allele in cases
column6: F_U, Frequency of this allele in controls
column7: A2, Major allele name
column8: P, Exact p-value for this test
column9: Estimated odds ratio (for A1)

第八列为fisher test的检验P值。同样是选取P值大于基因组或染色体显著水平的SNP，作为关联分析得到的显著位点。

3、model test

包括：Cochran-Armitage trend test(1df)、 Genotypic test(2df)、 Dominant gene action test(1df)、 Recessive gene action test(1df)。

plink --file sheep_qc --model --out result3
#输出result3.model文件
#也可以用fisher test的全模型检验
plink --file sheep_qc --model --fisher

每列信息如下：

column1: CHR, chromosome
column2: SNP, SNP ID
column3: TEST, type of test
column4: AFF, Genotypes/alleles in cases
column5: UNAFF, Genotypes/alleles in controls
column6: CHISQ, Chi-squated statistic
column7: DF, Degrees of freedom for test
column8: P, Asymptotic p-value for this test

4、logistic models

plink --file sheep_qc --logistic

5、Covariate adjust

#先生成MDS文件，plink.mds
plink --file sheep_qc --cluster --mds-plot 10
#应用logistic回归，加上MDS的C1和C2组分
plink --file sheep_qc --logistic --covar plink.mds --covar-name C1-C2 --out result4

6、计算Inflation λ 引入--adjust命令时，log文件中会给出λ值。

plink --file sheep_qc --assoc --adjust
#log文件中：
#Genomic inflation factor (based on median chi-squared) is ***
#Mean chi-squared statistic is ***

λ=1比较合适，如果λ＞1，需要慎重考虑一下群体结构的问题，而不仅仅是对全体P值除以λ。 #case-control设计的GWAS小结不考虑群体结构问题，一般上述1、2、3、4命令选一执行。 考虑群体结果，执行--cluster --mds-plot 10得到群体结构的mds文件，按5命令执行。 另外，对λ偏大的情况，需要慎重考虑一下群体结构问题对关联分析结果的假阳性和假阴性的影响。

GWAS中群体结构问题比较复杂，以后另写文章详述。