使用EMMAX做GWAS分析

EMMAX是一款非常高效和使用简单的GWAS分析软件，使用者可以进行模型设置，比如一般线性模型和混合线性模型，操作简单，只需要准备几个必须的文件即可。本博客记录本人在使用改软件的一些步骤。由于使用的基因型数据是已经过滤好的数据，我们在此就不在进行数据指控方面的操作；GWAS基因型质控的一款非常强大的软件是PLINK，可以设置多个参数来对样本和SNP进行过滤。使用EMMAX需注意几个方面：（1）基因型数据要求是PLINK的转置格式，在群体结构控制的文件中，需要多加一列1，表示截距值。群体结构控制的文件可以采用STRUCTURE或者ADMIXTURE，PCA软件的输出值，值得注意的是，在用STRUCTURE和ADMIXTURE软件输出的文件中，协变量的要减少最后一个变量，比如STRUCTURE分析后，得到群体的亚群个数是3，那么每个样品对会对应3个相应的系数，以此来表示这个个体属于这三个假定的祖先亚群的概率；3个系数的总和是1；在用这个协变量文件的时候，要删除最后一列，并在文件的第三列写1，表示截距。这样就可以进行EMMAX分析了。

（1） 准备基因型数据，这部分需要在PLINK软件中操作，本分析中，我从916分样本的数据中，提取本项目分析的320个个体的基因型数据，--keep参数后边接需要提取的样品名字，格式如下：-9表示样本的family ID是缺失的，后边的是每个材料的名字。最后输出的是三个带有后缀的文件，如：genotype_millet320.tfam， genotype_millet320.tped，genotype_millet320.map

-9	Ci007
-9	Ci012
-9	Ci013
-9	Ci021
-9	Ci026
-9	Ci032
-9	Ci046

命令如下：

plink --noweb --bfile ~/GWAS_foxmilt/plink/foxtail916 --keep sample_names.txt --nonfounders --recode12 --output-missing-genotype 0 --transpose --out genotype_millet320  # --transpose 表示输出数据转置，这个是PLINK软件的一种特殊的格式。

（2）计算样品的亲缘关系矩阵，这个工作可以直接用EMMAX软件自带的来计算，根据PLINK输出的基因型文件进行，文件输入的后缀名就是前一个步骤输出的后缀名，emmax-kin计算后会生成一个大的矩阵文件，用于后续MLM分析。

# 使用EMMAX自带的软件计算亲缘关系矩阵
~/software/emmax-beta-07Mar2010/emmax-kin -v -h -s -d 10 genotype_millet320

（3）计算协变量文件，可以用STRUCTURE分析后的协变量文件，也可以用主成分分析的协变量，这里我们用PLINK 的PCA来进行分析，并用前三个主成分作为协变量，计算命令如下：

# 使用PLINK软件计算群体结构的协变量，这里选择前三个主成分作为协变量，具体需要多少个协变量需要根据你的群体结构来进行判断，在本项目中，基因型数据的群体分层并不明显，选择前三个主成分就可以控制群体结构了
plink --noweb --tfile genotype_millet320 --pca 3 --out pca_320

协变量文件可根据上面的PCA分析结果进行准备，用AWK方便操作

#输出的协变量文件进行格式整理，这里第三列是截距值，前两列是family ID和 individual ID， 后边的三列是三个主成分
awk '{print $1"\t"$2"\t"1"\t"$3"\t"$4"\t"$5"\t"}' pca_320.eigenvec > pca_structure_320.txt # 展示如下：
-9	Ci007	1	-0.0540318	0.0158467	-0.0271305	
-9	Ci012	1	-0.0457195	0.0141429	-0.126416	
-9	Ci013	1	-0.0502624	0.00835785	-0.133798	
-9	Ci021	1	-0.055962	-0.00897986	0.0241683	
-9	Ci026	1	-0.0641019	-0.0368972	0.0350536	
-9	Ci032	1	-0.0508938	0.0761431	0.0161454	
-9	Ci046	1	-0.0531549	0.0864745	0.0476195	
-9	Ci053	1	-0.0547984	0.0798894	0.0657235	

