GWAS - 基因型与表型的关联分析流程

这里介绍的SNP数据与单一表型的相关分析流程

一、准备plink文件

ped和map文件转换为bed、fam、bim

plink --file mydata --out mydata --make-bed

mydata为ped和map的文件名，不需要加后缀
make-bed表示生成plink可以处理的二进制文件

二、准备表型文件

表型文件.pheno 前两列是FID and IID，第三列是表型Phenotype。
如果准备了表型文件，那么运行plink的同时会忽略原来 .fam 里的表型。
假如有多种表型，第一列和第二列还是Family ID、Individual ID，第三列及以后的每列都是表型，例如Phenotype A、Phenotype B、Phenotype C …

三、准备协变量

提前将协变量信息写入到.covar文件，前两列是FID and IID，后面的列是协变量，主要有年龄、身高、体重等。

四、plink做基因型与表型的关联分析

基因型可以是SNP数据，也可以是连续性变量
表型可以是连续性变量，也可以是二分类变量。

• 如果表型是二分类变量（0、1这种类型），用plink - -logistic
• 如果表型是连续的数量表型，用plink–linear，指的是用的连续型线性回归

plink --bfile mydata --linear --pheno pheno.txt --mpheno 1 --covar covar.txt --covar-number 1,2,3 --out mydata –noweb

• mydata为准备的基因型plink文件名，不需要后缀
• pheno.txt 为准备的表型
• mpheno 1指的是表型文件的第三列（即第一个表型）
• covar-number 1,2,3指的是协变量文件的第三列、第四列、第五列（即第一个、第二个、第三个协变量）
• out输出文件名

plink --bfile dopgene_range_snp --noweb --keep 980ID.txt --recode --make-bed --out 980_snp

生成的文件为.assoc.logistic 或者 .assoc.linear ，文件内容见下图。
主要需要的数据就是P值，此值越小代表SNP与表型关系越显著。

可以关注一下p值最小的位点信息，参考https://www.jianshu.com/p/4b03064e9147

五：结果可视化：做曼哈顿图

Manhattan plot是把GWAS分析之后所有SNP位点的p-value在整个基因组上从左到右依次画出来。

更厉害的是！为了可以更加直观地表达结果，通常都会将p-value转换为-log10(p-value)。这样的话，基因位点-log10(p-value)在Y轴的高度就对应了与表型性状或者疾病的关联程度，关联度越强，p-value越低，在图中的位置就越高。

GWAS研究中，p-value阈值一般要在10的-6次方甚至10-8次方以下，也就说曼哈顿图中Y轴大于6甚至大于8的那些SNP位点才是比较值得研究的，不过也可以根据数据情况调整阈值。

由于连锁不平衡（LD）关系的原因，那些在强关联位点周围的SNP也会跟着显示出类似的信号强度，并依次往两边递减。由于这个原因，我们在曼哈顿图上就会看到一个个整齐的信号峰（如上图红色部分）。而这些峰所处的位置一般也是整个研究中真正关心的地方。

• 怎么画曼哈顿图呢？ 当然是用R了！

library(qqman)
bleeding <- read.table("../data/bleeding.assoc.logistic", head=TRUE)#读取之前生成的相关分析文件
manhattan(bleeding,chr="CHR",bp="BP",p="P",snp="SNP", main = "Manhattan plot: logistic")
jpeg("bleeding_manhattan.jpeg")

六：结果可视化：QQ-plot

QQ-plot的纵轴是SNP位点的p-value值（这是实际得到的结果，observed），与曼哈顿图一样也是表示为 -log10(p-value)；横轴是则是均匀分布的概率值（这是Expecte的结果），同样也是换算为-log10。
得到QQ图之后，就可以判断应该用是否与表型性状相关，也就是GWAS结果好还是不好。P值越小的地方偏离直线越多，区别于预设的模型（随机分布），才证明我们挑出的关联的SNP越可靠。

参考：
https://www.jianshu.com/p/987859ae503c