全基因组关联分析（Genome-wide association study，GWAS）是应用基因组中数以百万级别的单核苷酸多样性（Single Nucleotide Ploymorphism，SNP）作为分子遗传标记，进行全基因组水平上表型与marker的相关性分析，从而发现影响复杂性状的基因突变的一种策略。
使用的软件为：
- plink
- EMMAX
- R
平台为Linux
00软件安装
Plink作为一款老牌软件，它的文件格式是多个计算软件所需要的。Plink自身也可以进行全基因组的关联分析，但是速度较慢不推荐。我们本次使用Plink的目的是利用它对SNP数据进行提取和过滤。Plink的安装也十分简单，可以直接用conda安装：
conda install -y plink
EMMAX是2010年Kang开发的EMMA升级版，它的计算速度特别快，并且在棉花，大豆，水稻甚至椰子等复杂性状探究中广泛应用。每一个SNP对复杂性状的解释率都很低，每个组件的方差在运算中都只计算一次，从而达到简化计算的目的。
软件下载地址：
EMMAX的下载地址非常友好，不需要翻墙。
EMMAX下载
我用的是老版本，应该没什么区别，看介绍只有编译上的差别，然后我们安装
tar xvf emmax-interl-binary-20120210.tar.gz
#然后加到环境变量
emmax-interl64
R语言大家应该熟知，不需要过度介绍
01phenotype和genotype文件准备
phenotype文件
emmax识别的phenotype文件需要一个表型一个单独的文本，我给大家截图展示一下：
phenotype
phenotype格式
文件的格式也非常简单，只有三列。前两列为品种编号和genotype中对应，第三列为表型。
genotype文件
plink格式的SNP文件需要用plink进行抽取和转置：
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 二进制plink文件
--keep 需要提取的品种集
--make-bed 生成二进制的plink文件
--maf 过滤最大等基因频率
--geno 过滤缺失率
--out 输出的文件名
02关联分析
亲缘关系矩阵
emmax有子命令可以计算亲缘关系矩阵，也可以用plink去计算。
emmax-kin-intel64 emmax_in -v -h -s -d 10 new_all_chr.filter
最后的new_all_chr.filter是genotype名
群体结构
群体结构对于群体分化不明显的数据，我觉得应该是可以用加的。群体结构我们使用Admixture软件进行计算，这个很多网文都有就不过多介绍，用conda是可以直接安装的。计算之后得到最小的CV值，然后把群体结构整理成协变量的格式。这里我使用的数据是亲缘关系较近的群体，就没有加群体结构。
计算关联性
计算关联性之前需要将我们得到的plink文件转置，emmax接受的是转置后的tped文件，代码如下：
plink --bfile gwas_all_chr.filter --recode12 --output-missing-genotype 0 --transpose --out new_all_chr.filter
接着是计算关联分析的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是我的表型文件的文件名，为了循环跑每一个sample
-t 是tped的文件名
-p 是phenotype文件名
-k 是刚才得到的kinship文件，如果要用群体结构的话可以加一个协方差的参数
-o 是输出文件的文件名
-v 我记得是输出当前进度的设置
-d 我记得是给予的线程数（ps：这两个参数记不太清了，反正没影响）
03结果文件解析
跑完上面的程序后，我们会得到三个文件分别是log，reml和ps文件。其中ps文件使我们需要的结果文件。目前手上的项目都已经结束，所以搬运一张ps文件的格式：
ps文件格式
这个文件四列分别为SNP_ID，回归系数，标准差和p值。
根据这个文件，我们把数据整理成SNP，CHR，BP（snp所在染色体的物理距离）和P（关联分析的P值）四列，然后用于下一步画图。
04结果画图
比较简单的曼哈顿图就是使用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是当前文件夹下所有的txt文件，这个地方我已经把ps文件进行抽提整理成上述所说的那四列的形式，在这里是为了批量画图。
最终结果
这里解释一下，为什么0-1的位置没有东西。这些点都是snp，0-1之间的snp的p值为0.1-1，也就是完全不显著的snp，为了节省画图时间和图片空间，我在ps到txt的文件整理过程就给过滤了。第二点，为什么有两个阈值线，这里分别代表显著和极显著。
好了，这就是GWAS工作的全部内容，有补充的小伙伴请在评论区留言。谢谢~
网友评论