作為一個統(tǒng)計遺傳學實驗室里的學生，怎么可以不會GWAS分析，雖然學的是生物信息學，但是每天聽師兄師姐在那里討論這個模型，那個矩陣啥的，多多少少有點印象，雖然不會推導公式吧，用用軟件總應該學會，所以我決定考試學習GWAS分析。

這個過程我要倒著來，假如說我已經(jīng)拿到了每個snp位點的P值，下一步就是畫曼哈頓圖，還記得第一次看到曼哈頓圖，感覺很是高大上。 后來師弟和我說，只要一個包就可以畫出來，這個R包就是qqman.

第一步，在R中安裝qqman這個R包：

install.packages("qqman")

第二步，查看學習手冊

??qqman

基本自學能力強的人看這個學習手冊就能學會，這個包還是不難的

help

建議這個學習步驟走一遍，就會了

第三步，仿照腳本，看里面的注釋內(nèi)容自己理解

setwd('/Users/mac/Desktop/123')#  設置工作目錄
library(qqman) # 載入包
data <- read.table("5filter_result.assoc.linear",header = TRUE) #讀取數(shù)據(jù)
data1 <- data[,c(1,2,3,9)] #按照規(guī)則截取列
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() 顯示當前圖像參數(shù)


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)

有點迷的顏色

但是似乎這個顏色有點丑哇，所以建議使用我?guī)煹芙o我的參數(shù)

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里面千萬不要加入中文的逗號，不然程序運行有問題，你還發(fā)現(xiàn)不了