作為一個(gè)統(tǒng)計(jì)遺傳學(xué)實(shí)驗(yàn)室里的學(xué)生，怎么可以不會(huì)GWAS分析，雖然學(xué)的是生物信息學(xué)，但是每天聽(tīng)?zhēng)熜謳熃阍谀抢镉懻撨@個(gè)模型，那個(gè)矩陣啥的，多多少少有點(diǎn)印象，雖然不會(huì)推導(dǎo)公式吧，用用軟件總應(yīng)該學(xué)會(huì)，所以我決定考試學(xué)習(xí)GWAS分析。

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

第一步，在R中安裝qqman這個(gè)R包：

install.packages("qqman")

第二步，查看學(xué)習(xí)手冊(cè)

??qqman

基本自學(xué)能力強(qiáng)的人看這個(gè)學(xué)習(xí)手冊(cè)就能學(xué)會(huì)，這個(gè)包還是不難的

help

建議這個(gè)學(xué)習(xí)步驟走一遍，就會(huì)了

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

setwd('/Users/mac/Desktop/123')#  設(shè)置工作目錄
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) #設(shè)置點(diǎn)的大小
color_set <- rainbow(9) # 設(shè)置顏色集合 建議c("#8DA0CB","#E78AC3","#A6D854","#FFD92F","#E5C494","#66C2A5","#FC8D62")

svg(file="manpic.svg", width=12, height=8)# 保存svg格式的圖片 設(shè)置名字 
#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() 顯示當(dāng)前圖像參數(shù)


str(gwasResults)  #zscore beita 值除以standard error 這個(gè)值越大 P越小
head(gwasResults) # 看前面幾行
tail(data2) #看后面幾行
as.data.frame(table(gwasResults$CHR))# 這個(gè)是沒(méi)根染色體上有多少SNP
as.data.frame(table(data2$CHR)) # 這個(gè)是沒(méi)根染色體上有多少SNP

qq(gwasResults$P) # 畫(huà)qq圖
qq(data2$P) # 畫(huà)qq圖

manhattan(gwasResults, annotatePval = 0.01) # 這個(gè)可以對(duì)每根染色體上最高的那個(gè)點(diǎn)注釋出來(lái)

腳本里面記錄我覺(jué)得比較重要的幾條命令

我這里來(lái)詳細(xì)介紹一個(gè)
如果我們啥也不設(shè)置，我感覺(jué)圖片有點(diǎn)難看的

manhattan(gwasResults)

raw pic

我們來(lái)加點(diǎn)彩色的

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

有點(diǎn)迷的顏色

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

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

感覺(jué)柔和了一點(diǎn)

然后我標(biāo)記一下最高點(diǎn)

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

注釋一下

好啦，差不多啦，功能夠用就行了，如果要再個(gè)性化，建議看一個(gè)這個(gè)R包的源碼

注： R里面千萬(wàn)不要加入中文的逗號(hào)，不然程序運(yùn)行有問(wèn)題，你還發(fā)現(xiàn)不了