群体遗传学生物信息

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

2020-06-21  本文已影响0人  SnowPye

这里介绍的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数据,也可以是连续性变量
表型可以是连续性变量,也可以是二分类变量。

plink --bfile mydata --linear --pheno pheno.txt --mpheno 1 --covar covar.txt --covar-number 1,2,3 --out mydata –noweb
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也会跟着显示出类似的信号强度,并依次往两边递减。由于这个原因,我们在曼哈顿图上就会看到一个个整齐的信号峰(如上图红色部分)。而这些峰所处的位置一般也是整个研究中真正关心的地方。

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越可靠

附:keep命令提取表型

举例,如果自己的bfile中有1000个样本,然而自己做关联分析只需要其中的980个,另外20个被剔除,这个时候就需要用keep命令提取我们所需要的数据。
plink中的介绍为:


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

上一篇下一篇

猜你喜欢

热点阅读