这里介绍的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值最小的位点信息,参考//www.greatytc.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越可靠。
附:keep命令提取表型
举例,如果自己的bfile中有1000个样本,然而自己做关联分析只需要其中的980个,另外20个被剔除,这个时候就需要用keep命令提取我们所需要的数据。
plink中的介绍为:
- bfile为提取数据的源头文件
- keep为要保留的ID,这里的文件需要自己准备:
--keep accepts a space/tab-delimited text file with family IDs in the first column and within-family IDs in the second column
说明这个文件必须是要有两行:第一行family IDs,第2行within-family IDs(与fam文件类似) - out定义输出文件名