作为一个统计遗传学实验室里的学生，怎么可以不会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里面千万不要加入中文的逗号，不然程序运行有问题，你还发现不了
网友评论