（4）准备表型文件。EMMAX软件要求每一个表型一个单独的文件，所以这里对每个表型都准备一个单独的文件，如果表型很多，就可以在R中一次性操作；表型文件有三列，分别是family ID，individual ID和phenotype，三列都不需要表头，如下：

-9	Ci007	0.44
-9	Ci012	0.626666666666667
-9	Ci013	0.44
-9	Ci021	0.62
-9	Ci026	0.56
-9	Ci032	0.753333333333333
-9	Ci046	0.62
-9	Ci053	0.753333333333333
-9	Ci054	0.793333333333333
-9	Ci055	0.866666666666667

（5）利用MLM模型进行全基因组关联分析。根据上面准备的文件，这里就可以使用混合线性模型进行分析了，代码如下：

~/software/emmax-beta-07Mar2010/emmax -v -d 10 -t genotype_millet320   -p final_data/phenotype_seed_germination_308.txt -k genotype_millet320.hIBS.kinf  -c pca_structure_320.txt  -o seed_germination_320

如果表型文件非常多，就可以写一个shell循环，如下边的分析中有68个表型文件：

ls *.txt | while read id; do trait_name=${id%.*}; ~/software/emmax-beta-07Mar2010/emmax -v -d 10 -t  ../genotype_499_transpose  -p $id -k ../genotype_499_transpose.hIBS.kinf  -c ../pca_structure_pca_genotype_499_transpose.txt  -o $trait_name;done

（6）结果查看和数据可视化。EMMAX输出的文件有三个，分别是一个日志文件，一个关联结果的p值文件(.ps 后缀)和一个包含似然值文件，里边包含样品的遗传力。三个文件中，关联结果文件一般有三列或者四列构成，不同的操作系统下EMMAX软件输出会有差异，如我在Linux中，采用旧版本的就输出三列，分别是标记名称，标准差以及p值，标准差是对标记效应的标准差。需要根据前面基因型的map文件，展示每个标记的分布以及具体的p值，可以用qqman这个软件包进行操作。为了方便，我写了一个R脚本 plot_emmax.R：

#!/usr/bin/Rscript 
library(optparse)
library(qqman)
option_list = list(make_option(c("-m","--map"), help="map file",type='character'),make_option(c("-e","--emmax"), help = "the emmax results",type='character'),make_option(c("-o","--output"), help = "the output name", type='character'),make_option(c("-p","--significant"), help = "p value for significant signals", default = -log10(1e-8), type = 'numeric' )
)opts = parse_args(OptionParser(option_list=option_list))
map_df = read.table(file=opts$map, header = F,sep="\t")
names(map_df) = c('CHR','SNP',"Genetic_distance","BP")
emmax_df = read.table(file=opts$emmax, header = F,sep="\t")
names(emmax_df) = c('SNP','SE',"P")
merge_map_emmax = merge(map_df, emmax_df, by="SNP")
merge_map_emmax_sort = merge_map_emmax[order(merge_map_emmax$CHR, merge_map_emmax$BP),]
out_manhattan = paste(opts$output,'.manhattan.pdf')
pdf(out_manhattan)
manhattan(merge_map_emmax_sort,col = c("blue4", "orange3"),main = "Manhattan Plot", suggestiveline = FALSE, genomewideline = opts$significant)
dev.off()
qq_name = paste(opts$output,'.qq.pdf')
pdf(qq_name)
qq(merge_map_emmax_sort$P)
dev.off()

最后通过一个循环运行每个结果文件，输出曼哈顿图和QQ图：

ls *.ps | while read id; do sample_name=${id%.*};Rscript plot_emmax.R -m ../genotype_499.plk.map  -e $id -o $sample_name  -p 7.16; done
# -p 是显著性p值的-log10(p threshold)的数,p值显著性可以用FDR来控制

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！