這里介紹的SNP數(shù)據(jù)與單一表型的相關(guān)分析流程
一、準(zhǔn)備plink文件
ped和map文件轉(zhuǎn)換為bed青灼、fam暴心、bim
plink --file mydata --out mydata --make-bed
mydata為ped和map的文件名，不需要加后綴
make-bed表示生成plink可以處理的二進(jìn)制文件
二杂拨、準(zhǔn)備表型文件
表型文件.pheno 前兩列是FID and IID专普，第三列是表型Phenotype。
如果準(zhǔn)備了表型文件弹沽，那么運(yùn)行plink的同時(shí)會(huì)忽略原來(lái) .fam 里的表型脆诉。
假如有多種表型，第一列和第二列還是Family ID贷币、Individual ID击胜，第三列及以后的每列都是表型，例如Phenotype A役纹、Phenotype B偶摔、Phenotype C ...
三、準(zhǔn)備協(xié)變量
提前將協(xié)變量信息寫入到.covar文件促脉，前兩列是FID and IID辰斋，后面的列是協(xié)變量，主要有年齡瘸味、身高宫仗、體重等。
四旁仿、plink做基因型與表型的關(guān)聯(lián)分析
基因型可以是SNP數(shù)據(jù)藕夫，也可以是連續(xù)性變量
表型可以是連續(xù)性變量，也可以是二分類變量。
- 如果表型是二分類變量（0毅贮、1這種類型）办悟，用plink - -logistic
- 如果表型是連續(xù)的數(shù)量表型，用plink--linear滩褥，指的是用的連續(xù)型線性回歸
plink --bfile mydata --linear --pheno pheno.txt --mpheno 1 --covar covar.txt --covar-number 1,2,3 --out mydata –noweb
- mydata為準(zhǔn)備的基因型plink文件名病蛉，不需要后綴
- pheno.txt 為準(zhǔn)備的表型
- mpheno 1指的是表型文件的第三列（即第一個(gè)表型）
- covar-number 1,2,3指的是協(xié)變量文件的第三列、第四列瑰煎、第五列（即第一個(gè)铺然、第二個(gè)、第三個(gè)協(xié)變量）
- out輸出文件名
plink --bfile dopgene_range_snp --noweb --keep 980ID.txt --recode --make-bed --out 980_snp
生成的文件為.assoc.logistic 或者 .assoc.linear 酒甸，文件內(nèi)容見下圖探熔。
主要需要的數(shù)據(jù)就是P值，此值越小代表SNP與表型關(guān)系越顯著烘挫。
可以關(guān)注一下p值最小的位點(diǎn)信息，參考http://www.reibang.com/p/4b03064e9147
五：結(jié)果可視化：做曼哈頓圖
Manhattan plot是把GWAS分析之后所有SNP位點(diǎn)的p-value在整個(gè)基因組上從左到右依次畫出來(lái)柬甥。
更厲害的是饮六！為了可以更加直觀地表達(dá)結(jié)果，通常都會(huì)將p-value轉(zhuǎn)換為-log10(p-value)苛蒲。這樣的話卤橄，基因位點(diǎn)-log10(p-value)在Y軸的高度就對(duì)應(yīng)了與表型性狀或者疾病的關(guān)聯(lián)程度，關(guān)聯(lián)度越強(qiáng)臂外，p-value越低窟扑，在圖中的位置就越高。
GWAS研究中漏健，p-value閾值一般要在10的-6次方甚至10-8次方以下嚎货，也就說(shuō)曼哈頓圖中Y軸大于6甚至大于8的那些SNP位點(diǎn)才是比較值得研究的，不過(guò)也可以根據(jù)數(shù)據(jù)情況調(diào)整閾值蔫浆。
由于連鎖不平衡（LD）關(guān)系的原因殖属，那些在強(qiáng)關(guān)聯(lián)位點(diǎn)周圍的SNP也會(huì)跟著顯示出類似的信號(hào)強(qiáng)度，并依次往兩邊遞減瓦盛。由于這個(gè)原因洗显，我們?cè)诼D圖上就會(huì)看到一個(gè)個(gè)整齊的信號(hào)峰（如上圖紅色部分）。而這些峰所處的位置一般也是整個(gè)研究中真正關(guān)心的地方原环。
- 怎么畫曼哈頓圖呢挠唆？ 當(dāng)然是用R了！
library(qqman)
bleeding <- read.table("../data/bleeding.assoc.logistic", head=TRUE)#讀取之前生成的相關(guān)分析文件
manhattan(bleeding,chr="CHR",bp="BP",p="P",snp="SNP", main = "Manhattan plot: logistic")
jpeg("bleeding_manhattan.jpeg")
六：結(jié)果可視化：QQ-plot
QQ-plot的縱軸是SNP位點(diǎn)的p-value值（這是實(shí)際得到的結(jié)果嘱吗，observed）玄组，與曼哈頓圖一樣也是表示為 -log10(p-value)；橫軸是則是均勻分布的概率值（這是Expecte的結(jié)果），同樣也是換算為-log10巧勤。
得到QQ圖之后嵌灰，就可以判斷應(yīng)該用是否與表型性狀相關(guān)，也就是GWAS結(jié)果好還是不好颅悉。P值越小的地方偏離直線越多沽瞭，區(qū)別于預(yù)設(shè)的模型（隨機(jī)分布），才證明我們挑出的關(guān)聯(lián)的SNP越可靠剩瓶。
附：keep命令提取表型
舉例驹溃，如果自己的bfile中有1000個(gè)樣本，然而自己做關(guān)聯(lián)分析只需要其中的980個(gè)延曙，另外20個(gè)被剔除豌鹤，這個(gè)時(shí)候就需要用keep命令提取我們所需要的數(shù)據(jù)。
plink中的介紹為：
- bfile為提取數(shù)據(jù)的源頭文件
- keep為要保留的ID枝缔，這里的文件需要自己準(zhǔn)備：
--keep accepts a space/tab-delimited text file with family IDs in the first column and within-family IDs in the second column
說(shuō)明這個(gè)文件必須是要有兩行：第一行family IDs布疙，第2行within-family IDs（與fam文件類似） - out定義輸出文件名
