寫(xiě)在前面：本文為自己在學(xué)習(xí)過(guò)程中看過(guò)各類資料后整理做的筆記，純屬方便自己學(xué)習(xí)以及后期回顧，也希望有幸能幫助到和我一樣的小白們，各位大神如果發(fā)現(xiàn)有錯(cuò)誤即使糾正涂滴，我們共同學(xué)習(xí)！特別感謝所有我閱讀過(guò)的各類推文的作者晴音，我會(huì)附上來(lái)源柔纵，以示版權(quán)及尊重！大家快樂(lè)地學(xué)習(xí)吧锤躁！
全基因組關(guān)聯(lián)分析（Genome-Wide Association Study搁料，GWAS）流程
一、準(zhǔn)備plink文件
1、準(zhǔn)備PED文件
PED文件有六列郭计，六列內(nèi)容如下：
Family ID
Individual ID
Paternal ID
Maternal ID
Sex (1=male; 2=female; other=unknown)
Phenotype
PED
文件是空格（空格或制表符）分隔的文件霸琴。
PED文件長(zhǎng)這個(gè)樣：
2、準(zhǔn)備MAP文件
MAP文件有四列昭伸，四列內(nèi)容如下：
chromosome (1-22, X, Y or 0 if unplaced)
rs# or snp identifier
Genetic distance (morgans)
Base-pair position (bp units)
MAP文件長(zhǎng)這個(gè)樣：
3梧乘、生成bed、fam庐杨、bim宋下、文件
輸入命令：
plink --file mydata --out mydata --make-bed
注：plink指的是plink軟件，如果軟件安裝在某個(gè)指定的路徑的話辑莫，前面還要加上路徑，比如安裝在路徑為/your/pathway/的文件夾下罩引，則命令應(yīng)該為“/your/pathway/plink --file mydata --out mydata --make-bed”
mydata指的是1和2生成的PED和MAP文件名各吨，不需要寫(xiě).ped和.map后綴
bed：
bim：
fam：
注：有時(shí)fam文件最后一列diagnosis會(huì)為 -9 或 NA，都表示缺失值
插一步：質(zhì)控步驟（QC）
可以參照這篇推文袁铐，http://www.reibang.com/p/57c2dbda8a86
填補(bǔ)后質(zhì)控（Post-imputation quality control）
千人基因組大概有83 million變異位點(diǎn)揭蜒，經(jīng)過(guò)填補(bǔ)后有許多質(zhì)量不好的位點(diǎn)，需要過(guò)濾掉剔桨。
去除MAF = 0的位點(diǎn)
去除MAF<0.01 和 info>0.3的位點(diǎn)屉更。info值用來(lái)衡量填充位點(diǎn)的質(zhì)量，一般較差的位點(diǎn)info <0.15洒缀，較好的位點(diǎn)info >0.85瑰谜。所以過(guò)濾閾值一般在0.15-0.85之間。對(duì)于同一個(gè)位點(diǎn)來(lái)說(shuō)树绩，MAF值越小萨脑，info值也越小〗确梗可以將MFA值和info值畫(huà)出柱狀圖渤早，找到一個(gè)比較好的閾值進(jìn)行過(guò)濾。
去除缺失率過(guò)多的位點(diǎn)（98%以上）
質(zhì)控后去除被試命令：?
```
plink --bfile test_2kw -remove remove.txt --out 939_test --make-bed --noweb
```
二瘫俊、準(zhǔn)備表型文件（Alternate phenotype files）
一般表型文件為txt格式鹊杖，表型文件有三列，分別為：
Family ID
Individual ID
Phenotype
假如有多種表型扛芽，第一列和第二列還是Family ID骂蓖、Individual ID，第三列及以后的每列都是表型胸哥，例如以下：
Family ID
Individual ID
Phenotype A
Phenotype B
Phenotype C
……
表型文件長(zhǎng)這樣：
三涯竟、準(zhǔn)備協(xié)變量文件（Covariate files）
協(xié)變量文件同表型文件類似，第一列和第二列是Family ID、Individual ID庐船，第三列及以后的每列都是協(xié)變量
Family ID
Individual ID
Covariate A
Covariate B
Covariate C
……
協(xié)變量文件長(zhǎng)這個(gè)樣（這里有三個(gè)協(xié)變量银酬，分別為Sex,Age,temperature）：
四、plink進(jìn)行表型和基因型以及協(xié)變量的關(guān)聯(lián)分析
命令如下：
plink --bfile mydata --linear --pheno pheno.txt --mpheno 1 --covar
covar.txt --covar-number 1,2,3 --out mydata –noweb
生成的文件為mydata.assoc.linear
注：“mydata”mydata文件不需要后綴筐钟，“--mpheno 1”指的是表型文件的第三列（即第一個(gè)表型）
“--covar-number 1,2,3”指的是協(xié)變量文件的第三列揩瞪、第四列、第五列（即第一個(gè)篓冲、第二個(gè)李破、第三個(gè)協(xié)變量）
“--linear”指的是用的連續(xù)型線性回歸，如果表型為二項(xiàng)式（即0壹将、1）類型嗤攻，則用“--logistic”
五、畫(huà)曼哈頓圖
安裝R語(yǔ)言的qqman包诽俯，其中的manhattan()妇菱，即可畫(huà)曼哈頓圖
https://mp.weixin.qq.com/s/3KN0WHccKv5EbLQUCQVqmw
R包畫(huà)圖
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)根染色體上有多
SNPas.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)
還可以使用haploview畫(huà)曼哈頓圖http://www.reibang.com/p/609149db6fab
自己入的小坑： assoc.linear文件可以用notepad++打開(kāi)，但是? .linear不是制表符暴区，需要linux下轉(zhuǎn)換制表符闯团，sed 's/[ ][ ]*/,/g' file.linear > out.csv，sed -i 's/^,//g' out.csv仙粱，加了個(gè)逗號(hào)房交，最后轉(zhuǎn)成.csv文件，還進(jìn)行了cut伐割，cut -f 1,2,3,9 file_name > new_file候味，提取做曼哈頓圖時(shí)需要的那幾列，得到了最終R語(yǔ)言中要做曼哈頓圖的文件隔心。我滴個(gè)天负溪！哎~
最后需要的文件長(zhǎng)這樣，好有成就感
其實(shí)......不用轉(zhuǎn)格式也可以济炎，后來(lái)找到了這個(gè)http://www.reibang.com/p/e914ecb99fcc
直接加載 .linear格式
終于做完了川抡，放一張我自己做的圖，啊须尚，雖然還部分有細(xì)節(jié)問(wèn)題崖堤，但感覺(jué)擁有了全世界，啊~
部分還待解決的問(wèn)題：
1.多個(gè)pheno的問(wèn)題：與單個(gè)pheno跑出來(lái)的數(shù)據(jù)大小是一樣的耐床，格式也完全一樣密幔，p值也只有一列（我以為是多列，每個(gè)pheno一個(gè)p值）撩轰，但他們的p值不一樣胯甩，所以問(wèn)題來(lái)了昧廷，計(jì)算多表型association時(shí)得到的p值是代表什么？多表型一起相關(guān)的結(jié)果偎箫？
2.協(xié)變量問(wèn)題：如果有多個(gè)協(xié)變量木柬，為什么每個(gè)協(xié)變量都會(huì)跑出來(lái)一個(gè)p值？我有50萬(wàn)個(gè)snp淹办，8個(gè)協(xié)變量眉枕，這樣就會(huì)出來(lái)幾百萬(wàn)個(gè)p值，都要加載進(jìn)去做曼哈頓圖嗎怜森？
3.? .assoc.linear文件中最后幾行是什么呀速挑，沒(méi)搞懂
解決問(wèn)題2：
plink manual中是這樣給出的，線性模型的問(wèn)題副硅，我數(shù)學(xué)知識(shí)有限姥宝，講不太清楚，有大趴制＃看到的話可以同通俗的語(yǔ)言給我留言供大家學(xué)習(xí)伶授。總之呢流纹，我們關(guān)心的是第一個(gè)ADD，所以违诗，manual給出了下面的
1.https://www.cnblogs.com/leezx/p/9013615.html
2.http://www.cnblogs.com/chenwenyan/p/6095531.html
3.https://mp.weixin.qq.com/mp/profile_ext?action=home&__biz=MzIwMzQyMTA2NA==&scene=124#wechat_redirect
4.https://mp.weixin.qq.com/s/2Uy9TXbsV267CCnRCa9vQQ
5.GWAS+公共數(shù)據(jù)庫(kù)上nature
https://mp.weixin.qq.com/s/TOtYKabcEkkwDVOorB_cYQ
6.找到的視頻漱凝，不過(guò)沒(méi)有質(zhì)控
