全基因組關(guān)聯(lián)分析（Genome-wide association study，GWAS）是應(yīng)用基因組中數(shù)以百萬級別的單核苷酸多樣性（Single Nucleotide Ploymorphism，SNP）作為分子遺傳標(biāo)記，進(jìn)行全基因組水平上表型與marker的相關(guān)性分析，從而發(fā)現(xiàn)影響復(fù)雜性狀的基因突變的一種策略。

使用的軟件為：

• plink
• EMMAX
• R
  平臺為Linux

00軟件安裝

Plink作為一款老牌軟件，它的文件格式是多個(gè)計(jì)算軟件所需要的。Plink自身也可以進(jìn)行全基因組的關(guān)聯(lián)分析，但是速度較慢不推薦。我們本次使用Plink的目的是利用它對SNP數(shù)據(jù)進(jìn)行提取和過濾。Plink的安裝也十分簡單，可以直接用conda安裝：
conda install -y plink
EMMAX是2010年Kang開發(fā)的EMMA升級版，它的計(jì)算速度特別快，并且在棉花，大豆，水稻甚至椰子等復(fù)雜性狀探究中廣泛應(yīng)用。每一個(gè)SNP對復(fù)雜性狀的解釋率都很低，每個(gè)組件的方差在運(yùn)算中都只計(jì)算一次，從而達(dá)到簡化計(jì)算的目的。
軟件下載地址：
EMMAX的下載地址非常友好，不需要翻墻。

EMMAX下載

tar xvf emmax-interl-binary-20120210.tar.gz
#然后加到環(huán)境變量
emmax-interl64

R語言大家應(yīng)該熟知，不需要過度介紹

01phenotype和genotype文件準(zhǔn)備

phenotype文件

phenotype

phenotype格式

文件的格式也非常簡單，只有三列。前兩列為品種編號和genotype中對應(yīng)，第三列為表型。

genotype文件

plink格式的SNP文件需要用plink進(jìn)行抽取和轉(zhuǎn)置：
plink --noweb --bfile /public/home/hguo/public/data/SNP/GWAS_Analysis/Rice_GWAS/gwas_all_chr.filter --keep keep_individual.txt --make-bed --maf 0.05 --geno 0.1 --out gwas_all_chr.filter
--bfile 二進(jìn)制plink文件
--keep 需要提取的品種集
--make-bed 生成二進(jìn)制的plink文件
--maf 過濾最大等基因頻率
--geno 過濾缺失率
--out 輸出的文件名

02關(guān)聯(lián)分析

親緣關(guān)系矩陣

emmax有子命令可以計(jì)算親緣關(guān)系矩陣，也可以用plink去計(jì)算。
emmax-kin-intel64 emmax_in -v -h -s -d 10 new_all_chr.filter
最后的new_all_chr.filter是genotype名

群體結(jié)構(gòu)

群體結(jié)構(gòu)對于群體分化不明顯的數(shù)據(jù)，我覺得應(yīng)該是可以用加的。群體結(jié)構(gòu)我們使用Admixture軟件進(jìn)行計(jì)算，這個(gè)很多網(wǎng)文都有就不過多介紹，用conda是可以直接安裝的。計(jì)算之后得到最小的CV值，然后把群體結(jié)構(gòu)整理成協(xié)變量的格式。這里我使用的數(shù)據(jù)是親緣關(guān)系較近的群體，就沒有加群體結(jié)構(gòu)。

計(jì)算關(guān)聯(lián)性

計(jì)算關(guān)聯(lián)性之前需要將我們得到的plink文件轉(zhuǎn)置，emmax接受的是轉(zhuǎn)置后的tped文件，代碼如下：

plink --bfile gwas_all_chr.filter --recode12 --output-missing-genotype 0 --transpose --out new_all_chr.filter

接著是計(jì)算關(guān)聯(lián)分析的p值。這一步很簡單，直接上代碼：

for line in $(cat name.txt)
do
 emmax -v -d 10 -t new_all_chr.filter -p pheno/$line -k new_all_chr.filter.hIBS.kinf  -o result/$line
done

解釋一下，name.txt是我的表型文件的文件名，為了循環(huán)跑每一個(gè)sample
-t 是tped的文件名
-p 是phenotype文件名
-k 是剛才得到的kinship文件，如果要用群體結(jié)構(gòu)的話可以加一個(gè)協(xié)方差的參數(shù)
-o 是輸出文件的文件名
-v 我記得是輸出當(dāng)前進(jìn)度的設(shè)置
-d 我記得是給予的線程數(shù)（ps：這兩個(gè)參數(shù)記不太清了，反正沒影響）

03結(jié)果文件解析

ps文件格式

這個(gè)文件四列分別為SNP_ID，回歸系數(shù)，標(biāo)準(zhǔn)差和p值。
根據(jù)這個(gè)文件，我們把數(shù)據(jù)整理成SNP，CHR，BP（snp所在染色體的物理距離）和P（關(guān)聯(lián)分析的P值）四列，然后用于下一步畫圖。

04結(jié)果畫圖

比較簡單的曼哈頓圖就是使用R包qqman，只要把文件格式整理成我上述描述的那樣可以直接使用下面的代碼：

temp<-list.files(pattern = "*.txt")
for (i in 1:length(temp)) {
library("qqman")
mydata<-read.table(temp[i],header=TRUE,sep="\t")
mydata$SNP<-as.character(mydata$SNP)
mydata$CHR<-as.numeric(mydata$CHR)
mydata1<-na.omit(mydata)

outfile<-paste("MM/",temp[i],".png",sep = "")
png(outfile,height = 900,width = 1200,pointsize = 12)
manhattan(mydata1,main="Manhattan Plot",cex=1.5,cex.axis=1,col=c("red","blue"),suggestiveline =F,
genomewideline=F,,chrlabs = c("1:12"))
abline(h=-log10(1.33E-6), col="orange",lty=2,lwd=2)
abline(h=-log10(6.66E-8), col="purple",lty=2,lwd=2)
dev.off()
}

解釋一下，*txt是當(dāng)前文件夾下所有的txt文件，這個(gè)地方我已經(jīng)把ps文件進(jìn)行抽提整理成上述所說的那四列的形式，在這里是為了批量畫圖。

最終結(jié)果

好了，這就是GWAS工作的全部內(nèi)容，有補(bǔ)充的小伙伴請?jiān)谠u論區(qū)留言。謝謝~