2018-10-25 GWAS实战（一） qqman绘制曼哈顿图

作为一个统计遗传学实验室里的学生，怎么可以不会GWAS分析，虽然学的是生物信息学，但是每天听师兄师姐在那里讨论这个模型，那个矩阵啥的，多多少少有点印象，虽然不会推导公式吧，用用软件总应该学会，所以我决定考试学习GWAS分析。

这个过程我要倒着来，假如说我已经拿到了每个snp位点的P值，下一步就是画曼哈顿图，还记得第一次看到曼哈顿图，感觉很是高大上。后来师弟和我说，只要一个包就可以画出来，这个R包就是qqman.

第一步，在R中安装qqman这个R包：

install.packages("qqman")

第二步，查看学习手册

??qqman

基本自学能力强的人看这个学习手册就能学会，这个包还是不难的

help

建议这个学习步骤走一遍，就会了

第三步，仿照脚本，看里面的注释内容自己理解

setwd('/Users/mac/Desktop/123')#  设置工作目录
library(qqman) # 载入包
data <- read.table("5filter_result.assoc.linear",header = TRUE) #读取数据
data1 <- data[,c(1,2,3,9)] #按照规则截取列
data2 <- na.omit(data1) # 删除含有NA的整行
par(cex=0.8) #设置点的大小
color_set <- rainbow(9) # 设置颜色集合 建议c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")

svg(file="manpic.svg", width=12, height=8)# 保存svg格式的图片 设置名字 
#manhattan(data2,main="Manhattan Plot",col = color_set) #suggestiveline = FALSE 更加显著
manhattan(data2,main="Manhattan Plot",col = c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62"),suggestiveline = FALSE,annotatePval = 0.01) #suggestiveline = FALSE 更加显著
dev.off() # 保存图片

#par() 显示当前图像参数


str(gwasResults)  #zscore beita 值除以standard error 这个值越大 P越小
head(gwasResults) # 看前面几行
tail(data2) #看后面几行
as.data.frame(table(gwasResults$CHR))# 这个是没根染色体上有多少SNP
as.data.frame(table(data2$CHR)) # 这个是没根染色体上有多少SNP

qq(gwasResults$P) # 画qq图
qq(data2$P) # 画qq图

manhattan(gwasResults, annotatePval = 0.01) # 这个可以对每根染色体上最高的那个点注释出来

脚本里面记录我觉得比较重要的几条命令

我这里来详细介绍一个
如果我们啥也不设置，我感觉图片有点难看的

manhattan(gwasResults)

raw pic

我们来加点彩色的

color_set <- rainbow(22)
manhattan(gwasResults,col = color_set)

有点迷的颜色

但是似乎这个颜色有点丑哇，所以建议使用我师弟给我的参数

color_set <- c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")
par(cex=0.8)
manhattan(gwasResults,col = color_set)

感觉柔和了一点

然后我标记一下最高点

color_set <- c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")
par(cex=0.8)
manhattan(gwasResults,col = color_set,annotatePval = 0.01)

注释一下

好啦，差不多啦，功能够用就行了，如果要再个性化，建议看一个这个R包的源码

注： R里面千万不要加入中文的逗号，不然程序运行有问题，你还发现不了

最后编辑于：2018.10.25 09:27:54