介绍如何使用SNP数据进行GWAS分析，软件选择为emmax和tassel。

软件选择

emmax
plink
qqman
CMplot

数据准备

vcf文件：all_snp.vcf
表型数据文件：sample.table

sample.table

参考脚本

vcf格式转tped格式

plink --vcf ./all_snp.vcf   \  # 输入文件
--maf 0.05 --geno 0.1  \ #过滤
--recode12   \#输出格式
--output-missing-genotype 0  \ #指定缺失数据用0表示
--transpose   \ #指定输出的tped格式
--out snp_filter  \  #输出文件前缀
--set-missing-var-ids @:#  \
--allow-extra-chr   

结果文件

表型数据准备
表型数据的样品id要求有两列，且id的排列顺序与tfam文件的id排列顺序一致

计算亲缘关系矩阵

/home/software/emmax/emmax-kin-intel64 -v  \#打开日志模式
-d 10  \#指定精度
-o ./pop.kinship  \#指定输出文件名称
snp_filter   #指定输入SNP数据前缀

gwas分析

/home/software/emmax/emmax-intel64 -v -d 10  \ #
-t snp_filter  \ # 输出snp数据的前缀
-p sample.table  \ # 表型数据
-k pop.kinship   \ # 亲缘关系矩阵
-o emmax.out  \ # 输出文件前缀
1> emmax.log 2>emmax.error

结果文件为emmax.out.ps，最后一列为pvalue

emmax.out.ps

gwas结果绘图

#生成绘图的输入文件
paste  snp_filter.map  emmax.out.ps | awk  'BEGIN{print "SNP\tCHR\tBP\tP"}{if($2==$5){print $2"\t"$1"\t"$4"\t"$NF}}'  > emmax.out.ps.manht_input

输入文件需要有4列，snp名称、染色体、位置、pvalue

emmax.out.ps.manht_input

Rscript ./manhattan_cmplot.R  emmax.out.ps.manht_input  emmax.out.ps.manht_figure

曼哈顿图
QQ图

欢迎关注Bioinfor生信云！