作為一個(gè)統(tǒng)計(jì)遺傳學(xué)實(shí)驗(yàn)室里的學(xué)生，怎么可以不會(huì)GWAS分析，雖然學(xué)的是生物信息學(xué)窑眯，但是每天聽師兄師姐在那里討論這個(gè)模型砚作，那個(gè)矩陣啥的，多多少少有點(diǎn)印象，雖然不會(huì)推導(dǎo)公式吧，用用軟件總應(yīng)該學(xué)會(huì)，所以我決定考試學(xué)習(xí)GWAS分析黍氮。
這個(gè)過程我要倒著來，假如說我已經(jīng)拿到了每個(gè)snp位點(diǎn)的P值，下一步就是畫曼哈頓圖沫浆，還記得第一次看到曼哈頓圖觉壶，感覺很是高大上。 后來師弟和我說件缸，只要一個(gè)包就可以畫出來铜靶，這個(gè)R包就是qqman.
第一步，在R中安裝qqman這個(gè)R包：
install.packages("qqman")
第二步他炊，查看學(xué)習(xí)手冊(cè)
??qqman
基本自學(xué)能力強(qiáng)的人看這個(gè)學(xué)習(xí)手冊(cè)就能學(xué)會(huì)争剿，這個(gè)包還是不難的
建議這個(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è)是沒根染色體上有多少SNP
as.data.frame(table(data2$CHR)) # 這個(gè)是沒根染色體上有多少SNP
qq(gwasResults$P) # 畫qq圖
qq(data2$P) # 畫qq圖
manhattan(gwasResults, annotatePval = 0.01) # 這個(gè)可以對(duì)每根染色體上最高的那個(gè)點(diǎn)注釋出來
腳本里面記錄我覺得比較重要的幾條命令
我這里來詳細(xì)介紹一個(gè)
如果我們啥也不設(shè)置，我感覺圖片有點(diǎn)難看的
manhattan(gwasResults)
我們來加點(diǎn)彩色的
color_set <- rainbow(22)
manhattan(gwasResults,col = color_set)
但是似乎這個(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)
然后我標(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里面千萬不要加入中文的逗號(hào)，不然程序運(yùn)行有問題炒刁，你還發(fā)現(xiàn)不了
